text:
export PathTracker,
    track!,
    track,
    PathTrackerStatus,
    is_at_infinity,
    is_success,
    is_tracking,
    is_failed,
    is_terminated_callback,
    is_invalid_startvalue,
    pathtracker,
    pathtracker_startsolutions,
    PathResult,
    solution,
    start_solution,
    accuracy,
    residual,
    winding_number,
    multiplicity,
    condition_jacobian,
    path_number,
    is_finite,
    is_singular,
    is_nonsingular,
    is_real,
    is_projective,
    is_affine

const pathtracker_supported_keywords = [
    :at_infinity_check,
    :endgame_start,
    :endgame_start_callback,
    :max_winding_number,
    :min_accuracy,
    :min_cond_eg,
    :min_step_size_before_eg,
    :min_step_size_eg,
    :s_always_consider_valuation,
    :samples_per_loop,
    :precision_strategy,
]

############
## STATUS ##
############

module PathTrackerStatus
"""
    enum PathTrackerStatus

The possible states a `PathTracker` can be in:

* `PathTrackerStatus.tracking`
* `PathTrackerStatus.success`
* `PathTrackerStatus.at_infinity`
* `PathTrackerStatus.excess_solution`
* `PathTrackerStatus.post_check_failed`
* `PathTrackerStatus.terminated_accuracy_limit`
* `PathTrackerStatus.terminated_ill_conditioned`
* `PathTrackerStatus.terminated_invalid_startvalue`
* `PathTrackerStatus.terminated_max_winding_number`
* `PathTrackerStatus.terminated_max_iters`
* `PathTrackerStatus.terminated_step_size_too_small`
"""
@enum states begin
    tracking
    success
    at_infinity
    terminated_accuracy_limit
    terminated_invalid_startvalue
    terminated_ill_conditioned
    terminated_max_iters
    terminated_max_winding_number
    terminated_step_size_too_small
    post_check_failed
    excess_solution
end
end

"""
    path_tracker_status(code::CoreTrackerStatus.states)

Construct a [`PathTrackerStatus.states`](@ref) from a [`CoreTrackerStatus.states`](@ref).
"""
function path_tracker_status(code::CoreTrackerStatus.states)
    if code == CoreTrackerStatus.success
        return PathTrackerStatus.success
    elseif code == CoreTrackerStatus.terminated_invalid_startvalue
        return PathTrackerStatus.terminated_invalid_startvalue
    elseif code == CoreTrackerStatus.terminated_maximal_iterations
        return PathTrackerStatus.terminated_max_iters
    elseif code == CoreTrackerStatus.terminated_step_size_too_small
        return PathTrackerStatus.terminated_step_size_too_small
    elseif code == CoreTrackerStatus.terminated_ill_conditioned
        return PathTrackerStatus.terminated_ill_conditioned
    else
        return PathTrackerStatus.tracking
    end
end

"""
    is_success(status::PathTrackerStatus.states)

Returns `true` if `status` indicates a success in tracking.
"""
is_success(status::PathTrackerStatus.states) = status == PathTrackerStatus.success

"""
    is_success(status::PathTrackerStatus.states)

Returns `true` if `status` indicates that a path diverged towards infinity.
"""
is_at_infinity(status::PathTrackerStatus.states) = status == PathTrackerStatus.at_infinity

"""
    is_tracking(status::PathTrackerStatus.states)

Returns `true` if `status` indicates the tracking is not going on.
"""
is_tracking(status::PathTrackerStatus.states) = status == PathTrackerStatus.tracking

"""
    is_invalid_startvalue(status::PathTrackerStatus.states)

Returns `true` if the provided start value was not valid.
"""
is_invalid_startvalue(status::PathTrackerStatus.states) =
    status == PathTrackerStatus.terminated_invalid_startvalue

"""
    is_terminated_callback(status::PathTrackerStatus.states)

Returns `true` if the provided callback indicated a termination of the path.
"""
is_terminated_callback(status::PathTrackerStatus.states) =
    status == PathTrackerStatus.terminated_callback

#############
## Options ##
#############
@enum PrecisionStrategy begin
    PREC_STRATEGY_NEVER
    PREC_STRATEGY_FINITE
    PREC_STRATEGY_ALWAYS
end

mutable struct PathTrackerOptions
    endgame_start::Float64
    at_infinity_check::Bool
    min_accuracy::Float64
    min_cond_eg::Float64
    min_step_size_before_eg::Float64
    min_step_size_eg::Float64
    precision_strategy::PrecisionStrategy
    s_always_consider_valuation::Float64
end

Base.show(io::IO, opts::PathTrackerOptions) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::PathTrackerOptions) = opts

###########
## STATE ##
###########

mutable struct PathTrackerState{AV<:AbstractVector}
    status::PathTrackerStatus.states
    valuation::Valuation
    prediction::AV
    solution::AV
    intermediate_solution::AV
    solution_accuracy::Float64
    solution_cond::Float64
    solution_residual::Float64
    winding_number::Union{Nothing,Int}
    max_winding_number_hit::Bool
    s::Float64
    eg_started::Bool
end

function PathTrackerState(x::AbstractVector; at_infinity_check::Bool = true)
    status = PathTrackerStatus.tracking
    valuation = Valuation(x; affine = at_infinity_check)
    prediction = copy(x)
    solution = copy(x)
    intermediate_solution = copy(x)
    solution_accuracy = NaN
    solution_cond = NaN
    solution_residual = NaN
    winding_number = nothing
    max_winding_number_hit = false
    s = 0.0
    eg_started = false
    PathTrackerState(
        status,
        valuation,
        prediction,
        solution,
        intermediate_solution,
        solution_accuracy,
        solution_cond,
        solution_residual,
        winding_number,
        max_winding_number_hit,
        s,
        eg_started,
    )
end

Base.show(io::IO, S::PathTrackerState) = print_fieldnames(io, S)
Base.show(io::IO, ::MIME"application/prs.juno.inline", S::PathTrackerState) = S

function init!(state::PathTrackerState, s::Float64)
    state.status = PathTrackerStatus.tracking
    init!(state.valuation)
    state.prediction .= 0.0
    state.solution .= 0.0
    state.intermediate_solution .= NaN
    state.solution_accuracy = NaN
    state.solution_cond = NaN
    state.solution_residual = NaN
    state.winding_number = nothing
    state.max_winding_number_hit = false
    state.s = s
    state.eg_started = false
end

##################
## PATH TRACKER ##
##################

abstract type AbstractPathTracker end
Base.broadcastable(T::AbstractPathTracker) = Ref(T)

"""
    PathTracker

The `PathTracker` combines the path tracking (with [`CoreTracker`](@ref)) with an *endgame*
routine. The *endgame*  is the name for special algorithms the end of the path tracking.
These enable to detect if a path is diverging or if it ends in a singular solution.
The `PathTracker` is more opinionated than the `CoreTracker` and implements additional logic
to handle numerical difficult solutions. In particular it always reparametrizes a solution
path to use a logarithmic time scale, i.e., ``x(t) → x(e^{-s})`` and we track for ``s`` from
``0`` to ``∞``.

In order to construct a `PathTracker` it is recommended to use the [`pathtracker`](@ref) and
[`pathtracker_startsolutions`](@ref) helper functions.
With a `PathTracker` constructed you can track a single path using the [`track`](@ref) method.
The result of this will be a [`PathResult`](@ref).

The `PathTracker` assumes that the provided homotopy `H` is defined in such a way that
`H(x,1)` is the start system and `H(x,0)` the target system.
During the path tracking an approximation of the valuation of a Puiseux series expansion of
the solution is computed. This is used to decide whether a path is diverging.
Before a certain treshold (`s_always_consider_valuation`) this approximation is only trusted
if a path gives some sign of ill-conditioning.
To compute singular solutions the *Cauchy endgame* is used which is based on
Cauchy's integral formula. For this we have to track solutions along a loop around the origin.
The number of loops necessary to arrive back at the start point is called the *winding number*.

## Options

The path tracker accepts all options which [`CoreTracker`](@ref) accepts. Furthermore the
following options are accepted:

* `at_infinity_check`: This is true if the provided start system is an affine polynomial system.
* `endgame_start` (default `2.0`): The value of `s` where the endgame is started earliest.
* `max_winding_number` (default `12`): The maximal winding number tried in the Cauchy endgame.
* `min_accuracy` (default `1e-5`): The `PathTracker` automatically lowers the desired
  accuracy automatically to `min_accuracy` if path tracking would fail othwerwise (since
  the desired cannot be reached)
* `min_cond_eg` (default `1e5`): The minimal condition number before the Cauchy endgame is
  started or a path is cut off.
* `min_step_size_before_eg` (default `exp2(-40)`): The minimal allowed step size before the
  endgame starts.
* `min_step_size_eg` (default `exp2(-120)`): The minimal allowed step size during the endgame.
  This is also control what is considered `∞` for the path tracking.
* `s_always_consider_valuation` (default `-log(1e-16)`) A threshold after which we always consider
  the valuation.
* `samples_per_loop` (default `8`): The number of samples used during the endgame.
* `precision_strategy` (default `:adaptive_finite`): This controls whether `H(x,t)` is
  possibly evaluated with higher than machine precision during the endgame
  (the Jacobian is always computed with machine precision). The `:adaptive_finite` allows
  this only if we are optimistic that we can still obtain a finite solution.
  Other options are `:adaptive_never` where this is never allowed and `:adaptive_always`
  where it is always enabled.
"""
struct PathTracker{
    AV<:AbstractVector{Complex{Float64}},
    Prob<:AbstractProblem,
    CT<:CoreTracker{Float64,AV},
} <: AbstractPathTracker
    problem::Prob
    core_tracker::CT
    endgame::CauchyEndgame{AV}
    state::PathTrackerState{AV}
    options::PathTrackerOptions
    # We modify options of the core tracker dynamically. Therefore we need to be able
    # so start fresh for a new path
    default_ct_options::CoreTrackerOptions
end

function PathTracker(
    prob::AbstractProblem,
    x::AbstractVector{<:Number};
    min_step_size = 1e-30,
    kwargs...,
)

    core_tracker_kwargs, options = splitkwargs(kwargs, coretracker_supported_keywords)
    core_tracker = CoreTracker(
        prob,
        x;
        log_transform = true,
        predictor = MixPredictor(Pade21(), Heun()),
        min_step_size = min_step_size,
        core_tracker_kwargs...,
    )
    PathTracker(prob, core_tracker; options...)
end

function PathTracker(
    prob::Problem,
    core_tracker::CoreTracker;
    at_infinity_check = default_at_infinity_check(prob),
    endgame_start::Float64 = -log(0.1),
    max_winding_number::Int = 12,
    min_accuracy::Float64 = max(core_tracker.options.accuracy, 1e-5),
    min_cond_endgame::Float64 = 1e4,
    min_step_size_before_eg::Float64 = exp2(-40),
    min_step_size_eg = exp2(-120),
    s_always_consider_valuation::Float64 = -log(1e-16),
    samples_per_loop::Int = 8,
    precision_strategy::Symbol = :adaptive_finite,
)
    state = PathTrackerState(current_x(core_tracker); at_infinity_check = at_infinity_check)
    endgame = CauchyEndgame(
        current_x(core_tracker);
        samples_per_loop = samples_per_loop, max_winding_number = max_winding_number,
    )
    default_ct_options = copy(core_tracker.options)

    options = PathTrackerOptions(
        endgame_start,
        at_infinity_check,
        min_accuracy,
        min_cond_endgame,
        min_step_size_before_eg,
        min_step_size_eg,
        make_precision_strategy(precision_strategy),
        s_always_consider_valuation,
    )

    PathTracker(prob, core_tracker, endgame, state, options, default_ct_options)
end

default_at_infinity_check(prob::Problem{AffineTracking}) = true
default_at_infinity_check(prob::Problem{ProjectiveTracking}) = homvars(prob) !== nothing

always_true(_, _) = true

function make_precision_strategy(precision_strategy::Symbol)
    if precision_strategy == :adaptive_finite
        PREC_STRATEGY_FINITE
    elseif precision_strategy == :adaptive_never
        PREC_STRATEGY_NEVER
    elseif precision_strategy == :adaptive_always
        PREC_STRATEGY_ALWAYS
    else
        ArgumentError(
            "Unknown argument `precision_strategy = $(precision_strategy)`." *
            "Possible values are `:adaptive_finite`, `:adaptive_never`" *
            " or `:adaptive_always`.",
        ) |> throw
    end
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", PT::PathTracker) = PT
function Base.show(io::IO, S::PathTracker{AV}) where {AV}
    println(io, "PathTracker with solution type $AV")
end

seed(PT::PathTracker) = seed(PT.problem)

