text:
export AbstractPathTracker,
    Tracker,
    TrackerResult,
    TrackerOptions,
    TrackerParameters,
    TrackerCode,
    DEFAULT_TRACKER_PARAMETERS,
    FAST_TRACKER_PARAMETERS,
    CONSERVATIVE_TRACKER_PARAMETERS,
    track,
    init!,
    track!,
    step!,
    start_parameters!,
    target_parameters!,
    status,
    state,
    is_success,
    is_terminated,
    is_invalid_startvalue,
    is_tracking,
    solution,
    steps,
    accepted_steps,
    rejected_steps,
    iterator

###
### Options and Parameters
###

"""
    TrackerParameters

Parameters that control the performance and robustness characteristics of the path tracking
algorithm. See [^Tim20] for an explanation and derivation of the parameters.
We provide three sets of parameters for common use cases:

* [`DEFAULT_TRACKER_PARAMETERS`](@ref)
* [`FAST_TRACKER_PARAMETERS`](@ref)
* [`CONSERVATIVE_TRACKER_PARAMETERS`](@ref)

[^Tim20]: Timme, S. "Mixed Precision Path Tracking for Polynomial Homotopy Continuation". arXiv:1902.02968 (2020)
"""
Base.@kwdef mutable struct TrackerParameters
    a::Float64 = 0.125
    β_a::Float64 = 1.0
    # compared to β_ω in the article, β_ω_p goes into the formula by β_ω_p^p, i.e.
    # the relative reduction in step size is independent of the order p
    β_ω_p::Float64 = 3.0
    β_τ::Float64 = 0.4
    strict_β_τ::Float64 = min(0.75β_τ, 0.4)
    min_newton_iters::Int = 2
end
Base.show(io::IO, TP::TrackerParameters) = print_fieldnames(io, TP)

"The default [`TrackerParameters`](@ref) which have a good balance between robustness and efficiency."
const DEFAULT_TRACKER_PARAMETERS = TrackerParameters()
"[`TrackerParameters`](@ref) which trade speed against a higher chance of path jumping."
const FAST_TRACKER_PARAMETERS = TrackerParameters(β_τ = 0.75, β_ω_p = 2.0)
"[`TrackerParameters`](@ref) which trade robustness against some speed."
const CONSERVATIVE_TRACKER_PARAMETERS = TrackerParameters(β_ω_p = 4.0, β_τ = 0.25)

"""
    TrackerOptions(; options...)

The set of options for a [`Tracker`](@ref).

## Options

* `automatic_differentiation = 1`: The value `automatic_differentiation` determines
  up to which order the derivative is computed using automatic differentiation.
  Otherwise numerical differentiation is used. The automatic differentiation results
  in additional compilation time, however for numerically challenging paths it is strongly
  recommended to use `automatic_differentiation = 3`.
* `max_steps = 10_000`: The maximal number of steps a tracker attempts
* `max_step_size = Inf`: The maximal size of a step
* `max_initial_step_size = Inf`: The maximal size of the first step
* `min_step_size = 1e-48`: The minimal step size. If a smaller step size would
  be necessary, then the tracking gets terminated.
* `extended_precision = true`: Whether to allow for the use of extended precision,
  if necessary, in some computations. This can greatly improve the ability to track
  numerically difficult paths.
* `terminate_cond = 1e13`: If the relative component-wise condition number
  `cond(H_x, ẋ)` is larger than `terminate_cond` then the path is terminated as too
  ill-conditioned.
* `parameters::Union{Symbol,TrackerParameters} = :default` Set the
  [`TrackerParameters`](@ref) to control the performance of the path tracking algorithm.
  The values `:default`, `:conservative` and `:fast` are shorthands for using
  [`DEFAULT_TRACKER_PARAMETERS`](@ref), [`CONSERVATIVE_TRACKER_PARAMETERS`](@ref) resp.
  [`FAST_TRACKER_PARAMETERS`](@ref).
"""
mutable struct TrackerOptions
    automatic_differentiation::Int
    max_steps::Int
    max_step_size::Float64
    max_initial_step_size::Float64
    extended_precision::Bool
    min_step_size::Float64
    min_rel_step_size::Float64
    terminate_cond::Float64
    parameters::TrackerParameters
end
function TrackerOptions(;
    automatic_differentiation::Int = 1,
    max_steps::Int = 10_000,
    max_step_size::Float64 = Inf,
    max_initial_step_size::Float64 = Inf,
    extended_precision::Bool = true,
    min_step_size::Float64 = 1e-48,
    min_rel_step_size::Float64 = 0.0,
    terminate_cond::Float64 = 1e13,
    parameters::Union{Symbol,TrackerParameters} = :default,
)

    if parameters isa Symbol
        if parameters == :default
            parameters = DEFAULT_TRACKER_PARAMETERS
        elseif parameters == :fast
            parameters = FAST_TRACKER_PARAMETERS
        elseif parameters == :conservative
            parameters = CONSERVATIVE_TRACKER_PARAMETERS
        else
            throw(ArgumentError("Unsupported `parameters` value $parameters"))
        end
    end

    TrackerOptions(
        automatic_differentiation,
        max_steps,
        max_step_size,
        max_initial_step_size,
        extended_precision,
        min_step_size,
        min_rel_step_size,
        terminate_cond,
        parameters,
    )
end

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

###
### TrackerCode
###

@doc """
    TrackerCode

The possible states a `CoreTracker` can have are of type `TrackerCode.codes` and can be

* `TrackerCode.success`: Indicates a successfull tracking.
* `TrackerCode.tracking`: The tracking is still in progress.
* `TrackerCode.terminated_accuracy_limit`: Tracking terminaed since the accuracy was insufficient.
* `TrackerCode.terminated_invalid_startvalue`: Tracking terminated since the provided start value was invalid.
* `TrackerCode.terminated_ill_conditioned`: Tracking terminated since the path was too ill-conditioned.
* `TrackerCode.terminated_max_steps`: Tracking terminated since maximal number of steps is reached.
* `TrackerCode.terminated_step_size_too_small`: Trackint terminated since the step size was too small.
* `TrackerCode.terminated_unknown`: An unintended error occured. Please consider reporting an issue.
"""
module TrackerCode

@enum codes begin
    tracking
    success
    terminated_max_steps
    terminated_accuracy_limit
    terminated_ill_conditioned
    terminated_invalid_startvalue
    terminated_invalid_startvalue_singular_jacobian
    terminated_step_size_too_small
    terminated_unknown
end

end

"""
    is_success(code::TrackerCode.codes)

Returns `true` if `code` indicates a success in the path tracking.
"""
is_success(S::TrackerCode.codes) = S == TrackerCode.success

"""
    is_terminated(code::TrackerCode.codes)

Returns `true` if `code` indicates that the path tracking got terminated.
"""
is_terminated(S::TrackerCode.codes) = S ≠ TrackerCode.tracking && S ≠ TrackerCode.success

"""
    is_invalid_startvalue(code::TrackerCode.codes)

Returns `true` if `code` indicates that the path tracking got terminated since the start
value was not a regular zero.
"""
function is_invalid_startvalue(S::TrackerCode.codes)
    S == TrackerCode.terminated_invalid_startvalue ||
        S == TrackerCode.terminated_invalid_startvalue_singular_jacobian
end

"""
    is_tracking(code::TrackerCode.codes)

Returns `true` if `code` indicates that the path tracking is not yet finished.
"""
is_tracking(S::TrackerCode.codes) = S == TrackerCode.tracking

###
### RESULT
###

"""
    TrackerResult

Containing the result of tracking a path with a [`Tracker`](@ref).

## Fields

* `return_code::Symbol`: A code indicating whether the tracking was successfull (`:success`).
  See [`TrackerCode`](@ref) for all possible values.
* `solution::V`: The solution when the tracking stopped.
* `t::ComplexF64`: The value of `t` when the tracking stopped.
* `accuracy::Float64`: Estimate of the relative accuracy of the `solution`.
* `accepted_steps::Int`: Number of steps that got accepted.
* `rejected_steps::Int`: Number of steps that got rejected.
* `extended_precision::Bool`: Indicate whether extended precision is necessary to achieve
  the accuracy of the `solution`.
* `extended_precision_used::Bool`: This is `true` if during the tracking at any point
  extended precision was used.
"""
struct TrackerResult
    return_code::Symbol
    solution::Vector{ComplexF64}
    t::ComplexF64
    accuracy::Float64
    accepted_steps::Int
    rejected_steps::Int
    extended_precision::Bool
    extended_precision_used::Bool
    ω::Float64
    μ::Float64
    τ::Float64
end

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

"""
    is_success(result::TrackerResult)

Returns `true` if the path tracking was successfull.
"""
is_success(R::TrackerResult) = R.return_code == :success

"""
    is_invalid_startvalue(result::TrackerResult)

Returns `true` if the path tracking failed since the start value was invalid.
You can inspect `result.return_code` to get the exact return code. Possible values
if `is_invalid_startvalue(result) == true` are
* `:terminated_invalid_startvalue_singular_jacobian` indicates that the Jacobian of the homotopy at
  the provided start value is singular, i.e., if it has not full-column rank.
* `:terminated_invalid_startvalue` indicates that the the provided start value is not sufficiently
  close to a solution of the homotopy.
"""
function is_invalid_startvalue(R::TrackerResult)
    R.return_code == :terminated_invalid_startvalue ||
        R.return_code == :terminated_invalid_startvalue_singular_jacobian
end

"""
    solution(result::TrackerResult)

Returns the solutions obtained by the `Tracker`.
"""
solution(result::TrackerResult) = result.solution

"""
    steps(result::TrackerResult)

Returns the number of steps done.
"""
steps(result::TrackerResult) = accepted_steps(result) + rejected_steps(result)

"""
    accepted_steps(result::TrackerResult)

Returns the number of accepted steps.
"""
accepted_steps(result::TrackerResult) = result.accepted_steps

"""
    rejected_steps(result::TrackerResult)

Returns the number of rejected_steps steps.
"""
rejected_steps(result::TrackerResult) = result.rejected_steps

###
### STATE
###

mutable struct TrackerState{M<:AbstractMatrix{ComplexF64}}
    x::Vector{ComplexF64} # current x
    x̂::Vector{ComplexF64} # last prediction
    x̄::Vector{ComplexF64} # candidate for new x
    # internal step size
    segment_stepper::SegmentStepper
    Δs_prev::Float64 # previous step size
    # path tracking algorithm
    accuracy::Float64 # norm(x - x(t))
    ω::Float64 # liptschitz constant estimate, see arxiv:1902.02968
    ω_prev::Float64
    μ::Float64 # limit accuracy
    τ::Float64 # trust region size
    norm_Δx₀::Float64 # debug info only
    extended_prec::Bool
    used_extended_prec::Bool
    refined_extended_prec::Bool
    keep_extended_prec::Bool
    norm::WeightedNorm{InfNorm}
    use_strict_β_τ::Bool

    jacobian::Jacobian{M}
    cond_J_ẋ::Float64 # estimate of cond(H(x(t),t), ẋ(t))
    code::TrackerCode.codes

    # statistics
    accepted_steps::Int
    rejected_steps::Int
    last_steps_failed::Int
    ext_accepted_steps::Int
    ext_rejected_steps::Int
end

function TrackerState(H, x₁::AbstractVector, norm::WeightedNorm{InfNorm})
    x = convert(Vector{ComplexF64}, x₁)
    x̂ = zero(x)
    x̄ = zero(x)

    segment_stepper = SegmentStepper(1.0, 0.0)
    Δs_prev = 0.0
    accuracy = 0.0
    μ = eps()
    ω = 1.0
    ω_prev = 1.0
    τ = Inf
    norm_Δx₀ = NaN
    use_strict_β_τ = false
    used_extended_prec = extended_prec = keep_extended_prec = refined_extended_prec = false
    jacobian = Jacobian(zeros(ComplexF64, size(H)))
    cond_J_ẋ = NaN
    code = TrackerCode.tracking
    u = zeros(ComplexF64, size(H, 1))
    accepted_steps = rejected_steps = ext_accepted_steps = ext_rejected_steps = 0
    last_steps_failed = 0

    TrackerState(
        x,
        x̂,
        x̄,
        segment_stepper,
        Δs_prev,
        accuracy,
        ω,
        ω_prev,
        μ,
        τ,
        norm_Δx₀,
        extended_prec,
        used_extended_prec,
        refined_extended_prec,
        keep_extended_prec,
        norm,
        use_strict_β_τ,
        jacobian,
        cond_J_ẋ,
        code,
        accepted_steps,
        rejected_steps,
        last_steps_failed,
        ext_accepted_steps,
        ext_rejected_steps,
    )
end

Base.show(io::IO, state::TrackerState) = print_fieldnames(io, state)
Base.show(io::IO, ::MIME"application/prs.juno.inline", state::TrackerState) = state
function Base.getproperty(state::TrackerState, sym::Symbol)
    if sym === :t
        return getfield(state, :segment_stepper).t
    elseif sym == :Δt
        return getfield(state, :segment_stepper).Δt
    elseif sym == :t′
        return getfield(state, :segment_stepper).t′
    elseif sym == :last_step_failed
        return getfield(state, :last_steps_failed) > 0
    else # fallback to getfield
        return getfield(state, sym)
    end
