text:
"""
    HMCSampler

A `AbstractMCMC.AbstractSampler` for kernels in AdvancedHMC.jl.

# Fields

$(FIELDS)

# Notes

Note that all the fields have the prefix `initial_` to indicate
that these will not necessarily correspond to the `kernel`, `metric`,
and `adaptor` after sampling.

To access the updated fields use the resulting [`HMCState`](@ref).
"""
struct HMCSampler{K, M, A} <: AbstractMCMC.AbstractSampler
    "Initial [`AbstractMCMCKernel`](@ref)."
    initial_kernel::K
    "Initial [`AbstractMetric`](@ref)."
    initial_metric::M
    "Initial [`AbstractAdaptor`](@ref)."
    initial_adaptor::A
end
HMCSampler(kernel, metric) = HMCSampler(kernel, metric, Adaptation.NoAdaptation())

"""
    DifferentiableDensityModel(ℓπ, ∂ℓπ∂θ)
    DifferentiableDensityModel(ℓπ, m::Module)

A `AbstractMCMC.AbstractMCMCModel` representing a differentiable log-density.

If a module `m` is given as the second argument, then `m` is assumed to be an
automatic-differentiation package and this will be used to compute the gradients.

Note that the module `m` must be imported before usage, e.g.
```julia
using Zygote: Zygote
model = DifferentiableDensityModel(ℓπ, Zygote)
```
results in a `model` which will use Zygote.jl as its AD-backend.

# Fields
$(FIELDS)
"""
struct DifferentiableDensityModel{Tlogπ, T∂logπ∂θ} <: AbstractMCMC.AbstractModel
    "Log-density. Maps `AbstractArray` to value of the log-density."
    ℓπ::Tlogπ
    "Gradient of log-density. Returns a tuple of `ℓπ` and the gradient evaluated at the given point."
    ∂ℓπ∂θ::T∂logπ∂θ
end

struct DummyMetric <: AbstractMetric end
function DifferentiableDensityModel(ℓπ, m::Module)
    h = Hamiltonian(DummyMetric(), ℓπ, m)
    return DifferentiableDensityModel(h.ℓπ, h.∂ℓπ∂θ)
end

"""
    HMCState

Represents the state of a [`HMCSampler`](@ref).

# Fields

$(FIELDS)

"""
struct HMCState{
    TTrans<:Transition,
    TMetric<:AbstractMetric,
    TKernel<:AbstractMCMCKernel,
    TAdapt<:Adaptation.AbstractAdaptor
}
    "Index of current iteration."
    i::Int
    "Current [`Transition`](@ref)."
    transition::TTrans
    "Current [`AbstractMetric`](@ref), possibly adapted."
    metric::TMetric
    "Current [`AbstractMCMCKernel`](@ref)."
    κ::TKernel
    "Current [`AbstractAdaptor`](@ref)."
    adaptor::TAdapt
end

"""
    $(TYPEDSIGNATURES)

A convenient wrapper around `AbstractMCMC.sample` avoiding explicit construction of [`HMCSampler`](@ref).
"""
function AbstractMCMC.sample(
    model::DifferentiableDensityModel,
    kernel::AbstractMCMCKernel,
    metric::AbstractMetric,
    adaptor::AbstractAdaptor,
    N::Integer;
    kwargs...
)
    return AbstractMCMC.sample(Random.GLOBAL_RNG, model, kernel, metric, adaptor, N; kwargs...)
end

function AbstractMCMC.sample(
    rng::Random.AbstractRNG,
    model::DifferentiableDensityModel,
    kernel::AbstractMCMCKernel,
    metric::AbstractMetric,
    adaptor::AbstractAdaptor,
    N::Integer;
    progress = true,
    verbose = false,
    callback = nothing,
    kwargs...
)
    sampler = HMCSampler(kernel, metric, adaptor)
    if callback === nothing
        callback = HMCProgressCallback(N, progress = progress, verbose = verbose)
        progress = false # don't use AMCMC's progress-funtionality
    end

    return AbstractMCMC.mcmcsample(
        rng, model, sampler, N;
        progress = progress,
        verbose = verbose,
        callback = callback,
        kwargs...
    )
end

function AbstractMCMC.sample(
    model::DifferentiableDensityModel,
    kernel::AbstractMCMCKernel,
    metric::AbstractMetric,
    adaptor::AbstractAdaptor,
    parallel::AbstractMCMC.AbstractMCMCEnsemble,
    N::Integer,
    nchains::Integer;
    kwargs...
)
    return AbstractMCMC.sample(
        Random.GLOBAL_RNG, model, kernel, metric, adaptor, N, nchains;
        kwargs...
    )
end

function AbstractMCMC.sample(
    rng::Random.AbstractRNG,
    model::DifferentiableDensityModel,
    kernel::AbstractMCMCKernel,
    metric::AbstractMetric,
    adaptor::AbstractAdaptor,
    parallel::AbstractMCMC.AbstractMCMCEnsemble,
    N::Integer,
    nchains::Integer;
    progress = true,
    verbose = false,
    callback = nothing,
    kwargs...
)
    sampler = HMCSampler(kernel, metric, adaptor)
    if callback === nothing
        callback = HMCProgressCallback(N, progress = progress, verbose = verbose)
        progress = false # don't use AMCMC's progress-funtionality
    end

    return AbstractMCMC.mcmcsample(
        rng, model, sampler, parallel, N, nchains;
        progress = progress,
        verbose = verbose,
        callback = callback,
        kwargs...
    )
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::DifferentiableDensityModel,
    spl::HMCSampler;
    init_params = nothing,
    kwargs...
)
    metric = spl.initial_metric
    κ = spl.initial_kernel
    adaptor = spl.initial_adaptor

    if init_params === nothing
        init_params = randn(rng, size(metric, 1))
    end

    # Construct the hamiltonian using the initial metric
    hamiltonian = Hamiltonian(metric, model.ℓπ, model.∂ℓπ∂θ)

    # Get an initial sample.
    h, t = AdvancedHMC.sample_init(rng, hamiltonian, init_params)

    # Compute next transition and state.
    state = HMCState(0, t, h.metric, κ, adaptor)

    # Take actual first step.
    return AbstractMCMC.step(rng, model, spl, state; kwargs...)
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::DifferentiableDensityModel,
    spl::HMCSampler,
    state::HMCState;
    nadapts::Int = 0,
    kwargs...
)
    # Get step size
    @debug "current ϵ" getstepsize(spl, state)

    # Compute transition.
    i = state.i + 1
    t_old = state.transition
    adaptor = state.adaptor
    κ = state.κ
    metric = state.metric

    # Reconstruct hamiltonian.
    h = Hamiltonian(metric, model.ℓπ, model.∂ℓπ∂θ)

    # Make new transition.
    t = transition(rng, h, κ, t_old.z)

    # Adapt h and spl.
    tstat = stat(t)
    h, κ, isadapted = adapt!(h, κ, adaptor, i, nadapts, t.z.θ, tstat.acceptance_rate)
    tstat = merge(tstat, (is_adapt=isadapted,))

    # Compute next transition and state.
    newstate = HMCState(i, t, h.metric, κ, adaptor)

    # Return `Transition` with additional stats added.
    return Transition(t.z, tstat), newstate
end

################
### Callback ###
################
"""
    HMCProgressCallback

A callback to be used with AbstractMCMC.jl's interface, replicating the
logging behavior of the non-AbstractMCMC [`sample`](@ref).

# Fields
$(FIELDS)
"""
struct HMCProgressCallback{P}
    "`Progress` meter from ProgressMeters.jl."
    pm::P
    "Specifies whether or not to use display a progress bar."
    progress::Bool
    "If `progress` is not specified and this is `true` some information will be logged upon completion of adaptation."
    verbose::Bool
end

function HMCProgressCallback(n_samples; progress=true, verbose=false)
    pm = progress ? ProgressMeter.Progress(n_samples, desc="Sampling", barlen=31) : nothing
    HMCProgressCallback(pm, progress, verbose)
end

function (cb::HMCProgressCallback)(
    rng, model, spl, t, state, i;
    nadapts = 0,
    kwargs...
)
    progress = cb.progress
    verbose = cb.verbose
    pm = cb.pm

    metric = state.metric
    adaptor = state.adaptor
    κ = state.κ
    tstat = t.stat
    isadapted = tstat.is_adapt

    # Update progress meter
    if progress
        # Do include current iteration and mass matrix
        pm_next!(
            pm,
            (iterations=i, tstat..., mass_matrix=metric)
        )
        # Report finish of adapation
    elseif verbose && isadapted && i == nadapts
        @info "Finished $nadapts adapation steps" adaptor κ.τ.integrator metric
    end
end


text:
function selection!(context::NamedTuple,selections::NamedTuple,names::NamedTuple)
    request = context.request
	try
		response = HTTP.Response(200)
        label = JSON.parse(String(request.body))
        t = @elapsed selection!(label,selections,names)

		@info """$(request.method) $(replace(request.target, r"\?.+" => "")) | $t seconds"""
		return response

	catch exception
		printstyled("ERROR $(replace(request.target, r"\?.+" => "")) | $exception\n",color=:red)
		return HTTP.Response(500)
	end
end

function selection!(label::Dict,selections::NamedTuple,names::NamedTuple)

    if label["name"] ∈ selections.names.keys
	    selections.names[label["name"]] = label["selected"]

    elseif label["name"] == "populations"
        for name ∈ names.populations selections.names[name]=label["selected"] end

    elseif label["name"] == "conditions"
        for name ∈ names.conditions selections.names[name]=label["selected"] end

    elseif label["name"] == "groups"
        for name ∈ names.groups selections.names[name]=label["selected"] end

    elseif typeof(label["name"]) <: Dict
		for key ∈ filter(key->key≠:channels,keys(names))

			for _ ∈ 1:length(names[key]) pop!(names[key]) end
			append!(names[key],label["name"][String(key)])
		end
    end

    selectionCodes = selections.codes .& encode(selections.names.vals)
	rows = @. selectionCodes ≠ 0
	for key ∈ filter( key -> (key≠:channels) & ~isempty(names[key]), keys(names) )

		groupCodes = encode(map( x -> x ∈ names[key], selections.names.keys ))
		@. rows &= ( selectionCodes & groupCodes ≠ 0 )
	end

    @. selections.rows = rows
	return selections.rows
end

function encode( input; T::DataType=Int )
	v,encodedInt = 1,zero(T)
	iter = Iterators.reverse(input)

	for i ∈ iter
		encodedInt += v * i
		v <<= 1
	end
	return encodedInt
end

function decode( code::Int, codeMap::Vector{<:Any} )
	return codeMap[[ bit=='1' for bit ∈ bitstring(code)[end-length(codeMap)+1:end] ]]
end

ext: jl

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
module Day2

export parse_input
export part_1
export part_2

function parse_input(input::String)
	ids = split(input)
	return ids
end

function check_id(id::AbstractString)
	has_two = false
	has_three = false
	for chr in Set{Char}(id)
		if has_two & has_three
			break
		end
		n = 0
		for other in id
			if chr == other
				n += 1
			end
		end
		has_two = (n == 2) | has_two
		has_three = (n == 3) | has_three
	end

	return has_two, has_three
end

function part_1(ids::Array{<:AbstractString})
	n_two = 0
	n_three = 0
	for id in ids
		has_two, has_three = check_id(id)
		n_two += has_two
		n_three += has_three
	end

	return n_two * n_three
end

function is_similar(id1::AbstractString, id2::AbstractString)
	n_different = 0
	for (a, b) in zip(id1, id2)
		if a != b
			n_different += 1
		end
		if n_different > 1
			return false
		end
	end
	return return true
end

function get_common(id1::AbstractString, id2::AbstractString)
	common = ""
	for (a, b) in zip(id1, id2)
		if a == b
			common *= a
		end
	end
	return common
end

function part_2(ids::Array{<:AbstractString})
	ids_to_check = copy(ids)
	while length(ids_to_check) > 0
		id_to_check = pop!(ids_to_check)
		for other in ids_to_check
			if is_similar(id_to_check, other)
				return get_common(id_to_check, other)
			end
		end
	end
end

end


text:
# MIT License

# Copyright (c) 2021 Fran Urbano
# Copyright (c) 2021 Val Lyashov

# 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.

using Test
using Forecast
import Forecast: hmaSymmetricWeights