"""
    init!(tracker::PathTracker, x)

Prepare the `PathTracker`` `tracker` to track a path with start solution `x`.
"""
function init!(
    tracker::PathTracker,
    x;
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
)
    options, ct_options = tracker.options, tracker.core_tracker.options
    copy!(ct_options, tracker.default_ct_options)
    if accuracy !== nothing
        ct_options.accuracy = accuracy
    end
    if max_corrector_iters !== nothing
        ct_options.max_corrector_iters = max_corrector_iters
    end
    if options.precision_strategy == PREC_STRATEGY_ALWAYS ||
       options.precision_strategy == PREC_STRATEGY_FINITE
        ct_options.precision = PRECISION_ADAPTIVE
    end
    if tracker.core_tracker.predictor isa MixPredictorCache
        tracker.core_tracker.predictor.predictor1_active[] = true
    end

    init!(tracker.core_tracker, x, 0.0, options.endgame_start)
    init!(tracker.state, 0.0)

    # before endgame don't track the condition number
    # and we only update the limiting accuracy etc after a step failed
    ct_options.track_cond = false
    ct_options.steps_jacobian_info_update = -1
    ct_options.min_step_size = options.min_step_size_before_eg

    tracker
end

"""
    step!(tracker::PathTracker)

Perform a single step of the `PathTracker`.
"""
function step!(tracker::PathTracker)
    @unpack core_tracker, state, options, endgame = tracker
    step!(core_tracker)

    state.s = s = real(current_t(core_tracker))
    ct_status = status(core_tracker)

    # If s > 18.42 we have exp(-18.42) ≈ 1e-8 and we get to the range where numerical
    # differentation with the Padé predictor can become tricky.
    # Therefore we switch to the other predictor
    s > 18.42 && activate_predictor2!(core_tracker)

    if is_tracking(ct_status)
        state.eg_started || return state.status
        # If we didn't move forward there is nothing to do
        core_tracker.state.last_step_failed && return state.status

        # update the valuation
        update!(state.valuation, core_tracker)

        # We only care about the valuation if the path is starting to get worse
        # or we passed a certain threshold.
        cond_bad = cond(core_tracker) > options.min_cond_eg
        consider_always = s > options.s_always_consider_valuation
        if !(consider_always || cond_bad)
            return state.status
        end

        # Judge the current valuation to determine how to proceed
        verdict = judge(state.valuation; tol = 1e-2, tol_at_infinity = 1e-4)
        if verdict == VAL_AT_INFINITY &&
           (consider_always || cond_bad) && options.at_infinity_check
            state.status = PathTrackerStatus.at_infinity
        # If we expect a finite value and there is some ill conditioning let's do the
        # Cauchy endgame
        elseif verdict == VAL_FINITE && cond_bad
            # Perform endgame to estimate singular solution
            @label run_cauchy_eg
            activate_predictor2!(core_tracker)
            retcode, m, p_accuracy = predict!(state.prediction, core_tracker, endgame)

            if retcode == CAUCHY_SUCCESS
                # We only accept a result of the Cauchy endgame only if two consecutive
                # loops resulted in the same number of loop
                m′ = state.winding_number
                if m′ === nothing
                    state.winding_number = m
                    state.solution .= state.prediction
                elseif m′ == m &&
                       norm(core_tracker)(state.solution, state.prediction) ≤ 4p_accuracy
                    @label cauchy_eg_success
                    state.status = PathTrackerStatus.success
                    state.solution_accuracy = p_accuracy
                    converged, s_accuracy, s_cond, s_res = check_converged!(
                        state.solution,
                        core_tracker,
                        state.prediction,
                        Inf;
                        tol = tracker.default_ct_options.accuracy,
                    )
                    state.solution_cond = s_cond
                    state.solution_residual = s_res
                    if converged && s_cond < 1e8 && s_accuracy < p_accuracy
                        state.solution_accuracy = s_accuracy
                    else
                        state.solution .= state.prediction
                    end
                else
                    state.solution .= state.prediction
                    state.winding_number = m
                end

            elseif retcode == CAUCHY_TERMINATED_ACCURACY_LIMIT
                if core_tracker.options.precision == PRECISION_FIXED_64
                    core_tracker.options.precision = PRECISION_ADAPTIVE
                    @goto run_cauchy_eg
                end

            elseif options.precision_strategy == PREC_STRATEGY_NEVER
                state.status = PathTrackerStatus.terminated_accuracy_limit
                # If we are here, we have a precision strategy which allows to use
                # higher precision. Therefore switch to it if this not yet happened
                # and retry to cauchy eg

            elseif retcode == CAUCHY_TERMINATED_MAX_WINDING_NUMBER
                # Terminate if we hit twice max winding number and make one step pause
                if state.max_winding_number_hit
                    state.status = PathTrackerStatus.terminated_max_winding_number
                else
                    state.max_winding_number_hit = true
                end
            elseif retcode == CAUCHY_TERMINATED_MAX_ITERS
                state.status = PathTrackerStatus.terminated_max_iters
            else
                state.status = PathTrackerStatus.terminated_ill_conditioned
            end
        end

    # We split the tracking in two parts. Pre endgame and endgame. Therefore we have to
    # handle the success case twice
    elseif is_success(ct_status) && !state.eg_started
        state.intermediate_solution .= current_x(core_tracker)
        state.eg_started = true
        # set min_step size
        core_tracker.options.min_step_size = options.min_step_size_eg
        # setup path tracker to continue tracking
        t₀ = -log(4 * options.min_step_size_eg)
        init!(
            tracker.core_tracker,
            current_x(core_tracker),
            options.endgame_start,
            t₀;
            loop = true,
        )
        # start keeping track of the condition number during the endgame
        tracker.core_tracker.options.track_cond = true
        # update the limiting accuracy and cond every other step
        tracker.core_tracker.options.steps_jacobian_info_update = 2
        # update the valuation
        update!(state.valuation, core_tracker)

        # If our strategy is to allow adaptive precision only for finite values
        # then we start with not allowing adaptive precision during the endgame
        # and wait to hit the accuracy limit.
        if options.precision_strategy == PREC_STRATEGY_FINITE
            tracker.core_tracker.options.precision = PRECISION_FIXED_64
        end

    elseif is_success(ct_status) ||
           ct_status == CoreTrackerStatus.terminated_step_size_too_small
        converged, s_acc, s_cond, s_res = check_converged!(
            state.solution,
            core_tracker,
            current_x(core_tracker),
            Inf;
            tol = tracker.default_ct_options.accuracy,
        )
        if converged
            state.status = PathTrackerStatus.success
            state.solution_accuracy = s_acc
            state.solution_cond = s_cond
            state.solution_residual = s_res
            state.winding_number = nothing
            state.s = Inf
        elseif ct_status == CoreTrackerStatus.terminated_step_size_too_small
            state.status = PathTrackerStatus.terminated_step_size_too_small
        else
            state.status = PathTrackerStatus.post_check_failed
        end

    elseif ct_status == CoreTrackerStatus.terminated_accuracy_limit
        # First update the valuation  and check whether we are good
        update!(state.valuation, core_tracker)
        # we put less requirement on tol_at_infinity
        verdict = judge(state.valuation; tol = 1e-3, tol_at_infinity = 1e-4)
        # We also terminate at this point if the valuation indicates that we are
        # substantially less than 0
        if verdict == VAL_AT_INFINITY && options.at_infinity_check
            state.status = PathTrackerStatus.at_infinity
        # Now, we have to differentiate 3 different cases:
        # 1) Current accuracy is larger than the defined minimal accuracy
        #     -> decrease accuracy to minimal accuracy
        # 2) We are already at minimal accuracy and allow adaptive precision for finite
        #    values
        #  a) it looks like we could end up with something finite
        #     -> enable adaptive precision
        #  b) else
        #     -> terminate
        # 3) We are already at minimal accuracy and don't allow adaptive precision
        #     -> terminate
        # 1)
        elseif core_tracker.options.accuracy > options.min_accuracy
            # TODO: Currently we never increase this again
            core_tracker.options.accuracy = options.min_accuracy
            core_tracker.state.status = CoreTrackerStatus.tracking
        # 2)
        elseif options.precision_strategy == PREC_STRATEGY_FINITE
            # 2.a) We say that a path could still become finite if it is not classified
            #      as going to infinity for a loose tolerance
            #      Also if the current valuation is less than -1 (so significantly less than 0)
            #      we do not continue the tracking
            loose_verdict = judge(state.valuation; tol = 1e-2, tol_at_infinity = 1e-2)
            if loose_verdict != VAL_AT_INFINITY
                core_tracker.options.precision = PRECISION_ADAPTIVE
                core_tracker.state.status = CoreTrackerStatus.tracking
            # 2.b)
            else
                state.status = PathTrackerStatus.terminated_accuracy_limit
            end
        # 3)
        else
            state.status = PathTrackerStatus.terminated_accuracy_limit
        end
    elseif state.winding_number !== nothing
        p_accuracy = 1e-5
        @goto cauchy_eg_success
    else
        state.status = path_tracker_status(ct_status)
    end

    if !is_success(state.status) && !is_tracking(state.status)
        state.solution .= current_x(core_tracker)
        state.solution_cond = LA.cond(core_tracker)
    end

    state.status
end

"Update the valuation with the current state."
update!(val::Valuation, T::CoreTracker) =
    update!(val, T.state.x, T.state.ẋ, real(current_t(T)), T.predictor)

function activate_predictor2!(core_tracker::CoreTracker)
    if core_tracker.predictor isa MixPredictorCache
        core_tracker.predictor.predictor1_active[] = false
    end
    nothing
end

"""
    check_converged!(y, tracker::CoreTracker, x::AbstractVector, t)

Check that the CoreTracker is converged also converged at `(x,t)`. Returns a named tuple
`(converged, accuracy, cond)` and stores a (possibly) refined solution into `y`.
"""
function check_converged!(
    y,
    tracker::CoreTracker,
    x::AbstractVector,
    t::Number;
    tol::Float64 = 1e-6,
)
    init!(tracker.state.norm, x)
    # Make sure to evaluate the Jacobian only *once*. Otherwise it can happen at singuar
    # solutions that we bounce away from a good solution.
    result =
        correct!(y, tracker, x, t; tol = tol, min_iters = 2, max_iters = 3, full_steps = 1)
    if is_converged(result)
        tracker.state.residual .= abs.(tracker.corrector.r)
    else
        evaluate!(tracker.corrector.r, tracker.homotopy, x, t)
        tracker.state.residual .= abs.(tracker.corrector.r)
    end
    res = LA.norm(tracker.state.residual, InfNorm())
    # Fix the residual to the norm otherwise we can scale away zero rows.
    tracker.state.residual .= res
    # If the accuracy is NaN we assume that we dived by 0 -> Jacobian singular
    if isnan(result.accuracy)
        cond_jac = Inf
    else
        cond_jac = cond!(tracker.state.jacobian, tracker.state.norm, tracker.state.residual)
    end
    (
        converged = is_converged(result),
        accuracy = result.accuracy,
        cond = cond_jac,
        res = res,
    )
end

function Base.iterate(tracker::PathTracker, state::Int = 0)
    if is_tracking(tracker.state.status)
        step!(tracker)
        tracker, state + 1
    else
        nothing
    end
end

####################
## QUERYING STATE ##
####################

status(tracker::PathTracker) = tracker.state.status

"""
    solution(tracker::PathTracker)

Obtain the solution computed by the `PathTracker`.
"""
solution(tracker::PathTracker) = pull_back(tracker.problem, tracker.state.solution)

"""
    norm(tracker::PathTracker)

Obtain the norm used by the `PathTracker`.
"""
LA.norm(tracker::PathTracker) = LA.norm(tracker.core_tracker)

"""
    winding_number(tracker::PathTracker)

Obtain the estimate of the winding number computed by the `PathTracker`. Returns `nothing`
if the Cauchy endgame was not run.
"""
winding_number(tracker::PathTracker) = winding_number(tracker.state)
winding_number(state::PathTrackerState) = state.winding_number

###############################
## Convencience constructors ##
###############################

const pathtracker_startsolutions_supported_keywords = [
    problem_startsolutions_supported_keywords
    coretracker_supported_keywords
    pathtracker_supported_keywords
]

"""
    pathtracker_startsolutions(args...; kwargs...)

Construct a [`PathTracker`](@ref) and start solutions in the same way [`solve`](@ref) does it.
This also takes the same input arguments as `solve`. This is convenient if you want
to investigate single paths.
"""
function pathtracker_startsolutions(args...; kwargs...)
    invalid = invalid_kwargs(kwargs, pathtracker_startsolutions_supported_keywords)
    check_kwargs_empty(invalid, pathtracker_startsolutions_supported_keywords)
    supported, rest = splitkwargs(kwargs, problem_startsolutions_supported_keywords)
    prob, startsolutions = problem_startsolutions(args...; supported...)
    construct_tracker(prob, startsolutions; rest...), startsolutions
end

function construct_tracker(prob::Problem, startsolutions; kwargs...)
    PathTracker(prob, start_solution_sample(startsolutions); kwargs...)
end

"""
    pathtracker(args...; kwargs...)

Construct a [`PathTracker`](@ref) in the same way [`solve`](@ref) does it.
This also takes the same input arguments as `solve` with the exception that you do not need
to specify startsolutions.

## Examples

### Obtain single solution
We want to construct a path tracker to track a parametrized system `f` with parameters `p`
from the parameters `a` to `b`.
```julia
tracker = pathtracker(f, parameters=p, p₁=a, p₀=b)
```
You then can obtain a single solution at `b` by using
```julia
x_b = solution(track(tracker, x_a))
```
"""
pathtracker(args...; kwargs...) = first(pathtracker_startsolutions(args...; kwargs...))