end

steps(S::TrackerState) = S.accepted_steps + S.rejected_steps
ext_steps(S::TrackerState) = S.ext_accepted_steps + S.ext_rejected_steps

###
### TRACKER
###

"""
    Tracker(H::AbstractHomotopy;
            options = TrackerOptions(),
            weighted_norm_options = WeightedNormOptions())

Construct a tracker for the given homotopy `H`. The algorithm computes along the path ``x(t)``
the local derivatives up to order 4.
For `options` see also [`TrackerOptions`](@ref).
The algorithm uses as a weighted infinity norm to measure distances.
See also [`WeightedNorm`](@ref).

[^Tim20]: Timme, S. "Mixed Precision Path Tracking for Polynomial Homotopy Continuation". arXiv:1902.02968 (2020)

## Example

We want to solve the system
```julia
@var x y t
F = System([x^2 + y^2 - 3, 2x^2 + 0.5x*y + 3y^2 - 2])
```
using a total degree homotopy and `Tracker`.
```julia
# construct start system and homotopy
G = System(im * [x^2 - 1, y^2 - 1])
H = StraightLineHomotopy(G, F)
start_solutions = [[1,1], [-1,1], [1,-1], [-1,-1]]
# construct tracker
tracker = Tracker(H)
# track each start solution separetely
results = track.(tracker, start_solutions)
println("# successfull: ", count(is_success, results))
```
We see that we tracked all 4 paths successfully.
```
# successfull: 4
```
"""
struct Tracker{H<:AbstractHomotopy,M<:AbstractMatrix{ComplexF64}}
    homotopy::H
    predictor::Predictor
    corrector::NewtonCorrector
    # these are mutable
    state::TrackerState{M}
    options::TrackerOptions
end

Tracker(H::ModelKit.Homotopy; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[], kwargs...) =
    Tracker(fixed(H; compile = compile); kwargs...)
function Tracker(
    H::AbstractHomotopy,
    x::AbstractVector = zeros(size(H, 2));
    weighted_norm_options::WeightedNormOptions = WeightedNormOptions(),
    options = TrackerOptions(),
)
    if !isa(options, TrackerOptions)
        options = TrackerOptions(; options...)
    else
        options = deepcopy(options)
    end
    norm = WeightedNorm(ones(size(H, 2)), InfNorm(), weighted_norm_options)
    state = TrackerState(H, x, norm)
    predictor = Predictor(H)
    corrector = NewtonCorrector(options.parameters.a, state.x, size(H, 1))
    Tracker(H, predictor, corrector, state, options)
end

Base.show(io::IO, C::Tracker) = print(io, typeof(C), "()")
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::Tracker) = x
Base.broadcastable(C::Tracker) = Ref(C)

"""
    state(tracker::Tracker)

Return the state of the tracker.
"""
state(tracker::Tracker) = tracker.state

"""
    status(tracker::Tracker)

Get the current [`TrackerCode`](@ref) of `tracker`.
"""

"""
    solution(tracker::Tracker)

Get the current solution.
"""
solution(T::Tracker) = get_solution(T.homotopy, T.state.x, T.state.t)

status(tracker::Tracker) = tracker.state.code

LA.cond(tracker::Tracker) = tracker.state.cond_J_ẋ

function LA.cond(tracker::Tracker, x, t, d_l = nothing, d_r = nothing)
    J = tracker.state.jacobian
    evaluate_and_jacobian!(tracker.corrector.r, matrix(J), tracker.homotopy, x, t)
    updated!(J)
    LA.cond(J, d_l, d_r)
end
# Step Size

_h(a) = 2a * (√(4 * a^2 + 1) - 2a)

# intial step size
function initial_step_size(
    state::TrackerState,
    predictor::Predictor,
    options::TrackerOptions,
)
    a = options.parameters.β_a * options.parameters.a
    p = order(predictor)
    ω = state.ω
    e = local_error(predictor)
    if isinf(e)
        # don't have anything we can use, so just use a conservative number
        # (in all well-behaved cases)
        e = 1e5
    end
    τ = trust_region(predictor)
    Δs₁ = nthroot((√(1 + 2 * _h(a)) - 1) / (ω * e), p) / options.parameters.β_ω_p
    Δs₂ = options.parameters.β_τ * τ
    Δs = nanmin(Δs₁, Δs₂)
    min(Δs, options.max_step_size, options.max_initial_step_size)
end

function update_stepsize!(
    state::TrackerState,
    result::NewtonCorrectorResult,
    options::TrackerOptions,
    predictor::Predictor;
    ad_for_error_estimate::Bool = true,
)
    a = options.parameters.β_a * options.parameters.a
    p = order(predictor)
    ω = clamp(state.ω + 2(state.ω - state.ω_prev), state.ω, 8state.ω)
    τ = state.τ #trust_region(predictor)
    if is_converged(result)
        e = local_error(predictor)
        Δs₁ = nthroot((√(1 + 2 * _h(a)) - 1) / (ω * e), p) / options.parameters.β_ω_p
        Δs₂ = options.parameters.β_τ * τ
        if state.use_strict_β_τ || dist_to_target(state.segment_stepper) < Δs₂
            Δs₂ = options.parameters.strict_β_τ * τ
        end
        Δs = min(nanmin(Δs₁, Δs₂), options.max_step_size)
        if state.use_strict_β_τ && dist_to_target(state.segment_stepper) < Δs
            Δs *= options.parameters.strict_β_τ
        end
        # increase step size by a factor of at most 10 in one step
        Δs = min(Δs, 10 * state.Δs_prev)
        if state.last_step_failed
            Δs = min(Δs, state.Δs_prev)
        end
    else
        j = result.iters - 2
        Θ_j = nthroot(result.θ, 1 << j)
        h_Θ_j = _h(Θ_j)
        h_a = _h(0.5a)
        if isnan(Θ_j) ||
           result.return_code == NEWT_SINGULARITY ||
           isnan(result.accuracy) ||
           result.iters == 1 ||
           h_Θ_j < h_a

            Δs = 0.25 * state.segment_stepper.Δs
        else
            Δs =
                nthroot((√(1 + 2 * _h(0.5a)) - 1) / (√(1 + 2 * _h(Θ_j)) - 1), p) *
                state.segment_stepper.Δs
        end
    end
    propose_step!(state.segment_stepper, Δs)
    nothing
end

function check_terminated!(state::TrackerState, options::TrackerOptions)
    if state.extended_prec || !options.extended_precision
        @unpack a, min_newton_iters = options.parameters
        tol_acc = a^(2^min_newton_iters - 1) * _h(a)
    else
        tol_acc = Inf
    end

    t′ = state.segment_stepper.t′
    t = state.segment_stepper.t
    if is_done(state.segment_stepper)
        state.code = TrackerCode.success
    elseif steps(state) ≥ options.max_steps
        state.code = TrackerCode.terminated_max_steps
    elseif state.ω * state.μ > tol_acc
        state.code = TrackerCode.terminated_accuracy_limit
        # elseif state.last_steps_failed ≥ 3 && state.cond_J_ẋ > options.terminate_cond
        #     state.code = TrackerCode.terminated_ill_conditioned
    elseif state.segment_stepper.Δs < options.min_step_size
        state.code = TrackerCode.terminated_step_size_too_small
    elseif fast_abs(t′ - t) ≤ 2eps(fast_abs(t))
        state.code = TrackerCode.terminated_step_size_too_small
    elseif options.min_rel_step_size > 0 &&
           t′ != state.segment_stepper.target &&
           fast_abs(t′ - t) < fast_abs(t) * options.min_rel_step_size
        state.code = TrackerCode.terminated_step_size_too_small
    end
    nothing
end

function update_predictor!(tracker::Tracker, x̂ = nothing, Δs = nothing, t_prev = NaN)
    @unpack predictor, homotopy, state = tracker
    update!(predictor, homotopy, state.x, state.t, state.jacobian, state.norm, x̂)
end

"""
    init!(tracker::Tracker, x₁, t₁, t₀)

Setup `tracker` to track `x₁` from `t₁` to `t₀`.

    init!(tracker::Tracker, t₀)

Setup `tracker` to continue tracking the current solution to `t₀`.
This keeps the current state.
"""
init!(tracker::Tracker, r::TrackerResult, t₁::Number = 1.0, t₀::Number = 0.0) =
    init!(tracker, solution(r), t₁, t₀; ω = r.ω, μ = r.μ)

function init!(
    tracker::Tracker,
    x₁::AbstractVector,
    t₁::Number = 1.0,
    t₀::Number = 0.0;
    ω::Float64 = NaN,
    μ::Float64 = NaN,
    τ::Float64 = Inf,
    max_initial_step_size::Float64 = Inf,
    keep_steps::Bool = false,
    extended_precision::Bool = false,
)
    @unpack state, predictor, corrector, homotopy, options = tracker
    @unpack x, x̄, norm, jacobian = state

    # intialize state
    set_solution!(x, homotopy, x₁, t₁)
    init!(state.segment_stepper, t₁, t₀)
    state.Δs_prev = 0.0
    state.accuracy = eps()
    state.ω = 1.0
    state.keep_extended_prec = false
    state.use_strict_β_τ = false
    init!(norm, x)
    init!(jacobian)
    state.code = TrackerCode.tracking
    if !keep_steps
        state.accepted_steps = state.rejected_steps = 0
        state.ext_accepted_steps = state.ext_rejected_steps = 0
    end
    state.last_steps_failed = 0

    # compute ω and limit accuracy μ for the start value
    t = state.t
    if isnan(ω) || isnan(μ)
        a = options.parameters.a
        valid, ω, μ = init_newton!(
            x̄,
            corrector,
            homotopy,
            x,
            t,
            jacobian,
            norm;
            a = a,
            extended_precision = extended_precision,
        )
        if !valid && !extended_precision
            extended_precision = true
            valid, ω, μ = init_newton!(
                x̄,
                corrector,
                homotopy,
                x,
                t,
                jacobian,
                norm;
                a = a,
                extended_precision = true,
            )
        end
    else
        valid = true
    end
    state.used_extended_prec = state.extended_prec = extended_precision

    if !isnan(ω)
        state.ω = ω
    end
    if valid
        state.accuracy = μ
        state.μ = max(μ, eps())
    else
        evaluate_and_jacobian!(corrector.r, workspace(jacobian), homotopy, x, t)
        J = matrix(workspace(jacobian))
        if J isa StructArrays.StructArray
            corank = size(J, 2) - LA.rank(Matrix(J), rtol = 1e-14)
        else
            corank = size(J, 2) - LA.rank(J, rtol = 1e-14)
        end
        if corank > 0
            state.code = TrackerCode.terminated_invalid_startvalue_singular_jacobian
        else
            state.code = TrackerCode.terminated_invalid_startvalue
        end

        return false
    end

    # initialize the predictor
    state.τ = τ
    evaluate_and_jacobian!(corrector.r, workspace(jacobian), homotopy, state.x, t)
    updated!(jacobian)
    init!(tracker.predictor)
    update_predictor!(tracker)
    state.τ = trust_region(predictor)
    # compute initial step size
    Δs = initial_step_size(state, predictor, tracker.options)
    Δs = min(Δs, max_initial_step_size)
    propose_step!(state.segment_stepper, Δs)
    state.ω_prev = state.ω

    is_tracking(state.code)
end

function init!(tracker::Tracker, t₀::Number; max_initial_step_size::Float64 = Inf)
    @unpack state, predictor, options = tracker
    state.code = TrackerCode.tracking
    init!(state.segment_stepper, state.t, t₀)
    Δs = initial_step_size(state, predictor, tracker.options)
    Δs = min(Δs, max_initial_step_size)
    propose_step!(state.segment_stepper, Δs)
    state.Δs_prev = 0.0

    tracker
end

function update_precision!(tracker::Tracker, μ_low)
    @unpack homotopy, corrector, state, options = tracker
    @unpack μ, ω, x, t, jacobian, norm = state
    @unpack a = options.parameters

    options.extended_precision || return false

    if state.extended_prec && !state.keep_extended_prec && !isnan(μ_low) && μ_low > μ
        # check if we can go low again
        if μ_low * ω < a^7 * _h(a)
            state.extended_prec = false
            state.μ = μ_low
        end
    elseif μ * ω > a^5 * _h(a)
        use_extended_precision!(tracker)
    end

    state.extended_prec
end

function use_extended_precision!(tracker::Tracker)
    @unpack homotopy, corrector, state, options = tracker
    @unpack μ, ω, x, t, jacobian, norm = state
    @unpack a = options.parameters

    options.extended_precision || return state.μ
    !state.extended_prec || return state.μ

    state.extended_prec = true
    state.used_extended_prec = true
    # do two refinement steps
    for i = 1:2
        μ = extended_prec_refinement_step!(
            x,
            corrector,
            homotopy,
            x,
            t,
            jacobian,
            norm;
            simple_newton_step = false,
        )
    end
    state.μ = max(μ, eps())
    state.μ