@testset "hma" begin

    # weights
    w = hmaSymmetricWeights(8)
    @test ismissing(w[end])

    w = hmaSymmetricWeights(7)
    @test sum(w) ≈ 1.0

    w = hmaSymmetricWeights(13)
    @test sum(w) ≈ 1.0
    # (https://www.mathworks.com/help/econ/seasonal-adjustment-using-snxd7m-seasonal-filters.html)
    @test length(w) == 13
    @test round.(w,digits=3) ≈ [-0.019, -0.028, 0.0, 0.065, 0.147, 0.214, 0.24,
                                 0.214, 0.147, 0.065, 0.0, -0.028, -0.019]

    # 13-term moving average
    x = hma(sin.(1:100), 13)
    @test length(x) == 100
    @test first(x)  ≈ 0.5636212576875559
    @test last(x)   ≈ -0.6658500507547161
    @test x[50]     ≈ -0.043655947941143455

end


text:
using .Fields: F2GNB, F2PB, FP, Field, PrimeField, BinaryField, tobits, value, reducer
using .Curves: AbstractPoint, ECPoint, AffinePoint, Weierstrass, BinaryCurve, gx, gy, field, eq
using .Specs: MODP, Koblitz, ECP, EC2N, Spec, PB, GNB

specialize(::Type{FP}, p::Integer) = FP{static(p)}

specialize(::Type{BinaryCurve}, a::BitVector, b::BitVector) = BinaryCurve{StaticBitVector(a), StaticBitVector(b)}
specialize(::Type{BinaryCurve}, a::F, b::F) where F <: BinaryField = specialize(BinaryCurve, tobits(a), tobits(b))

static_int(x::Integer) = x
static_int(x::BigInt) = StaticBigInt(x)

specialize(::Type{Weierstrass}, a::Integer, b::Integer) = Weierstrass{static_int(a), static_int(b)}

specialize(::Type{Weierstrass}, a::BitVector, b::BitVector) = Weierstrass{StaticBitVector(a), StaticBitVector(b)}
specialize(::Type{Weierstrass}, a::F, b::F) where F <: BinaryField = specialize(Weierstrass, tobits(a), tobits(b))

function specialize(::Type{ECPoint{P}}, curve::Spec; name=nothing) where P <: AbstractPoint

    Q = specialize(P, curve)

    _order = order(curve)
    _cofactor = 1 # Need to update this one

    R = ECPoint{Q}(_order, _cofactor; name)

    return R
end

specialize(::Type{ECPoint}, spec::Spec; name = nothing) = specialize(ECPoint{AffinePoint}, spec; name)

##################### Macro for defining curve as also a group ###################

# macro def(constant_name, type, group_spec)

#     name_str = string(constant_name)

#     M = @__MODULE__

#     return esc(quote
#         const $constant_name = $M.specialize($type, $group_spec)

#         $M.name(::Type{$constant_name}) = $name_str

#         local ORDER = $group_spec.n
#         $M.order(::Type{$constant_name}) = ORDER

#         local GENERATOR = $constant_name($M.generator($group_spec)...)
#         $M.generator(::Type{$constant_name}) = GENERATOR
#     end)
# end

macro def(constant_name, type, group_spec)

    M = @__MODULE__

    return esc(quote
        const $constant_name = $M.specialize($type, $group_spec; name = $(QuoteNode(constant_name)))
    end)
end

function specialize(::Type{F2GNB}, N::Int)

    @assert div(N, 8) != 0 "Out of X9.62 spec"
    T = gn_basis_representation_rule(m)

    return F2GNB{N, T}
end

function F2GNB(x::BitVector)

    N = length(x)
    F = specialize(F2GNB, N)

    return F(x)
end

function specialize(::Type{AffinePoint{Weierstrass, F}}, curve::ECP) where F <: PrimeField

    (; p, a, b) = curve

    P = AffinePoint{specialize(Weierstrass, a, b), specialize(F, p)}

    return P
end

specialize(::Type{AffinePoint}, spec::ECP) = specialize(AffinePoint{Weierstrass, FP}, spec)

function specialize(::Type{F}, basis::PB) where F <: BinaryField
    (; f) = basis
    return F(f)
end

function specialize(::Type{F}, basis::GNB) where F <: BinaryField
    (; m, T) = basis
    return F{m, T}
end

specialize(::Type{F}, basis::Integer) where F <: FP = F{static(basis)}

specialize(::Type{BinaryCurve}, curve::EC2N) = specialize(BinaryCurve, a(curve), b(curve))

function specialize(::Type{AffinePoint{BinaryCurve, F}}, curve) where F <: BinaryField
    P = AffinePoint{specialize(BinaryCurve, curve), specialize(F, curve.basis)}
    return P
end

specialize(::Type{AffinePoint}, curve::EC2N{GNB}) = specialize(AffinePoint{BinaryCurve, F2GNB}, curve)
specialize(::Type{AffinePoint}, curve::EC2N{PB}) = specialize(AffinePoint{BinaryCurve, F2PB}, curve)

specialize(::Type{AffinePoint{BinaryCurve, F}}, curve::Koblitz) where F <: BinaryField = specialize(AffinePoint{BinaryCurve, F}, curve.bec)

function spec(::Type{P}; kwargs...) where P <: AbstractPoint

    n = hasmethod(order, Tuple{Type{P}}) ? order(P) : nothing
    h = hasmethod(cofactor, Tuple{Type{P}}) ? cofactor(P) : nothing

    EQ = eq(P)
    F = field(P)

    return spec(EQ, F; n, h, kwargs...)
end

spec(p::P) where P <: AbstractPoint = spec(P; Gx=value(gx(p)), Gy=value(gy(p)))

function spec(::Type{EQ}, ::Type{F}; n=nothing, h=nothing, Gx=nothing, Gy=nothing) where {EQ <: Weierstrass, F <: PrimeField}

    _a = a(EQ)
    _b = b(EQ)

    p = modulus(F)

    return ECP(p, n, _a, _b, Gx, Gy) # I will need to add H in the end here
end

spec(::Type{F}) where F <: F2PB = PB(reducer(F))
spec(::Type{F2GNB{N, T}}) where {N, T} = GNB(N, T)

function spec(::Type{EQ}, ::Type{F}; n=nothing, h=nothing, Gx=nothing, Gy=nothing) where {EQ <: BinaryCurve, F <: BinaryField}

    _a = convert(BitVector, a(EQ))
    _b = convert(BitVector, b(EQ))  ### Perhaps I would be better to do conversion at acessor methods!

    basis = spec(F)

    #p = modulus(F)
    if isnothing(Gx) && isnothing(Gy)
        return EC2N(basis, n, _a, _b) # I will need to add H in the end here
    else
        return EC2N(basis, n, _a, _b, Gx, Gy)
    end
end

specialize(::Type{ECGroup{P}}, spec::Spec; name = nothing) where P <: ECPoint = ECGroup{specialize(P, spec; name)}
specialize(::Type{ECGroup}, spec::Spec; name = nothing) = ECGroup{specialize(ECPoint, spec; name)}

specialize(::Type{PGroup}, spec::MODP; name = nothing) = PGroup(spec.p, spec.q; name)

function spec(x::Symbol)
    if x == :P_192
        return Specs.Curve_P_192
    elseif x == :P_244
        return Specs.Curve_P_244
    elseif x == :P_256
        return Specs.Curve_P_256
    elseif x == :P_384
        return Specs.Curve_P_384
    elseif x == :P_521
        return Specs.Curve_P_521
    else
        error("$x not implemented")
    end
end

specialize(::Type{PGroup}, p, q) = PGroup(p, q)

spec(::Type{ECGroup{P}}) where P = spec(P)
spec(g::ECGroup) = spec(g.x)

spec(::Type{G}) where G <: PGroup = MODP(; p = modulus(G), q = order(G))


text:
export qmr

"""
 Quasi Minimal Residual (QMR) method applied to the complex symmetric system of
 linear equations Ax=b. The algorithm employed here was modified from Algorithm 8.1
 in the paper:
    Freund R W and Nachtigal N M. An implementation of the QMR method based on
      coupled two-term recurrences. SIAM J Sci., Comput., 1994, 15 (2): 313-337.

 - by Bo Han, May 2016, hanbo1735@163.com, China University of Geosciences.

Input:
    A       :: Function   - computing A*x.
    b       :: Vector     - right hand side vector.
    tol     :: Real       - error tolerance.
    maxIter :: Int        - maximum number of iterations.
    M       :: Function   - preconditioners, either matrices or functions computing M\\x.
    x0      :: Vector     - starting guess.
    out     :: Int        - flag for output (-1: no output, 0: only errors,
                                              1: final status, 2: relres at each iteration).

Output:
    x       :: Vector     - solution.
    flag    :: Int        - exit flag ( 0: desired tolerance achieved
                                       -1: maxIter reached without converging
                                       -2: delta equal to zero
                                       -3: epsilon equal to zero
                                       -9: right hand side was zero ).
    err     :: Real       - residual, i.e., norm(A*x-b)/norm(b).
    iter    :: Int        - number of iterations.
    resvec  :: Vector     - residual at each iteration.

"""
function qmr(A::SparseMatrixCSC{T,Int},b::Array{T,1}; kwargs...) where {T}
	Ax = zeros(T,size(A,2))                  # pre-allocate
	return qmr(x -> mul!(Ax,A,x,1.0,0.0),b;kwargs...)
end

qmr(A,b::Vector;kwargs...) = qmr(x -> A*x ,b;kwargs...)