#################
## PATH RESULT ##
#################

"""
    PathResult{V<:AbstractVector}

A `PathResult` is the result of tracking of a path with [`track`](@ref) using a
[`PathTracker`](@ref).

# Fields
* `return_code`: See the list of return codes below.
* `solution::V`: The solution vector.
* `t::Float64`: The value of `t` at which `solution` was computed. Note that if
  `return_code` is `:at_infinity`, then `t` is the value when this was decided.
* `accuracy::Float64`: An approximation of ``||x-x̄||`` where ``x`` is the computed solution
  and ``x̄`` is the true solution.
* `residual::Float64`: The value of the infinity-norm of `H(solution, 0)`.
* `multiplicity::Union{Nothing, Int}` is the multiplicity of the `solution`. This is
  only assigned by [`singular`](@ref).
* `condition_jacobian::Union{Nothing, Float64}`: This is the condition number of the
  row-equilibrated Jacobian at the solution. A high condition number indicates a singularity.
* `winding_number:Union{Nothing, Int}`: The estimated winding number. This is a lower bound
  on the multiplicity of the solution.
* `start_solution::Union{Nothing, Int}`: The start solution of the path.
* `accepted_steps::Int`: The number of accepted steps during the path tracking.
* `rejected_steps::Int`: The number of rejected steps during the path tracking.
* `valuation::Union{Nothing, Vector{Float64}}`: An approximation of the valuation of the
  Puiseux series expansion of `x(t)`.
* `valuation_accuracy::Union{Nothing, Vector{Float64}}`: An estimate of the accuracy of the
  valuation of the Puiseux series expansion of `x(t)`.

# Return codes

These is the list of possible return codes:

* `:success`: The `PathTracker` obtained a solution.
* `:at_infinity`: The `PathTracker` stopped the tracking of the path since it determined
  that that path is diverging towards infinity.
* `:terminated_callback`: One of the optional `PathTracker` callbacks terminated the tracking.
* `:terminated_max_iters`: The `PathTracker` terminated since it reached the limit accuracy.
* `:terminated_invalid_startvalue`: The `PathTracker` terminated since the provided start
  value is invalid.
* `:terminated_step_size_too_small`: The `PathTracker` terminated since the step size
  became smaller than the provided threshold.
* `:terminated_accuracy_limit`: The `PathTracker` terminated since the problem was too
  ill-conditioned to be tracked further with the desired minimal accuracy.
* `:terminated_ill_conditioned`: The `PathTracker` terminated since the Jacobian of the
  homotopy was too ill-conditioned.
* `:post_check_failed`: The verification of a non-singular solution failed.

# Constructors

     PathResult(tracker::PathTracker,
                start_solution=nothing,
                path_number::Union{Nothing,Int}=nothing;
                details=:default)

Construct a `PathResult` using the current state of the `PathTracker`.
Possible values for `details` are `:minimal` (minimal details), `:default` (default) and
`:extensive` (all information possible).
"""
struct PathResult{V<:AbstractVector}
    return_code::Symbol
    solution::V
    t::Float64
    accuracy::Float64
    residual::Float64
    condition_jacobian::Float64
    winding_number::Union{Nothing,Int}
    path_number::Base.RefValue{Union{Nothing,Int}}
    start_solution::Union{Nothing,V} # level 1+
    intermediate_solution::Union{Nothing,V}
    multiplicity::Base.RefValue{Union{Nothing,Int}} # only assigned by singular(Result)
    # performance stats
    accepted_steps::Int
    rejected_steps::Int
    valuation::Union{Nothing,Vector{Float64}} # level 2+
    valuation_accuracy::Union{Nothing,Vector{Float64}} # level 2+
end

function PathResult(
    tracker::PathTracker,
    start_solution::Union{Nothing,AbstractVector} = nothing,
    path_number::Union{Nothing,Int} = nothing;
    details::Symbol = :default,
)
    @unpack state, core_tracker = tracker
    details_level = detailslevel(details)

    return_code = Symbol(state.status)
    x = solution(tracker)
    t = return_code == :success ? 0.0 : exp(-state.s)
    accuracy = state.solution_accuracy
    residual = state.solution_residual
    condition_jac = state.solution_cond
    windingnumber = state.winding_number
    # this needs to be assigned manually
    multiplicity = Base.RefValue{Union{Nothing,Int}}(nothing)

    startsolution = nothing
    if details_level ≥ 1 && start_solution !== nothing
        # mimic the behaviour in track! to get a start solution of the same type as x
        embed!(core_tracker.state.x, tracker.problem, start_solution)
        startsolution = pull_back(tracker.problem, core_tracker.state.x; regauge = false)
    end
    intermediate_solution = nothing
    if !isnan(state.intermediate_solution[1])
        intermediate_solution = pull_back(tracker.problem, state.intermediate_solution)
    end

    accepted_steps = core_tracker.state.accepted_steps
    rejected_steps = core_tracker.state.rejected_steps

    valuation = valuation_accuracy = nothing
    if details_level == 2
        valuation = copy(tracker.state.valuation.ν)
        valuation_accuracy = abs.(tracker.state.valuation.ν̇)
    end

    PathResult(
        return_code,
        x,
        t,
        accuracy,
        residual,
        condition_jac,
        windingnumber,
        Base.RefValue{Union{Nothing,Int}}(path_number),
        startsolution,
        intermediate_solution,
        multiplicity,
        accepted_steps,
        rejected_steps,
        valuation,
        valuation_accuracy,
    )
end

function detailslevel(details::Symbol)
    if details == :minimal
        return 0
    elseif details == :extensive
        return 2
    else
        return 1
    end
end

"""
    result_type(tracker::PathTracker)

Returns the type of result `track` will return.
"""
result_type(tracker::PathTracker) = PathResult{typeof(solution(tracker))}

function Base.show(io::IO, r::PathResult{AV}) where {AV}
    iscompact = get(io, :compact, false)
    if iscompact || haskey(io, :typeinfo)
        println(io, "• return_code → :$(r.return_code)")
        if r.return_code ≠ :success
            println(io, " • t → $(r.t)")
        end
        println(io, " • solution → ", r.solution)
        r.accuracy !== nothing && println(
            io,
            " • accuracy → $(Printf.@sprintf "%.3e" r.accuracy)",
        )
        winding_number(
            r,
        ) !== nothing && println(io, " • winding_number → $(winding_number(r))")
        if multiplicity(r) !== nothing
            println(io, " • multiplicity → $(multiplicity(r))")
        end
        path_number(r) !== nothing && println(io, " • path_number: ", path_number(r))
    else
        header = "PathResult{$AV}"
        compact_io = IOContext(io, :compact => true)
        println(io, "PathResult{$AV}")
        println(io, "="^length(header))
        path_number(r) !== nothing && println(io, " • path_number → ", path_number(r))
        println(io, " • return_code → :$(r.return_code)")
        if r.return_code ≠ :success
            println(compact_io, " • t → $(r.t)")
        end
        println(compact_io, " • solution → ", r.solution)
        r.accuracy !== nothing && println(compact_io, " • accuracy → ", r.accuracy)
        r.residual !== nothing && println(compact_io, " • residual → ", r.residual)
        if r.condition_jacobian !== nothing
            println(compact_io, " • condition_jacobian → ", r.condition_jacobian)
        end

        if multiplicity(r) !== nothing
            println(io, " • multiplicity → $(multiplicity(r))")
        end
        if winding_number(r) !== nothing
            println(io, " • winding_number → $(winding_number(r))")
        end
        if r.valuation !== nothing
            println(compact_io, " • valuation → ", r.valuation)
        end
        if r.valuation_accuracy !== nothing
            println(compact_io, " • valuation_accuracy → ", r.valuation_accuracy)
        end
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::PathResult) = x

"""
    solution(r::PathResult)

Get the solution of the path.
"""
solution(r::PathResult) = r.solution

"""
    start_solution(r::PathResult)

Get the start solution of the path.
"""
start_solution(r::PathResult) = r.start_solution

"""
    accuracy(r::PathResult)

Get the accuracy of the solution. This is an estimate of the (relative) distance to the
true solution.
"""
accuracy(r::PathResult) = r.accuracy

"""
    residual(r::PathResult)

Get the residual of the solution ``x`` of the path, i.e., ``||H(x, 0)||₂``.
"""
residual(r::PathResult) = r.residual

"""
    winding_number(r::PathResult)

Get the winding number of the solution of the path. Returns `nothing` if it wasn't computed.
"""
winding_number(r::PathResult) = r.winding_number

"""
    multiplicity(P::PathResult)

Returns the multiplicity of `P`.
"""
multiplicity(P::PathResult) = P.multiplicity[]

"""
    condition_jacobian(r::PathResult)

Return the condition number of the Jacobian of the result.
"""
condition_jacobian(r::PathResult) = r.condition_jacobian

"""
    cond(r::PathResult)

Return the condition number of the Jacobian of the result.
"""
LA.cond(r::PathResult) = r.condition_jacobian

"""
    path_number(r::PathResult)

The number of the path.
"""
path_number(r::PathResult) = r.path_number[]

"""
    is_success(r::PathResult)

Checks whether the path is successfull.
"""
is_success(r::PathResult) = r.return_code == :success

"""
    is_failed(r::PathResult)

Checks whether the path failed.
"""
is_failed(r::PathResult) = !(r.return_code == :at_infinity || r.return_code == :success)

"""
    is_at_infinity(r::PathResult)

Checks whether the path goes to infinity.
"""
is_at_infinity(r::PathResult) = r.return_code == :at_infinity

"""
    is_finite(r::PathResult)

Checks whether the path result is finite.
"""
is_finite(r::PathResult) = r.return_code == :success

# Base fallback
Base.isfinite(r::PathResult) = is_finite(r)

"""
    is_singular(r::PathResult; tol=1e10)

Checks whether the path result is singular. This is true if
the multiplicity is larger than  1 or if the condition number of the Jacobian
is larger than `tol`.
"""
is_singular(r::PathResult; tol = 1e10) = is_singular(r, tol)
function is_singular(r::PathResult, tol::Real)
    (
        unpack(r.condition_jacobian, 1.0) > tol ||
        unpack(multiplicity(r), 1) > 1 || unpack(winding_number(r), 1) > 1
    ) && is_success(r)
end

"""
    is_nonsingular(r::PathResult; tol=1e10)

Checks whether the path result is non-singular. This is true if
it is not singular.
"""
is_nonsingular(r::PathResult; kwargs...) = !is_singular(r; kwargs...) && is_success(r)
is_nonsingular(r::PathResult, tol::Real) = !is_singular(r, tol) && is_success(r)

"""
    is_real(r::PathResult; tol=1e-6)

We consider a result as `real` if the 2-norm of the imaginary part of the solution is at most `tol`.
"""
is_real(r::PathResult; tol = 1e-6) = is_real(r, tol)
is_real(r::PathResult, tol::Real) = is_real_vector(r.solution, tol)
# provide fallback since this in in Base
Base.isreal(r::PathResult, tol) = is_real(r, tol)
Base.isreal(r::PathResult; kwargs...) = is_real(r; kwargs...)

"""
    is_projective(r::PathResult)

Return`s true if the solution is a projective vector.
"""
is_projective(r::PathResult) = false
is_projective(r::PathResult{<:PVector}) = true

"""
    is_affine(pathresult)

Return`s true if the solution is an affine vector.
"""
is_affine(r::PathResult) = !is_projective(r)

###########
## track ##
###########

"""
    track(tracker::PathTracker, x₁;
            path_number=nothing,
            details::Symbol=:default,
            start_parameters = nothing,
            target_parameters = nothing)

Track the path `x(t)` with start solution `x₁` from ``1`` towards ``0``.
Returns a [`PathResult`](@ref).

The `details` options controls the level of details of the information available
in the [`PathResult`](@ref).
If `tracker` uses a parameter homotopy you can set the start and target parameters
by setting the corresponding fields.
To investigate the behaviour of a particular take a look at [`path_info`](@ref).
"""
function track(
    tracker::PathTracker,
    x,
    t₁ = nothing,
    t₀ = nothing;
    path_number::Union{Int,Nothing} = nothing,
    details::Symbol = :default,
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
    start_parameters = nothing,
    target_parameters = nothing,
)
    if t₁ !== nothing
        @warn(
            "track(tracker, x, t₁, t₀) is deprecated. `t₁` is always assumed to be `1` " *
            " and `t₀` set automatically.",
        )
    end

    start_parameters !== nothing && start_parameters!(tracker, start_parameters)
    target_parameters !== nothing && target_parameters!(tracker, target_parameters)

    track!(tracker, x; accuracy = accuracy, max_corrector_iters = max_corrector_iters)
    PathResult(tracker, x, path_number; details = details)
end