end

function refine_current_solution!(tracker; min_tol::Float64 = 4 * eps())
    @unpack homotopy, corrector, state, options = tracker
    @unpack x, x̄, t, jacobian, norm = state

    μ = state.accuracy
    μ̄ = extended_prec_refinement_step!(
        x̄,
        corrector,
        homotopy,
        x,
        t,
        jacobian,
        norm;
        simple_newton_step = false,
    )
    if μ̄ < μ
        x .= x̄
        μ = μ̄
    end
    k = 1
    while (μ > min_tol && k ≤ 3)
        μ̄ = extended_prec_refinement_step!(x̄, corrector, homotopy, x, t, jacobian, norm)
        if μ̄ < μ
            x .= x̄
            μ = μ̄
        end
        k += 1
    end
    μ
end

"""
    step!(tracker::Tracker, debug::Bool = false)

Perform a single tracking step. Returns `true` if the step was accepted.
"""
function step!(tracker::Tracker, debug::Bool = false)
    @unpack homotopy, corrector, predictor, state, options = tracker
    @unpack t, Δt, t′, x, x̂, x̄, jacobian, norm = state

    debug && printstyled(
        "\nt: ",
        round(t, sigdigits = 5),
        " Δt: ",
        round(Δt; sigdigits = 5),
        "\n";
        color = state.extended_prec ? :blue : :yellow,
    )
    # Use the current approximation of x(t) to obtain estimate
    # x̂ ≈ x(t + Δt) using the predictor
    predict!(x̂, predictor, homotopy, x, t, Δt)
    update!(state.norm, x̂)

    # Correct the predicted value x̂ to obtain x̄.
    # If the correction is successfull we have x̄ ≈ x(t+Δt).
    result = newton!(
        x̄,
        corrector,
        homotopy,
        x̂,
        t′,
        jacobian,
        norm;
        ω = state.ω,
        μ = state.μ,
        extended_precision = state.extended_prec,
        first_correction = state.accepted_steps == 0,
    )

    if debug
        println("Prediction method: ", predictor.method)
        println("τ: ", predictor.trust_region)
        printstyled(result, "\n"; color = is_converged(result) ? :green : :red)
    end

    if is_converged(result)
        # move forward
        x .= x̄
        state.Δs_prev = state.segment_stepper.Δs
        step_success!(state.segment_stepper)
        state.accuracy = result.accuracy
        state.μ = max(result.accuracy, eps())
        state.ω_prev = state.ω
        state.ω = max(result.ω, 0.5 * state.ω, 0.1)
        update_precision!(tracker, result.μ_low)

        if is_done(state.segment_stepper) &&
           options.extended_precision &&
           state.accuracy > 1e-14
            state.accuracy = refine_current_solution!(tracker; min_tol = 1e-14)
            state.refined_extended_prec = true
        end
        update_predictor!(tracker, x̂, state.Δs_prev, t)
        state.τ = trust_region(predictor)

        state.accepted_steps += 1
        state.ext_accepted_steps += state.extended_prec
        state.last_steps_failed = 0
    else
        # Step failed, so we have to try with a new (smaller) step size
        state.rejected_steps += 1
        state.ext_rejected_steps += state.extended_prec
        state.last_steps_failed += 1
    end
    state.norm_Δx₀ = result.norm_Δx₀
    update_stepsize!(state, result, options, predictor; ad_for_error_estimate = false)

    check_terminated!(state, options)

    !state.last_step_failed
end

"""
    track!(tracker::Tracker, x, t₁ = 1.0, t₀ = 0.0; debug::Bool = false)

The same as [`track`](@ref) but only returns the final [`TrackerCode`](@ref).

    track!(tracker::Tracker, t₀; debug::Bool = false)

Track with `tracker` the current solution to `t₀`. This keeps the current state.
"""
function track!(
    tracker::Tracker,
    x::AbstractVector,
    t₁ = 1.0,
    t₀ = 0.0;
    ω::Float64 = NaN,
    μ::Float64 = NaN,
    extended_precision::Bool = false,
    τ::Float64 = Inf,
    keep_steps::Bool = false,
    max_initial_step_size::Float64 = Inf,
    debug::Bool = false,
)
    init!(
        tracker,
        x,
        t₁,
        t₀;
        ω = ω,
        μ = μ,
        extended_precision = extended_precision,
        τ = τ,
        keep_steps = keep_steps,
        max_initial_step_size = max_initial_step_size,
    )

    while is_tracking(tracker.state.code)
        step!(tracker, debug)
    end

    tracker.state.code
end

function track!(tracker::Tracker, r::TrackerResult, t₁ = 1.0, t₀ = 0.0; debug::Bool = false)
    track!(
        tracker,
        solution(r),
        t₁,
        t₀;
        debug = debug,
        ω = r.ω,
        μ = r.μ,
        τ = r.τ,
        extended_precision = r.extended_precision,
    )
end
function track!(
    tracker::Tracker,
    t₀;
    debug::Bool = false,
    max_initial_step_size::Float64 = Inf,
)
    init!(tracker, t₀; max_initial_step_size = max_initial_step_size)

    while is_tracking(tracker.state.code)
        step!(tracker, debug)
    end

    tracker.state.code
end

function TrackerResult(H::AbstractHomotopy, state::TrackerState)
    TrackerResult(
        Symbol(state.code),
        get_solution(H, state.x, state.t),
        state.t,
        state.accuracy,
        state.accepted_steps,
        state.rejected_steps,
        state.extended_prec || state.refined_extended_prec,
        state.used_extended_prec || state.refined_extended_prec,
        state.ω,
        state.μ,
        state.τ,
    )
end

"""
     track(tracker::Tracker, x::AbstractVector, t₁ = 1.0, t₀ = 0.0; debug::Bool = false)

Track the given solution `x` at `t₁` using `tracker` to a solution at `t₀`.

    track(tracker::Tracker, r::TrackerResult, t₁ = 1.0, t₀ = 0.0; debug::Bool = false)

Track the solution of the result `r` from `t₁` to `t₀`.
"""
@inline function track(tracker::Tracker, x, t₁ = 1.0, t₀ = 0.0; kwargs...)
    track!(tracker, x, t₁, t₀; kwargs...)
    TrackerResult(tracker.homotopy, tracker.state)
end

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

Set the start parameters of the homotopy of the tracker.
"""
start_parameters!(T::Tracker, p) = (start_parameters!(T.homotopy, p); T)

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

Set the target parameters of the homotopy of the tracker.
"""
target_parameters!(T::Tracker, p) = (target_parameters!(T.homotopy, p); T)

parameters!(T::Tracker, p, q) = (parameters!(T.homotopy, p, q); T)

# PathIterator #
struct PathIterator{T<:Tracker}
    tracker::T
    t_real::Bool
end
Base.IteratorSize(::Type{<:PathIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:PathIterator}) = Base.HasEltype()

"""
    iterator(tracker::Tracker, x₁, t₁=1.0, t₀=0.0)

Prepare a tracker to make it usable as a (stateful) iterator. Use this if you want to inspect a specific
path. In each iteration the tuple `(x,t)` is returned.

## Example

Assume you have `Tracker` `tracker` and you wan to track `x₁` from 1.0 to 0.25:
```julia
for (x,t) in iterator(tracker, x₁, 1.0, 0.25)
    println("x at t=\$t:")
    println(x)
end
```

Note that this is a stateful iterator. You can still introspect the state of the tracker.
For example to check whether the tracker was successfull
(and did not terminate early due to some problem) you can do
```julia
println("Success: ", is_success(status(tracker)))
```
"""
function iterator(tracker::Tracker, x₁, t₁ = 1.0, t₀ = 0.0; kwargs...)
    init!(tracker, x₁, t₁, t₀; kwargs...)
    PathIterator(tracker, typeof(t₁ - t₀) <: Real)
end

function current_x_t(iter::PathIterator)
    @unpack x, t = iter.tracker.state
    (copy(x), iter.t_real ? real(t) : t)
end

function Base.iterate(iter::PathIterator, state = nothing)
    state === nothing && return current_x_t(iter), 1
    iter.tracker.state.code != TrackerCode.tracking && return nothing

    while is_tracking(iter.tracker.state.code)
        step_failed = !step!(iter.tracker)
        step_failed || break
    end
    current_x_t(iter), state + 1
end


text:
"""
All the terminals are under my command. `Comonicon` is a CLI (Command Line Interface) generator
that features light-weight dependency (optional to have zero dependency),
fast start-up time and easy to use. See the [website](https://comonicon.org)
for more info.
"""
module Comonicon

using Pkg
using Markdown
using ExproniconLite

export @cast, @main

include("compat.jl")
include("options.jl")
include("ast/ast.jl")
include("codegen/julia.jl")
include("codegen/zsh.jl")
include("tools.jl")

include("frontend/frontend.jl")
include("builder/builder.jl")

end # module

ext: jl


text:
my_func(x,y) = 2x + y

derivative_of_my_func(x,y) = ForwardDiff.derivative(x->my_func(x,y),x)


text:
module Reexport

macro reexport(ex)
    isa(ex, Expr) && (ex.head == :module ||
                      ex.head == :using ||
                      (ex.head == :toplevel &&
                       all(e->isa(e, Expr) && e.head == :using, ex.args))) ||
    error("@reexport: syntax error")

#    if ex.head == :module
#        modules = {ex.args[2]}
#        ex = Expr(:toplevel, ex, Expr(:using, :., ex.args[2]))
#    elseif ex.head == :using
#        modules = {ex.args[end]}
#    else
#        modules = {e.args[end] for e in ex.args}
#        end

    if ex.head == :module
        modules = Any[ex.args[2]]
        ex = Expr(:toplevel, ex, Expr(:using, :., ex.args[2]))
    elseif ex.head == :using
        modules = Any[ex.args[end]]
    else
        modules = Any[e.args[end] for e in ex.args]
    end

    esc(Expr(:toplevel, ex,
             [:(eval(Expr(:export, names($(mod))...))) for mod in modules]...))
end

export @reexport

end # module


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Pkg25

const DEFAULT_META = "https://github.com/KristofferC/METADATA.jl"

const META_BRANCH = "metadata-v3"

isdevmetadata() = Dir.getmetabranch() == META_BRANCH

type Pkg25Error <: Exception
    msg::AbstractString
    ex::Nullable{Exception}
end
Pkg25Error(msg::AbstractString) = Pkg25Error(msg, Nullable{Exception}())
function Base.showerror(io::IO, pkgerr::Pkg25Error)
    print(io, pkgerr.msg)
    if !isnull(pkgerr.ex)
        pkgex = get(pkgerr.ex)
        if isa(pkgex, CompositeException)
            for cex in pkgex
                print(io, "\n=> ")
                showerror(io, cex)
            end
        else
            print(io, "\n")
            showerror(io, pkgex)
        end
    end
end

for file in split("dir types reqs cache read query resolve write entry")
    include("$file.jl")
end
const cd = Dir.cd

dir(path...) = Dir.path(path...)

"""
    dir() -> AbstractString

Returns the absolute path of the package directory. This defaults to
`joinpath(homedir(),".julia","v\$(VERSION.major).\$(VERSION.minor)")` on all platforms (i.e.
`~/.julia/v$(VERSION.major).$(VERSION.minor)` in UNIX shell syntax). If the `JULIA_PKGDIR`
environment variable is set, then that path is used in the returned value as
`joinpath(ENV["JULIA_PKGDIR"],"v\$(VERSION.major).\$(VERSION.minor)")`. If `JULIA_PKGDIR` is
a relative path, it is interpreted relative to whatever the current working directory is.
"""
dir()

"""
    dir(names...) -> AbstractString

Equivalent to `normpath(Pkg25.dir(),names...)` – i.e. it appends path components to the
package directory and normalizes the resulting path. In particular, `Pkg25.dir(pkg)` returns
the path to the package `pkg`.
"""
dir(names...)

"""
    init(meta::AbstractString=DEFAULT_META, branch::AbstractString=META_BRANCH)

Initialize `Pkg25.dir()` as a package directory. This will be done automatically when the
`JULIA_PKGDIR` is not set and `Pkg25.dir()` uses its default value. As part of this process,
clones a local METADATA git repository from the site and branch specified by its arguments,
which are typically not provided. Explicit (non-default) arguments can be used to support a
custom METADATA setup.
"""
init(meta::AbstractString=DEFAULT_META, branch::AbstractString=META_BRANCH) = Dir.init(meta,branch)

function __init__()
    vers = "v$(VERSION.major).$(VERSION.minor)"
    unshift!(Base.LOAD_CACHE_PATH, abspath(Dir._pkgroot(), "lib", vers))
end

"""
    edit()

Opens `Pkg25.dir("REQUIRE")` in the editor specified by the `VISUAL` or `EDITOR` environment
variables; when the editor command returns, it runs `Pkg25.resolve()` to determine and install
a new optimal set of installed package versions.
"""
edit() = cd(Entry.edit)