function qmr(A::Function, b::Vector; tol::Real=1e-6, maxIter::Int=100, M=x->copy(x),
             x::Vector=[],out::Int=0)
	println("WARNING: QMR: this (new) funtion does not seem to work very well at the moment. Need to check algorithm correctness. ")
    n  = length(b)
    if norm(b)==0; return zeros(eltype(b),n),-9,0,0,zeros(1); end
    Mf = isa(M,Function) ? M : x -> M\x

    # allocate vector: p,d,s
    p = zeros(eltype(b),n)
    d = zeros(eltype(b),n)
    s = zeros(eltype(b),n)

    # specify scalar
    c_n      =  one(eltype(b))
    eps      =  one(eltype(b))
    theta_n  =  zero(eltype(b))
    eta      =  -one(eltype(b))

    if isempty(x)
        x = zeros(eltype(b),n)
        r  = copy(b)
    else
		r  = b - A(x)
    end

	if eltype(b) <: Complex
       x = complex(x)
    end

    resvec = zeros(maxIter+1)
    bnrm2 = norm( b )
    err   = norm( r ) / bnrm2; resvec[1] = err

    rho_np = norm( r )
    v = r / rho_np
    # v = BLAS.scal(n,one(eltype(r))/rho_np,r,1)

    if out==2
        println("=== qmr ===")
        println(@sprintf("%4s\t%7s","iter","relres"))
    end

    iter = 1
    flag = -1

    # start the iteration
    while iter <= maxIter
        c_nm     = c_n
        rho_n    = rho_np
        theta_nm = theta_n

        mv = Mf(v)
        #delta = At_mul_B(v, mv)    # delta = v.' * mv
        delta = transpose(v) * mv
        abs(delta) < 1e-40 && (flag=-2; break)

        # Update p:   p = M \ v - p * (rho*delta/eps)
        sctmp = -rho_n * delta / eps
        BLAS.scal!(n,sctmp,p,1)                 #  - p * (rho*delta/eps)
        BLAS.axpy!(n,one(eltype(p)),mv,1,p,1)   # p = p + mv

        ap   = A(p)
        # eps  = At_mul_B(p, ap)    # eps  = p.' * ap
        eps  = transpose(p) * ap
        abs(eps) < 1e-40 && (flag=-3; break)
        beta = eps / delta

        # Update v:   v = A * p - beta * v,   v=v/norm(v)
        BLAS.scal!(n,-beta,v,1)                  #  -beta * v
        BLAS.axpy!(n,one(eltype(v)),ap,1,v,1)    # v = v + ap
        rho_np = norm( v )
        BLAS.scal!(n,one(eltype(v))/rho_np,v,1)  # v = v / rho_np

        theta_n = rho_np^2 / c_nm / (real(beta)^2+imag(beta)^2)    # theta_n^2
        c_n     = 1.0 / (1.0 + theta_n)                            # c_n^2
        eta     = -eta * rho_n * c_n / (beta * c_nm)
        k2      = theta_nm * c_n

        # Update d:   d = eta * p + k2 * d
        BLAS.scal!(n,one(eltype(b))*k2,d,1)     # k2 * d
        BLAS.axpy!(n,eta,p,1,d,1)               # d = d + eta * p

        # Update x:  x = x + d
        BLAS.axpy!(n,one(eltype(d)),d,1,x,1)

        # Update s:  s = eta * ap + k2 * s
        BLAS.scal!(n,one(eltype(s))*k2,s,1)     # k2 * s
        BLAS.axpy!(n,eta,ap,1,s,1)              # s = s + eta * ap

        # Update r:  r = r - s
        BLAS.axpy!(n,-one(eltype(r)),s,1,r,1)

        err = norm( r ) / bnrm2
        resvec[iter+1] = err

        out==2 && println(@sprintf("%3d\t%1.2e",iter,resvec[iter+1]))
        iter += 1

        err <= tol && (flag=0; break)
    end  # iter = 1:maxIter
    iter = min(iter, maxIter)

    if out>=0
      if flag==-1
        println(@sprintf("qmr iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
						 maxIter,resvec[iter],tol))
      elseif flag==-2
        println(@sprintf("qmr: delta equal to zero at iteration %d. Returned residual has norm %1.2e.", iter,resvec[iter]))
      elseif flag==-3
        println(@sprintf("qmr: epsilon equal to zero at iteration %d. Returned residual has norm %1.2e.", iter,resvec[iter]))
      elseif out>=1
	    println(@sprintf("qmr achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,resvec[iter+1]))
	  end
    end

	return x, flag,resvec[iter+1],iter,resvec[1:iter+1]

end


text:
@testset "generic" begin
    k = SqExponentialKernel()
    @test length(k) == 1
    @test iterate(k) == (k, nothing)
    @test iterate(k, 1) == nothing
end


text:
module TBLIS

using Libdl
using LinearAlgebra
using Hwloc_jll
using tblis_jll

import LinearAlgebra: mul!

export TTensor
export mul!
export add!

global tblis = C_NULL
global tci = C_NULL
global hwloc = C_NULL

function __init__()
    if Sys.isapple()
        global tblis = dlopen(joinpath(dirname(pathof(TBLIS)),"libtblis"))
    else
        global tci = dlopen(tblis_jll.tci)
        global hwloc = dlopen(Hwloc_jll.libhwloc)
        global tblis = dlopen(tblis_jll.tblis)
    end
end

function set_num_threads(n)
    # Previous number of Threads to be returned
    original = get_num_threads()

    tblis_set_num_threads = dlsym(tblis, :tblis_set_num_threads)
    ccall(tblis_set_num_threads, Cvoid, (Cuint,), n)
    return original
end

function get_num_threads()
    tblis_get_num_threads = dlsym(tblis, :tblis_get_num_threads)
    return ccall(tblis_get_num_threads, Cint, ())
end

include("TTensor.jl")
include("Ops.jl")

end #module


text:
"""
    downloadchapter(chapter, dir)

Downloads and unzips a complete chapter from Boletin Estadístico of Bank of Spain

Creates a subfolder with the name of the chapter if it does not exists already

# Example

    downloadchapter(2, "C:\\Users\\Weyland\\Documents\\")
"""

function downloadchapter(capitulo::Int, directory::String)

    n = lpad(capitulo,2,"0")

    #TODO: sustituir por Downloads en base con 1.6
    #Descarga de ficheros desde la página del Banco de España. Cambiar por
    file = HTTP.get("https://www.bde.es/webbde/es/estadis/infoest/series/be"*n*".zip").body
    r = ZipFile.Reader(IOBuffer(file));

    isdir(joinpath(directory,"be"*n)) || mkdir(joinpath(directory,"be"*n)) # Si no hay directorio lo creamos

    for f in r.files
        println("Filename: $(f.name)")
        io = open(joinpath(directory,"be"*n,f.name),"w") #Creación de fichero con nombre igual al del ZIP
        write(io,read(f))
        close(io)
    end
    close(r)

end #function

"""
downloadchapter(f, chapter, dir)

Downloads and unzips a complete chapter from Boletin Estadístico of Bank of Spain.
Then applies a function that acts on the files on the same directory

# Example

    downloadchapter(CSVtoDataFrame, 2, "C:\\Users\\Weyland\\Documents\\")

"""
function downloadchapter(f::Function, capitulo::Int, directory::String)

    downloadchapter(capitulo::Int, directory::String)

    f(joinpath(directory,"be"*lpad(capitulo,2,"0")))

end #function


text:
using OrdinaryDiffEq, DiffEqFlux, GalacticOptim, Optim, Plots, Flux

function lotka_volterra!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

# Initial condition
u0 = [1.0, 1.0]

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

# LV equation parameter. p = [α, β, δ, γ]
p = [1.5, 1.0, 3.0, 1.0]

# Setup the ODE problem, then solve
prob_ode = ODEProblem(lotka_volterra!, u0, tspan, p)
sol_ode = solve(prob_ode, Tsit5())

function predict_adjoint(p)
    return Array(solve(prob_ode, Tsit5(), p=p, saveat = tsteps))
end

function loss_adjoint(p)
    prediction = predict_adjoint(p)
    loss = sum(abs2, x-1 for x in prediction)
    return loss, prediction
end

list_plots = []
iter = 0
callback = function (p, l, pred)
  global iter, list_plots

  if iter == 0
    list_plots = []
  end
  iter += 1

  display(l)

  # using `remake` to re-create our `prob` with current parameters `p`
  remade_solution = solve(remake(prob_ode, p = p), Tsit5(), saveat = tsteps)
  plt = plot(remade_solution, ylim = (0, 6))

  push!(list_plots, plt)
  display(plt)

  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  return false
end

optprob = OptimizationFunction((x,p) -> loss_adjoint(x), GalacticOptim.AutoForwardDiff())

prob = GalacticOptim.OptimizationProblem(optprob, p)

result_ode = GalacticOptim.solve(prob,
                                 BFGS(initial_stepnorm = 0.0001),
                                 cb = callback)

u0 = Float32[2.0; 0.0]
datasize = 30
tspan = (0.0f0, 1.5f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))

dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, tanh),
                  FastDense(50, 2))
prob_neuralode = NeuralODE(dudt2, tspan, Tsit5(), saveat = tsteps)

dudt2 = Chain(x -> x.^3,
              Dense(2, 50, tanh),
              Dense(50, 2))

function predict_neuralode(p)
  Array(prob_neuralode(u0, p))
end

function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end

list_plots = []
iter = 0
callback = function (p, l, pred; doplot = false)
  global list_plots, iter

  if iter == 0
    list_plots = []
  end
  iter += 1

  display(l)

  # plot current prediction against data
  plt = scatter(tsteps, ode_data[1,:], label = "data")
  scatter!(plt, tsteps, pred[1,:], label = "prediction")
  push!(list_plots, plt)
  if doplot
    display(plot(plt))
  end

  return false
end

optprob = OptimizationFunction( (p,x) -> loss_neuralode(p), GalacticOptim.AutoForwardDiff())

prob = GalacticOptim.OptimizationProblem(optprob, prob_neuralode.p)

result_neuralode = GalacticOptim.solve(prob,
                                ADAM(), cb = callback,
                                maxiters = 300)

prob2 = remake(prob,u0=result_neuralode.minimizer)
result_neuralode2 = GalacticOptim.solve(prob2,
                                        BFGS(initial_stepnorm=0.0001),
                                        cb = callback,
                                        maxiters = 100)
@test result_neuralode2.minimum < 10


text:
begin
    struct OldStructSyntax; end
    register(OldStructSyntax, Deprecation(
        "The type-definition keywords (type, immutable, abstract) where changed in Julia 0.6",
        "julia",
        v"0.6.0",
        v"0.7.0-DEV.198",
        typemax(VersionNumber)
    ))
    # Despite having different syntax, all of these parse the same,
    # so macros can't distinguish
    applies_in_macrocall(dep::OldStructSyntax, context) = true

    # Special purpose formatter to unindent multi-line arglists
    function format_paramlist(orig_tree, tree, matches)
        if isexpr(tree, CSTParser.Struct)
            isexpr(children(tree)[2], CSTParser.Curly) || return tree
            children(tree)[2] = format_addindent_body(children(tree)[2], -3, nothing)
        else
            isexpr(children(tree)[3], CSTParser.Curly) || return tree
            children(tree)[3] = format_addindent_body(children(tree)[3], 10, nothing)
        end
        tree
    end

    match(OldStructSyntax,
        "immutable \$name\n\$BODY...\nend",
        "struct \$name\n\$BODY!...\nend",
        format_paramlist
    )
    match(OldStructSyntax,
        "type \$name\n\$BODY...\nend",
        "mutable struct \$name\n\$BODY!...\nend",
        format_paramlist
    )
    match(OldStructSyntax,
        "abstract \$name",
        "abstract type \$name end"
    )
    match(OldStructSyntax,
        "bitstype \$size \$name",
        "primitive type \$name \$size! end"
    )
end

begin
    struct OldTypeAliasSyntax; end
    register(OldTypeAliasSyntax, Deprecation(
        "The `typealias` keyword is deprecated in Julia 0.6",
        "julia",
        v"0.6.0",
        v"0.6.0",
        typemax(VersionNumber)
    ))

    match(OldTypeAliasSyntax,
        "typealias \$X{\$T...} \$B",
        "\$X{\$T...} = \$B"
    )
    match(OldTypeAliasSyntax,
        "typealias \$A \$B",
        "const \$A = \$B"
    )
end

begin
    struct OldStyleConstructor
        # In 0.6/0.7 F{T}() still means a parametric function
        # definition, so we can't perform this upgrade, if there's
        # no `where` clause.
        non_where_curly::Bool
        # In 0.6/0.7 F() inside a struct still refers to the old inner constructor syntax
        # so we can't perform this upgrade
        inner_constructor::Bool
    end
    OldStyleConstructor() = OldStyleConstructor(false, false)

    function dep_for_vers(::Type{OldStyleConstructor}, vers)
        is10 = all(interval->(v"1.0-DEV" <= interval.lower), vers["julia"].intervals)
        OldStyleConstructor(is10, is10)
    end

    register(OldStyleConstructor, Deprecation(
        "This constructor syntax is no longer required",
        "julia",
        v"0.6.0", v"0.6.0", typemax(VersionNumber)
    ))
    function format_old_constructor(orig_tree, tree, matches)
        ok, call′ = format_new_call_expr(tree, orig_tree, tree, orig_tree)
        ok || return tree
        return call′
    end
    function filter_non_where_curly(S, dep, tree, matches)
        dep.non_where_curly && return true
        (isexpr(parent(tree), BinarySyntaxOpCall) || isexpr(parent(tree), CSTParser.WhereOpCall)) || return false
        isexpr(children(parent(tree))[2], OPERATOR, Tokens.WHERE) || return false
        return true
    end
    match(OldStyleConstructor,
        "(::Type{\$NAME{\$T...}})(\$ARGS...)",
        "\$NAME{\$T...}(\$ARGS...)",
        format_old_constructor,
        filter = filter_non_where_curly
    )
    function filter_in_struct(S, dep, tree, matches)
        dep.inner_constructor && return true
        return get_struct_parent(tree.parent) === nothing
    end
    function filter_params(S, dep, tree, matches)
        # Handled by the above
        name = first(matches[:NAME][2])
        isexpr(name, Curly) && return false
        # Filter out (::Type{<:Foo})
        if isexpr(name, CSTParser.UnarySyntaxOpCall) &&
            isexpr(children(name)[1], CSTParser.OPERATOR, Tokens.ISSUBTYPE)
            return false
        end
        p = parent(tree)
        isexpr(p, BinarySyntaxOpCall)  || isexpr(p, CSTParser.WhereOpCall)  || return true
        isexpr(children(p)[2], OPERATOR, Tokens.WHERE) || return true
        # Get all the parameter names
        names = extract_identifiers(children(p)[3:end])
        !(Expr(name) in names)
    end
    match(OldStyleConstructor,
        "(::Type{\$NAME})(\$ARGS...)",
        "\$NAME(\$ARGS...)",
        format_old_constructor,
        filter = (args...)->filter_params(args...) && filter_in_struct(args...)
    )
end

begin
    struct ObsoleteCompatString; end
    register(ObsoleteCompatString, Deprecation(
        "This Compat definition is deprecated",
        "julia",
        v"0.5.0", v"0.6.0", typemax(VersionNumber)
    ))

    match(ObsoleteCompatString,
        "Compat.UTF8String",
        "String",
    )

    match(ObsoleteCompatString,
        "Compat.ASCIIString",
        "String",
    )
end

begin
    struct ReadString; end
    register(ReadString, Deprecation(
        "readstring(x) is deprecated to read(x, String)",
        "julia",
        v"0.7.0-DEV.1053", v"1.0", typemax(VersionNumber)
    ))

    match(ReadString,
        "readstring(\$A)",
        "read(\$A, String)"
    )
end