"""
    start_parameters!(tracker::PathTracker, p)

Set the start parameters of the homotopy in in `tracker` to `p`.
"""
start_parameters!(T::PathTracker, p) = start_parameters!(T.core_tracker, p)

"""
    target_parameters!(tracker::PathTracker, p)

Set the target parameters of the homotopy in in `tracker` to `p`.
"""
target_parameters!(T::PathTracker, p) = target_parameters!(T.core_tracker, p)

"""
    track!(tracker::PathTracker, x₁)::PathTrackerStatus.states

Track the path `x(t)` with start solution `x₁` from ``1`` towards ``0``.
Returns a [`PathTrackerStatus.states`](@ref).
"""
function track!(
    tracker::PathTracker,
    x;
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
)
    init!(tracker, x; accuracy = accuracy, max_corrector_iters = max_corrector_iters)
    while is_tracking(tracker.state.status)
        step!(tracker)
    end
    tracker.state.status
end


text:
include("../test_exception.jl")
nparts = (2,2,2)
prun(test_exception,mpi,nparts)

ext: jl


text:
@testset "Examples" begin
    using Markdown: @md_str
    include("examples/quicktutorial.jl")
    include("examples/searchspace.jl")
    include("examples/mutation.jl")
    include("examples/crossover.jl")
    include("examples/fitness.jl")
    include("examples/candidate.jl")
    include("examples/evolution.jl")
    include("examples/iterators.jl")
end


text:
using CuArrays
CuArrays.allowscalar(false)

include("../models/AttentionEmbed/Transformer.jl")

TimeSteps = 100
EmbDims = 128
Vocab = Integer(1e4)
Layers = 8
Stacks = 1

X = map(a->Int32.(a), [[3,7,3,5], [8,6,4,5,2,1], [10,7,5], [11,3]]);
Y = map(a->Int32.(a), [[3,7,3,5], [8,6,4,5,2,1], [10,7,5], [11,3]]);

Input = SequenceFeed(EmbDims, Vocab, TimeSteps)
Encode = Encoding(EmbDims, TimeSteps, Layers, Stacks)
Decode = Decoding(EmbDims, TimeSteps, Layers, Stacks)
Output = DenseProjection(Vocab)

function Model(x::AbstractArray{T, 1}, y::AbstractArray{T, 1}) where T
    EncSeqLens = SequenceLengths(x)
    x = Input(x)
    x, EncMask, context = Encode(x, EncSeqLens)

    DecSeqLens = SequenceLengths(y)
    y = Input(y)
    y = Decode(x, y, DecSeqLens, context, EncMask, Encode.Attn)

    y = Output(y, Input.Emb.E)
    return y
end

@time output = Model(X, Y);
# println(output)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ b129ba7c-953a-11ea-3379-17adae34924c
md"# _Welcome to Pluto!_
We're here to help you tell stories with your code. Let's go! 🚀"

# ╔═╡ 4d88b926-9543-11ea-293a-1379b1b5ae64
md"## Cats
Let's say you're like my grandma, and you have a lot of cats. Our story will be about them."

# ╔═╡ 3e8e381e-953f-11ea-3d3e-71d0fea52560
cat = "Katrin"

# ╔═╡ aeb3a6bc-9540-11ea-0b8f-6d37412bfe68
if cat == "Ks*;lj"
	md"Oh no! Someone messed with my pretty introduction. Change the code above to give our cat a proper name!"
else
	md"**Well done, your cat is called $cat now.** This text gets updated every time you change the name. To see how the magic works, click on the 👁 to the left of this text."
end

# ╔═╡ 611c28fa-9542-11ea-1751-fbdedcfb7690
md"To edit any code, just click on it. When you're done programming, press the ▶ in the lower-right to run the code. You can also use `Shift+Enter` if you are in a hurry."

# ╔═╡ 7a81b97e-fae7-11ea-2ef2-9768a7d4362f
friend = "Ryan"

# ╔═╡ 6f7eecec-9543-11ea-1284-dd52fce3ecca
md"I feel like our cat needs a friend. Let's call them $(friend)."

# ╔═╡ a1a20314-9543-11ea-17de-0b658da18992
if !@isdefined friend
	md"Uh oh, what is this? I forgot to add a cell defining our friend. Can you do it for me?"
else
	md"**Well done!** $cat and $friend are both happy with your performance."
end

# ╔═╡ 90708aaa-fae7-11ea-31ee-d3d84717608a
.1 + .2 ≈ 0.3

# ╔═╡ 0508c34e-fae8-11ea-17f8-332503fbd513
√4.2

# ╔═╡ 72c00416-fae9-11ea-24eb-7dee7217d4eb
2+2

# ╔═╡ f112b662-9543-11ea-3dcb-2906a99b2188
md"A cell is a container for code & output. To add one, click on the ➕ above or below another cell. You can do it wherever you like. After you're done writing code in your cell, remember to run it!"

# ╔═╡ e0642f42-9545-11ea-14ee-fde52cb54ccc
md"## Feeding neighbors
Our cats have some neighbors. Let's involve them in the story too!"

# ╔═╡ 3653b1ac-9546-11ea-2a44-ddd3054636fe
neighbors = [cat, friend, "Smerfetka", "Latte", "Fozzy", "Mongo", "Mango", "Dweeba", "Miffy", "Spiffy", "Biff"]

# ╔═╡ 19ff8d36-9547-11ea-0e08-e5cdd8338673
md"Now, if you're like my grandma, you're feeding the entire neigborhood by yourself. Let's see how many cans of cat food you need to prepare."

# ╔═╡ 270ac49e-9549-11ea-3ffd-71ddaee9f134
md"But what does `confusing_function` do? If you ever need help, click on 📚 **Live docs** in the lower right, and then place your cursor on the code you need help with. It... usually works, we're polishing it right now!"

# ╔═╡ eac62fea-954e-11ea-2768-39ce6f4059ab
begin
	cans_in_stock = "🥫🥫🥫🥫🥫🥫🥫🥫🥫🥫"
	cans_in_stock = "🥫🥫🥫🥫"
end

# ╔═╡ 745a4584-954a-11ea-028e-59011f268ec6
cans_in_stock

# ╔═╡ 55ad7422-954e-11ea-0a33-a3b03febb56e
if @isdefined cans_in_stock
	md"Actually, I have a hunch there will be another cat coming. Uncomment the code below (remove the #) to add one more can. Remember to run it after making the change!"
else
	md"**Whoopsie!** Because Pluto figures out execution order for you, it doesn't really make sense to assign to the same variable twice. A smarter way to plan ahead is to write `cans_in_stock = consumption` - Pluto will take care of updating everything."
end

# ╔═╡ 6c8e2108-9550-11ea-014d-235770ed4771
md"## Saving cats and notebooks
Alright, we have a neighborhood full of well-fed cats. But oh no, here comes..."

# ╔═╡ bdd5d268-9550-11ea-1046-31efedc36872
if @isdefined scary_dog
	md" $scary_dog is terrorizing the neighborhood! We must do something about it!"
else
	md"You saved the neighborhood! Referencing `scary_dog` leads to an `UndefVarError`, as if it never even existed."
end

# ╔═╡ 36cd006a-9551-11ea-3c0c-df8b7f2843c1
md"To delete a cell like the one defining $scary_dog, click on the ❌ to the right of its code."

# ╔═╡ fb4e471c-9551-11ea-1ab5-41bbd5de76b8
md"Speaking of saving, this notebook is autosaved whenever you change something. The default location for new notebooks is $(tempdir()), a special directory that gets _emptied_ when you restart your computer. To change the save location to something more permanent, scroll to the top - you'll see it next to the Pluto logo."

# ╔═╡ 9d3af596-9554-11ea-21bd-bf427c91c424
md"## ⚡ Pluto power ⚡
Oof, that dog situation was intense. Let's give our cats superpowers to make sure it never happens again!"

# ╔═╡ 3150bf1a-9555-11ea-306f-0fd4d9229a51
md"Remember learning HTML in junior high? Now you can use it for something! Pluto lets you `@bind` variables to HTML elements. As always, every time you change something, Pluto knows what to update!"

# ╔═╡ f2c79746-9554-11ea-39ca-298fd09248ad
@bind power_level html"<input type='range'>"

# ╔═╡ 0b094cea-9556-11ea-268e-0d270fd04d56
md"The power level is $power_level, but we should do more than just say it - let's equip our cats with $power_level emoji!"

# ╔═╡ 1908f9f2-9557-11ea-2abd-dd52f8d776f4
begin
	power_emoji = "⚡"
	power = repeat(power_emoji, power_level)
end

# ╔═╡ 784b1774-9557-11ea-315e-d1ea277ce0fd
if !@isdefined power
	md"Uh oh! Pluto doesn't support multiple expressions per cell. This is a conscious choice - this restriction helps you write less buggy code once you get used to it. To fix the code, you can either split the above cell, or wrap it in a `begin ... end` block. Both work."
else
	md"**Well done!** Your cats have powers now."
end

# ╔═╡ 5edadcd2-9554-11ea-1714-b5b7692c4797
md"We're almost done! It's time to share your amazing story. Press `Ctrl+P` to \"print\" this page to pdf - or you can always share this notebook's save file. (The file is pure Julia, by the way, and it's runnable! You'll learn more about this in the advanced introduction.)"

# ╔═╡ 4634c856-9553-11ea-008d-3539195970ea
md"## Final notes"

# ╔═╡ 653646bc-faea-11ea-301b-3b8ef032f870
@bind nickname html"<input type='text'>"

# ╔═╡ 4d0ebb46-9553-11ea-3431-2d203f594815
md"If anything about this introduction left you confused, something doesn't work, or you have a cool new idea - don't hesitate to contact us! You can do it right from this screen, using the `Instant feedback` form in the bottom right."

# ╔═╡ d736e096-9553-11ea-3ba5-277afde1afe7
md"Also, if you were wondering where `confusing_function` came from, here you go! Remember that you, too, can place code wherever you like."

# ╔═╡ 7366f1b6-954c-11ea-3f68-b12444c902c3
"""
confusing_function(text::String, array::Array)

Repeats the `text` as many times as there are elements in `array`.
"""
confusing_function(text::String, array::Array) = repeat(text, length(array))

# ╔═╡ a4a60262-9547-11ea-3a81-5bf7f9ee5d16
consumption = confusing_function("🥫", neighbors)

# ╔═╡ e11e1660-9549-11ea-22f6-8bb53dc045fe
md"Now we know to prepare $(length(consumption)) cans. Let's stock up!"

# ╔═╡ f27f90c2-954f-11ea-3f93-17acb2ce4280
md"We have $(length(cans_in_stock)) cans of cat food, and our cats need $(length(consumption)). Try adding another cat to the neighborhood to see what changes!"

# ╔═╡ 42446826-faea-11ea-2db4-e12b4de8d601
@bind x html"<input type='range'>"

# ╔═╡ 18d86874-faeb-11ea-0edb-8959e9340df3
f(x) = x^2

# ╔═╡ feba9bd8-faea-11ea-3b03-4194b0343d24
y = f(x)

# ╔═╡ 70fca2b4-faea-11ea-09f9-23d0ee01d59f
repeat(nickname, y)

# ╔═╡ 1deaaf36-9554-11ea-3dae-85851f73dbc6
md"_Have fun using Pluto!_

_~ Fons van der Plas & Nicholas Bochenski_"

# ╔═╡ Cell order:
# ╟─b129ba7c-953a-11ea-3379-17adae34924c
# ╟─4d88b926-9543-11ea-293a-1379b1b5ae64
# ╠═3e8e381e-953f-11ea-3d3e-71d0fea52560
# ╟─aeb3a6bc-9540-11ea-0b8f-6d37412bfe68
# ╟─611c28fa-9542-11ea-1751-fbdedcfb7690
# ╠═7a81b97e-fae7-11ea-2ef2-9768a7d4362f
# ╠═6f7eecec-9543-11ea-1284-dd52fce3ecca
# ╟─a1a20314-9543-11ea-17de-0b658da18992
# ╠═90708aaa-fae7-11ea-31ee-d3d84717608a
# ╠═0508c34e-fae8-11ea-17f8-332503fbd513
# ╠═72c00416-fae9-11ea-24eb-7dee7217d4eb
# ╟─f112b662-9543-11ea-3dcb-2906a99b2188
# ╟─e0642f42-9545-11ea-14ee-fde52cb54ccc
# ╠═3653b1ac-9546-11ea-2a44-ddd3054636fe
# ╟─19ff8d36-9547-11ea-0e08-e5cdd8338673
# ╠═a4a60262-9547-11ea-3a81-5bf7f9ee5d16
# ╟─270ac49e-9549-11ea-3ffd-71ddaee9f134
# ╟─e11e1660-9549-11ea-22f6-8bb53dc045fe
# ╠═745a4584-954a-11ea-028e-59011f268ec6
# ╟─55ad7422-954e-11ea-0a33-a3b03febb56e
# ╠═eac62fea-954e-11ea-2768-39ce6f4059ab
# ╟─f27f90c2-954f-11ea-3f93-17acb2ce4280
# ╟─6c8e2108-9550-11ea-014d-235770ed4771
# ╟─bdd5d268-9550-11ea-1046-31efedc36872
# ╟─36cd006a-9551-11ea-3c0c-df8b7f2843c1
# ╟─fb4e471c-9551-11ea-1ab5-41bbd5de76b8
# ╟─9d3af596-9554-11ea-21bd-bf427c91c424
# ╟─3150bf1a-9555-11ea-306f-0fd4d9229a51
# ╠═f2c79746-9554-11ea-39ca-298fd09248ad
# ╟─0b094cea-9556-11ea-268e-0d270fd04d56
# ╠═1908f9f2-9557-11ea-2abd-dd52f8d776f4
# ╟─784b1774-9557-11ea-315e-d1ea277ce0fd
# ╟─5edadcd2-9554-11ea-1714-b5b7692c4797
# ╟─4634c856-9553-11ea-008d-3539195970ea
# ╠═653646bc-faea-11ea-301b-3b8ef032f870
# ╟─4d0ebb46-9553-11ea-3431-2d203f594815
# ╟─d736e096-9553-11ea-3ba5-277afde1afe7
# ╠═7366f1b6-954c-11ea-3f68-b12444c902c3
# ╠═42446826-faea-11ea-2db4-e12b4de8d601
# ╠═70fca2b4-faea-11ea-09f9-23d0ee01d59f
# ╠═18d86874-faeb-11ea-0edb-8959e9340df3
# ╠═feba9bd8-faea-11ea-3b03-4194b0343d24
# ╟─1deaaf36-9554-11ea-3dae-85851f73dbc6

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
export CuModule