"""
    rm(pkg)

Remove all requirement entries for `pkg` from `Pkg25.dir("REQUIRE")` and call `Pkg25.resolve()`.
"""
rm(pkg::AbstractString) = cd(Entry.rm,pkg)

"""
    add(pkg, vers...)

Add a requirement entry for `pkg` to `Pkg25.dir("REQUIRE")` and call `Pkg25.resolve()`. If
`vers` are given, they must be `VersionNumber` objects and they specify acceptable version
intervals for `pkg`.
"""
add(pkg::AbstractString, vers::VersionNumber...) = cd(Entry.add,pkg,vers...)

"""
    available() -> Vector{String}

Returns the names of available packages.
"""
available() = cd(Entry.available)

"""
    available(pkg) -> Vector{VersionNumber}

Returns the version numbers available for package `pkg`.
"""
available(pkg::AbstractString) = cd(Entry.available,pkg)

"""
    installed() -> Dict{String,VersionNumber}

Returns a dictionary mapping installed package names to the installed version number of each
package.
"""
installed() = cd(Entry.installed)

"""
    installed(pkg) -> Void | VersionNumber

If `pkg` is installed, return the installed version number, otherwise return `nothing`.
"""
installed(pkg::AbstractString) = cd(Entry.installed,pkg)

"""
    status()

Prints out a summary of what packages are installed and what version and state they're in.
"""
status(io::IO=STDOUT) = cd(Entry.status,io)

"""
    status(pkg)

Prints out a summary of what version and state `pkg`, specifically, is in.
"""
status(pkg::AbstractString, io::IO=STDOUT) = cd(Entry.status,io,pkg)

"""
    clone(pkg)

If `pkg` has a URL registered in `Pkg25.dir("METADATA")`, clone it from that URL on the
default branch. The package does not need to have any registered versions.
"""
clone(url_or_pkg::AbstractString) = cd(Entry.clone,url_or_pkg)

"""
    clone(url, [pkg])

Clone a package directly from the git URL `url`. The package does not need to be registered
in `Pkg25.dir("METADATA")`. The package repo is cloned by the name `pkg` if provided; if not
provided, `pkg` is determined automatically from `url`.
"""
clone(url::AbstractString, pkg::AbstractString) = cd(Entry.clone,url,pkg)

"""
    checkout(pkg, [branch="master"]; merge=true, pull=true)

Checkout the `Pkg25.dir(pkg)` repo to the branch `branch`. Defaults to checking out the
"master" branch. To go back to using the newest compatible released version, use
`Pkg25.free(pkg)`. Changes are merged (fast-forward only) if the keyword argument `merge ==
true`, and the latest version is pulled from the upstream repo if `pull == true`.
"""
checkout(pkg::AbstractString, branch::AbstractString="master"; merge::Bool=true, pull::Bool=true) =
    cd(Entry.checkout,pkg,branch,merge,pull)

"""
    free(pkg)

Free the package `pkg` to be managed by the package manager again. It calls `Pkg25.resolve()`
to determine optimal package versions after. This is an inverse for both `Pkg25.checkout` and
`Pkg25.pin`.

You can also supply an iterable collection of package names, e.g., `Pkg25.free(("Pkg251",
"Pkg252"))` to free multiple packages at once.
"""
free(pkg) = cd(Entry.free,pkg)

"""
    pin(pkg)

Pin `pkg` at the current version. To go back to using the newest compatible released
version, use `Pkg25.free(pkg)`
"""
pin(pkg::AbstractString) = cd(Entry.pin,pkg)

"""
    pin(pkg, version)

Pin `pkg` at registered version `version`.
"""
pin(pkg::AbstractString, ver::VersionNumber) = cd(Entry.pin,pkg,ver)

"""
    update(pkgs...)

Update the metadata repo – kept in `Pkg25.dir("METADATA")` – then update any fixed packages
that can safely be pulled from their origin; then call `Pkg25.resolve()` to determine a new
optimal set of packages versions.

Without arguments, updates all installed packages. When one or more package names are provided as
arguments, only those packages and their dependencies are updated.
"""
update(upkgs::AbstractString...) = cd(Entry.update,Dir.getmetabranch(),Set{String}([upkgs...]))

"""
    resolve()

Determines an optimal, consistent set of package versions to install or upgrade to. The
optimal set of package versions is based on the contents of `Pkg25.dir("REQUIRE")` and the
state of installed packages in `Pkg25.dir()`, Packages that are no longer required are moved
into `Pkg25.dir(".trash")`.
"""
resolve() = cd(Entry.resolve)

"""
    build()

Run the build scripts for all installed packages in depth-first recursive order.
"""
build() = cd(Entry.build)

"""
    build(pkgs...)

Run the build script in `deps/build.jl` for each package in `pkgs` and all of their
dependencies in depth-first recursive order. This is called automatically by `Pkg25.resolve()`
on all installed or updated packages.
"""
build(pkgs::AbstractString...) = cd(Entry.build,[pkgs...])

"""
    test(; coverage=false)

Run the tests for all installed packages ensuring that each package's test dependencies are
installed for the duration of the test. A package is tested by running its
`test/runtests.jl` file and test dependencies are specified in `test/REQUIRE`.
Coverage statistics for the packages may be generated by passing `coverage=true`.
The default behavior is not to run coverage.
"""
test(;coverage::Bool=false) = cd(Entry.test; coverage=coverage)

"""
    test(pkgs...; coverage=false)

Run the tests for each package in `pkgs` ensuring that each package's test dependencies are
installed for the duration of the test. A package is tested by running its
`test/runtests.jl` file and test dependencies are specified in `test/REQUIRE`.
Coverage statistics for the packages may be generated by passing `coverage=true`.
The default behavior is not to run coverage.
"""
test(pkgs::AbstractString...; coverage::Bool=false) = cd(Entry.test,AbstractString[pkgs...]; coverage=coverage)

"""
    dependents(packagename)

List the packages that have `packagename` as a dependency.
"""
dependents(packagename::AbstractString) = Reqs.dependents(packagename)

"""
    setprotocol!(proto)

Set the protocol used to access GitHub-hosted packages. Defaults to 'https', with a blank
`proto` delegating the choice to the package developer.
"""
setprotocol!(proto::AbstractString) = Cache.setprotocol!(proto)

# point users to Pkg25Dev
register(args...) =
    error("Pkg25.register(pkg,[url]) has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

tag(pkg, ver=nothing, commit=nothing) =
    error("Pkg25.tag(pkg, [ver, [commit]]) has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

publish() =
    error("Pkg25.publish() has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

generate(pkg, license) =
    error("Pkg25.generate(pkg, license) has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

license(lic=nothing) =
    error("Pkg25.license([lic]) has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

submit(pkg, commit=nothing) =
    error("Pkg25.submit(pkg[, commit]) has been moved to the package Pkg25Dev.jl.\n",
          "Run Pkg25.add(\"Pkg25Dev\") to install Pkg25Dev on Julia v0.5-")

end # module


text:
include(joinpath("docs", "make.jl"))

cd(@__DIR__)

# fix all links to prepend them with `JuDocTemplates.jl/`
for (root, _, files) ∈ walkdir(joinpath(@__DIR__, "docs", "build"))
    for file ∈ files
        endswith(file, ".html") || continue
        path = joinpath(root, file)
        html = read(path, String)
        html = replace(html, "href=\"/" => "href=\"/JuDocTemplates.jl/")
        html = replace(html, "src=\"/" => "src=\"/JuDocTemplates.jl/")
        write(path, html)
    end
end

const JS_GHP = """
    var ghpages = require('gh-pages');
    ghpages.publish('docs/build', function(err) {});
    """
run(`node -e $JS_GHP`)


text:
using Gillespie
using Base.Test

include("../examples/lotka.jl")
@test isequal(gillespie_data[end,2],741)
@test isequal(gillespie_data[end,3],1048)

include("../examples/logistic_growth.jl")
@test isequal(data[end,2],1008)

include("../examples/decaying_dimer.jl")
@test isequal(data[end,2],312)
@test isequal(data[end,3],714)
@test isequal(data[end,4],856)

include("../examples/sir.jl")
@test isequal(data[end,2],1)
@test isequal(data[end,3],181)
@test isequal(data[end,4],818)

include("../examples/sir2.jl")
@test isequal(data[end,2],1)
@test isequal(data[end,3],181)
@test isequal(data[end,4],818)


text:
function Qfunc_UICM(x)
    u1 = x + 2.0
    if(u1) < -30
        u1 = -30
        wn1 = 1e-10
    elseif u1 <= 0.6
        x1 = exp(u1-1);
	    wn1 = x1 * (1.0 + 4.0 * x1 / 3.0) / (1.0 + x1 * (7 / 3 + 5 * x1 / 6));
    else
        wn1 = (u1-1) - 24 * (((u1-1) + 2) * (u1-1) - 3) / ((7 * (u1-1) + 58) * (u1-1) + 127);
    end
    zn1 = (u1-1) - wn1 - log(wn1);
    temp3 = 1 + wn1;
    y2 = 2 * temp3 * (temp3 + 2 * zn1/3) - zn1;
    en1 = zn1 * y2 / (temp3 * (y2 - zn1));
    wn1 = wn1 * (1 + en1)
    return wn1 + 1.0
end

function ϕsa_UICM(Vg,Vfb,Na,tox,T)
    G = γ(Na,tox)
    phit = kb*T
    #f(x) = (Vg-Vfb-x)^2/G^2-phit*(exp(-x/phit)-1)-x
    #return fzero(f,Vg-Vfb,Roots.Order0(),atol=1e-6, rtol=1e-6)
    alpha=-G*sqrt(phit);
    phisa=0;
    ug=Vg-Vfb
    gammasquare = G*G
    phisa = 0.0
    if (ug<alpha)
        phisa = - 2 * phit * log( ug / alpha )
    end
    if (ug > phit)
        auxphisat=1+4*(ug-phit)/(gammasquare)
        phisa=ug+(gammasquare/2)*(1-sqrt(auxphisat))
    end
    if (ug>0)
        phisa = ug<phisa ? ug : phisa

	else
        phisa = phisa>ug ? phisa : ug
    end
    f(x) = (Vg-Vfb-x)^2/G^2-(phit*(exp(-x/phit)-1)+x)
    return fzero(f,phisa,Roots.Order0(),atol=1e-12,rtol=1e-6,maxevals=10000)

end

function n_UICM(Vg,Vfb,Na,tox,T)
    G = γ(Na,tox)
    phit = kb*T
    phisa =ϕsa_UICM(Vg,Vfb,Na,tox,T)
    signphisa = sign(phisa)
    phisa2= abs(phisa);
    safesqrt(x) = x> 0.0 ? sqrt(x) : 0.0
	if (phisa2<=0.00008)
	    phisa2=-0.00002008;
	    n=1+(G*(1-exp(-(phisa2)/phit)))/(-2*safesqrt((phisa2)+phit*(exp(-(phisa2)/phit))-phit));
	else
	    n=1+(G*(1-exp(-(phisa)/phit)))/(2*signphisa*safesqrt((phisa)+phit*(exp(-(phisa)/phit))-phit));
    end
    return n
end

function Vp_UICM(Vg,Vfb,Na,tox,T,phif,sigma =0.0 )
    G = γ(Na,tox)
    phit = kb*T
    phisa = ϕsa_UICM(Vg,Vfb,Na,tox,T)
    η = n_UICM(Vg,Vfb,Na,tox,T)
	vp = phisa - 2 * phif - phit*(1+log(η/(η-1)))
    #vp = vp + sigma / n * (vd + vs);
	return vp
end

function Vth_UICM(Vfb,Na,tox,T,phif,sigma=0.0)
    phit = kb*T
    G = γ(Na,tox)
    k = 0#log(n/(n-1))
    #uprintln(n,",",phit,",",G,",",k,",",Vfb,",",Na,",",tox,",",T,",",phif)
    Vt0 = Vfb+2*phif+phit*k+G*sqrt(2*phif+phit*k)
    n = n_UICM(Vt0,Vfb,Na,tox,T)
    k = log(n/(n-1))
    return Vfb+2*phif+phit*k+G*sqrt(2*phif+phit*k)
end

function id_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    Cox = ϵ₀*3.9/tox
    return Cox/2.0*n*phit^2*((qis^2-2qis)-(qid^2-2*qid))
end

function Cgs_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    #n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    Cox = ϵ₀*3.9/tox
    α = (qid+1)/(qis+1)
    return (2/3.0)*Cox*(1+2*α)/(1+α)/(1+α) * qis/(1+qis)
end

function Cgd_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    #n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    Cox = ϵ₀*3.9/tox
    α = (qid+1)/(qis+1)
    return (2/3.0)*Cox*(α*α+2*α)/(1+α)/(1+α) * qid/(1+qid)
end

function Cbs_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    η = n_UICM(Vg,Vfb,Na,tox,T)
    return (η-1)*Cgs_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma)
end

function Cds_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    η = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    Cox = ϵ₀*3.9/tox
    α = (qid+1)/(qis+1)
    return (-4.0/15.0)*η*Cox*(α+3*α+α*α)/(1+α)/(1+α)/(1+α) *qis/(1+qis)
end

function Cbd_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    η = n_UICM(Vg,Vfb,Na,tox,T)
    return (η-1.0)*Cgd_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma)
end
function Csd_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    Cds_UICM(Vg,Vd,Vs,Vfb,Na,tox,T,phif)
end

function Cgb_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    η = n_UICM(Vg,Vfb,Na,tox,T)
    Cox = ϵ₀*3.9/tox
    return (η-1)/η*(Cox-Cgs_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif)-Cgd_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif))
end

function α_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    return  (qid+1)/(qis+1)
end

function ic_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif,sigma=0.0)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    #n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    #Cox = ϵ₀*3.9/tox
    IF = (qis-1)^2-1.0
    IR = (qid-1)^2-1.0

    return IF-IR