begin
    struct ConditionalWhitespace; end
    register(ConditionalWhitespace, Deprecation(
        "The ternary operator now requires whitespace on both sides of the punctuation.",
        "julia",
        v"0.4.0", v"0.7.0-DEV.797", typemax(VersionNumber)
    ))

    function add_ws_lead_trail!(ret, idx, exprs)
        any = false
        expr = exprs[2]
        allws = string(trailing_trivia(exprs[1]), leading_trivia(expr))
        if isempty(allws) || !isspace(allws[end])
            expr = children(ret)[idx] = TriviaReplacementNode(ret, expr, string(leading_trivia(expr), " "), trailing_trivia(expr))
            any = true
        end
        allws = string(trailing_trivia(expr), leading_trivia(exprs[3]))
        if isempty(allws) || !isspace(allws[end])
            expr = children(ret)[idx] = TriviaReplacementNode(ret, expr, leading_trivia(expr), string(" ", trailing_trivia(expr)))
            any = true
        end
        any
    end

    match(ConditionalWhitespace, CSTParser.ConditionalOpCall) do x
        dep, expr, resolutions, context = x
        ret = ChildReplacementNode(nothing, collect(children(expr)), expr)
        changed = add_ws_lead_trail!(ret, 2, children(expr)[1:3])
        changed |= add_ws_lead_trail!(ret, 4, children(expr)[3:5])
        if changed
            buf = IOBuffer()
            print_replacement(buf, ret, false, false)
            push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
        end
    end
end

begin
    struct GeneratorWhitespace; end
    register(GeneratorWhitespace, Deprecation(
        "Generators and comprehensions now require whitespace before the `for`",
        "julia",
        v"0.4.0", v"0.7.0-DEV.797", typemax(VersionNumber)
    ))

    match(GeneratorWhitespace, CSTParser.Generator) do x
        dep, expr, resolutions, context = x
        ret = ChildReplacementNode(nothing, collect(children(expr)), expr)
        body, fornode, iterand = children(expr)
        allws = string(trailing_trivia(body), leading_trivia(fornode))
        if isempty(allws) || !isspace(allws[end])
            children(ret)[2] = TriviaReplacementNode(ret, fornode, string(leading_trivia(fornode), " "), trailing_trivia(fornode))
            buf = IOBuffer()
            print_replacement(buf, ret, false, false)
            repl = String(take!(buf))
            push!(resolutions, TextReplacement(dep, expr.span, repl))
        end
    end
end

# This doesn't work becuase the `...` what we want to match against
# is also used as a wildcard sentinel
begin
    struct TupleSplat; end

    register(TupleSplat, Deprecation(
        "Splat of single value in tuple needs trailing comma",
        "julia",
        v"0.7.0-DEV.2559", v"1.0", typemax(VersionNumber)
    ))

    function filter_funcdef_broadcast_multiarg(S, dep, expr, matches)
        length(children(expr)) == 3     || return false  # 3 = '(' + ')' + 1 arg
        isexpr(children(expr)[2], UnarySyntaxOpCall) || return false
        CSTParser.has_sig(parent(expr)) && return false
        # AST Interpolation $(x...)
        if isexpr(CSTParser.parent(expr), CSTParser.UnarySyntaxOpCall) &&
            isexpr(children(CSTParser.parent(expr))[1], CSTParser.OPERATOR, Tokens.EX_OR)
            return false
        # Broadcasting
        elseif isexpr(CSTParser.parent(expr), CSTParser.BinarySyntaxOpCall) &&
            isexpr(children(CSTParser.parent(expr))[2], CSTParser.OPERATOR, Tokens.DOT)
            return false
        end
        return true
    end

    match(TupleSplat,
          "(\$ID...)",
          "(\$ID...,)",
          filter = filter_funcdef_broadcast_multiarg
    )
end

begin
    struct Void2Nothing; end
    register(Void2Nothing, Deprecation(
        "The type `Void` is renamed to `Nothing` (and a synonym `Cvoid` is added)",
        "julia",
        v"0.7.0-DEV.3137",
        v"1.0",
        typemax(VersionNumber)
    ))

    match(Void2Nothing,
        "Void",
        "Nothing",
    )

    match(Void2Nothing,
        "Ptr{Void}",
        "Ptr{Cvoid}",
    )

    match(Void2Nothing, CSTParser.Call) do x
        # Match ccall only
        dep, expr, resolutions, context = x
        is_identifier(children(expr)[1], "ccall") || return
        _replace_void_cvoid!(dep, resolutions, expr)
    end

    function _replace_void_cvoid!(dep, resolutions, expr)
        if expr.expr isa CSTParser.IDENTIFIER && expr.expr.val == "Void"
            push!(resolutions, TextReplacement(dep, expr.span, "Cvoid"))
        end
        isempty(children(expr)) && return
        for c in children(expr)
            _replace_void_cvoid!(dep, resolutions, c)
        end
        return
    end
end

begin
    struct IssubtypeToInfix; end
    register(IssubtypeToInfix, Deprecation(
        "issubtype is deprecated for `<:`",
        "julia",
        v"0.7.0-DEV.1162", v"1.0", typemax(VersionNumber)
    ))

    match(IssubtypeToInfix,
        "issubtype(\$A, \$B)",
        "\$A <: \$B!",
    )
end

function is_scope_Base(scope)
    while scope !== nothing
        if scope.t == "Module" && scope.namespace == "Base"
            return true
        end
        scope = scope.parent
    end
    return false
end

function filter_base_id(analysis, expr)
    S, file_scope = analysis
    @assert isa(expr, OverlayNode{CSTParser.IDENTIFIER})
    isexpr(parent(expr), CSTParser.Quotenode) && return false
    # Not on the LHS of a keyword arg
    if isexpr(parent(expr), CSTParser.Kw)
        expr == children(parent(expr))[1] && return false
    end
    scope = CSTAnalyzer.find_scope(file_scope, expr.span)
    binding = CSTAnalyzer.find_ref(span_text(expr), scope, S)
    !isa(binding, CSTAnalyzer.Binding) && return false
    return binding.t == "BaseCore" || is_scope_Base(scope)
end

struct SimpleRename
    min_ver::VersionNumber
end
SimpleRename() = SimpleRename(v"0.6.0")
function dep_for_vers(::Type{SimpleRename}, vers)
    SimpleRename(minimum(map(interval->interval.lower, vers["julia"].intervals)))
end
register(SimpleRename, Deprecation(
    "a number of identifiers were renamed",
    "julia",
    v"0.7.0-DEV.1", v"0.7.0-DEV.1", typemax(VersionNumber)
))

macro add_rename(pattern, replacement, ver, hascompat = false)
    :(($ver, $(string(pattern))=>$(string(replacement)), $hascompat))
end

match(SimpleRename, CSTParser.IDENTIFIER) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    for rename in renames
        if rename[1] <= dep.min_ver && is_identifier(expr, rename[2][1])
            filter_base_id(analysis, expr) || return
            push!(resolutions, TextReplacement(dep, expr.span, rename[2][2]))
        end
    end
end

const renames = [
    @add_rename ipermute!      invpermute!      v"0.7.0-DEV.3173" true
    @add_rename unshift!       pushfirst!       v"0.7.0-DEV.3155" true
    @add_rename shift!         popfirst!        v"0.7.0-DEV.3155" true
    @add_rename JULIA_HOME     Sys.BINDIR       v"0.7.0-DEV.3073" true
    @add_rename CartesianRange CartesianIndices v"0.7.0-DEV.3025" true
    @add_rename Display        AbstractDisplay  v"0.7.0-DEV.2695" true
    @add_rename strwidth       textwidth        v"0.7.0-DEV.1930" true
    @add_rename charwidth      textwidth        v"0.7.0-DEV.1930" true
    @add_rename sqrtm          sqrt             v"0.7.0-DEV.1599"
    @add_rename logm           log              v"0.7.0-DEV.1597"
    @add_rename expm           exp              v"0.7.0-DEV.1486"
    @add_rename Complex32      ComplexF16       v"0.7.0-DEV.2919" true
    @add_rename Complex64      ComplexF32       v"0.7.0-DEV.2919" true
    @add_rename Complex128     ComplexF64       v"0.7.0-DEV.2919" true
    @add_rename EnvHash        EnvDict          v"0.7.0-DEV.2265"
    @add_rename ctranspose     adjoint          v"0.7.0-DEV.1415" true
    @add_rename ctranspose!    adjoint!         v"0.7.0-DEV.1415" true
    @add_rename writecsv       writedlm         v"0.7.0-DEV.1737"
    @add_rename readcsv        readdlm          v"0.7.0-DEV.1740"
    @add_rename is_linux       Sys.islinux      v"0.7.0-DEV.914" true
    @add_rename is_bsd         Sys.isbsd        v"0.7.0-DEV.914" true
    @add_rename is_apple       Sys.isapple      v"0.7.0-DEV.914" true
    @add_rename is_unix        Sys.isunix       v"0.7.0-DEV.914" true
    @add_rename is_windows     Sys.iswindows    v"0.7.0-DEV.914" true
    @add_rename AbstractIOBuffer GenericIOBuffer v"0.7.0-DEV.961"
    @add_rename select         partialsort      v"0.7.0-DEV.1535" true
    @add_rename select!        partialsort!     v"0.7.0-DEV.1535" true
    @add_rename selectperm     partialsortperm  v"0.7.0-DEV.1535" true
    @add_rename selectperm!    partialsortperm! v"0.7.0-DEV.1535" true
    @add_rename Range          AbstractRange    v"0.7.0-DEV.1721" true
    @add_rename isleaftype     isconcretetype   v"0.7.0-DEV.1775" true
    @add_rename isnumber       isnumeric        v"0.7.0-DEV.1775" true
    @add_rename Associative    AbstractDict     v"0.7.0-DEV.2951" true
    @add_rename find           findall          v"0.7.0-DEV.3415" true
    @add_rename module_parent   parentmodule    v"0.7.0-DEV.3460" true
    @add_rename datatype_module parentmodule    v"0.7.0-DEV.3460" true
    @add_rename function_module parentmodule    v"0.7.0-DEV.3460" true
    @add_rename indmin         argmin           v"0.7.0-DEV.3516" true
    @add_rename indmax         argmax           v"0.7.0-DEV.3516" true
    @add_rename module_name    nameof           v"0.7.0-DEV.3539" true
    @add_rename function_name  nameof           v"0.7.0-DEV.3539" true
    @add_rename datatype_name  nameof           v"0.7.0-DEV.3539" true
    @add_rename method_exist   hasmethod        v"0.7.0-DEV.3455" true
    @add_rename object_id      objectid         v"0.7.0-DEV.3455" true
    @add_rename DevNull        devnull          v"0.7.0-DEV.4068" true
    @add_rename STDIN          stdin            v"0.7.0-DEV.4068" true
    @add_rename STDOUT         stdout           v"0.7.0-DEV.4068" true
    @add_rename STDERR         stderr           v"0.7.0-DEV.4068" true
    @add_rename reprmime       repr             v"0.7.0-DEV.4010"
    @add_rename isabstract     isabstracttype   v"0.7.0-DEV.1775" true
    @add_rename iteratorsize   IteratorSize     v"0.7.0-DEV.3309" true
    @add_rename iteratoreltype IteratorEltype   v"0.7.0-DEV.3309" true
    @add_rename nb_available   bytesavailable   v"0.7.0-DEV.3477" true
    @add_rename broadcast_indices broadcast_axes  v"0.7.0-DEV.4936"
    @add_rename atan2          atan             v"0.7.0-alpha.44"
    @add_rename srand          Random.seed!     v"0.7.0-beta2.171"
    @add_rename realmin        floatmin         v"0.7.0-beta2.169"
    @add_rename realmax        floatmax         v"0.7.0-beta2.169"
]

for (ver, repl, hascompat) in renames
    if hascompat
        match(SimpleRename, "Compat.$(repl.second)", repl.second)
    end
end

function filter_no_kw_already(S, dep, tree, matches)
    @assert isexpr(tree, CSTParser.Call)
    for c in children(tree)
        isexpr(c, CSTParser.Kw) && return false
    end
    return true
end

function filter_not_def(S, dep, tree, matches)
    p = parent(tree)
    isexpr(p, CSTParser.WhereOpCall) && (p = parent(p))
    isexpr(p, CSTParser.Macro) && return false
    CSTParser.defines_function(p) && return false
    return true
end

function resolve_qualified_expr(id, analysis)
    S, file_scope = analysis
    if isexpr(id, CSTParser.BinarySyntaxOpCall) && isexpr(children(id)[2], CSTParser.OPERATOR, Tokens.DOT)
        isexpr(children(id)[3], CSTParser.Quotenode) || return nothing
        cid = children(children(id)[3])[1]
        isexpr(cid, CSTParser.IDENTIFIER) || return nothing
        return (is_identifier(children(id)[1], "Base") ||
                is_identifier(children(id)[1], "Core")) ? (cid, "BaseCore") : nothing
    else
        isexpr(id, CSTParser.IDENTIFIER) || return nothing
        binding = CSTAnalyzer.resolve(S, file_scope, id)
        !isa(binding, CSTAnalyzer.Binding) && return nothing
        return (id, binding.t)
    end