type CuModule
  ptr::CUmodule

  function CuModule(ptr::CUmodule)
    m = new(ptr)
    finalizer(m, cuModuleUnload)
    m
  end
end

function CuModule(filename::String)
  p = CUmodule[0]
  cuModuleLoad(p, filename)
  CuModule(p[1])
end

function CuModule(image::Vector{UInt8})
  p = CUmodule[0]
  cuModuleLoadData(p, image)
  #cuModuleLoadDataEx(p, image, 0, CUjit_option[], Ptr{Void}[])
  CuModule(p[1])
end

Base.unsafe_convert(::Type{CUmodule}, m::CuModule) = m.ptr


text:
"""
    slippage!(slipx, slipy, height, velx, vely, δ, μ)

Fluid substrate interaction that effectively mimics a velocity boundary condition at ``h=0``.

# Arguments

- `slipx :: Array{<:Number,2}`: The x-component of the force due the velocity boundary condition
- `slipy :: Array{<:Number,2}`: The y-component of the force due the velocity boundary condition
- `height::Array{<:Number,2}`: Height field `` h(\\mathbf{x},t)``
- `velx::Array{<:Number,2}`: x-component of the macroscopic velocity vector
- `vely::Array{<:Number,2}`: y-component of the macroscopic velocity vector
- `δ <: Number`: Extrapolation length into the substrate where the **no-slip** is met
- `μ <: Number`: Kinematic viscosity of the simulation, dependent on the value of **τ**

# Mathematics

With the velocity boundary condition at the fluid substrate interface we build the second main descriptor between our model and the thin film equation.
One well studied assumption is that the fluid velocity vanishes at `` h(\\mathbf{x}) = 0 `` which is called **no-slip** condition.
In terms of thin film mobility `M(h)` this means

`` M(h) = \\frac{h^3}{3\\mu}. ``

This assumption, however, can be relaxed to allow for some slippage with a further parameter δ which determines the slip length.
The here presented model matches the thin film equation by modification of the momentum equation of the shallow water model.
Therefore the force we add to the momentum equation to compensate for the substrate fluid friction is given as

`` F_{fric} = -\\nu \\alpha_{\\delta}(h)\\mathbf{u},  ``

where ν is the kinematic viscosity and α

`` \\alpha_{\\delta}(h) = \\frac{6h}{2h^2 + 6h\\delta +3\\delta^2}, ``

and δ being the slip length, or the distance inside the substrate where the fluids velocity vanishes.

# References

- [Zitz, Scagliarini and Harting](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.033313)
- [Münch, Wagner und Witelski](https://link.springer.com/article/10.1007/s10665-005-9020-3)
- [Oron, Davis and Bankoff](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.69.931)
"""
function slippage!(slipx, slipy, height, velx, vely, δ, μ)
    slipx .= (6μ .* height .* velx) ./ (2 .* height.^2 .+ 6δ .* height .+ 3δ^2 )
    slipy .= (6μ .* height .* vely) ./ (2 .* height.^2 .+ 6δ .* height .+ 3δ^2 )
    return nothing
end
# with state struct
function slippage!(state::LBM_state_2D, sys::SysConst)
    state.slipx .= (6*sys.μ .* state.height .* state.velx) ./ (2 .* state.height.^2 .+ 6sys.δ .* state.height .+ 3*sys.δ^2 )
    state.slipy .= (6*sys.μ .* state.height .* state.vely) ./ (2 .* state.height.^2 .+ 6sys.δ .* state.height .+ 3*sys.δ^2 )
    return nothing
end

function slippage!(slip, height, vel, δ, μ)
    slip .= (6μ .* height .* vel) ./ (2 .* height.^2 .+ 6δ .* height .+ 3δ^2 )
    return nothing
end
# with state struct
function slippage!(state::State_1D, sys::SysConst_1D)
    state.slip .= (6*sys.μ .* state.height .* state.vel) ./ (2 .* state.height.^2 .+ 6*sys.δ .* state.height .+ 3*sys.δ^2 )
    return nothing
end