end

function Id_UICM(Vg,Vs,Vd,Vfb,Na,tox,T,phif;sigma=0.0,Kp=5e-2,β=-1.5)
    vp = Vp_UICM(Vg,Vfb,Na,tox,T,phif)
    phit = kb*T
    n = n_UICM(Vg,Vfb,Na,tox,T)
    qis = Qfunc_UICM((vp-Vs)/phit)
    qid = Qfunc_UICM((vp-Vd)/phit)
    Cox = ϵ₀*3.9/tox
    return Kp*(T/300)^(β)*Cox/2.0*n*phit^2*((qis^2-2qis)-(qid^2-2*qid))
end

mutable struct ACMModel <: TransistorModel
    Vfb
    Nb
    tox
    T0
    ϕB
    sigma::Number
    W::Number
    L::Number
    μ_0
    β
    mos::Union{MOSStructure,Nothing}
end

ACMModel(m::MOSStructure,W=1.0,L=1.0;sigma=0.0,μ_0=10,β=-1.5) = ACMModel(VFB(m),abs(C(m)),m.tox,Temperature(m),ϕb(m),sigma,W,L,μ_0,β,m)

id(vg,vd,vs,m::ACMModel) = id_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
ic(vg,vd,vs,m::ACMModel) = ic_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Id(vg,vd,vs,m::ACMModel) = Id_UICM(vg,vs,vd,m.Vfb,m.Nb,m.tox,m.T0,m.ϕB; sigma = m.sigma, β=m.β,Kp=m.W/m.L*m.μ_0)
IdNR(vg,vd,vs,m::ACMModel) = Id(vg,vd,vs,m)-gm(vg,vd,vs,m)*(vg-vs)-gds(vg,vd,vs,m)*(vd-vs)
gm(vg,vd,vs,m::ACMModel) = Zygote.ForwardDiff.derivative(x-> Id(x, vd, vs,m),vg)
gds(vg,vd,vs,m::ACMModel) = Zygote.ForwardDiff.derivative(x-> Id(vg, x, vs,m),vd)
Vp(Vg,m::ACMModel) = Vp_UICM(Vg,m.Vfb,m.Nb,m.tox,m.T0,m.ϕB)
nq(Vg,m::ACMModel) = n_UICM(Vg,m.Vfb,m.Nb,m.tox,m.T)
ϕsa(Vg,m::ACMModel) = ϕsa_UICM(Vg,m.Vfb,m.Nb,m.tox,m.T)
function If(Vg,Vd,m::ACMModel)
    vp = Vp(Vg,m)
    phit = kb*m.T0
    qid = Qfunc_UICM((vp-Vd)/phit)
    return (qid+1)^2-1.0
end
Vth(m::ACMModel) = Vth_UICM(m.Vfb,m.Nb,m.tox,m.T0,m.ϕB)

Cgs(vg,vd,vs,m::ACMModel) = Cgs_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Cgd(vg,vd,vs,m::ACMModel) = Cgd_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Cbs(vg,vd,vs,m::ACMModel) = Cbs_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Cgb(vg,vd,vs,m::ACMModel) = Cgb_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Csd(vg,vd,vs,m::ACMModel) = Csd_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Cds(vg,vd,vs,m::ACMModel) = Cds_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
Cbd(vg,vd,vs,m::ACMModel) = Cbd_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)
α(vg,vd,vs,m::ACMModel) = α_UICM(vg,vs,vd,m.Vfb,m.ϕB,m.tox,m.T0,m.ϕB,m.sigma)

setW!(Wn,m::ACMModel) = begin m.W = Wn end
setL!(Ln,m::ACMModel) = begin m.L = Ln end
function setT!(T,m::ACMModel)
    Sdist = m.mos.semiconductor.dist
    Ssemi = m.mos.semiconductor.material
    nMos = MOSStructure(m.mos.metal,m.mos.oxide,SemiconductorData(T,Sdist,Ssemi),m.mos.tox,m.mos.tsi)
    nModel = ACMModel(nMos,m.W,m.L;μ_0=m.μ_0)
    m.Vfb = nModel.Vfb
    m.Nb = nModel.Nb
    m.tox = nModel.tox
    m.T0 = nModel.T0
    m.ϕB = nModel.ϕB
    m.sigma = nModel.sigma
    m.W = nModel.W
    m.L = nModel.L
    m.μ_0 = nModel.μ_0
    m.β = nModel.β
    m.mos = nMos
end


text:
function my_f(x,y)
    return 2*x + 3*y
end

function derivative_of_my_f(x,y)
    return ForwardDiff.derivative(x -> my_f(x,y),x)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineSizeListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{VirtualMachineSize} } # spec name: value

    function VirtualMachineSizeListResult(;value=nothing)
        o = new()
        validate_property(VirtualMachineSizeListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type VirtualMachineSizeListResult

const _property_map_VirtualMachineSizeListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_VirtualMachineSizeListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{VirtualMachineSize}")
Base.propertynames(::Type{ VirtualMachineSizeListResult }) = collect(keys(_property_map_VirtualMachineSizeListResult))
Swagger.property_type(::Type{ VirtualMachineSizeListResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineSizeListResult[name]))}
Swagger.field_name(::Type{ VirtualMachineSizeListResult }, property_name::Symbol) =  _property_map_VirtualMachineSizeListResult[property_name]

function check_required(o::VirtualMachineSizeListResult)
    true
end

function validate_property(::Type{ VirtualMachineSizeListResult }, name::Symbol, val)
end


text:
function show()
    myfunc = (x, y) -> (x + y)
    println(join([myfunc(1, 2)], " "))
end

function main()
    show()
end

main()


text:
using Pkg
Pkg.update()
Pkg.add(PackageSpec(url="https://github.com/essenciary/Genie.jl"))
Pkg.add(PackageSpec(url="https://github.com/JuliaDB/DBI.jl"))
Pkg.add(PackageSpec(url="https://github.com/JuliaDB/PostgreSQL.jl"))
Pkg.add(PackageSpec(url="https://github.com/essenciary/Flax.jl"))
Pkg.add(PackageSpec(url="https://github.com/essenciary/SearchLight.jl"))
Pkg.add("JuMP")
Pkg.add("Convex")
Pkg.add("MathProgBase")
Pkg.add("MathOptInterface")
Pkg.add("Clp")
Pkg.add("Cbc")
Pkg.add("GLPK")
Pkg.add("GLPKMathProgInterface")
Pkg.add("Ipopt")
Pkg.add("Optim")

#RUN julia -e 'Pkg.clone("https://github.com/JuliaDB/DBI.jl.git")'

max_stars_repo_name: rpb/Pivotal-Optimization

max_issues_repo_name: rpb/Pivotal-Optimization

max_forks_repo_name: rpb/Pivotal-Optimization


text:
# -----------------------------------------------------------------------
# Cross entropy seach algorithm
# -----------------------------------------------------------------------
# Required packages
using Statistics, DataFrames
include("random_init.jl")

"""
    apply_expansion(v,γ)

## Description
Function to apply bound expansion to the domain
## Arguments
- v: value
- γ: expansion factor
"""
#apply_expansion(v,γ) = Int(round(v*γ))

"""
    find_newdomain(data::DataFrame, nbest::Int)

## Description
Find domain range for the next iteration (CE method)
## Arguments
- data: dataframe with coordinates x and y and the values measured (works only for OneCylType by now)
- nbest: the number of best values that should be taken
"""
function find_newdomain(data::DataFrame, nbest::Int)

    sort!(data, order(:pm, rev=true))
    best = data[1:nbest,:]
    ncoord = ncol(data) - 1

    range_coords = [minimum(best[!,i]):maximum(best[!,i]) for i in 1:ncoord]

    return range_coords
end

"""
    crossentropy_criteria(new_mean::AbstractVector, new_std::AbstractVector, previous_mean::AbstractVector; mean_tol::AbstractVector, std_tol::AbstractVector)

## Description
Criteria for cross entropy search method
## Arguments
- new_mean: array with current mean coordinate values
- new_std: array with current standard deviation coordinate values
- previous_means: array with previous mean coordinate values
- mean_tol: tolerance criteria for mean variation (for each direction)
- std_tol: tolerance criteria for standard deviation (for each direction)
"""
function crossentropy_criteria(new_mean::AbstractVector, new_std::AbstractVector, previous_mean::AbstractVector;
                                mean_tol::AbstractVector, std_tol::AbstractVector)

    Δmean = abs.(new_mean - previous_mean)
    if all(Δmean .< mean_tol) && all(new_std .< std_tol)
        return true
    else
        return false
    end
end

"""
    crossentropy(df::AbstractVector, nbest::Int; mean_tol::AbstractVector, std_tol::AbstractVector)

## Description
Implementation of one iteration for the cross entropy seach method
## Arguments
- df: vector with all the dataframes obtained during the process (check experiments.jl)
- nbest: number of values which should be selected for the CE
- mean_tol: tolerance criteria for mean variation (for each direction)
- std_tol: tolerance criteria for standard deviation (for each direction)
## Output
Returns the bounds for a new domain.
"""
function crossentropy(df::AbstractVector, nbest::Int; mean_tol::AbstractVector, std_tol::AbstractVector)

    n = length(df)
    i = n - 1

    if i == 0
        converge = false
    else
        ncoord = ncol(df[end]) - 1
        previous_mean = [mean(df[i][:,q]) for q in 1:ncoord]
        new_mean = [mean(df[end][:,q]) for q in 1:ncoord]
        new_std = [std(df[end][:,q]) for q in 1:ncoord]
        converge = crossentropy_criteria(new_mean, new_std, previous_mean, mean_tol=mean_tol, std_tol=std_tol)
    end

    lims = find_newdomain(df[end], nbest)

    return lims, converge
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    # need to turn on comments to skip the "fold along ..." lines at the bottom
    dotsMtx = readdlm(fname, ',', Int, '\n', comments=true, comment_char='f')
    dots = [dotsMtx[i,:] for i=1:size(dotsMtx)[1]]

    # read in folds
    folds = []
    open(fname, "r") do f
        readuntil(f, "f")
        line = readline(f)
        push!(folds, (line[11], parse(Int, SubString(line,13))))
        for line in eachline(f)
            push!(folds, (line[12], parse(Int, SubString(line,14))))
        end
    end

    return dots, folds
end

function Fold!(dots::Array{Array{Int,1},1}, fold)
    ndx = fold[1] == 'x' ? 1 : 2
    newDots = Array{Int,1}[]

    for i=1:length(dots)
        if fold[2] < dots[i][ndx] <= 2*fold[2]
            dots[i][ndx] = 2*fold[2] - dots[i][ndx]
            if count(x -> x == dots[i], newDots) == 0
                push!(newDots, dots[i])
            end
        elseif dots[i][ndx] < fold[2]
            if count(x -> x == dots[i], newDots) == 0
                push!(newDots, dots[i])
            end
        end
    end

    return newDots
end

function Part1(dots::Array{Array{Int,1},1}, folds::Array)
    return length(Fold!(dots, folds[1]))
end

function Part2(dots::Array{Array{Int,1},1}, folds::Array)
    for fold in folds
        dots = Fold!(dots, fold)
    end

    grid = fill(" ", max([x[1] for x in dots]...)+1, max([x[2] for x in dots]...)+1)

    foreach(x -> grid[x[1]+1, x[2]+1] = "#", dots)

    foreach(x -> println(join(grid[:,x], "")), 1:size(grid)[2])
end

#dataFile = "../data/test13.txt"
dataFile = "../data/day13.txt"

data, folds = ReadData(dataFile)
println(Part1(data, folds))

data, folds = ReadData(dataFile)
Part2(data, folds)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthenticationV1TokenRequestSpec <: SwaggerModel
    audiences::Any # spec type: Union{ Nothing, Vector{String} } # spec name: audiences
    boundObjectRef::Any # spec type: Union{ Nothing, IoK8sApiAuthenticationV1BoundObjectReference } # spec name: boundObjectRef
    expirationSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: expirationSeconds

    function IoK8sApiAuthenticationV1TokenRequestSpec(;audiences=nothing, boundObjectRef=nothing, expirationSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("audiences"), audiences)
        setfield!(o, Symbol("audiences"), audiences)
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("boundObjectRef"), boundObjectRef)
        setfield!(o, Symbol("boundObjectRef"), boundObjectRef)
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("expirationSeconds"), expirationSeconds)
        setfield!(o, Symbol("expirationSeconds"), expirationSeconds)
        o
    end