end

const BaseCoreNumberTypes = [
    # Core
    :AbstractFloat, :Bool, :Float16, :Float32, :Float64, :Int, :Int128,
    :Int16, :Int32, :Int64, :Int8, :Integer, :Number, :Real, :Signed,
    :UInt, :UInt128, :UInt16, :UInt32, :UInt64, :UInt8, :Unsigned,
    # Base
    :BigFloat, :BigInt, :Cchar, :Cdouble, :Cfloat, :Cint, :Cintmax_t,
    :Clong, :Clonglong, :Complex, :Complex128, :Complex32, :Complex64,
    :Cptrdiff_t, :Cshort, :Csize_t, :Cssize_t, :Cuchar, :Cuint, :Cuintmax_t,
    :Culong, :Culonglong, :Cushort, :Cwchar_t, :Irrational, :Rational
]

# Ideally we'd do the opposite and find things that are likely ::Number, but
# that's harder
function filter_not_likely_type(analysis, dep, tree, matches)
    @assert isexpr(tree, CSTParser.Call)
    x = resolve_qualified_expr(children(tree)[3], analysis)
    x === nothing && return true
    (id, t) = x
    t == :DataType && return false
    t == "BaseCore" || return true
    (Symbol(id_name(id)) in BaseCoreNumberTypes) && return false
    return true
end

const RoundingMode = [
    :RoundNearest, :RoundToZero, :RoundUp, :RoundDown,
    :RoundFromZero, :RoundNearestTiesAway, :RoundNearestTiesUp
]

function filter_not_rounding_mode(analysis, dep, tree, matches)
    @assert isexpr(tree, CSTParser.Call)
    x = resolve_qualified_expr(children(tree)[end - 1], analysis)
    x === nothing && return true
    (id, t) = x
    t == "BaseCore" || return true
    (Symbol(id_name(id)) in RoundingMode) && return false
    return true
end

keyword_default_filter(args...) = filter_no_kw_already(args...) && filter_not_def(args...)

struct KeywordsUnlocked; end
begin
    register(KeywordsUnlocked, Deprecation(
        "a bunch of functions got keyword arguments instead of positional arguments",
        "julia",
        v"0.7.0-DEV.1", v"1.0", typemax(VersionNumber)
    ))
    match(KeywordsUnlocked,
        "Timer(\$timeout, \$repeat)",
        "Timer(\$timeout, interval=\$repeat!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "Timer(\$callback, \$delay, \$repeat)",
        "Timer(\$callback, \$delay!, interval=\$repeat!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "names(\$m, \$all)",
        "names(\$m, all=\$all!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "names(\$m, \$all, \$imported)",
        "names(\$m, all=\$all!, imported=\$imported!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "code_native(\$io, \$f, \$types, \$syntax)",
        "code_native(\$io, \$f!, \$types!, syntax=\$syntax!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "code_native(\$f, \$types, \$syntax)",
        "code_native(\$f, \$types!, syntax=\$syntax!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "eachmatch(\$re, \$str, \$overlap)",
        "eachmatch(\$re, \$str!, overlap=\$overlap!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "matchall(\$re, \$str, \$overlap)",
        "matchall(\$re, \$str!, overlap=\$overlap!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "chop(\$s, \$head)",
        "chop(\$s, head=\$head!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "chop(\$s, \$head, \$tail)",
        "chop(\$s, head=\$head!, tail=\$tail!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "tryparse(\$T, \$s, \$base)",
        "tryparse(\$T, \$s!, base=\$base!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "mkdir(\$path, \$mode)",
        "mkdir(\$path, mode=\$mode!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "mkpath(\$path, \$mode)",
        "mkpath(\$path, mode=\$mode!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "countlines(\$x, \$eol)",
        "countlines(\$x, eol=\$eol!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "PipeBuffer(\$data, \$maxsize)",
        "PipeBuffer(\$data, maxsize=\$maxsize!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "unsafe_wrap(\$T, \$pointer, \$dims, \$own)",
        "unsafe_wrap(\$T, \$pointer!, \$dims!, own=\$own!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "linspace(\$start, \$stop)",
        "range(\$start, stop=\$stop, length=50)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "logspace(\$start, \$stop)",
        "exp10.(range(\$start!, stop=\$stop!, length=50))",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "linspace(\$start, \$stop, \$length)",
        "range(\$start!, stop=\$stop!, length=\$length)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "logspace(\$start, \$stop, \$length)",
        "exp10.(range(\$start, stop=\$stop, length=\$length))",
        filter = keyword_default_filter
    )
    for f in ("trunc", "floor", "ceil", "round")
        match(KeywordsUnlocked,
            "$f(\$x, \$digits)",
            "$f(\$x!, digits=\$digits!)",
            filter = (args...)->filter_no_kw_already(args...) &&
                                filter_not_likely_type(args...) &&
                                filter_not_def(args...) &&
                                filter_not_rounding_mode(args...)
        )
        match(KeywordsUnlocked,
            "$f(\$x, \$digits, \$base)",
            "$f(\$x!, digits=\$digits!, base=\$base!)",
            filter = (args...)->filter_no_kw_already(args...) &&
                                filter_not_likely_type(args...) &&
                                filter_not_def(args...) &&
                                filter_not_rounding_mode(args...)
        )
    end
    match(KeywordsUnlocked,
        "signif(\$x, \$digits)",
        "round(\$x!, sigdigits=\$digits!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "signif(\$x, \$digits, \$base)",
        "round(\$x!, sigdigits=\$digits!, base=\$base!)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "method_exists(\$f, \$t, \$world)",
        "hasmethod(\$f, \$t, world=\$world)",
        filter = keyword_default_filter
    )
    match(KeywordsUnlocked,
        "print_with_color(\$color, \$args...)",
        "printstyled(\$args!..., color=\$color)",
        filter = keyword_default_filter
    )
end

function filter_func_base_id(analysis, dep, tree, matches)
    @assert isexpr(tree, CSTParser.Call)
    id = first(children(tree))
    return filter_base_id(analysis, id)
end

struct Misc07; end
begin
    register(Misc07, Deprecation(
        "a number of function APIs were changed",
        "julia",
        v"0.7.0-DEV.1", v"1.0", typemax(VersionNumber)
    ))
    match(Misc07,
        "ismatch(\$r, \$s)",
        "occursin(\$r, \$s)"
    )
    match(Misc07,
        "findin(\$a, \$b)",
        "findall(in(\$b), \$a)"
    )
    match(Misc07,
        "method_exists(\$f, \$t)",
        "hasmethod(\$f, \$t)"
    )
    match(Misc07,
        "parse(\$s)",
        "Meta.parse(\$s)",
        filter = filter_func_base_id
    )
    match(Misc07,
        "prevpow2(\$x)",
        "prevpow(2, \$x)"
    )
    match(Misc07,
        "nextpow2(\$x)",
        "nextpow(2, \$x)"
    )
end

struct RenamedKeyword
    min_ver::VersionNumber
end
RenamedKeyword() = RenamedKeyword(v"0.6.0")
function dep_for_vers(::Type{RenamedKeyword}, vers)
    RenamedKeyword(minimum(map(interval->interval.lower, vers["julia"].intervals)))
end
register(RenamedKeyword, Deprecation(
    "a number of keyword arguments were renamed",
    "julia",
    v"0.7.0-DEV.1", v"0.7.0-DEV.1", typemax(VersionNumber)
))

function process_renamed_keyword(x)
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    for rename in keyword_renames
        if rename[1] <= dep.min_ver && is_identifier(expr, rename[3][1])
            # Check if we're in keyword context
            p = parent(expr)
            isexpr(p, CSTParser.Kw) || return
            # Check if the call we're in is on the list
            pp = parent(p)
            isexpr(pp, CSTParser.Parameters) && (pp = parent(pp))
            isexpr(pp, CSTParser.Call) || return
            call_id = children(pp)[1]
            any(name->is_identifier(call_id, name), rename[2]) || return
            push!(resolutions, TextReplacement(dep, expr.span, rename[3][2]))
        end
    end
end

match(RenamedKeyword, CSTParser.IDENTIFIER) do x
    process_renamed_keyword(x)
end

const keyword_renames = [
    (v"0.7.0-DEV.3995", ["cp", "mv", "cptree"], "remove_destination"=>"force")
    (v"0.7.0-DEV.4724", ["split", "rsplit"], "keep"=>"keepempty")
]

# =======================
struct LocalConst; end
register(LocalConst, Deprecation(
    "The `const` keyword in local scope has been deprecated (it never worked anyway)",
    "julia",
    v"0.1.0", v"0.7.0-DEV.1", typemax(VersionNumber)
))

match(LocalConst, CSTParser.Const) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    S, file_scope = analysis
    isexpr(parent(expr), CSTParser.Global) && return
    isexpr(children(expr)[2], CSTParser.Global) && return
    scope = CSTAnalyzer.find_scope(file_scope, expr.span)
    if !(scope.t in ("__toplevel__", "Module")) && !(scope.t in ("Quote", "@eval"))
        buf = IOBuffer()
        print_replacement(buf, children(expr)[2], false, false)
        push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
    end
    return
end

# =======================
struct ArrayUndef; end
register(ArrayUndef, Deprecation(
    "Arrays now require an explicit `undef` initializer",
    "julia",
    v"0.7.0-alpha", v"0.7.0-DEV.1", typemax(VersionNumber)
))

match(ArrayUndef, CSTParser.Call) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    curly = children(expr)[1]
    isexpr(curly, CSTParser.Curly) || return
    isexpr(children(curly)[1], CSTParser.IDENTIFIER) || return
    id = Symbol(id_name(children(curly)[1]))
    (id == :Array || id == :Vector) || return
    args = children(expr)[3:end-1]
    is_vector = id == :Vector || (
        length(children(curly)) == 6 &&
        isexpr(children(curly)[5], CSTParser.LITERAL) &&
        Expr(children(curly)[5]) == 1
    )
    # Vector{T}(0) -> Vector{T}()
    if is_vector
        if length(args) == 1
            arg = args[1]
            if isexpr(arg, CSTParser.LITERAL) && Expr(arg.expr) == 0
                repl = ChildReplacementNode(nothing, [
                    children(expr)[1:2]...,
                    children(expr)[4:end]...
                ], expr)
                buf = IOBuffer()
                print_replacement(buf, repl, false, false)
                push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
                return
            end
        elseif length(args) == 0
            # This is ok
            return
        end
    end
    if all(args) do arg
                isexpr(arg, CSTParser.PUNCTUATION, Tokens.COMMA) && return true
                isexpr(arg, CSTParser.LITERAL) || return false
                isa(Expr(arg.expr), Integer) || return false
                return true
            end
        repl = ChildReplacementNode(nothing, [
            children(expr)[1:2]...,
            ReplacementNode("undef", "", ""),
            (isempty(args) ? () : (ReplacementNode(",","", " "),))...,
            children(expr)[3:end]...
        ], expr)
        buf = IOBuffer()
        print_replacement(buf, repl, false, false)
        push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
    end
end

# =======================
struct LoggingMacros; end
register(LoggingMacros, Deprecation(
    "Several logging functions became macros",
    "julia",
    v"0.7.0-DEV.2979", v"0.7.0-DEV.1", typemax(VersionNumber)
))

match(LoggingMacros, CSTParser.Call) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    fname = children(expr[1])
    (is_identifier(fname, "warn") || is_identifier(fname, "info")) || return
    filter_base_id(analysis, fname) || return
    # For now, only support a single string argument
    length(children(expr)) == 4 || return
    isexpr(children(expr)[3], CSTParser.LITERAL) || return
    repl = ChildReplacementNode(nothing, [
        ReplacementNode(string("@",id_name(fname)), leading_trivia(fname), trailing_trivia(fname)),
        children(expr)[2:end]...
    ], expr)
    buf = IOBuffer()
    print_replacement(buf, repl, false, false)
    push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
end

# =======================
struct SortSlices; end
register(SortSlices, Deprecation(
    "sortrows/sortcols have been deprecated in favor of sortslices",
    "julia",
    v"0.7.0-beta2.193", v"0.7.0-beta2.193", typemax(VersionNumber)
))