"""
    h∇p!(state)

Computation of the pressure gradient multiplied with the height.

# Mathematics

The pressure gradient (``\\nabla p``) is **the** driving force of the standard thin film equation `` \\partial_t h = (M(h)\\nabla p)``.
Our approach however does not solve the thin film equation directly,
Instead we have to add the pressure gradient as a force which is given as

`` F_{film} = -\\frac{1}{\\rho_0} h \\nabla p_{film}, ``

where the term ``p_{film}`` describes the film pressure

`` p_{film} = -\\gamma [\\Delta h - \\Pi(h)] .``

As such it is the combination of a laplacian term that minimizes the surface area as well as a interfacial potential between substrate and fluid.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> state = Swalbe.Sys(Swalbe.SysConst(Lx=10, Ly=10), "CPU");

julia> state.pressure .= reshape(collect(1:100),10,10);

julia> Swalbe.h∇p!(state)

julia> @test all(state.h∇px[1,:] .== -4.0) # at boundary
Test Passed
  Expression: all(state.h∇px[1, :] .== -4.0)

julia> @test all(state.h∇px[2,:] .== 1.0) # inside
Test Passed
  Expression: all(state.h∇px[2, :] .== 1.0)

```

# References

- [Zitz, Scagliarini and Harting](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.033313)
- [Craster, Matar](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.1131)
- [Oron, Davis and Bankoff](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.69.931)

See also: [`Swalbe.filmpressure!`](@ref)
"""
function h∇p!(state::LBM_state_2D)
    fip, fjp, fim, fjm, fipjp, fimjp, fimjm, fipjm = viewneighbors(state.dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(fip, state.pressure, (1,0))
    circshift!(fjp, state.pressure, (0,1))
    circshift!(fim, state.pressure, (-1,0))
    circshift!(fjm, state.pressure, (0,-1))
    # Diagonal elements
    circshift!(fipjp, state.pressure, (1,1))
    circshift!(fimjp, state.pressure, (-1,1))
    circshift!(fimjm, state.pressure, (-1,-1))
    circshift!(fipjm, state.pressure, (1,-1))
    # In the end it is just a weighted sum...
    state.h∇px .= state.height .* (-1/3 .* (fip .- fim) .- 1/12 .* (fipjp .- fimjp .- fimjm .+ fipjm))
    state.h∇py .= state.height .* (-1/3 .* (fjp .- fjm) .- 1/12 .* (fipjp .+ fimjp .- fimjm .- fipjm))

    return nothing
end
# One dimensional implementation
function h∇p!(state::State_1D)
    fip, fim = viewneighbors_1D(state.dgrad)
    # One dim case, central differences
    circshift!(fip, state.pressure, 1)
    circshift!(fim, state.pressure, -1)

    # In the end it is just a weighted sum...
    state.h∇p .= state.height .* -0.5 .* (fip .- fim)

    return nothing
end

"""
    thermal!(fx, fy, height, kᵦT, μ, δ)

Computations of force due to thermal fluctuations.

# Arguments

- `fx :: Array{<:Number,2}`: x-component of the force due to the fluctuations
- `fy :: Array{<:Number,2}`: y-component of the force due to the fluctuations
- `height :: Array{<:Number,2}`: Height field ``h(\\mathbf{x},t) ``
- `kᵦT <: Number`: Strenght of thermal fluctuations in lattice units
- `μ <: Number`: The kinetic viscosity
- `δ <: Number`: Slip length, needed to normalize

# Mathematics

The classical thin film equation is an equation without thermal noise which is defined as

``\\partial_t h = \\nabla\\cdot(M(h)\\nabla p).``

Classically thermal excitations are neglected for thin film flows.
It is of course possible to introduce for example an envolving surface tension or viscosity but both do not account for e.g. the spectrum of thermal capillary waves.
The reason for this short coming is the complex from the thin film equation takes if derived from the Landau-Lifshitz Navier-Stokes equation.
One addition term due to the stochastic stress tensor makes it somewhat impossible to solve the eqaution self-consistently.
However Grün et al. showed that it is not nessary to use the full stochastic stress tensor, but simple a multiplicative noise term

``\\partial_t h = \\nabla\\cdot[M(h)(\\nabla p + \\sigma \\mathcal{N}], ``

where ``\\sigma`` and ``\\mathcal{N}`` are a dimensionless temperture and a gaussian white noise.
Similar to the pressure gradient the addition of this term is introduced as force in our model

`` F_{fluc} = \\frac{1}{\\rho_0}\\sqrt{2k_BT\\mu\\alpha_{\\delta}(h)}\\mathcal{N},``

with ``\\alpha_{\\delta}(h)`` being the force generated due to substrate slip.

# Examples
```jldoctest; output = false
julia> using Swalbe, Statistics, Test, Random

julia> Random.seed!(1234); # Set a seed

julia> x = ones(50,50); y = ones(50,50); h = ones(50,50);

julia> Swalbe.thermal!(x, y, h, 0.1, 1/6, 1)

julia> @test mean(x) ≈ 0.0 atol=1e-2
Test Passed
  Expression: ≈(mean(x), 0.0, atol = 0.01)
   Evaluated: ≈(-0.004867491968048309, 0.0; atol = 0.01)

julia> @test mean(y) ≈ 0.0 atol=1e-2
Test Passed
  Expression: ≈(mean(y), 0.0, atol = 0.01)
   Evaluated: ≈(-0.000805723303937491, 0.0; atol = 0.01)

julia> @test var(x) ≈ 2*0.1/11 atol=(2*0.1/11)/10 # var = 2kbt*6*μ/slip
Test Passed
  Expression: ≈(var(x), (2 * 0.1) / 11, atol = ((2 * 0.1) / 11) / 10)
   Evaluated: ≈(0.01919664721873892, 0.018181818181818184; atol = 0.0018181818181818184)
```

# References

- [Grün, Mecke and Rauscher](https://link.springer.com/article/10.1007/s10955-006-9028-8)
- [Mecke, Rauscher](https://iopscience.iop.org/article/10.1088/0953-8984/17/45/042/meta)
- [Davidovitch, Moro and Stone](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.95.244505)

See also: [`Swalbe.slippage!`](@ref)
"""
function thermal!(fluc_x, fluc_y, height, kᵦT, μ, δ)
    randn!(fluc_x)
    randn!(fluc_y)
    fluc_x .*= sqrt.(2 .* kᵦT .* μ .* 6 .* height ./
                    (2 .* height.^2 .+
                     6 .* height .* δ .+
                     3 .* δ^2))
    fluc_y .*= sqrt.(2 .* kᵦT .* μ .* 6 .* height ./
                    (2 .* height.^2 .+
                     6 .* height .* δ .+
                     3 .* δ^2))
    return nothing
end

function thermal!(fluc, height, kᵦT, μ, δ)
    randn!(fluc)
    fluc .*= sqrt.(2 .* kᵦT .* μ .* 6 .* height ./
                  (2 .* height.^2 .+
                   6 .* height .* δ .+
                   3 .* δ^2))

    return nothing
end

function thermal!(state::State_thermal, sys::SysConst)
    randn!(state.kbtx)
    randn!(state.kbty)
    state.kbtx .*= sqrt.(2 .* sys.kbt .* sys.μ .* 6 .* state.height ./
                    (2 .* state.height.^2 .+
                     6 .* state.height .* sys.δ .+
                     3 .* sys.δ^2))
    state.kbty .*= sqrt.(2 .* sys.kbt .* sys.μ .* 6 .* state.height ./
                    (2 .* state.height.^2 .+
                     6 .* state.height .* sys.δ .+
                     3 .* sys.δ^2))
    return nothing
end
# With thermal state
function thermal!(state::State_thermal_1D, sys::SysConst_1D)
    randn!(state.kbt)
    state.kbt .*= sqrt.(2 .* sys.kbt .* sys.μ .* 6 .* state.height ./
                  (2 .* state.height.^2 .+
                   6 .* state.height .* sys.δ .+
                   3 .* sys.δ^2))

    return nothing
end

"""
   inclination!(α, state)

Force that mimics the effect of an inclined plate.

Simple model to add a body force on the fluid that should mimic the effect of an inclined plate.
The model includes a smoothing `tanh` function to absorb shocks which might occure.

# Arguments

-`α :: Vector`: Force vector that both hosts direction of the force as well as strength
-`state::State`: Lattice Boltzmann state of the fluid, here we need the `state.Fx`, `state.Fy` fields
-`t::Int`: Time step, used for the smoothing `tanh` factor
-`tstart::Int`: Time delay at which the `tanh` becomes positve
-`tsmooth::Int`: Time interval over which the `tanh` is smeared

# Mathematics

This body force is simply force strength time the mass of the fluid or even simpler the height (assuming ρ=1).
Thus the force becomes

`` \\mathbf{F} = \\mathbf{\\alpha} h \\tanh\\bigg(\\frac{t-t_0}{t_s}\\bigg),``

with ``t_0`` being the time lag at which the `tanh` changes sign and ``t_s`` is width of interval between -1 and 1.

See also: [Swalbe.run_dropletforced](@ref)
"""
function inclination!(α::Vector, state::LBM_state_2D; t=1000, tstart=0, tsmooth=1)
    state.Fx .+= state.height .* α[1] .* (0.5 .+ 0.5 .* tanh((t - tstart)/tsmooth))
    state.Fy .+= state.height .* α[2] .* (0.5 .+ 0.5 .* tanh((t - tstart)/tsmooth))

    return nothing
end

function inclination!(α::Vector, state::State_1D; t=1000, tstart=0, tsmooth=1)
    state.F .+= state.height .* α[1] .* (0.5 .+ 0.5 .* tanh((t - tstart)/tsmooth))

    return nothing
end
"""
    update_rho()

Time evolution of the `active` field rho.

TODO: @Tilman!
"""
function update_rho!(rho, rho_int, height, dgrad, differentials; D=1.0, M=0.0)
    lap_rho, grad_rho, lap_h, grad_h = view_four(differentials)
    ∇²f!(lap_rho, rho, dgrad)
    ∇f!(grad_rho, rho, dgrad)
    ∇²f!(lap_h, height, dgrad)
    ∇f!(grad_h, height, dgrad)
    rho_int .= (D .* lap_rho .-
               M .* (grad_rho.^2 .+ rho .* lap_rho) .-
               D .* (grad_rho .* grad_h ./ height + rho .* (lap_h ./ height .- (grad_h ./ height).^2 )))

    rho .+= rho_int

    return nothing
end

"""
    view_four()

Splits a chuck of memory in four equivalent chucks
"""
function view_four(f)
    f1 = view(f, :, 1)
    f2 = view(f, :, 2)
    f3 = view(f, :, 3)
    f4 = view(f, :, 4)

    return f1, f2, f3, f4
end


text:
using Atum
using Atum.EulerGravity

using StaticArrays: SVector
using WriteVTK

import Atum: boundarystate
function boundarystate(law::EulerGravityLaw, n⃗, x⃗, q⁻, _)
  ρ⁻, ρu⃗⁻, ρe⁻ = EulerGravity.unpackstate(law, q⁻)
  ρ⁺, ρe⁺ = ρ⁻, ρe⁻
  ρu⃗⁺ = ρu⃗⁻ - 2 * (n⃗' * ρu⃗⁻) * n⃗
  SVector(ρ⁺, ρu⃗⁺..., ρe⁺)
end

function risingbubble(law, x⃗, add_perturbation=true)
  FT = eltype(law)
  x, z = x⃗

  Φ = grav(law) * z

  cv_d = FT(719)
  cp_d = γ(law) * cv_d
  R_d = cp_d - cv_d

  θref = FT(300)
  p0 = FT(1e5)
  xc = FT(500)
  zc = FT(350)
  rc = FT(250)
  δθc = FT(1 / 2)

  r = sqrt((x - xc) ^ 2 + (z - zc) ^ 2)
  δθ = r <= rc ? δθc : zero(FT)

  θ = θref
  if add_perturbation
    θ += δθ * (1 + cos(π * r / rc)) / 2
  end
  π_exner = 1 - grav(law) / (cp_d * θ) * z
  ρ = p0 / (R_d * θ) * π_exner ^ (cv_d / R_d)

  ρu = FT(0)
  ρv = FT(0)

  T = θ * π_exner
  ρe = ρ * (cv_d * T + Φ)

  SVector(ρ, ρu, ρv, ρe)
end

function run(A, FT, N, K; outputvtk=true)
  Nq = N + 1

  law = EulerGravityLaw{FT, 2}()

  cell = LobattoCell{FT, A}(Nq, Nq)
  vx = range(FT(0), stop=FT(1e3), length=K+1)
  vz = range(FT(0), stop=FT(1e3), length=K+1)
  grid = brickgrid(cell, (vx, vz); periodic = (true, false))

  dg = DGSEM(; law, cell, grid, numericalflux = RoeFlux())

  cfl = FT(1 // 3)
  dt = cfl * step(vz) / N / 330
  timeend = FT(500)

  q = risingbubble.(Ref(law), points(grid))
  qref = risingbubble.(Ref(law), points(grid), false)

  if outputvtk
    vtkdir = joinpath("output", "euler_gravity", "risingbubble")
    mkpath(vtkdir)
    pvd = paraview_collection(joinpath(vtkdir, "timesteps"))
  end

  do_output = function(step, time, q)
    if outputvtk && step % ceil(Int, timeend / 100 / dt) == 0
      filename = "step$(lpad(step, 6, '0'))"
      vtkfile = vtk_grid(joinpath(vtkdir, filename), grid)
      P = Bennu.toequallyspaced(cell)
      ρ, ρu, ρv, ρe = components(q)
      ρ_ref, ρu_ref, ρv_ref, ρe_ref = components(qref)
      vtkfile["ρ"] = vec(Array(P * (ρ - ρ_ref)))
      vtkfile["ρu"] = vec(Array(P * (ρu - ρu_ref)))
      vtkfile["ρv"] = vec(Array(P * (ρv - ρv_ref)))
      vtkfile["ρe"] = vec(Array(P * (ρe - ρe_ref)))
      vtk_save(vtkfile)
      pvd[time] = vtkfile
    end
  end

  odesolver = LSRK54(dg, q, dt)

  outputvtk && do_output(0, FT(0), q)
  solve!(q, timeend, odesolver; after_step=do_output)
  outputvtk && vtk_save(pvd)
end

let
  A = Array
  FT = Float64
  N = 4
  K = 10
  run(A, FT, N, K)
end


text:
module EditorConfig

###############
## Parameter ##
###############

struct Parameter{T}
    value::T
    validate::Union{Function, Nothing}
    desc::String # Used when calling help <param name>
end

validate(p::Parameter) = p.validate == nothing ? true : p.validate(p.value)

function set(p::Parameter, x)
    old_val = p.value

    # Correct type?
    try
        p.value = x
    catch Exception
        p.value = old_val
        throw(ArgumentError("Invalid parameter assignment: $sym, $x"))
    end

    # Valid?
    if !validate(p)
        p.value = old_val
        throw(ArgumentError("Invalid parameter assignmnt: $sym, $x"))
    end
end

############
## CONFIG ##
############

const CONFIG = Dict{Symbol, Parameter}()

function configSet(sym::Symbol, x)
    # Check if parameter exists
    if !(sym in keys(CONFIG))
        throw(ArgumentError("No parameter named $sym"))
    end

    p = CONFIG[sym]
    set(p, x)

    CONFIG[sym] = p
end

configGet(sym::Symbol) = CONFIG[sym].value
configDesc(sym::Symbol) = CONFIG[sym].desc
configIsParam(sym::Symbol) = sym in keys(CONFIG)

##################
## KEY BINDINGS ##
##################

const KEY_BINDINGS = Dict{UInt32, String}()

"""Remove a keybinding"""
function rmKeyBinding(c::Char)
    delete!(KEY_BINDINGS, UInt32(c) & 0x1f)
end

"""Set a keybinding"""
function setKeyBinding(c::Char, s::String)
    KEY_BINDINGS[UInt32(c) & 0x1f] = s
end

"""Get command from keybinding"""
function getKeyBinding(c::Char) ::String
    get(KEY_BINDINGS, UInt32(c) & 0x1f, "")
end

"""Return true if the given key is bound to a command"""
isKeyBound(c::Char) = (UInt32(c) & 0x1f) in keys(KEY_BINDINGS)

########################
## DEFAULT PARAMETERS ##
########################

CONFIG[:tab_stop] = Parameter{Int}(4, n-> n > 0 && n <= 16, "visual size of a tab in number of spaces")
CONFIG[:expandtab] = Parameter{Bool}(false, nothing, "if true, use spaces instead of tabs when pressing <tab>")
CONFIG[:status_fullpath] = Parameter{Bool}(false, nothing, "show full path to current file")

##########################
## DEFAULT KEY BINDINGS ##
##########################

setKeyBinding('s', "save")
setKeyBinding('o', "open")
setKeyBinding('f', "find")
setKeyBinding('q', "quit")

export
    configGet,
    configSet,
    configIsParam,
    configDesc,
    rmKeyBinding,
    setKeyBinding,
    getKeyBinding,
    isKeyBound

end #module


text:
"""
Implement supremum and infimum functions as extensions of maximum and minimum.
Cover the corner case of empty collections and partially ordered types.
"""

module SupremumInfimum

export supremum, infimum, sup, inf, leaf_type
export is_ordered, is_partially_ordered, is_totally_ordered, is_iterable
export IntSup, IntInf, typesup, typeinf

include("supinf.jl")
include("check.jl")
include("limits.jl")

end # module


text:
using Base.MathConstants

# barebones polynomial type supporting +, -, *
struct SimplePoly{T}
    coeffs::Vector{T}
end

Base.:*(p::SimplePoly, c) = SimplePoly(p.coeffs * c)
Base.:*(c, p::SimplePoly) = p * c

function Base.:+(p::SimplePoly{S}, q::SimplePoly{T}) where {S, T}
    n, m = length(p.coeffs), length(q.coeffs)
    ext1, ext2 = max(0, m - n), max(0, n - m)
    c1 = [p.coeffs ; zeros(S, ext1)]
    c2 = [q.coeffs ; zeros(T, ext2)]
    return SimplePoly(c1 + c2)
end

Base.:-(p::SimplePoly) = -1 * p
Base.:-(p::SimplePoly, q::SimplePoly) = p + (-q)

function Base.:*(p::SimplePoly{S}, q::SimplePoly{T}) where {S, T}
    n, m = length(p.coeffs) - 1, length(q.coeffs) - 1
    c = zeros(promote_type(S, T), m + n + 1)
    for i in 0:n, j in 0:m
        c[i + j + 1] += p.coeffs[i + 1] * q.coeffs[j + 1]
    end
    return SimplePoly(c)
end

function E₁_cfpoly_approx(n::Integer, ::Type{T}=BigInt) where {T<:Real}
    q = SimplePoly(T[1])
    p = x = SimplePoly(T[0,1])
    for i = n:-1:1
        p, q = x*p+(1+i)*q, p # from cf = x + (1+i)/cf = x + (1+i)*q/p
        p, q = p + i*q, p     # from cf = 1 + i/cf = 1 + i*q/p
    end
    # do final 1/(x + inv(cf)) = 1/(x + q/p) = p/(x*p + q)
    return p, x*p + q
end

macro E₁_cf64(x, n::Integer)
    # consider using BigFloat?
    p, q = E₁_cfpoly_approx(n, Float64)
    xesc = esc(x)

    num_expr =  :(@evalpoly $xesc)
    append!(num_expr.args, Float64.(p.coeffs))
    den_expr = :(@evalpoly $xesc)
    append!(den_expr.args, Float64.(q.coeffs))
    :( exp(-$xesc) * $num_expr / $den_expr )
end

function E₁_taylor_coefficients(::Type{T}, n::Integer) where {T<:Number}
    n < 0 && throw(ArgumentError("$n ≥ 0 is required"))
    n == 0 && return T[]
    n == 1 && return T[-eulergamma]
    # iteratively compute the terms in the series, starting with k=1
    term::T = 1
    terms = T[-eulergamma, term]
    for k=2:n
        term = -term * (k-1) / (k * k)
        push!(terms, term)
    end
    return terms
end

# inline the Taylor expansion for a given order n, in double precision
macro E₁_taylor64(z, n::Integer)
    c = E₁_taylor_coefficients(Float64, n)
    zesc = esc(z)
    taylor = :(@evalpoly $zesc)
    append!(taylor.args, c)
    :( $taylor - log($zesc) )
end

# adapted from Steven G Johnson's initial implementation: issue #19
function expint_opt(x::Float64)
    x < 0 && throw(DomainError(x, "negative argument, convert to complex first"))
    x == 0 && return Inf
    if x > 2.15
        # minimax rational approximations to E_1(x)/exp(-x):
        x < 4.0     && return exp(-x) *
            @evalpoly(x, 3.600530862438501481559423277418128014798, 28.73031134165011013783185685393062481126, 46.04314409968065653003548224846863877481, 21.47189493062368074985000918414086604187, 2.719957622921613321844755385973197500235, 1.508750885580864752293599048121678438568e-6) /
            @evalpoly(x, 1.0, 18.06743589038646654075831055159865459831, 61.19456872238615922515377354442679999566, 64.81772518730116701207299231777089576859, 24.19034591054828198408354214931112970741, 2.720026796991567940556850921390829046015)
        x < 10.0    && return exp(-x) *
            @evalpoly(x, 3.149019890512432117647119992448352099575, 14.77395058090815888166486507990452627136, 14.78214309058953358717796744960600201013, 4.559401130686434886620102186841739864936, 0.4027858394909585103775445204576054721422, 2.302781920509468929446800686773538387432e-9) /
            @evalpoly(x, 1.0, 11.65960373479520446458792926669115987821, 26.20023773503894444535165299118172674849, 18.93899893550582921168134979000319186841, 4.962178168140565906794561070524079194193, 0.4027860481050182948737116109665534358805)
        x < 20.0    && return exp(-x) *
            @evalpoly(x, 2.564801308922428705318296668924369454617, 5.482252510134574167659359513298970499768, 2.379528224853089764405551768869103662657, 0.2523431276282591480166881146593592650031, 1.444719769329975045925053905197199934930e-9, -8.977332061106791470409502623202677045468e-12) /
            @evalpoly(x, 1.0, 6.421214405318108272004472721910023284626, 7.609584052290707052877352911548283916247, 2.631866613851763364839413823973711355890, 0.2523432345539146248733539983749722854603)
        x < 200.0 && return @E₁_cf64(x, 8) # fallback continued fraction
        return x < 740.0 ? @E₁_cf64(x, 4) : 0.0 # underflow
    else
        # crossover point to taylor should be tuned more
        return x ≤ 0.6 ? (x ≤ 0.053 ? (x ≤ 4.4e-3 ? @E₁_taylor64(x,4) :
                                                       @E₁_taylor64(x,8)) :
                                       @E₁_taylor64(x,15)) :
                          @E₁_taylor64(x,37)
    end
end

function expint_opt(z::Complex{Float64})
    x² = real(z)^2
    y² = imag(z)^2
    if x² + 0.233*y² ≥ 7.84 # use cf expansion, ≤ 30 terms
        if (x² ≥ 546121) & (real(z) > 0) # underflow
            return zero(z)
        elseif x² + 0.401*y² ≥ 58.0 # ≤ 15 terms
            if x² + 0.649*y² ≥ 540.0 # ≤ 8 terms
                x² + y² ≥ 4e4 && return @E₁_cf64(z, 4)
                return @E₁_cf64(z, 8)
            end
            return @E₁_cf64(z, 15)
        end
        return @E₁_cf64(z, 30)
    else # use Taylor expansion, ≤ 37 terms
        r² = x² + y²
        return r² ≤ 0.36 ? (r² ≤ 2.8e-3 ? (r² ≤ 2e-7 ? @E₁_taylor64(z,4) :
                                                       @E₁_taylor64(z,8)) :
                                         @E₁_taylor64(z,15)) :
                          @E₁_taylor64(z,37)
    end
end

function expint(z::Complex{Float64})
    if real(z) < 0
        return expint(1, z)
    else
        return expint_opt(z)
    end
end
expint(x::Float64) = expint_opt(x)

function expint(x::BigFloat)
    iszero(x) && return Inf
    signbit(x) && throw(DomainError(x, "negative argument, convert to complex first"))
    return -expinti(-x)
end

expint(z::Union{T,Complex{T},Rational{T},Complex{Rational{T}}}) where {T<:Integer} = expint(float(z))
expint(x::Number) = expint(1, x)
expint(z::Float32) = Float32(expint(Float64(z)))
expint(z::ComplexF32) = ComplexF32(expint(ComplexF64(z)))
expint(z::Float16) = Float16(expint(Float64(z)))
expint(z::ComplexF16) = ComplexF16(expint(ComplexF64(z)))

# Continued fraction for En(ν, z) that doesn't use a term with
# the gamma function: https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/10/0001/
function En_cf_nogamma(ν::Number, z::Number, n::Int=1000)
    B = float(z + ν)
    Bprev::typeof(B) = z
    A::typeof(B) = 1
    Aprev::typeof(B) = 1
    ϵ = 10*eps(real(B))
    scale = sqrt(floatmax(real(A)))

    # two recurrence steps / loop
    iters = 0
    for i = 2:n
        iters += 1

        A′ = A
        A = z*A + (i-1) * Aprev
        Aprev = A′
        B′ = B
        B = z*B + (i-1) * Bprev
        Bprev = B′

        A′ = A
        A = A + (ν+i-1) * Aprev
        Aprev = A′
        B′ = B
        B = B + (ν+i-1) * Bprev
        Bprev = B′

        conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
        conv && i > 4 && break

        # rescale
        if fastabs(A) > scale
            A     /= scale
            Aprev /= scale
            B     /= scale
            Bprev /= scale
        end
    end

    exppart = exp(-z)
    if isinf(exppart)
        return exp(-z + log(A) - log(B)), iters
    else
        cfpart = A/B
        return cfpart * exppart, iters
    end
end

# Calculate Γ(1 - ν) * z^(ν-1) safely
En_safe_gamma_term(ν::Number, z::Number) = exp((ν - 1)*log(z) + loggamma(1 - oftype(z, ν)))

# continued fraction for En(ν, z) that uses the gamma function:
# https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/10/0005/
# returns the two terms from the above equation separately
function En_cf_gamma(ν::Number, z::Number, n::Int=1000)
    A = float(1 - ν)
    B::typeof(A) = 1
    Bprev::typeof(A) = 0
    Aprev::typeof(A) = 1
    ϵ = 10*eps(real(B))
    scale = sqrt(floatmax(real(A)))

    iters = 0
    j = 1
    for i = 2:n
        iters += 1

        A′ = A
        term = iseven(i) ? -(i÷2 - 1 - ν)*z : z*(i÷2)
        A = (i - ν)*A + term * Aprev
        Aprev = A′
        B′ = B
        B = (i - ν)*B + term * Bprev
        Bprev = B′

        conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
        conv && break

        if max(abs(real(A)), abs(imag(A))) > scale
            A     /= scale
            Aprev /= scale
            B     /= scale
            Bprev /= scale
        end
    end

    gammapart = En_safe_gamma_term(ν, z)
    cfpart = exp(-z)
    if isinf(cfpart)
        factor = sign(real(cfpart)) + sign(imag(cfpart))*im
        cfpart = exp(-z + log(B) - log(A))
    else
        cfpart *= B/A
    end
    return gammapart, -cfpart, iters
end

# picks between continued fraction representations in
# En_cf_nogamma and En_cf_gamma
# returns (evaluated result, # iterations used, whether En_cf_gamma was chosen)
function En_cf(ν::Number, z::Number, niter::Int=1000)
    gammapart, cfpart, iters = En_cf_gamma(ν, z, niter)
    gammaabs, cfabs = abs(gammapart), abs(cfpart)
    if gammaabs != Inf && gammaabs > 1.0 && gammaabs > cfabs
        # significant gamma part, use this
        return gammapart + cfpart, iters, true
    else
        return En_cf_nogamma(ν, z, niter)..., false
    end
end

# Compute expint(ν, z₀+Δ) given start = expint(ν, z₀), as described by [Amos 1980].
# This is used to incrementally approach the negative real axis.
function En_taylor(ν::Number, start::Number, z₀::Number, Δ::Number)
    a = exp(z₀) * start
    k, iters = 0, 0
    asum = a
    Δ_prod_fact = -Δ
    ϵ = 10*eps(real(asum))

    for k = 0:100
        a_pre = Δ_prod_fact + a*Δ*(ν - k - 1)/(k + 1)
        a = a_pre / z₀
        asum_prev = asum
        asum += a

        if abs(asum_prev - asum) < ϵ
            break
        end

        Δ_prod_fact *= -Δ / (k + 2)
    end

    res = exp(-z₀) * asum
    return res
end

# series about origin, general ν
# https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/06/01/04/01/01/0003/
function En_expand_origin_general(ν::Number, z::Number)
    gammaterm = En_safe_gamma_term(ν, z)
    frac = one(real(z))
    sumterm = frac / (1 - ν)
    k, maxiter = 1, 100
    ϵ = 10*eps(real(sumterm))
    while k < maxiter
        frac *= -z / k
        prev = sumterm
        sumterm += frac / (k + 1 - ν)
        if abs(sumterm - prev) < ϵ
            break
        end
        k += 1
    end
    return gammaterm - sumterm
end

# series about the origin, special case for integer n > 0
# https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/06/01/04/01/02/0005/
function En_expand_origin_posint(n::Integer, z::Number)
    gammaterm = 1 # (-z)^(n-1) / (n-1)!
    for i = 1:n-1
        gammaterm *= -z / i
    end

    frac = one(real(z))
    gammaterm *= digamma(oftype(frac,n)) - log(z)
    sumterm = n == 1 ? zero(frac) : frac / (1 - n)
    k, maxiter = 1, 100
    ϵ = 10*eps(real(sumterm))
    while k < maxiter
        frac *= -z / k
        # skip term with zero denominator
        if k != n-1
            prev = sumterm
            sumterm += frac / (k + 1 - n)
            if abs(sumterm - prev) < ϵ
                break
            end
        end
        k += 1
    end
    return gammaterm - sumterm
end

function En_expand_origin(ν::Number, z::Number)
    if isinteger(ν) && real(ν) > 0
        return En_expand_origin_posint(Integer(ν), z)
    else
        return En_expand_origin_general(ν, z)
    end
end

# can find imaginary part of E_ν(x) for x on negative real axis analytically
# https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/04/05/01/0003/
function En_imagbranchcut(ν::Number, z::Number)
    a = real(z)
    e1 = exp(π*imag(ν))
    e2 = Complex(cospi(real(ν)), -sinpi(real(ν)))
    impart = π * im * e1 * e2 * exp((ν-1)*log(complex(a)) - loggamma(ν))
    impart *= signbit(imag(z)) ? -1 : 1
    return imag(impart) * im # get rid of any real error
end

"""
    expint(z)
    expint(ν, z)

Computes the exponential integral ``\\operatorname{E}_\\nu(z) = \\int_0^\\infty \\frac{e^{-zt}}{t^\\nu} dt``.
If ``\\nu`` is not specified, ``\\nu=1`` is used. Arbitrary complex ``\\nu`` and ``z`` are supported.

External links: [DLMF](https://dlmf.nist.gov/8.19), [Wikipedia](https://en.wikipedia.org/wiki/Exponential_integral)
"""
function expint(ν::Number, z::Number, niter::Int=1000)
    if abs(ν) > 50 && !(isreal(ν) && real(ν) > 0)
        throw(ArgumentError("Unsupported order |ν| > 50 off the positive real axis"))
    end

    z, = promote(float(z), ν)
    if isnan(ν) || isnan(z)
        return oftype(z, NaN) * z
    end

    if z isa Real && (isinteger(ν) ? (z < 0 && ν > 0) : (z < 0 || ν < 0))
        throw(DomainError(z, "En will only return a complex result if called with a complex argument"))
    end

    if z == 0
        return oftype(z, real(ν) > 0 ? 1/(ν-1) : Inf)
    end

    if ν == 0
        return exp(-z) / z
    elseif ν == 1 && real(z) > 0 && z isa Union{Float64, Complex{Float64}}
        return expint_opt(z)
    end
    # asymptotic test for underflow when Re z → ∞
    # https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/06/02/0003/
    if real(z) > -log(nextfloat(zero(real(z))))+1 # exp(-z) is zero
        return zero(z)
    end

    if abs2(z) < 9
        # use Taylor series about the origin for small z
        return En_expand_origin(ν, z)
    end

    if z isa Real
        return first(z > 0 ? En_cf(ν, z, niter) : En_cf_nogamma(ν, z, niter))
    else
        # Procedure for z near the negative real axis based on
        # (without looking at the accompanying source code):
        #   Amos, D. E. (1990). Computation of exponential integrals
        #   of a complex argument. ACM Transactions on Mathematical Software,
        #   16(2), 169–177. https://doi.org/10.1145/78928.78933
        doconj = imag(z) < 0
        rez, imz = real(z), abs(imag(z))
        z = doconj ? conj(z) : z
        ν = doconj ? conj(ν) : ν

        quick_niter, nmax = 50, 45
        # start with small imaginary part if exactly on negative real axis
        imstart = (imz == 0) ? abs(z)*sqrt(eps(typeof(real(z)))) : imz
        z₀ = rez + imstart*im
        E_start, i, _ = En_cf(ν, z₀, quick_niter)
        isnan(E_start) && return E_start
        if imz > 0 && i < nmax
            # didn't need to take any steps
            return doconj ? conj(E_start) : E_start
        end
        while i > nmax
            # double imaginary part until in region with fast convergence
            imstart *= 2
            z₀ = rez + imstart*im
            E_start, i, _ = En_cf(ν, z₀, quick_niter)
        end

        # nsteps chosen so |Δ| ≤ 0.5
        nsteps = ceil(2 * (imstart - imz))
        Δ = (imz - imstart)*im / nsteps

        for j = 1:nsteps
            # take Δ sized steps towards the desired z
            E_start = En_taylor(ν, E_start, z₀, Δ)
            z₀ += Δ
        end

        # more exact imaginary part available for non-integer ν
        if imz == 0
            E_start = real(E_start) + En_imagbranchcut(ν, z)
        end

        return doconj ? conj(E_start) : E_start
    end
    throw("unreachable")
end

##############################################################################
# expinti function Ei

"""
    expinti(x::Real)

Computes the exponential integral function ``\\operatorname{Ei}(x) = \\int_{-\\infty}^x \\frac{e^t}{t} dt``,
which is equivalent to ``-\\Re[\\operatorname{E}_1(-x)]`` where ``\\operatorname{E}_1`` is the `expint` function.
"""
expinti(x::Real) = x ≤ 0 ? -expint(-x) : -real(expint(complex(-x)))

# inline the Taylor expansion of Ei for a given order n, in double precision
macro Ei_taylor64(z, n::Integer)
    c = .- E₁_taylor_coefficients(Float64, n) .* (-1).^(0:n)
    zesc = esc(z)
    taylor = :(@evalpoly $zesc)
    append!(taylor.args, c)
    :( $taylor + log($zesc) )
end

# optimized double-precision version
function expinti(x::Float64)
    x ≤ 0 && return -expint(-x)

    if x > 2.15
        # minimax rational approximants:
        x < 4 && return @evalpoly(x,-2.43791466332154621,3.09402100064798205,9.35202477109609955,0.152659977028953397,0.0157273683896079142,0.0345566671015011426,-0.000421531433157416203) /
                        @evalpoly(x,1.0,4.28055563991564399,0.537599625698465573,-0.511064414527643313,0.0867748262262250088,-0.00623913330836521800,0.000172066498182538260)
        x < 10 && return exp(x) *
            @evalpoly(x,-1.58447018083420958,4.71806833998906997,-0.587691572500210206,0.125012472861504555,-0.00178055441724967428,0.000633648975971195928,0.0000147213934578379204,2.12391754244415544e-6) /
            @evalpoly(x,1.0,1.93297600031287800,0.660790440069106542,0.198322636197663277,0.0272447293513279631,0.00399501571688512611,0.000362510989191199243,0.0000182930089855534336,2.06800780072894204e-6)
        x < 20 && return exp(x) *
            @evalpoly(x,-1.77183291754640123,0.795659966861260409,-0.221223333413388642,0.0328877243243796815,-0.00331846947191676458,0.000180945604349930285,-5.97641401680304362e-6,2.42151808626299747e-11) /
            @evalpoly(x,1.0,-2.10926998628216150,0.933357955421497965,-0.245433884954174394,0.0356954809772243699,-0.00348034743685382360,0.000186615220328647350,-5.97232581353392099e-6)

        x > 710 && return Inf # overflow
        x⁻¹ = inv(x)

        # minimax rational approximant in x⁻¹ for x ∈ [20,200]
        x < 200 && return exp(x) *
            @evalpoly(x⁻¹, -5.29842699621003563e-14, +1.00000000004732488, -60.4361334939888359, +1327.83891720487710, -6810.63668974273961, -177755.383525765400,+3.00773484037048848e6, -1.53642380695372707e7, +2.08174653368702692e7) /
            @evalpoly(x⁻¹, 1.0, -61.4361334756161381, +1387.27504658395142, -8081.03888544858393, -172104.333927401741, +3.18903576285551101e6, -1.81873890267574206e7, +3.37312131843327704e7, -1.22198734384213631e7)

        # asymptotic series eˣ/x ∑k!/xᵏ for x ≥ 200:
        return exp(x)*x⁻¹ * @evalpoly(x⁻¹, 1,1,2,6,24,120,720,5040)
    end

    root = 0.37250741078136663446 # OEIS A091723
    dx = x - root
    if abs(dx) < 0.03
        # taylor series around real root of Ei
        return dx * @evalpoly(dx, 3.896215733907167310, -3.281607866398561671, 6.52237614543892570, -12.96969738353651704, 27.88629796294204998, -62.3788015289154187, 143.5349488096750988, -337.155827178746892, 804.531839982138251, -1943.79664572349884, 4743.76565040243084, -11673.46399116716364, 28926.9553054354509)
    else
        # TODO: crossover points and degrees could be tuned more
        return x ≤ 0.6 ? (x ≤ 0.053 ? (x ≤ 4.4e-3 ? @Ei_taylor64(x,4) :
                                                    @Ei_taylor64(x,8)) :
                                       @Ei_taylor64(x,15)) :
                          @Ei_taylor64(x,37)
    end
end

expinti(z::Union{T,Rational{T}}) where {T<:Integer} = expinti(float(z))
expinti(z::Float32) = Float32(expinti(Float64(z)))
expinti(z::Float16) = Float16(expinti(Float64(z)))

function expinti(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_eint, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, Base.MPFR.ROUNDING_MODE[])
    return z
end


text:
include("clusters.jl")

sim_param = setup_sim_param_model()
add_param_fixed(sim_param,"num_targets_sim_pass_one",150060)
add_param_fixed(sim_param,"max_incl_sys",80.0) #degrees; 0 (deg) for isotropic system inclinations; set closer to 90 (deg) for more transiting systems

##### To load a file with all the best active parameters from a set of optimization runs:

files_directory = "ACI/Model_Optimization/Clustered_P_R/All_params_random_weighted_targs150060_maxincl80/"
#files_directory = "ACI/Model_Optimization/Non_clustered/Some9_params1_random_weighted_targs150060_maxincl80/"
file_name = "Active_params_best_weighted_all.txt"

active_params_best_all = CSV.read(files_directory*file_name, delim=" ")
active_params_names = names(active_params_best_all)[2:end]

##### To simulate a catalog for each set of best active parameters:

for num_cat in 1:size(active_params_best_all)[1] #1:size(active_params_best_all)[1]
    println("Generating simulated catalog ", num_cat)
    tic()

    #To set up the model parameters:
    for (i,param_name) in enumerate(active_params_names)
        add_param_active(sim_param, string(param_name), active_params_best_all[num_cat, param_name])
    end

    save_name_end = string(active_params_best_all[num_cat, :run_number])*".out"

    #To generate the underlying systems:
    cat_phys = generate_kepler_physical_catalog(sim_param)

    #To save the underlying/true planets/systems:
    f = open(files_directory*"periods_all"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].orbit) > 0
            periods_sys = Array{Float64}(length(targ.sys[1].orbit))
            for (j,planet) in enumerate(targ.sys[1].orbit)
                periods_sys[j] = planet.P #days
            end
            println(f, periods_sys)
        end
    end
    close(f)

    f = open(files_directory*"eccentricities_all"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].orbit) > 0
            ecc_sys = Array{Float64}(length(targ.sys[1].orbit))
            for (j,planet) in enumerate(targ.sys[1].orbit)
                ecc_sys[j] = planet.ecc
            end
            println(f, ecc_sys)
        end
    end
    close(f)

    f = open(files_directory*"radii_all"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].planet) > 0
            radii_sys = Array{Float64}(length(targ.sys[1].planet))
            for (j,planet) in enumerate(targ.sys[1].planet)
                radii_sys[j] = planet.radius #solar radii
            end
            println(f, radii_sys)
        end
    end
    close(f)

    f = open(files_directory*"masses_all"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].planet) > 0
            masses_sys = Array{Float64}(length(targ.sys[1].planet))
            for (j,planet) in enumerate(targ.sys[1].planet)
                masses_sys[j] = planet.mass #solar masses
            end
            println(f, masses_sys)
        end
    end
    close(f)

    #To save the stellar properties of all the systems:
    f = open(files_directory*"stellar_masses_with_planets"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].planet) > 0
            star_mass = targ.sys[1].star.mass #solar masses
            println(f, star_mass)
        end
    end
    close(f)

    f = open(files_directory*"stellar_radii_with_planets"*save_name_end, "w")
    write_model_params(f, sim_param)
    for (i,targ) in enumerate(cat_phys.target)
        if length(targ.sys[1].planet) > 0
            star_radius = targ.sys[1].star.radius #solar radii
            println(f, star_radius)
        end
    end
    close(f)

    #To generate the observed systems:
    cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param)
    cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param)
    summary_stat = calc_summary_stats_model(cat_obs,sim_param)

    #To save the observed planets/systems:
    f = open(files_directory*"periods"*save_name_end, "w")
    write_model_params(f, sim_param)
    for num_pl_in_sys in 1:length(summary_stat.cache["idx_n_tranets"])
        num_targets = length(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
        period_array = Array{Float64}(num_pl_in_sys,num_targets)
        for (i,j) in enumerate(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
            period_array[:,i] = map(ExoplanetsSysSim.period, cat_obs.target[j].obs)[1:num_pl_in_sys]
        end
        println(f,"# Periods of systems with ", num_pl_in_sys, " detected planets.")
        if length(period_array) > 0
            println(f,period_array')
        end
    end
    close(f)

    f = open(files_directory*"depths"*save_name_end, "w")
    write_model_params(f, sim_param)
    for num_pl_in_sys in 1:length(summary_stat.cache["idx_n_tranets"])
        num_targets = length(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
        depths_array = Array{Float64}(num_pl_in_sys,num_targets)
        for (i,j) in enumerate(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
            depths_array[:,i] = map(ExoplanetsSysSim.depth, cat_obs.target[j].obs)[1:num_pl_in_sys]
        end
        println(f,"# Transit depths of systems with ", num_pl_in_sys, " detected planets.")
        if length(depths_array) > 0
            println(f,depths_array')
        end
    end
    close(f)

    f = open(files_directory*"durations"*save_name_end, "w")
    write_model_params(f, sim_param)
    for num_pl_in_sys in 1:length(summary_stat.cache["idx_n_tranets"])
        num_targets = length(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
        duration_array = Array{Float64}(num_pl_in_sys,num_targets)
        for (i,j) in enumerate(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
            duration_array[:,i] = map(ExoplanetsSysSim.duration, cat_obs.target[j].obs)[1:num_pl_in_sys]
        end
        println(f,"# Transit durations of systems with ", num_pl_in_sys, " detected planets.")
        if length(duration_array) > 0
            println(f,duration_array')
        end
    end
    close(f)

    #To save the stellar properties of the systems with observed planets:
    f = open(files_directory*"stellar_masses_obs"*save_name_end, "w")
    write_model_params(f, sim_param)
    for num_pl_in_sys in 1:length(summary_stat.cache["idx_n_tranets"])
        num_targets = length(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
        stellar_mass_array = Array{Float64}(num_targets)
        for (i,j) in enumerate(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
            stellar_mass_array[i] = cat_obs.target[j].star.mass
        end
        println(f,"# Stellar masses of systems with ", num_pl_in_sys, " detected planets.")
        if length(stellar_mass_array) > 0
            println(f, stellar_mass_array)
        end
    end
    close(f)

    f = open(files_directory*"stellar_radii_obs"*save_name_end, "w")
    write_model_params(f, sim_param)
    for num_pl_in_sys in 1:length(summary_stat.cache["idx_n_tranets"])
        num_targets = length(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
        stellar_radius_array = Array{Float64}(num_targets)
        for (i,j) in enumerate(summary_stat.cache["idx_n_tranets"][num_pl_in_sys])
            stellar_radius_array[i] = cat_obs.target[j].star.radius
        end
        println(f,"# Stellar radii of systems with ", num_pl_in_sys, " detected planets.")
        if length(stellar_radius_array) > 0
            println(f, stellar_radius_array)
        end
    end
    close(f)

    toc()
end


text:
"""
    omnivory(N::T) where {T <: UnipartiteNetwork}

This function measure an omnivory value for each species.

    length(omnivory(N))

This command gives you the number of omnivories species in the network.

"""

function omnivory(N::T) where {T <: UnipartiteNetwork}
    OI = Dict([s => 0.0 for s in species(N)])

    TL = fractional_trophic_level(N)
    k = degree(N; dims=1)

    for sp_i in species(N)

        # Species with no interaction have an omnivory index of 0
        k[sp_i] > 0 || continue

        # For every species, we set its initial omnivory to 0
        oi = 0.0
        for (j, sp_j) in enumerate(species(N))
            # Then for every species it consumes, we ha
            tl_diff = (TL[sp_j] - (TL[sp_i]-1.0)).^2.0
            corr = N[sp_i,sp_j]/k[sp_i]
            oi += tl_diff * corr
        end
        OI[sp_i] = oi
    end

    return OI
end