end # type IoK8sApiAuthenticationV1TokenRequestSpec

const _property_map_IoK8sApiAuthenticationV1TokenRequestSpec = Dict{Symbol,Symbol}(Symbol("audiences")=>Symbol("audiences"), Symbol("boundObjectRef")=>Symbol("boundObjectRef"), Symbol("expirationSeconds")=>Symbol("expirationSeconds"))
const _property_types_IoK8sApiAuthenticationV1TokenRequestSpec = Dict{Symbol,String}(Symbol("audiences")=>"Vector{String}", Symbol("boundObjectRef")=>"IoK8sApiAuthenticationV1BoundObjectReference", Symbol("expirationSeconds")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }) = collect(keys(_property_map_IoK8sApiAuthenticationV1TokenRequestSpec))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1TokenRequestSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1TokenRequestSpec[property_name]

function check_required(o::IoK8sApiAuthenticationV1TokenRequestSpec)
    (getproperty(o, Symbol("audiences")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, name::Symbol, val)
end


text:
using HTTP
using JSONSchema
using JSON
using Test
using ZipFile

const TEST_SUITE_URL = "https://github.com/json-schema-org/JSON-Schema-Test-Suite/archive/2.0.0.zip"

const SCHEMA_TEST_DIR = let
    dest_dir = mktempdir()
    dest_file = joinpath(dest_dir, "test-suite.zip")
    download(TEST_SUITE_URL, dest_file)
    for f in ZipFile.Reader(dest_file).files
        filename = joinpath(dest_dir, "test-suite", f.name)
        if endswith(filename, "/")
            mkpath(filename)
        else
            write(filename, read(f, String))
        end
    end
    joinpath(dest_dir, "test-suite", "JSON-Schema-Test-Suite-2.0.0", "tests")
end

const LOCAL_TEST_DIR = mktempdir(SCHEMA_TEST_DIR)

# Write test files for locally referenced schema files.
#
# These files have the same format as JSON Schema org test files. They are written
# to a sibling directory to JSON-Schema-Test-Suite-master/tests/draft* directories
# so they can be consumed the same way as the draft*/*.json test files.
# sibling directory for testing a relative path containing "../"
const REF_LOCAL_TEST_DIR = mktempdir(SCHEMA_TEST_DIR)

write(
    joinpath(REF_LOCAL_TEST_DIR, "localReferenceSchemaOne.json"),
    """{
    "type": "object",
    "properties": {"localRefOneResult": {"type": "string"}}
}""")

write(
    joinpath(REF_LOCAL_TEST_DIR, "localReferenceSchemaTwo.json"),
    """{
    "type": "object",
    "properties": {"localRefTwoResult": {"type": "number"}}
}""")

write(
    joinpath(REF_LOCAL_TEST_DIR, "nestedLocalReference.json"),
    """{
    "type": "object",
    "properties": {
        "result": {
            "\$ref": "file:localReferenceSchemaOne.json#/properties/localRefOneResult"
        }
    }
}""")

write(
    joinpath(LOCAL_TEST_DIR, "localReferenceTest.json"),
    """[{
    "description": "test locally referenced schemas",
    "schema": {
        "type": "object",
        "properties": {
            "result1": { "\$ref": "file:../$(basename(abspath(REF_LOCAL_TEST_DIR)))/localReferenceSchemaOne.json#/properties/localRefOneResult" },
            "result2": { "\$ref": "../$(basename(abspath(REF_LOCAL_TEST_DIR)))/localReferenceSchemaTwo.json#/properties/localRefTwoResult" }
        },
        "oneOf": [{
            "required": ["result1"]
        }, {
            "required": ["result2"]
        }]
    },
    "tests": [{
        "description": "reference only local schema 1",
        "data": {"result1": "some text"},
        "valid": true
    }, {
        "description": "reference only local schema 2",
        "data": {"result2": 1234},
        "valid": true
    }, {
        "description": "incorrect reference to local schema 1",
        "data": {"result1": true},
        "valid": false
    }, {
        "description": "reference neither local schemas",
        "data": {"result": true},
        "valid": false
    }, {
        "description": "reference both local schemas",
        "data": {"result1": "some text", "result2": 500},
        "valid": false
    }]
}]""")

write(
    joinpath(LOCAL_TEST_DIR, "nestedLocalReferenceTest.json"),
    """[{
    "description": "test locally referenced schemas",
    "schema": {
        "type": "object",
        "properties": {
            "result": {
                "\$ref": "file:../$(basename(abspath(REF_LOCAL_TEST_DIR)))/nestedLocalReference.json#/properties/result"
            }
        }
    },
    "tests": [{
        "description": "nested reference, correct type",
        "data": {"result": "some text"},
        "valid": true
    }, {
        "description": "nested reference, incorrect type",
        "data": {"result": 1234},
        "valid": false
    }]
}]""")

@testset "Draft 4/6" begin
    # Note(odow): I didn't want to use a mutable reference like this for the web-server.
    # The obvious thing to do is to start a new server for each value of `draft_folder`,
    # however, shutting down the webserver asynchronously doesn't play well with
    # testsets, and I spent far too long trying to figure out what was going on.
    # This is a simple hack until someone who knows more about this comes along...
    GLOBAL_TEST_DIR = Ref{String}("")
    server = HTTP.Sockets.listen(
        HTTP.Sockets.InetAddr(parse(HTTP.Sockets.IPAddr, "127.0.0.1"), 1234)
    )
    @async HTTP.listen("127.0.0.1", 1234; server = server, verbose = true) do http
        # Make sure to strip first character (`/`) from the target, otherwise it will
        # infer as a file in the root directory.
        file = joinpath(GLOBAL_TEST_DIR[], "../../remotes", http.message.target[2:end])
        HTTP.setstatus(http, 200)
        startwrite(http)
        write(http, read(file, String))
    end
    @testset "$(draft_folder)" for draft_folder in [
        "draft4", "draft6", basename(abspath(LOCAL_TEST_DIR))
    ]
        test_dir = joinpath(SCHEMA_TEST_DIR, draft_folder)
        GLOBAL_TEST_DIR[] = test_dir
        @testset "$(file)" for file in filter(n -> endswith(n, ".json"), readdir(test_dir))
            file_path = joinpath(test_dir, file)
            @testset "$(schema["description"])" for schema in JSON.parsefile(file_path)
                spec = Schema(
                    schema["schema"];
                    parent_dir =
                        schema["schema"] isa Bool ? abspath(".") : dirname(file_path)
                )
                @testset "$(test["description"])" for test in schema["tests"]
                    @test isvalid(test["data"], spec) == test["valid"]
                end
            end
        end
    end
    close(server)
end

@testset "Validate and diagnose" begin
    schema = Schema(Dict(
        "properties" => Dict(
            "foo" => Dict(),
            "bar" => Dict()
        ),
        "required" => ["foo"]
    ))
    data_pass = Dict("foo" => true)
    data_fail = Dict("bar" => 12.5)
    @test JSONSchema.validate(data_pass, schema) === nothing
    ret = JSONSchema.validate(data_fail, schema)
    fail_msg = """Validation failed:
    path:         top-level
    instance:     $(data_fail)
    schema key:   required
    schema value: ["foo"]
    """
    @test ret !== nothing
    @test sprint(show, ret) == fail_msg
    @test diagnose(data_pass, schema) === nothing
    @test diagnose(data_fail, schema) == fail_msg
end

@testset "parentFileDirectory deprecation" begin
    schema = Schema("{}"; parentFileDirectory = ".")
    @test typeof(schema) == Schema
end

@testset "Schemas" begin
    schema = Schema("""{
        \"properties\": {
        \"foo\": {},
        \"bar\": {}
        },
        \"required\": [\"foo\"]
    }""")
    @test typeof(schema) == Schema
    @test typeof(schema.data) == Dict{String,Any}

    schema_2 = Schema(false)
    @test typeof(schema_2) == Schema
    @test typeof(schema_2.data) == Bool
end

@testset "Base.show" begin
    schema = Schema("{}")
    @test sprint(show, schema) == "A JSONSchema"
end

@testset "errors" begin
    @test_throws(
        ErrorException("missing property 'Foo' in $(Dict{String,Any}())."),
        Schema("""{
            "type": "object",
            "properties": {"version": {"\$ref": "#/definitions/Foo"}},
            "definitions": {}
        }""")
    )

    @test_throws(
        ErrorException("unmanaged type in ref resolution $(Int64): 1."),
        Schema("""{
            "type": "object",
            "properties": {"version": {"\$ref": "#/definitions/Foo"}},
            "definitions": 1
        }""")
    )
    @test_throws(
        ErrorException("expected integer array index instead of 'Foo'."),
        Schema("""{
            "type": "object",
            "properties": {"version": {"\$ref": "#/definitions/Foo"}},
            "definitions": [1, 2]
        }""")
    )
    @test_throws(
        ErrorException("item index 3 is larger than array $(Any[1, 2])."),
        Schema("""{
            "type": "object",
            "properties": {"version": {"\$ref": "#/definitions/3"}},
            "definitions": [1, 2]
        }""")
    )
    @test_throws(
        ErrorException("cannot support circular references in schema."),
        validate(
            Dict("version" => 1),
            Schema("""{
                "type": "object",
                "properties": {
                    "version": {
                        "\$ref": "#/definitions/Foo"
                    }
                },
                "definitions": {
                    "Foo": {
                        "\$ref": "#/definitions/Foo"
                    }
                }
            }""")
        )
    )
end

@testset "_is_type" begin
    for (key, val) in Dict(
        :array => [1, 2],
        :boolean => true,
        :integer => 1,
        :number => 1.0,
        :null => nothing,
        :object => Dict(),
        :string => "string",
    )
        @test JSONSchema._is_type(val, Val(Symbol(key)))
        @test !JSONSchema._is_type(:not_a_json_type, Val(Symbol(key)))
    end
    @test !JSONSchema._is_type(true, Val(:number))
    @test !JSONSchema._is_type(true, Val(:integer))
end


text:
##
## L1 regularization
##
@with_kw mutable struct regularizer_group{Tr<:Regularizer} <: Regularizer
    groups::Vector{UnitRange{Int64}}
    regs::Vector{Tr}
end

function <(M1::ProductManifold, M2::ProductManifold)
    return M1.manifolds < M2.manifolds
end

function wholespace_manifold(reg::regularizer_group{Tr}, x) where {Tr}
    return ProductManifold([
        wholespace_manifold(reg.regs[i], x[group]) for (i, group) in enumerate(reg.groups)
    ]...)
end

## 0th order
function g(reg::regularizer_group, x)
    return sum(g(reg.regs[i], x[group]) for (i, group) in enumerate(reg.groups))
end

## 1st order
function prox_αg!(reg::regularizer_group, res, x, α)
    Ms = []
    for (i, group) in enumerate(reg.groups)
        push!(Ms, prox_αg!(reg.regs[i], view(res, group), view(x, group), α))
    end
    return ProductManifold(Ms...)
end

# ∇M_g!(reg::regularizer_groupdistball, ::l1Manifold, ∇M_g, x) = (@. ∇M_g = reg.λ * sign.(x))
# function ∇M_g(reg::regularizer_groupdistball, M::l1Manifold, x)
#     ∇M_g = zeros(representation_size(M))
#     return ∇M_g!(reg, M, ∇M_g, x)
# end

# ## 2nd order
# ∇²M_g_ξ(::regularizer_groupdistball, M::l1Manifold, x, ξ) = zeros(size(ξ))
# ∇²M_g_ξ!(::regularizer_groupdistball, M::l1Manifold, res, x, ξ) = (res .= 0)


text:
# Total ordering on tensors which defines the sweep direction
const TensorOrder = Base.By(λ->(λ.y,λ.x))
isless(A::Tensor, B::Tensor) = lt(TensorOrder, A, B)
sortperm(TN::TensorNetwork) = Base.sortperm(TN, order=TensorOrder)
issorted(TN::TensorNetwork) = Base.issorted(TN, order=TensorOrder)

# Sorts the tensor network according to TensorOrder and updates the adjacency labels
function sort!(TN::TensorNetwork)
    issorted(TN) && return TN
    σ = sortperm(TN)
    permute!(TN, σ)
    τ = invperm(σ)
    for t ∈ TN
        t.adj = τ[t.adj]
    end
    return TN
end

struct edge
    lo::Int
    hi::Int
end

# Returns *twice* the perpendicular distance, as the extra division is unnecessary
perpdist(TN::TensorNetwork, p::Int, e::edge) =
    (TN[e.hi].x-TN[e.lo].x)*(TN[e.hi].y-2*TN[p].y+TN[e.lo].y)-
    (TN[e.hi].y-TN[e.lo].y)*(TN[e.hi].x-2*TN[p].x+TN[e.lo].x)

# Is a point p to the 'left' of an edge e?
toleft(TN::TensorNetwork, p::Int, e::edge) = perpdist(TN,p,e)<0

# Do two edges non-trivially intersect?
intersect(TN::TensorNetwork, a::edge, b::edge) =
    a.lo!=b.lo && a.lo!=b.hi && a.hi!=b.lo && a.hi!=b.hi &&
    perpdist(TN,a.lo,b)*perpdist(TN,a.hi,b)<0 && perpdist(TN,b.lo,a)*perpdist(TN,b.hi,a)<0

# A (partial) order on edges relevant to a Bentley-Ottmann sweep
function edgelt(TN::TensorNetwork, a::edge, b::edge)
    TN[a.lo]>TN[b.lo] && return toleft(TN,a.lo,b)
    TN[a.lo]<TN[b.lo] && return !toleft(TN,b.lo,a)
    θa = atan(TN[a.hi].x-TN[a.lo].x, TN[a.hi].y-TN[a.lo].y)
    θb = atan(TN[b.hi].x-TN[b.lo].x, TN[b.hi].y-TN[b.lo].y)
    return θa<θb
end

# Adds a dim=1 edge to the network
function addedge!(TN::TensorNetwork, u::Int, v::Int)
    u∈TN[v].adj && return
    push!(TN[u].adj, v)
    TN[u].arr = reshape(TN[u].arr, (size(TN[u].arr)...,1))
    push!(TN[v].adj, u)
    TN[v].arr = reshape(TN[v].arr, (size(TN[v].arr)...,1))
    return nothing
end

# Given that an intersection exists, what is its position?
function intersectpos(TN::TensorNetwork, a::edge, b::edge)
    Aa = TN[a.hi].y - TN[a.lo].y
    Ba = TN[a.hi].x - TN[a.lo].x
    Ca = TN[a.hi].y*TN[a.lo].x - TN[a.lo].y*TN[a.hi].x

    Ab = TN[b.hi].y - TN[b.lo].y
    Bb = TN[b.hi].x - TN[b.lo].x
    Cb = TN[b.hi].y*TN[b.lo].x - TN[b.lo].y*TN[b.hi].x

    return (Bb*Ca-Ba*Cb, Ab*Ca-Aa*Cb) ./ (Aa*Bb-Ab*Ba)
end

# Add dim=1 edges corresponding to the left hull of the network. Together with connect! this
# allows the imposition of sweep-connectivity on an arbitrary network
function hull!(TN::TensorNetwork)
    sort!(TN)
    # Left hull
    H = Int[]
    for h ∈ eachindex(TN)
        push!(H, h)
        while length(H)>2 && toleft(TN, H[end], edge(H[end-2], H[end-1]))
            H = H[[1:end-2;end]]
        end
    end
    for i ∈ 2:length(H)
        addedge!(TN, H[i-1], H[i])
    end
    return TN
end

# Adds swap tensors as appropriate to planarise an arbitrary network. Can throw
# InvalidTNError if the graph is not appropriately connected. This is based off a
# Bentley-Ottmann Algorithm approach, such as can be seen in https://youtu.be/be5y0BVQ5kg
function planarise!(TN::TensorNetwork)
    # Check the intersection of neighbouring edges, labelled r-1 and r, in the AVLTree B. If
    # an intersection is found, also updates the network, Q and B
    function checkintersection(r::Int)
        # If there is no intersection then return
        1<r<=length(B) || return nothing
        L = B[r-1]
        R = B[r]
        intersect(TN, L, R) || return nothing
        # Otherwise we need to include a new swap tensor
        len += 1;
        # As we are constantly adding to TN, resize it only infrequently to avoid overheads
        # from memory management
        len>length(TN) && resize!(TN,2*length(TN))
        bL = size(TN[L.lo].arr, findfirst(isequal(L.hi), TN[L.lo].adj))
        bR = size(TN[R.lo].arr, findfirst(isequal(R.hi), TN[R.lo].adj))
        TN[len] = Tensor{Int}(
            [L.lo, R.lo, L.hi, R.hi],
            reshape(Matrix{Float64}(LinearAlgebra.I,bL*bR,bL*bR), (bL,bR,bL,bR)),
            intersectpos(TN, L, R)...
        )
        # Add the new tensor to the queue of points that need to be analysed
        DataStructures.enqueue!(Q, len, (TN[len].y,TN[len].x));
        # Disconnect all of the relevant tensors, connecting them through the swap
        replace!(TN[L.lo].adj, L.hi=>len)
        replace!(TN[L.hi].adj, L.lo=>len)
        replace!(TN[R.lo].adj, R.hi=>len)
        replace!(TN[R.hi].adj, R.lo=>len)
        # Remove the old edges from B and include the new ones
        delete!(B, L)
        delete!(B, R)
        push!(B, edge(L.lo, len))
        push!(B, edge(R.lo, len))
    end

    len = length(TN)

    # Self-balancing tree containing the edges cut by the current sweepline
    B::AVLTree{edge} = AVLTree{edge}(Lt((a,b)->edgelt(TN,a,b)))
    # Queue of points above the sweepline yet to be analysed
    Q = PriorityQueue{Int64,Tuple{Float64,Float64}}()
    # Populate the queue
    for l ∈ eachindex(TN)
        DataStructures.enqueue!(Q, l, (TN[l].y,TN[l].x))
    end
    # Loop over points in the queue until empty
    while !isempty(Q)
        q = DataStructures.dequeue!(Q)
        # Loop over all 'downward' edges of q
        for n ∈ TN[q].adj
            TN[n]>TN[q] && continue
            # For each such edge, remove it, and then check whether the now-neighbouring
            # edges in B are intersecting. If the tree throws an error this is because of a
            # failure of the partial ordering on edges, suggesting overlapping edges
            e = edge(n, q)
            try
                r = sorted_rank(B, e)
                delete!(B, e)
                checkintersection(r)
            catch _
                throw(InvalidTNError("overlapping edges"))
            end
        end
        # Loop over all 'upward' edges of q
        for n ∈ TN[q].adj
            TN[n]<TN[q] && continue
            e = edge(q, n)
            # Push the new edge into the tree, once again checking for any intersections
            # with its neighbours therein
            push!(B, e)
            try
                r = sorted_rank(B, e)
                checkintersection(r)
                checkintersection(r+1)
            catch _
                throw(InvalidTNError("overlapping edges"))
            end
        end
    end
    # Truncate off any unnecessary memory that was included to avoid memory overheads
    resize!(TN, len)
    return TN
end

# Sweep-connects a tensor network. This assumes that all of the edges along the left hull
# are included and that the graph is planar, as are imposed by !hull and !connect
# respectively. As with planarise! this algorithm is also based on a Bentley-Ottmann-style
# sweepline algorithm. sweep. Instead of looking for intersections, this looks for tensors
# with no 'downard' edges (other than the bottommost tensor), and connects them to a tensor
# at a Steiner point immediate to the left.
function connect!(TN::TensorNetwork)
    len = length(TN)

    sort!(TN)

    # Standard Bentley-Ottmann setup as in planarise!
    B::AVLTree{edge} = AVLTree{edge}(Lt((a,b)->edgelt(TN,a,b)))
    Q = PriorityQueue{Int64,Tuple{Float64,Float64}}()
    for l ∈ eachindex(TN)
        DataStructures.enqueue!(Q, l, (TN[l].y,TN[l].x));
    end
    while !isempty(Q)
        q = DataStructures.dequeue!(Q)
        # Consider only points where are not sweep-connected
        steiner = q>1
        for n ∈ TN[q].adj
            TN[n]>TN[q] && continue
            steiner = false
            delete!(B, edge(n, q))
        end
        if steiner
            # Find the edge L which is to the immediate left of the point q
            if isempty(TN[q].adj)
                e = edge(q, q)
                push!(B, e)
                r = sorted_rank(B, e)
                L = B[r-1]
                delete!(B,e)
            else
                e = edge(q, first(TN[q].adj))
                push!(B, e)
                r = sorted_rank(B, e)
                L = B[r-1]
            end
            if TN[L.lo].y == TN[q].y
                # If there is a tensor at the steiner point, directly connect them
                addedge!(FTN, L.lo, q)
            else
                # Otherwise include a new point
                Al = TN[L.hi].y - TN[L.lo].y
                Bl = TN[L.hi].x - TN[L.lo].x
                Cl = TN[L.lo].x*TN[L.hi].y - TN[L.lo].y*TN[L.hi].x
                x = (Cl + Bl * TN[q].y) / Al
                len += 1;
                len>length(TN) && resize!(TN, 2*length(TN))

                b = size(TN[L.lo].arr, findfirst(isequal(L.hi), TN[L.lo].adj))
                TN[len]=Tensor{Int}(
                    [L.lo, L.hi],
                    Matrix{Float64}(LinearAlgebra.I, b, b),
                    x, TN[q].y
                )
                addedge!(TN, len, q)
                replace!(TN[L.lo].adj, L.hi=>len)
                replace!(TN[L.hi].adj, L.lo=>len)
                delete!(B, L)
                push!(B, edge(len,L.hi))
            end
        end
        for n ∈ TN[q].adj
            TN[n]<TN[q] && continue
            push!(B, edge(q, n))
        end
    end
    resize!(TN, len)
    return TN
end


text:
"""
Python PyTorch torchvision wrapper library

MIT License

Copyright (c) 2019 amenoyoya

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

module TorchVision
    export torch, torchvision, transforms, models
    export Color, Train, Valid, Phase, typestr
    export seed_random!, predict, save, load, DataLoader, set_params_to_update!, train!
    export make_transformer, make_transformer_for_vgg16, make_transformer_for_training, make_transformer_for_vgg16_training

    using PyCall, Random, ProgressMeter

    const numpy = pyimport("numpy")
    const torch = pyimport("torch")
    const torchvision = pyimport("torchvision")
    const transforms = torchvision.transforms
    const models = torchvision.models

    # 色表現型
    const Color = Tuple{AbstractFloat,AbstractFloat,AbstractFloat}

    # 機械学習モード表現型
    abstract type Train end
    abstract type Valid end
    const Phase = Union{Type{Train}, Type{Valid}}

    # 型（CamelCase）を文字列（snake_case）に変換する
    typestr(T) = replace(replace(string(T), r"([A-Z])" => c -> "_" * lowercase(c)), r"^_" => s"")

    # 乱数初期化
    ## 値を破壊的に変更する（もしくはグローバル変数の状態に影響を与える）関数には慣例的に`!`をつける
    seed_random!(seed::Int) = begin
        Random.seed!(seed)
        numpy.random.seed(seed)
        torch.manual_seed(seed) # PyTorchの乱数初期化
    end

    # torch.Modelに予測させる
    ## Julia標準のArrayデータを入力として受け取り、Arrayで返す
    predict(model::PyObject, data::Array) = model(torch.tensor(data)).detach().numpy()

    # torch.Modelを保存
    save(model::PyObject, filepath::AbstractString) = torch.save(model.state_dict(), filepath)

    # torch.Modelを読み込み
    load(model::PyObject, filepath::AbstractString) = model.load_state_dict(torch.load(filepath))

    # PIL.Imageオブジェクトに対して変換処理を行う関数を生成
    ## @param resize: transforms.Resizeオブジェクト
    ## @param crop: transforms.Cropオブジェクト
    ## @param normalize: transforms.Normalizeオブジェクト
    ## @return image::PyObject（PIL.Image） -> data::Array{AbstractFloat, 3}（色数×高さ×幅）
    make_transformer(resize::PyObject, crop::PyObject, normalize::PyObject) = begin
        transformer = transforms.Compose([
            resize
            crop
            transforms.ToTensor() # torch.Tensor形式に変換
            normalize
        ])
        return (image::PyObject) -> transformer(image).numpy()
    end

    # PIL.ImageオブジェクトをVGG16ニューラルネットワークモデル用のデータ形式に変換する関数を生成
    # () -> ((img::PyObject) -> Array{Float32,3})
    make_transformer_for_vgg16() = begin
        transformer = make_transformer(
            transforms.Resize(224), # 短い方の辺の長さが224になるようにリサイズ
            transforms.CenterCrop(224), # 画像中央を 224 x 224 で切り取り
            # 色平均化
            ## 色チャンネルの平均値: RGB = (0.485, 0.456, 0.406)
            ## 色チャンネルの標準偏差: RGB = (0.229, 0.224, 0.225)
            transforms.Normalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225))
        )
        # 値を0〜1の範囲内に収めて画像データ配列を返す関数を返す
        ## clamp(x, low, high): if x < low then x = low, if x > high then x = high
        ## `.`演算子: 関数を行列に対して適用させる
        return (image::PyObject) -> clamp.(transformer(image), 0, 1)
    end

    # モデル訓練用に入力画像を変換する関数を生成
    # (resize, std, mean) -> ((img::PyObject, phase::Phase) -> Array{Float32,3})
    ## @param resize::Int = リサイズする大きさ
    ## @param std::Color = 各色チャンネルの平均値
    ## @param mean::Color = 各色チャンネルの標準偏差
    make_transformer_for_training(resize::Int, std::Color, mean::Color) = begin
        trasformers = Dict(
            Train => make_transformer(
                transforms.RandomResizedCrop(resize; scale=(0.5, 1.0)),
                transforms.RandomHorizontalFlip(),
                transforms.Normalize(mean, std)
            ),
            Valid => make_transformer(
                transforms.Resize(resize),
                transforms.CenterCrop(resize),
                transforms.Normalize(mean, std)
            )
        )
        return (img::PyObject, phase::Phase) -> trasformers[phase](img)
    end

    # VGG-16モデル訓練用に入力画像を変換する関数を生成
    # (resize, std, mean) -> ((img::PyObject, phase::Phase) -> Array{Float32,3})
    make_transformer_for_vgg16_training() = begin
        transformer = make_transformer_for_training(224, (0.485, 0.456, 0.406), (0.229, 0.224, 0.225))
        return (img::PyObject, phase::Phase) -> clamp.(transformer(img, phase), 0, 1)
    end

    # PyTorch規格の画像データセット型生成マクロ
    # 要 PyCall, ./Image
    # @param TypeName: 型名
    # @param make_dataset_list_function(phase::Phase)::Array{String,1}: データセットのファイルパスリストを返す関数
    # @param image_transform_function(img::PyObject, phase::Phase)::Array{Float32,3}: 画像をモデル入力用に変換する関数
    # @param labeling_function(img_path::String, phase::Phase)::Int: データセットのラベルインデックスを返す関数
    ## esc(SymbolExpression) を使うことで マクロ実行位置にそのままコードを埋め込むことができる
    ## ※ そのまま埋め込まれるため変数スコープ等は効かない
    macro image_dataset(TypeName, make_dataset_list_function, image_transform_function, labeling_function)
        esc(quote
            @pydef mutable struct $TypeName <: TorchVision.torch.utils.data.Dataset
                __init__(self, phase::TorchVision.Phase) = begin
                    pybuiltin(:super)($TypeName, self).__init__()
                    self.phase = phase
                    self.file_list = $make_dataset_list_function(phase)
                end

                __len__(self) = length(self.file_list)

                __getitem__(self, index::Int) = begin
                    # index番目の画像をロード
                    ## Juliaのindexは1〜なので +1 する
                    img_path = self.file_list[index + 1]
                    img = Image.open(img_path).convert("RGB") # グレースケール画像は強制的にRGB画像に変換
                    img_transformed = $image_transform_function(img, self.phase)
                    # ラベリング
                    label = $labeling_function(img_path, self.phase)
                    return img_transformed, label
                end
            end
        end)
    end

    # データローダーオブジェクトの作成
    # @param Dataset::PyObject = @TensorVision.image_dataset マクロで作成したデータセット型
    # @param batch_size::Int = (default: 8) ミニバッチサイズ
    # @param shuffle::Bool = (default: false) trueなら学習時にデータをシャッフルする
    # @return Dict{DataType(Phase), PyObject(DataLoader)}
    DataLoader(Dataset::PyObject; batch_size::Int=8, shuffle::Bool=false) = Dict(
        Train => torch.utils.data.DataLoader(
            Dataset(Train); batch_size=batch_size, shuffle=shuffle
        ),
        Valid => torch.utils.data.DataLoader(
            Dataset(Valid); batch_size=batch_size, shuffle=shuffle
        )
    )

    # 転移学習で学習させるパラメータを設定し、各パラメータの現在値を取得
    # @param model::PyObject = torch.Model
    # @param param_names::Array{String,1} = 学習させるパラメータ名
    # @return Array{PyObject,1} = [torch.Tensor]: アップデートされるパラメータの配列
    set_params_to_update!(model::PyObject, param_names::Array{String,1})::Array{PyObject,1} = begin
        params_to_update::Array{PyObject,1} = []
        # 学習させるパラメータ以外は勾配計算させない
        for (name, param) = model.named_parameters()
            if in(name, param_names)
                param.required_grad = true
                push!(params_to_update, param)
            else
                param.required_grad = false
            end
        end
        return params_to_update
    end

    # 学習実行
    # @param model::PyObject = torch.Model
    # @param phase::Phase = Train: 訓練モード | Valid: 検証モード
    # @param dataloader::Dict{DataType,PyObject} = DataLoader関数で生成したデータローダーオブジェクト
    # @param optimizer::PyObject = torch.Optimizer: 最適化関数オブジェクト
    # @param criterion::PyObject = torch.LossFunction: 損失関数オブジェクト
    # @return 損失和::AbstractFloat, 正解数::Int
    train!(model::PyObject, phase::Phase, dataloader::Dict{DataType,PyObject}, optimizer::PyObject, criterion::PyObject) = begin
        epoch_loss = 0.0 # epochの損失和
        epoch_corrects = 0 # epochの正解数
        # プログレスバー
        progress = Progress(pybuiltin("len")(dataloader[phase]))
        for (inputs, labels) in pybuiltin("iter")(dataloader[phase])
            # optimizer初期化
            optimizer.zero_grad()

            # 順伝搬計算
            torch.set_grad_enabled(phase == Train)
            outputs = model(inputs)
            loss = criterion(outputs, labels) # 損失計算
            (max, preds) = torch.max(outputs, 1) # ラベルを予測

            # 訓練時はバックプロパゲーション
            if phase == Train
                loss.backward()
                optimizer.step()
            end

            # イテレーション結果の計算
            epoch_loss += loss.item() * inputs.size(0)
            epoch_corrects += torch.sum(preds == labels.data)
            torch.set_grad_enabled(false)

            next!(progress)
        end
        return epoch_loss, epoch_corrects.numpy()
    end
end

max_stars_repo_name: amenoyoya/pytorch-learning

max_issues_repo_name: amenoyoya/pytorch-learning

max_forks_repo_name: amenoyoya/pytorch-learning


text:
@testset "cheby1 correct values origin" begin
    # legendre == jacobi(a,b=0,0); same function body
    for n in 1:10
        if !iseven(n)
            @test cheby1(n, [0.]) ≈ [0.]
        end
    end
end

@testset "cheby2 correct values origin" begin
    # legendre == jacobi(a,b=0,0); same function body
    for n in 1:10
        if !iseven(n)
            @test cheby2(n, [0.]) ≈ [0.]
        end
    end
end

@testset "cheby1 sum correct values origin" begin
    terms = [1,3,5,7,9]
    weights = ones(size(terms)...)
    @test cheby1_sum(terms, weights, [0.]) ≈ [0.]
end

@testset "cheby2 sum correct values origin" begin
    terms = [1,3,5,7,9]
    weights = ones(size(terms)...)
    @test cheby2_sum(terms, weights, [0.]) ≈ [0.]
end

@testset "cheby1 series correct value origin" begin
    terms = [1,3,5]
    out = [0., 0., 0.,]
    series = cheby1_series(terms, [0.,])
    @test series ≈ out'
end

@testset "cheby2 series correct value origin" begin
    terms = [1,3,5]
    out = [0., 0., 0.,]
    series = cheby2_series(terms, [0.,])
    @test series ≈ out'
end


text:
@testset "layer printing" begin # 2-arg show, defined with layes

  @test repr(Dense(2,3)) == "Dense(2 => 3)"
  @test repr(Chain(Dense(2,3))) == "Chain(Dense(2 => 3))"
  @test repr(Chain(lay=Dense(2,3))) == "Chain(lay = Dense(2 => 3))"

end
@testset "nested model printing" begin # 3-arg show, defined in show.jl

  # Dense -- has parameter count, but not when inside a matrix:

  toplevel_dense = repr("text/plain", Dense(2,3))
  @test occursin("Dense(2 => 3)", toplevel_dense)
  @test occursin("# 9 parameters", toplevel_dense)

  @test Meta.isexpr(Meta.parse(toplevel_dense), :call)  # comment is ignored

  vector_dense = repr("text/plain", [Dense(2,3), Dense(2,3)])
  @test occursin("Dense(2 => 3)", vector_dense)
  @test occursin("# 9 parameters", vector_dense)

  matrix_dense = repr("text/plain", fill(Dense(2,3), 3, 3))
  @test occursin("Dense(2 => 3)", matrix_dense)
  @test !occursin("# 9 parameters", matrix_dense)

  tuple_dense = repr("text/plain", tuple(Dense(2,3)))
  @test occursin("Dense(2 => 3)", tuple_dense)
  @test !occursin("# 9 parameters", tuple_dense)

  # Chain -- gets split over lines at top level only

  toplevel_chain = repr("text/plain", Chain(Dense(2,3)))
  @test occursin("Chain(\n  Dense(2 => 3)", toplevel_chain)
  @test occursin("# 9 parameters", toplevel_chain)
  @test !occursin("# Total:", toplevel_chain)

  vector_chain = repr("text/plain", [Chain(Dense(2,3)), Chain(Dense(2,3))])
  @test occursin("Chain(Dense(2 => 3))", vector_chain)
  @test occursin("# 9 parameters", vector_chain)
  @test !occursin("# Total:", vector_chain)

  matrix_chain = repr("text/plain", fill(Chain(Dense(2,3)), 3,3))
  @test occursin("Chain(Dense(2 => 3))", matrix_chain)
  @test !occursin("# 9 parameters", matrix_chain)
  @test !occursin("# Total:", matrix_chain)

  # ... and only long enough chains get a total at the end:

  longchain = Chain(Dense(2 => 3), Dense(3 => 4), Dense(4 => 5), softmax)

  toplevel_longchain = repr("text/plain", longchain)
  @test occursin("Chain(\n  Dense(2 => 3)", toplevel_longchain)
  @test occursin("# 9 parameters", toplevel_longchain)
  @test occursin("# Total: 6 arrays, 50 parameters", toplevel_longchain)

  vector_longchain = repr("text/plain", [longchain, longchain]) # pretty ugly in reality
  @test occursin("Chain(Dense(2 => 3)", vector_longchain)
  @test occursin("# 50 parameters", vector_longchain)
  @test !occursin("# 9 parameters", vector_longchain)
  @test !occursin("# Total:", vector_longchain)

  matrix_longchain = repr("text/plain", fill(longchain, 3,3))
  @test occursin("Chain(Dense(2 => 3)", matrix_longchain)
  @test !occursin("# 9 parameters", matrix_longchain)
  @test !occursin("# Total:", matrix_longchain)

  @test Meta.isexpr(Meta.parse(toplevel_longchain), :call)  # comments are ignored
  @test Meta.parse(toplevel_longchain).args[1] == :Chain

end


text:
# 2D FWI on Overthrust model using minConf library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
#

using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO, PyPlot

# Load starting model
n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ model0.m)
vmin = ones(Float32,model0.n) .* 1.3f0
vmax = ones(Float32,model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("../../data/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0)  # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)

############################### FWI ###########################################
F0 = judiModeling(deepcopy(model0), src_geometry, d_obs.geometry)

# Optimization parameters
niterations = 10
batchsize = 16
fhistory_SGD = zeros(Float32,niterations)

# Projection operator for bound constraints
proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2),model0.n)
ls = BackTracking(order=3, iterations=10, )

# Main loop
for j=1:niterations

    # get fwi objective function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0,q[i],d_obs[i])
    p = -gradient/norm(gradient, Inf)

    println("FWI iteration no: ",j,"; function value: ",fval)
    fhistory_SGD[j] = fval

    # linesearch
    function ϕ(α)
        F0.model.m .= proj(model0.m .+ α * p)
        misfit = .5*norm(F0[i]*q[i] - d_obs[i])^2
        return misfit
    end
    step, fval = ls(ϕ, 1f0, fval, dot(gradient, p))

    # Update model and bound projection
    model0.m .= proj(model0.m .+ step .* p)
end

figure(); imshow(sqrt.(1f0./adjoint(model0.m))); title("FWI with SGD")

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using Wavelets

wf = waveletfilter("db6")
N = 1024*32;
x0 = rand(N);
L = int(log2(N)-2)
tn = 1000
fw=true

ft=Float64
si = Array(ft, wf.n-1)       # tmp filter vector
if fw
    dcfilter = mirror(convert(Vector{ft},wf.qmf))  #mqmf
    scfilter = reverse(convert(Vector{ft},wf.qmf))  #rqmf
else
    scfilter = convert(Vector{ft},wf.qmf)  #qmf
    dcfilter = reverse(mirror(convert(Vector{ft},wf.qmf)))  #mrqmf
end
snew = Array(ft, length(x0)>>1)

y=x0*0.0

println("dwt! (N=",N,"), ", L, " levels")
#f(x0, L, wf) = for i = 1:tn; fwt(x0, L, wf); end
f(x0, L, wf) = for i = 1:tn; dwt!(y, x0, wf, L, fw, dcfilter, scfilter, si, snew); end
f(x0, L, wf);
@time f(x0, L, wf);