match(SortSlices, CSTParser.Call) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    fname = children(expr[1])
    (is_identifier(fname, "sortrows") || is_identifier(fname, "sortcols")) || return
    isrows = is_identifier(fname, "sortrows")
    repl = ChildReplacementNode(nothing, [
        ReplacementNode("sortslices", leading_trivia(fname), trailing_trivia(fname)),
        children(expr)[2:3]...,
        ReplacementNode(string(", dims=", isrows ? 1 : 2), "", ""),
        children(expr)[4:end]...
    ], expr)
    buf = IOBuffer()
    print_replacement(buf, repl, false, false)
    push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
end

# =======================
struct Precompile; end
register(Precompile, Deprecation(
    "__precompile__([true]) is now the default",
    "julia",
    v"0.7.0-beta2.199", v"0.7.0-beta2.199", typemax(VersionNumber)
))

match(Precompile, CSTParser.Call) do x
    dep, expr, resolutions, context, analysis = x
    context.in_macrocall && return
    fname = children(expr[1])
    (is_identifier(fname, "__precompile__")) || return
    if length(children(expr)) == 4
        arg = children(expr)[3]
        if !isexpr(arg, CSTParser.LITERAL) || Expr(arg) !== true
            return
        end
    end
    push!(resolutions, TextReplacement(dep, expr.fullspan, in_statement_position(expr) ? "" : "nothing"))
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto
#see LICENSE.txt
#this work was supported in part by DARPA Contract D15PC00135
#expwalk.jl
#returns a series of unums that are comprehensively tested for satisfying
#the boolean function bf  Uses a binary search method.
function expwalk(bf, ESS, FSS, rlow, rhigh, sign)
  res = Utype[]
  (esize, exponent) = encode_exp(rlow)
  low_exact = Unum{ESS,FSS}(z16, esize, sign, z64, exponent)
  low_ulp = unum_unsafe(low_exact, sign | UNUM_UBIT_MASK)
  (esize, exponent) = encode_exp(rhigh)
  mid_exact = Unum{ESS,FSS}(z16, esize, sign, t64, exponent)
  high_ulp = unum_unsafe(mid_exact, sign | UNUM_UBIT_MASK)

  if rlow == rhigh
    bf(low_exact) && (res = vcat(res, low_exact))
    bf(low_ulp) && (res = vcat(res, low_ulp))
    bf(mid_exact) && (res = vcat(res, mid_exact))
    !isnan(high_ulp) && bf(high_ulp) && (res = vcat(res, high_ulp))
  elseif rlow == rhigh - 1
    pres1 = expwalk(bf, ESS, FSS, rlow, rlow, sign)
    pres2 = expwalk(bf, ESS, FSS, rhigh, rhigh, sign)
    res = [res, pres1, pres2]
  else
    isnan(high_ulp) && (high_ulp = mmr(Unum{ESS,FSS}, sign))
    params = (sign != 0) ? (high_ulp, low_ulp) : (low_ulp, high_ulp)
    bres = bf(Ubound{ESS,FSS}(params...))
    if bres #then we have to try another round.
      mid = (rlow + rhigh) ÷ 2
      pres1 = expwalk(bf, ESS, FSS, rlow, mid, sign)
      pres2 = expwalk(bf, ESS, FSS, mid + 1, rhigh, sign)
      res = [res, pres1, pres2]
    end
  end
  res
end

function subnormalwalk(bf, ESS, FSS, sign)
  res = Utype[]
  low_ulp = Unum{ESS,FSS}(z16, z16, sign | UNUM_UBIT_MASK, z64, z64)
  exact_mid = Unum{ESS,FSS}(z16, z16, sign, t64, z64)
  high_ulp = unum_unsafe(exact_mid, sign | UNUM_UBIT_MASK)
  bf(low_ulp) && (res = vcat(res, low_ulp))
  bf(exact_mid) && (res = vcat(res, exact_mid))
  bf(high_ulp) && (res = vcat(res, high_ulp))
  res
end

function fullwalk(bf, ESS, FSS)
  res = Utype[]
  res = vcat(res, expwalk(bf, ESS, FSS, min_exponent(ESS), max_exponent(ESS), UNUM_SIGN_MASK))
  res = vcat(res, subnormalwalk(bf, ESS, FSS, UNUM_SIGN_MASK))
  bf(zero(Unum{ESS,FSS})) && (res = vcat(res, zero(Unum{ESS,FSS})))
  res = vcat(res, subnormalwalk(bf, ESS, FSS, z16))
  res = vcat(res, expwalk(bf, ESS, FSS, min_exponent(ESS), max_exponent(ESS), z16))
  res
end

export expwalk, subnormalwalk, fullwalk

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
struct TrapezoidalPoints{T <: AbstractFloat} <: AbstractVector{T}
    data::LinRange{T}
    function TrapezoidalPoints{T}(n) where {T}
        if n <= 1
            return new{T}(LinRange{T}(0,0,n))
        else
            return new{T}(LinRange{T}(-1,1,n))
        end
    end
end
TrapezoidalPoints(n) = TrapezoidalPoints{Float64}(n)
Base.size(x::TrapezoidalPoints) = size(x.data)
Base.getindex(x::TrapezoidalPoints, i::Int) = x.data[i]


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "machines.jl"))

Blank = " " # "□"
# Tape = 11111111111
# Tape = "101010" * Blank * "1"

global programme_counter = 0
programme_description(description::AbstractString) = (global programme_counter; programme_counter += 1; println("\t$programme_counter) $description"))

println("The following are the list of programmes to choose from:")

programme_description("Remove last character from string")
programme_description("Duplicate input string")
programme_description("Calculate the tally-code successor")
programme_description("Calculate the binary successor")
programme_description("Calculate the binary predecessor")
programme_description("Calculates the sum of two binary numbers (blank delimited)")
programme_description("x → x mod 3, where x is in tally-code")
programme_description("Halts if and only if the input has an equal number of zeros and ones")

println("\n")
println("Please choose a number as defined above.")

ChosenProgramme = parse(Int, readline())

println("\n")
println("Now please enter a tape.")
Tape = readline()
println("\n")

if Tape == ""
    throw(error("Don't mess with me.  The answer is probably a blank string, unless you want to compute the successor function in which case the answer is one."))
end

println("Do you want to show the process of the Turing Machine?")
Show = nothing
verbosity_input = readline()
if occursin(r"y"i, verbosity_input)
    Show = true
elseif occursin(r"n"i, verbosity_input)
    Show = false
else
    throw(error("You need to tell the programme whether or not you want to show the process of the Turing Machine."))
end
println("\n")

# define programmes
const programmes = [
    (TMProgramme("Turing Machine to Truncate String", "q0", "halt", Blank,
        [
            Rule("q0", "q0", "0", "0", Right),
            Rule("q0", "q0", "1", "1", Right),
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "0", "0", Right),
            Rule("q1", "q1", "1", "1", Right),
            Rule("q1", "q2", Blank, Blank, Left),
            Rule("q2", "q3", "0", Blank, Left),
            Rule("q2", "q3", "1", Blank, Left),
            Rule("q3", "q3", "0", "0", Left),
            Rule("q3", "q3", "1", "1", Left),
            Rule("q3", "q4", Blank, Blank, Left),
            Rule("q4", "q4", "0", "0", Left),
            Rule("q4", "q4", "1", "1", Left),
            Rule("q4", "halt", Blank, Blank, Right)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to Duplicate Input", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q2", "1", "A", Right),
            Rule("q2", "q2", "1", "1", Right),
            Rule("q2", "q2", "0", "0", Right),
            Rule("q2", "q2", "A", "A", Right),
            Rule("q2", "q2", "B", "B", Right),
            Rule("q2", "q3", Blank, "A", Left),
            Rule("q3", "q3", "A", "A", Left),
            Rule("q3", "q3", "B", "B", Left),
            Rule("q3", "q6", "1", "1", Left),
            Rule("q3", "q6", "0", "0", Left),
            Rule("q3", "q7", Blank, Blank, Right),
            Rule("q6", "q6", "1", "1", Left),
            Rule("q6", "q6", "0", "0", Left),
            Rule("q6", "q1", "A", "A", Right),
            Rule("q6", "q1", "B", "B", Right),
            Rule("q1", "q4", "0", "B", Left),
            Rule("q4", "q4", "1", "1", Right),
            Rule("q4", "q4", "0", "0", Right),
            Rule("q4", "q4", "A", "A", Right),
            Rule("q4", "q4", "B", "B", Right),
            Rule("q4", "q5", Blank, "B", Left),
            Rule("q5", "q5", "A", "A", Left),
            Rule("q5", "q5", "B", "B", Left),
            Rule("q5", "q6", "1", "1", Left),
            Rule("q5", "q6", "0", "0", Left),
            Rule("q5", "q7", Blank, Blank, Right),
            Rule("q7", "q7", "A", "1", Right),
            Rule("q7", "q7", "B", "0", Right),
            Rule("q7", "q8", Blank, Blank, Left),
            Rule("q8", "q8", "1", "1", Left),
            Rule("q8", "q8", "0", "0", Left),
            Rule("q8", "halt", Blank, Blank, Right)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to Calculate the Tally-code Successor", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "1", "1", Right),
            Rule("q1", "halt", Blank, "1", Stay)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to Calculate the Binary Successor", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "0", "0", Right),
            Rule("q1", "q1", "1", "1", Right),
            Rule("q1", "q2", Blank, Blank, Left),
            Rule("q2", "q2", "1", "0", Left),
            Rule("q2", "halt", "0", "1", Stay),
            Rule("q2", "halt", Blank, "1", Stay)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to Calculate the Binary Predecessor", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "0", "0", Right),
            Rule("q1", "q1", "1", "1", Right),
            Rule("q1", "q2", Blank, Blank, Left),
            Rule("q2", "q2", "0", "1", Left),
            Rule("q2", "halt", "1", "0", Stay)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to Add Two Binary Numbers", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "0", "0", Right),
            Rule("q1", "q1", "1", "1", Right),
            Rule("q1", "q1", "A", "A", Right),
            Rule("q1", "q1", "B", "B", Right),
            Rule("q2", "q1", "0", "A", Right),
            Rule("q2", "q1", "1", "B", Right),
            Rule("q2", "q2", "A", "A", Left),
            Rule("q2", "q2", "B", "B", Left),
            Rule("q3", "q2", Blank, Blank, Left),
            Rule("q3", "q3", "1", "1", Left),
            Rule("q3", "q3", "0", "0", Left),
            Rule("q1", "q4", Blank, Blank, Right),
            Rule("q4", "q4", "1", "1", Right),
            Rule("q4", "q4", "0", "0", Right),
            Rule("q4", "q5", Blank, Blank, Left),
            Rule("q5", "q8", "1", Blank, Left),
            Rule("q8", "q8", "1", "1", Left),
            Rule("q8", "q8", "0", "0", Left),
            Rule("q8", "q7", Blank, Blank, Left),
            Rule("q7", "q7", "A", "A", Left),
            Rule("q7", "q7", "B", "B", Left),
            Rule("q7", "q6", "1", "A", Left),
            Rule("q7", "q1", "0", "B", Right),
            Rule("q6", "q6", "1", "0", Left),
            Rule("q6", "q1", Blank, "1", Right),
            Rule("q6", "q1", "0", "1", Right),
            Rule("q5", "q3", "0", Blank, Left),
            Rule("q5", "q9", Blank, Blank, Left),
            Rule("q9", "q9", "B", "1", Left),
            Rule("q9", "q9", "A", "0", Left),
            Rule("q9", "halt", "0", "0", Stay),
            Rule("q9", "halt", "1", "1", Stay),
            Rule("q9", "halt", Blank, Blank, Stay)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine to salculate x mod 3, where x is in tallycode", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q2", "1", Blank, Right),
            Rule("q2", "q3", "1", Blank, Right),
            Rule("q3", "q4", Blank, "1", Left),
            Rule("q4", "halt", Blank, "1", Stay),
            Rule("q2", "halt", Blank, "1", Stay),
            Rule("q3", "q1", "1", Blank, Right),
            Rule("q1", "halt", Blank, Blank, Stay)
        ]),
        Tape, Show
    ),
    (TMProgramme("Turing Machine which accepts strings with equal numbers of zeros and ones", "q0", "halt", Blank,
        [
            Rule("q0", "q1", Blank, Blank, Right),
            Rule("q1", "q1", "B", "B", Right),
            Rule("q1", "q2", "1", "A", Right),
            Rule("q2", "q2", "B", "B", Right),
            Rule("q2", "q2", "1", "1", Right),
            Rule("q2", "q3", "0", "B", Left),
            Rule("q3", "q3", "0", "0", Left),
            Rule("q3", "q3", "1", "1", Left),
            Rule("q3", "q3", "B", "B", Left),
            Rule("q3", "q1", "A", "A", Right),
            Rule("q1", "q4", "0", "A", Right),
            Rule("q4", "q4", "B", "B", Right),
            Rule("q4", "q4", "0", "0", Right),
            Rule("q4", "q5", "1", "B", Left),
            Rule("q5", "q5", "B", "B", Left),
            Rule("q5", "q5", "0", "0", Left),
            Rule("q5", "q5", "1", "1", Left),
            Rule("q5", "q1", "A", "A", Right),
            Rule("q1", "halt", Blank, Blank, Right)
        ]),
        Tape, Show
    ),
]

# parse tape into tuple of strings.
# you can apply the string. function to this if you want to
# ensure that the elements are strings and not substrings
Tape = tuple(split(string(Tape), "")...)

function main()
    (prog, tape_tuple, verbose) = programmes[ChosenProgramme]

    # pad tape with blanks
    tape_tuple = (Blank, tape_tuple..., Blank)

    # construct dictionary from tape
    tape = Dict{Int,String}()
    entry_count = 0

    for i in tape_tuple
        entry_count += 1
        push!(tape, entry_count => string(i))
    end

    run_turing_machine(prog, tape, verbose)
    print("\n")
end

@time main()


text:
function main(r)
    x = go_and_get_num_steps(r, West)
    y = go_and_get_num_steps(r, Sud)

    go_up_and_right(r)

    start_mark(r, Sud)
    if isborder(r, Sud)
        mark_side(r, Nord)
    end
    mark_side(r, Sud)

    move_n_steps(r, Nord, y)
    move_n_steps(r, Ost, x)
end

function go_and_get_num_steps(r::Robot, side::HorizonSide) :: Int
    n = 0
    while !isborder(r, side)
        move!(r,side)
        n += 1
    end
    return n
end

function go_up_and_right(r::Robot)
    while !isborder(r, Nord)
        move!(r,Nord)
    end
    while !isborder(r, Ost)
        move!(r, Ost)
    end
end

function start_mark(r, side)
    while !isborder(r, West)
        mark_side(r, side)
        side = inverse(side)
        move!(r, West)
    end
end

function mark_side(r, side)
    while !isborder(r, side)
        putmarker!(r)
        move!(r, side)
    end
    putmarker!(r)
end

inverse(side) = HorizonSide(mod(Int(side)+2,4))

function move_n_steps(r, side, n)
    while n>0
        move!(r, side)
        n = n - 1
    end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
# A KDNode stores the information needed in each non leaf node
# to make the needed distance computations
struct KDNode{T}
    lo::T           # The low boundary for the hyper rect in this dimension
    hi::T           # The high boundary for the hyper rect in this dimension
    split_val::T    # The value the hyper rectangle was split at
    split_dim::Int  # The dimension the hyper rectangle was split at
end

struct KDTree{V <: AbstractVector,M <: MinkowskiMetric,T} <: NNTree{V,M}
    data::Vector{V}
    hyper_rec::HyperRectangle{T}
    indices::Vector{Int}
    metric::M
    nodes::Vector{KDNode{T}}
    tree_data::TreeData
    reordered::Bool
end

"""
    KDTree(data [, metric = Euclidean(), leafsize = 10]) -> kdtree

Creates a `KDTree` from the data using the given `metric` and `leafsize`.
The `metric` must be a `MinkowskiMetric`.
"""
function KDTree(data::AbstractVector{V},
                metric::M = Euclidean();
                leafsize::Int = 10,
                storedata::Bool = true,
                reorder::Bool = true,
                reorderbuffer::Vector{V} = Vector{V}()) where {V <: AbstractArray, M <: MinkowskiMetric}
    reorder = !isempty(reorderbuffer) || (storedata ? reorder : false)

    tree_data = TreeData(data, leafsize)
    n_d = length(V)
    n_p = length(data)

    indices = collect(1:n_p)
    nodes = Vector{KDNode{eltype(V)}}(undef, tree_data.n_internal_nodes)

    if reorder
        indices_reordered = Vector{Int}(undef, n_p)
        if isempty(reorderbuffer)
            data_reordered = Vector{V}(undef, n_p)
        else
            data_reordered = reorderbuffer
        end
    else
        # Dummy variables
        indices_reordered = Vector{Int}()
        data_reordered = Vector{V}()
    end

    # Create first bounding hyper rectangle that bounds all the input points
    hyper_rec = compute_bbox(data)

    # Call the recursive KDTree builder
    build_KDTree(1, data, data_reordered, hyper_rec, nodes, indices, indices_reordered,
                 1, length(data), tree_data, reorder)
    if reorder
        data = data_reordered
        indices = indices_reordered
    end

    if metric isa Distances.UnionMetrics
        p = parameters(metric)
        if p !== nothing && length(p) != length(V)
            throw(ArgumentError(
                "dimension of input points:$(length(V)) and metric parameter:$(length(p)) must agree"))
        end
    end

    KDTree(storedata ? data : similar(data, 0), hyper_rec, indices, metric, nodes, tree_data, reorder)
end

 function KDTree(data::AbstractMatrix{T},
                metric::M = Euclidean();
                leafsize::Int = 10,
                storedata::Bool = true,
                reorder::Bool = true,
                reorderbuffer::Matrix{T} = Matrix{T}(undef, 0, 0)) where {T <: AbstractFloat, M <: MinkowskiMetric}
    dim = size(data, 1)
    npoints = size(data, 2)
    points = copy_svec(T, data, Val(dim))
    if isempty(reorderbuffer)
        reorderbuffer_points = Vector{SVector{dim,T}}()
    else
        reorderbuffer_points = copy_svec(T, reorderbuffer, Val(dim))
    end
    KDTree(points, metric, leafsize = leafsize, storedata = storedata, reorder = reorder,
            reorderbuffer = reorderbuffer_points)
end

function build_KDTree(index::Int,
                      data::AbstractVector{V},
                      data_reordered::Vector{V},
                      hyper_rec::HyperRectangle,
                      nodes::Vector{KDNode{T}},
                      indices::Vector{Int},
                      indices_reordered::Vector{Int},
                      low::Int,
                      high::Int,
                      tree_data::TreeData,
                      reorder::Bool) where {V <: AbstractVector, T}
    n_p = high - low + 1 # Points left
    if n_p <= tree_data.leafsize
        if reorder
            reorder_data!(data_reordered, data, index, indices, indices_reordered, tree_data)
        end
        return
    end

    mid_idx = find_split(low, tree_data.leafsize, n_p)

    split_dim = 1
    max_spread = zero(T)
    # Find dimension and spread where the spread is maximal
    for d in 1:length(V)
        spread = hyper_rec.maxes[d] - hyper_rec.mins[d]
        if spread > max_spread
            max_spread = spread
            split_dim = d
        end
    end

    select_spec!(indices, mid_idx, low, high, data, split_dim)

    split_val = data[indices[mid_idx]][split_dim]

    lo = hyper_rec.mins[split_dim]
    hi = hyper_rec.maxes[split_dim]

    nodes[index] = KDNode{T}(lo, hi, split_val, split_dim)

    # Call the left sub tree with an updated hyper rectangle
    hyper_rec.maxes[split_dim] = split_val
    build_KDTree(getleft(index), data, data_reordered, hyper_rec, nodes,
                  indices, indices_reordered, low, mid_idx - 1, tree_data, reorder)
    hyper_rec.maxes[split_dim] = hi # Restore the hyper rectangle

    # Call the right sub tree with an updated hyper rectangle
    hyper_rec.mins[split_dim] = split_val
    build_KDTree(getright(index), data, data_reordered, hyper_rec, nodes,
                  indices, indices_reordered, mid_idx, high, tree_data, reorder)
    # Restore the hyper rectangle
    hyper_rec.mins[split_dim] = lo
end

function _knn(tree::KDTree,
              point::AbstractVector,
              best_idxs::Vector{Int},
              best_dists::Vector,
              skip::Function)

    init_min = get_min_distance(tree.hyper_rec, point)
    knn_kernel!(tree, 1, point, best_idxs, best_dists, init_min, skip)
    @simd for i in eachindex(best_dists)
        @inbounds best_dists[i] = eval_end(tree.metric, best_dists[i])
    end
    return
end

function knn_kernel!(tree::KDTree{V},
                        index::Int,
                        point::AbstractVector,
                        best_idxs::Vector{Int},
                        best_dists::Vector,
                        min_dist,
                        skip::F) where {V, F}
    # At a leaf node. Go through all points in node and add those in range
    if isleaf(tree.tree_data.n_internal_nodes, index)
        add_points_knn!(best_dists, best_idxs, tree, index, point, false, skip)
        return
    end

    node = tree.nodes[index]
    p_dim = point[node.split_dim]
    split_val = node.split_val
    lo = node.lo
    hi = node.hi
    split_diff = p_dim - split_val
    M = tree.metric
    # Point is to the right of the split value
    if split_diff > 0
        close = getright(index)
        far = getleft(index)
        ddiff = max(zero(eltype(V)), p_dim - hi)
    else
        close = getleft(index)
        far = getright(index)
        ddiff = max(zero(eltype(V)), lo - p_dim)
    end
    # Always call closer sub tree
    knn_kernel!(tree, close, point, best_idxs, best_dists, min_dist, skip)

    split_diff_pow = eval_pow(M, split_diff)
    ddiff_pow = eval_pow(M, ddiff)
    diff_tot = eval_diff(M, split_diff_pow, ddiff_pow)
    new_min = eval_reduce(M, min_dist, diff_tot)
    if new_min < best_dists[1]
        knn_kernel!(tree, far, point, best_idxs, best_dists, new_min, skip)
    end
    return
end

function _inrange(tree::KDTree,
                  point::AbstractVector,
                  radius::Number,
                  idx_in_ball = Int[])
    init_min = get_min_distance(tree.hyper_rec, point)
    inrange_kernel!(tree, 1, point, eval_op(tree.metric, radius, zero(init_min)), idx_in_ball,
                   init_min)
    return
end

# Explicitly check the distance between leaf node and point while traversing
function inrange_kernel!(tree::KDTree,
                         index::Int,
                         point::AbstractVector,
                         r::Number,
                         idx_in_ball::Vector{Int},
                         min_dist)
    # Point is outside hyper rectangle, skip the whole sub tree
    if min_dist > r
        return
    end

    # At a leaf node. Go through all points in node and add those in range
    if isleaf(tree.tree_data.n_internal_nodes, index)
        add_points_inrange!(idx_in_ball, tree, index, point, r, false)
        return
    end

    node = tree.nodes[index]
    split_val = node.split_val
    lo = node.lo
    hi = node.hi
    p_dim = point[node.split_dim]
    split_diff = p_dim - split_val
    M = tree.metric

    if split_diff > 0 # Point is to the right of the split value
        close = getright(index)
        far = getleft(index)
        ddiff = max(zero(p_dim - hi), p_dim - hi)
    else # Point is to the left of the split value
        close = getleft(index)
        far = getright(index)
        ddiff = max(zero(lo - p_dim), lo - p_dim)
    end
    # Call closer sub tree
    inrange_kernel!(tree, close, point, r, idx_in_ball, min_dist)

    # TODO: We could potentially also keep track of the max distance
    # between the point and the hyper rectangle and add the whole sub tree
    # in case of the max distance being <= r similarly to the BallTree inrange method.
    # It would be interesting to benchmark this on some different data sets.

    # Call further sub tree with the new min distance
    split_diff_pow = eval_pow(M, split_diff)
    ddiff_pow = eval_pow(M, ddiff)
    diff_tot = eval_diff(M, split_diff_pow, ddiff_pow)
    new_min = eval_reduce(M, min_dist, diff_tot)
    inrange_kernel!(tree, far, point, r, idx_in_ball, new_min)
end


text:
export newton, lagrange

newton(f, x₀::Real; kwargs...) = newton_impl(f, x₀, derivative, abs; kwargs...)
newton(f, x₀::SReal; kwargs...) = newton_impl(f, x₀, jacobian, norm; kwargs...)

function newton_impl(f, x₀,
                     jfun,  # jacobian function
                     nfun;  # norm function
                     rtol::Real=τᵣ,
                     atol::Real=τₐ,
                     maxit=100,  # maximum iterations
                     maxit_ls=20)  # maximum iterations inside the line search
    # Solve f(x) = 0 using the Newton-Armijo method.
    # x₀: initial guess
    # f: function of x

    isconverged = true  # true if solution converged; false otherwise
    α = 1e-4
    ∆max = 1 / rtol  # maximum step size
    perturb_ls = atol^0.75  # ≈ 1e-12, between sqrt(eps) and eps; some perturbation allowed in line search

    # Initialize.
    n = 0
    xₙ = x₀
    fₙ = f(xₙ)
    τf = rtol*nfun(fₙ) + atol
    rx₀ = nfun(x₀)

    # Perform the Newton method.
    while nfun(fₙ) > τf
        λ = 1.
        n_ls = 0  # line search iteration counter

        f′ₙ = jfun(f, xₙ)

        local ∆x  # see https://discourse.julialang.org/t/cumbersome-scoping-rules-for-try-catch-finally/4582
        try
            ∆x = -f′ₙ \ fₙ
        catch err
            isa(err,SingularException) && @error "For xₙ = $xₙ at n = $n, f′ₙ = $f′ₙ is singular."
            break
        end

        # Avoid too large Newton steps.
        nfun(∆x) ≤ ∆max || (isconverged = false; break)
        xₙ₊₁ = xₙ + λ*∆x
        fₙ₊₁ = f(xₙ₊₁)

        # Perform the line search to determine λ.  The stopping criterion does not
        # have perturb_ls ≈ 1e-12 on the RHS, but I guess this kind of perturbation
        # allows update in xₙ even in the situation where line search is supposed
        # to fail.
        while nfun(fₙ₊₁) ≥ (1 - α*λ) * nfun(fₙ) + perturb_ls
            λ /= 2
            xₙ₊₁ = xₙ + λ*∆x
            fₙ₊₁ = f(xₙ₊₁)
            n_ls += 1

            # Too many iteration steps in line search
            n_ls ≤ maxit_ls || (isconverged = false; break)
        end

        # Step accepted; continue the Newton method.
        xₙ = xₙ₊₁
        n += 1

        # Too many iteration steps in Newton's method.
        n ≤ maxit || (isconverged = false; break)
        # n ≤ maxit || throw(ErrorException("Newton method fails to converge in $n iteration steps."))

        fₙ = f(xₙ)
    end
    isnan(nfun(fₙ)) && @error "For xₙ = $xₙ at n = $n, f(xₙ) = $fₙ has NaN."

    return (sol=xₙ, converged=isconverged)  # named tuple
end

# function newton_impl(f, x₀,
#                      jfun,  # jacobian function
#                      nfun;  # norm function
#                      rtol::Real=τᵣ,
#                      atol::Real=τₐ,
#                      maxit=100,  # maximum iterations
#                      maxit_ls=20)  # maximum iterations inside the line search
#     # Solve f(x) = 0 using the Newton-Armijo method.
#     # x₀: initial guess
#     # f: function of x
#     K = Size(x₀)[1]
#
#     isconverged = true  # true if solution converged; false otherwise
#     α = 1e-4
#     ∆max = 1 / rtol  # maximum step size
#     perturb_ls = atol^0.75  # ≈ 1e-12, between sqrt(eps) and eps; some perturbation allowed in line search
#
#     # Initialize.
#     n = 0
#     xₙ = x₀
#     fₙ = f(xₙ)
#     τf = rtol*nfun(fₙ) + atol
#     rx₀ = nfun(x₀)
#
#     # Perform the Newton method.
#     while nfun(fₙ) > τf
#         λ = 1.
#         n_ls = 0  # line search iteration counter
#
#         f′ₙ = jfun(f, xₙ)
#         bool_nz = SVec(ntuple(k->!iszero(f′ₙ[:,k]), Val(K)))
#         M = sum(bool_nz)
#         ind_nz = SVec{M}((1:K)[bool_nz])
#
#         local ∆x  # see https://discourse.julialang.org/t/cumbersome-scoping-rules-for-try-catch-finally/4582
#         try
#             ∆x = -f′ₙ[:,ind_nz] \ fₙ
#         catch err
#             isa(err,SingularException) && @error "For xₙ = $xₙ at n = $n, f′ₙ = $f′ₙ is singular."
#             break
#         end
#
#         # Avoid too large Newton steps.
#         nfun(∆x) ≤ ∆max || (isconverged = false; break)
#         # xₙ₊₁ = xₙ + λ*∆x
#         xₙ₊₁′ = MVec(xₙ)
#         xₙ₊₁′[ind_nz] .+= λ .* ∆x
#         xₙ₊₁ = SVec(xₙ₊₁′)
#         fₙ₊₁ = f(xₙ₊₁)
#
#         # # Perform the line search to determine λ.  The stopping criterion does not
#         # # have perturb_ls ≈ 1e-12 on the RHS, but I guess this kind of perturbation
#         # # allows update in xₙ even in the situation where line search is supposed
#         # # to fail.
#         # while nfun(fₙ₊₁) ≥ (1 - α*λ) * nfun(fₙ) + perturb_ls
#         #     λ /= 2
#         #     xₙ₊₁ = xₙ + λ*∆x
#         #     fₙ₊₁ = f(xₙ₊₁)
#         #     n_ls += 1
#         #
#         #     # Too many iteration steps in line search
#         #     n_ls ≤ maxit_ls || (isconverged = false; break)
#         # end
#
#         # Step accepted; continue the Newton method.
#         xₙ = xₙ₊₁
#         n += 1
#
#         # Too many iteration steps in Newton's method.
#         n ≤ maxit || (isconverged = false; break)
#         # n ≤ maxit || throw(ErrorException("Newton method fails to converge in $n iteration steps."))
#
#         fₙ = f(xₙ)
#     end
#     isnan(nfun(fₙ)) && @error "For xₙ = $xₙ at n = $n, f(xₙ) = $fₙ has NaN."
#
#     return (sol=xₙ, converged=isconverged)  # named tuple
# end

# Need to use different initial guesses for finding the left and right boundaries.
function lagrange(f_obj, f_con, x₀::SReal{K}; kwargs...) where {K}
    ind_x = SVec(ntuple(identity, Val(K)))

    y₀ = SVec(x₀.data..., 1)
    L(y) = f_obj(y[ind_x]) - y[K+1] * f_con(y[ind_x])
    ∇L(y) = gradient(L, y)

    res = newton(∇L, y₀; kwargs...)
    yₛ = res.sol
    converged = res.converged

    xₛ = yₛ[ind_x]
    λₛ = yₛ[K+1]
    Lₛ = L(yₛ)

    return (sol=xₛ, λ=λₛ, L=Lₛ, converged=converged)
end

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
@testset "Linear trend" begin
    @testset "Linear trend model" begin
        y = collect(1:15.0)

        unimodel1 = linear_trend(y)
        @test isa(unimodel1, StateSpaceModel)
        @test unimodel1.mode == "time-invariant"

        ss1 = statespace(unimodel1)

        unimodel2 = linear_trend(y)
        @test ss1.filter_type <: KalmanFilter
        @test isa(ss1, StateSpace)
        @test ss1.smoother.alpha[:, 2] ≈ ones(15) rtol = 1e-4
        compare_forecast_simulation(ss1, 20, 1000, 1e-3)

        ss2 = statespace(unimodel2; filter_type = SquareRootFilter{Float64})

        @test ss2.filter_type <: SquareRootFilter
        @test isa(ss2, StateSpace)
        @test ss2.smoother.alpha[:, 2] ≈ ones(15) rtol = 1e-4
        compare_forecast_simulation(ss2, 20, 1000, 1e-3)

        unimodel3 = linear_trend(y)

        ss3 = statespace(unimodel3; filter_type = UnivariateKalmanFilter{Float64})

        @test ss3.filter_type <: UnivariateKalmanFilter
        @test isa(ss3, StateSpace)
        @test ss3.smoother.alpha[:, 2] ≈ ones(15) rtol = 1e-4
        compare_forecast_simulation(ss3, 20, 1000, 1e-3)

        @test ss2.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-3
        @test ss3.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-3
        @test ss2.filter.a ≈ ss1.filter.a rtol = 1e-3
        @test ss3.filter.a ≈ ss1.filter.a rtol = 1e-3

        diags1 = diagnostics(ss1)
        diags2 = diagnostics(ss2)
        diags3 = diagnostics(ss3)

        @test diags1.p_jarquebera ≈ diags2.p_jarquebera atol = 1e-2
        @test diags1.p_jarquebera ≈ diags3.p_jarquebera atol = 1e-2
    end

    @testset "Linear trend model with missing values" begin
        y = collect(1:30.0)
        y[4:8] .= NaN
        model = linear_trend(y)
        ss = statespace(model)

        @test y[9:end] ≈ ss.smoother.alpha[9:end, 1] rtol = 1e-3
        @test sum(ss.filter.P[:, :, 8]) > sum(ss.filter.P[:, :, 7]) > sum(ss.filter.P[:, :, 6]) >
                sum(ss.filter.P[:, :, 5]) > sum(ss.filter.P[:, :, 4])
    end
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 UpgradePolicy <: SwaggerModel
    mode::Any # spec type: Union{ Nothing, String } # spec name: mode
    rollingUpgradePolicy::Any # spec type: Union{ Nothing, RollingUpgradePolicy } # spec name: rollingUpgradePolicy
    automaticOSUpgradePolicy::Any # spec type: Union{ Nothing, AutomaticOSUpgradePolicy } # spec name: automaticOSUpgradePolicy

    function UpgradePolicy(;mode=nothing, rollingUpgradePolicy=nothing, automaticOSUpgradePolicy=nothing)
        o = new()
        validate_property(UpgradePolicy, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(UpgradePolicy, Symbol("rollingUpgradePolicy"), rollingUpgradePolicy)
        setfield!(o, Symbol("rollingUpgradePolicy"), rollingUpgradePolicy)
        validate_property(UpgradePolicy, Symbol("automaticOSUpgradePolicy"), automaticOSUpgradePolicy)
        setfield!(o, Symbol("automaticOSUpgradePolicy"), automaticOSUpgradePolicy)
        o
    end
end # type UpgradePolicy

const _property_map_UpgradePolicy = Dict{Symbol,Symbol}(Symbol("mode")=>Symbol("mode"), Symbol("rollingUpgradePolicy")=>Symbol("rollingUpgradePolicy"), Symbol("automaticOSUpgradePolicy")=>Symbol("automaticOSUpgradePolicy"))
const _property_types_UpgradePolicy = Dict{Symbol,String}(Symbol("mode")=>"String", Symbol("rollingUpgradePolicy")=>"RollingUpgradePolicy", Symbol("automaticOSUpgradePolicy")=>"AutomaticOSUpgradePolicy")
Base.propertynames(::Type{ UpgradePolicy }) = collect(keys(_property_map_UpgradePolicy))
Swagger.property_type(::Type{ UpgradePolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpgradePolicy[name]))}
Swagger.field_name(::Type{ UpgradePolicy }, property_name::Symbol) =  _property_map_UpgradePolicy[property_name]

const _allowed_UpgradePolicy_mode = ["Automatic", "Manual", "Rolling"]

function check_required(o::UpgradePolicy)
    true
end

function validate_property(::Type{ UpgradePolicy }, name::Symbol, val)
    if name === Symbol("mode")
        Swagger.validate_param(name, "UpgradePolicy", :enum, val, _allowed_UpgradePolicy_mode)
    end
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# simplelazy.jl
#
# Solve a simple problem using lazy constraints. See documentation for more
# details about callbacks.
#############################################################################

using JuMP
using Gurobi

# We will use Gurobi
m = Model(solver=GurobiSolver())

# Define our variables to be inside a box, and integer
@variable(m, 0 <= x <= 2, Int)
@variable(m, 0 <= y <= 2, Int)

@objective(m, Max, y)

# We now define our callback function that takes one argument,
# the callback handle. Note that we can access m, x, and y because
# this function is defined inside the same scope
function corners(cb)
    x_val = getvalue(x)
    y_val = getvalue(y)
    println("In callback function, x=$x_val, y=$y_val")

    # We have two constraints, one cutting off the top
    # left corner and one cutting off the top right corner, e.g.
    # (0,2) +---+---+ (2,2)
    #       |xx/ \xx|
    #       |x/   \x|
    #       |/     \|
    #       +       +
    #       |       |
    #       |       |
    #       |       |
    # (0,0) +---+---+ (2,0)

    # Allow for some impreciseness in the solution
    TOL = 1e-6

    # Check top left, allowing some tolerance
    if y_val - x_val > 1 + TOL
        # Cut off this solution
        println("Solution was in top left, cut it off")
        # Use the original variables
        @lazyconstraint(cb, y - x <= 1)
    # Check top right
    elseif y_val + x_val > 3 + TOL
        # Cut off this solution
        println("Solution was in top right, cut it off")
        # Use the original variables
        @lazyconstraint(cb, y + x <= 3)
    end
end  # End of callback function

# Tell JuMP/Gurobi to use our callback function
addlazycallback(m, corners)

# Solve the problem
solve(m)

# Print our final solution
println("Final solution: [ $(getvalue(x)), $(getvalue(y)) ]")