text:
# Explicit Runge-Kutta solvers
##############################
# (Hairer & Wanner 1992 p.134, p.165-169)

###########################################
# Tableaus for explicit Runge-Kutta methods
###########################################

struct TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T}
    order::(Tuple{Vararg{Int}}) # the order of the methods
    a::Matrix{T}
    # one or several row vectors.  First row is used for the step,
    # second for error calc.
    b::Matrix{T}
    c::Vector{T}
    function TableauRKExplicit{Name, S, T}(order,a,b,c) where {Name, S, T}
        @assert isa(S,Integer)
        @assert isa(Name,Symbol)
        @assert c[1]==0
        @assert istril(a)
        @assert S==length(c)==size(a,1)==size(a,2)==size(b,2)
        @assert size(b,1)==length(order)
        # consistency:
        if T<:AbstractFloat
            @assert norm(sum(a,dims=2)-c'',Inf) < 100*eps(T)
        else
            @assert norm(sum(a,dims=2)-c'',Inf) < 100*eps(Float64)
        end
        new{Name, S, T}(order,a,b,c)
    end
end
function TableauRKExplicit(name::Symbol, order::(Tuple{Vararg{Int}}),
                   a::Matrix{T}, b::Matrix{T}, c::Vector{T}) where T
    TableauRKExplicit{name,length(c),T}(order, a, b, c)
end
function TableauRKExplicit(name::Symbol, order::(Tuple{Vararg{Int}}), T::Type,
                   a::Matrix, b::Matrix, c::Vector)
    TableauRKExplicit{name,length(c),T}(order, convert(Matrix{T},a),
                                        convert(Matrix{T},b), convert(Vector{T},c) )
end
conv_field(D,a::Array{T,N}) where {T,N} = convert(Array{D,N}, a)
function Base.convert(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T}) where {Tnew<:Real,Name,S,T}
    # Converts the tableau coefficients to the new type Tnew
    newflds = ()
    for n in fieldnames(typeof(tab))
        fld = getfield(tab,n)
        if eltype(fld)==T
            newflds = tuple(newflds..., conv_field(Tnew, fld))
        else
            newflds = tuple(newflds..., fld)
        end
    end
    TableauRKExplicit{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way?
end

isexplicit(b::TableauRKExplicit) = istril(b.a) # Test whether it's an explicit method
isadaptive(b::TableauRKExplicit) = size(b.b, 1)==2

# First same as last.  Means ks[:,end]=ks_nextstep[:,1], c.f. H&W p.167
isFSAL(btab::TableauRKExplicit) = btab.a[end,:]==btab.b[1,:] && btab.c[end]==1 # the latter is not needed really

## Tableaus for explicit RK methods
# Fixed step:
const bt_feuler = TableauRKExplicit(:feuler,(1,), Rational{Int64},
                                    zeros(Int,1,1),
                                    [1][:,:],
                                    [0]
                                    )
const bt_midpoint = TableauRKExplicit(:midpoint,(2,), Rational{Int64},
                                      [0  0
                                       1//2  0],
                                      [0  1],
                                      [0, 1//2]
                                      )
const bt_heun = TableauRKExplicit(:heun,(2,), Rational{Int64},
                                  [0  0
                                   1  0],
                                  [1//2  1//2],
                                  [0, 1])

const bt_rk4 = TableauRKExplicit(:rk4,(4,),Rational{Int64},
                                 [0    0    0 0
                                  1//2 0    0 0
                                  0    1//2 0 0
                                  0    0    1 0],
                                 [1//6  1//3  1//3  1//6],
                                 [0, 1//2, 1//2, 1])

# Adaptive step:
# Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods
const bt_rk21 = TableauRKExplicit(:heun_euler,(2,1), Rational{Int64},
                                  [0     0
                                   1     0],
                                  [1//2  1//2
                                   1     0],
                                  [0,    1])

# Bogacki–Shampine coefficients
const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64},
                                  [0           0      0      0
                                   1//2         0      0      0
                                   0         3//4      0      0
                                   2//9       1//3     4//9     0],
                                  [7//24 1//4 1//3 1//8
                                   2//9 1//3 4//9 0],
                                  [0, 1//2, 3//4, 1]
                         )

# Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
const bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Rational{Int64},
                          [  0           0           0            0         0     0
                             1//4        0           0            0         0     0
                             3//32       9//32       0            0         0     0
                          1932//2197 -7200//2197  7296//2197      0         0     0
                           439//216     -8        3680//513    -845//4104   0     0
                            -8//27       2       -3544//2565   1859//4104 -11//40 0 ],
                           [25//216      0        1408//2565   2197//4104  -1//5  0
                            16//135      0        6656//12825 28561//56430 -9//50 2//55],
                            [0,          1//4,       3//8,       12//13,    1,    1//2])

# Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method
const bt_dopri5 = TableauRKExplicit(:dopri, (5,4), Rational{Int64},
                     [0           0            0                   0            0      0 0
                      1//5        0            0                   0            0      0 0
                      3//40       9//40        0                   0            0      0 0
                      44//45      -56//15      32//9               0            0      0 0
                      19372//6561 -25360//2187 64448//6561 -212//729            0      0 0
                      9017//3168  -355//33     46732//5247   49//176 -5103//18656      0 0
                      35//384     0            500//1113    125//192  -2187//6784 11//84 0],
                     [35//384         0     500//1113       125//192      -2187//6784         11//84      0
                      5179//57600     0     7571//16695     393//640     -92097//339200     187//2100     1//40],
                     [0, 1//5, 3//10, 4//5, 8//9, 1, 1]
                     )

# Fehlberg 7(8) coefficients
# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control".
# National Aeronautics and Space Administration.
const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64},
                            [     0       0      0       0         0          0        0        0      0       0     0 0 0
                                  2//27   0      0       0         0          0        0        0      0       0     0 0 0
                                  1//36   1//12  0       0         0          0        0        0      0       0     0 0 0
                                  1//24   0      1//8    0         0          0        0        0      0       0     0 0 0
                                  5//12   0    -25//16  25//16     0          0        0        0      0       0     0 0 0
                                  1//20   0      0       1//4      1//5       0        0        0      0       0     0 0 0
                                -25//108  0      0     125//108  -65//27    125//54    0        0      0       0     0 0 0
                                 31//300  0      0       0       61//225    -2//9     13//900   0      0       0     0 0 0
                                  2       0      0     -53//6    704//45   -107//9    67//90    3      0       0     0 0 0
                                -91//108  0      0      23//108 -976//135   311//54  -19//60   17//6  -1//12   0     0 0 0
                               2383//4100 0      0    -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0
                                  3//205  0      0       0        0        -6//41     -3//205  -3//41  3//41   6//41 0 0 0
                              -1777//4100 0      0    -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0],
                              [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0
                               0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0     41//840 41//840],
                               [0,    2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1]
                            )

################################
# Fixed step Runge-Kutta methods
################################

# TODO: iterator method
ode1(fn, y0, tspan; kwargs...) = oderk_fixed(fn, y0, tspan, bt_feuler; kwargs...)
ode2_midpoint(fn, y0, tspan; kwargs...) = oderk_fixed(fn, y0, tspan, bt_midpoint; kwargs...)
ode2_heun(fn, y0, tspan; kwargs...) = oderk_fixed(fn, y0, tspan, bt_heun; kwargs...)
ode4(fn, y0, tspan; kwargs...) = oderk_fixed(fn, y0, tspan, bt_rk4; kwargs...)

function oderk_fixed(fn, y0, tspan, btab::TableauRKExplicit; kwargs...)
    # Non-arrays y0 treat as scalar
    fn_(t, y) = [fn(t, y[1])]
    t,y = oderk_fixed(fn_, [y0], tspan, btab)
    return t, vcat_nosplat(y)
end
function oderk_fixed(fn, y0::AbstractVector, tspan,
                     btab_::TableauRKExplicit{N,S}, alias_u0=false) where {N,S}
    # TODO: instead of AbstractVector use a Holy-trait

    # Needed interface:
    # On components:
    # On y0 container: length, deepcopy, similar, setindex!
    # On time container: getindex, convert. length

    Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_)
    dof = length(y0)

    ys = Vector{Ty}(undef,length(tspan))
    allocate!(ys, y0, dof)
    if alias_u0
        ys[1] = y0
    else
        ys[1] = deepcopy(y0)
    end

    tspan = convert(Vector{Et}, tspan)
    # work arrays:
    ks = Vector{Ty}(undef,S)
    # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place
    ytmp = similar(y0, Eyf, dof)
    for i=1:length(tspan)-1
        dt = tspan[i+1]-tspan[i]
        ys[i+1][:] = ys[i]
        for s=1:S
            calc_next_k!(ks, ytmp, ys[i], s, fn, tspan[i], dt, dof, btab)
            for d=1:dof
                ys[i+1][d] += dt * btab.b[s]*ks[s][d]
            end
        end
    end
    return tspan, ys
end

##############################
# Adaptive Runge-Kutta methods
##############################

ode21(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk21; kwargs...)
ode23(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk23; kwargs...)
ode45_fe(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk45; kwargs...)
ode45_dp(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_dopri5; kwargs...)
# Use Dormand-Prince version of ode45 by default
ode45(fn, y0, tspan; kwargs...) = ode45_dp(fn, y0, tspan; kwargs...)
ode78(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_feh78; kwargs...)

function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...)
    # For y0 which don't support indexing.
    fn_ = (t, y) -> [fn(t, y[1])]
    t,y = oderk_adapt(fn_, [y0], tspan, btab; kwords...)
    return t, vcat_nosplat(y)
end
function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S};
                     reltol = 1.0e-5, abstol = 1.0e-8,
                     norm=LinearAlgebra.norm,
                     minstep=abs(tspan[end] - tspan[1])/1e18,
                     maxstep=abs(tspan[end] - tspan[1])/2.5,
                     initstep=0,
                     points=:all,
                     alias_u0=false
                     ) where {N,S}
    # Needed interface:
    # On components:
    #  - note that the type of the components might change!
    # On y0 container: length, similar, setindex!
    # On time container: getindex, convert, length

    # For y0 which support indexing.  Currently y0<:AbstractVector but
    # that could be relaxed with a Holy-trait.
    !isadaptive(btab_) && error("Can only use this solver with an adaptive RK Butcher table")
    Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_)
    # parameters
    order = minimum(btab.order)
    timeout_const = 5 # after step reduction do not increase step for
                      # timeout_const steps

    ## Initialization
    dof = length(y0)
    tspan = convert(Vector{Et}, tspan)
    t = tspan[1]
    tstart = tspan[1]
    tend = tspan[end]

    # work arrays:
    if alias_u0
        y = y0                     # y at time t
    else
        y  = similar(y0, Eyf, dof) # y at time t
    end
    y[:]   = y0
    ytrial = similar(y0, Eyf, dof) # trial solution at time t+dt
    yerr   = similar(y0, Eyf, dof) # error of trial solution
    ks = Vector{Ty}(undef,S)
    # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place
    ytmp   = similar(y0, Eyf, dof)

    # output ys
    nsteps_fixed = length(tspan) # these are always output
    ys = Vector{Ty}(undef,nsteps_fixed)
    allocate!(ys, y0, dof)
    ys[1] = y0

    # Option points determines where solution is returned:
    if points == :all
        tspan_fixed = tspan
        tspan = Et[tstart]
        iter_fixed = 2 # index into tspan_fixed
        sizehint!(tspan, nsteps_fixed)
    elseif points != :specified
        error("Unrecognized option points==$points")
    end
    # Time
    dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0
    if initstep!=0
        dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.")
    end
    # Diagnostics
    dts = Et[]
    errs = Float64[]
    steps = [0,0]  # [accepted, rejected]

    ## Integration loop
    islaststep = abs(t+dt-tend)<=eps(tend) ? true : false
    timeout = 0 # for step-control
    iter = 2 # the index into tspan and ys
    while true
        # do one step (assumes ks[1]==f0)
        rk_embedded_step!(ytrial, yerr, ks, ytmp, y, fn, t, dt, dof, btab)
        # Check error and find a new step size:
        err, newdt, timeout = stepsize_hw92!(dt, tdir, y, ytrial, yerr, order, timeout,
                                            dof, abstol, reltol, maxstep, norm)

        if err<=1 # accept step
            # diagnostics
            steps[1] +=1
            push!(dts, dt)
            push!(errs, err)

            # Output:
            f0 = ks[1]
            f1 = isFSAL(btab) ? ks[S] : fn(t+dt, ytrial)
            if points==:specified
                # interpolate onto given output points
                while iter-1<nsteps_fixed && (tdir*tspan[iter]<tdir*(t+dt) || islaststep) # output at all new times which are < t+dt
                    hermite_interp!(ys[iter], tspan[iter], t, dt, y, ytrial, f0, f1) # TODO: 3rd order only!
                    iter += 1
                end
            else
                # first interpolate onto given output points
                while iter_fixed-1<nsteps_fixed && tdir*t<tdir*tspan_fixed[iter_fixed]<tdir*(t+dt) # output at all new times which are < t+dt
                    yout = hermite_interp(tspan_fixed[iter_fixed], t, dt, y, ytrial, f0, f1)
                    index_or_push!(ys, iter, yout) # TODO: 3rd order only!
                    push!(tspan, tspan_fixed[iter_fixed])
                    iter_fixed += 1
                    iter += 1
                end
                # but also output every step taken
                index_or_push!(ys, iter, deepcopy(ytrial))
                push!(tspan, t+dt)
                iter += 1
            end
            ks[1] = f1 # load ks[1]==f0 for next step

            # Break if this was the last step:
            islaststep && break

            # Swap bindings of y and ytrial, avoids one copy
            y, ytrial = ytrial, y

            # Update t to the time at the end of current step:
            t += dt
            dt = newdt

            # Hit end point exactly if next step within 1% of end:
            if tdir*(t+dt+dt/100) >= tdir*tend
                dt = tend-t
                islaststep = true # next step is the last, if it succeeds
            end
        elseif abs(newdt)<minstep  # minimum step size reached, break
            println("Warning: dt < minstep.  Stopping.")
            break
        else # redo step with smaller dt
            islaststep = false
            steps[2] +=1
            dt = newdt
            timeout = timeout_const
        end
    end
    return tspan, ys
end

function rk_embedded_step!(ytrial, yerr, ks, ytmp, y, fn, t, dt, dof, btab::TableauRKExplicit{N,S}) where {N,S}
    # Does one embedded R-K step updating ytrial, yerr and ks.
    #
    # Assumes that ks[:,1] is already calculated!
    #
    # Modifies ytrial, yerr, ks, and ytmp

    # Needed interface:
    # On components: arithmetic, zero
    # On y0 container: fill!, setindex!, getindex

    fill!(ytrial, zero(eltype(ytrial)) )
    fill!(yerr, zero(eltype(ytrial)) )
    for d=1:dof
        ytrial[d] += btab.b[1,1]*ks[1][d]
        yerr[d]   += btab.b[2,1]*ks[1][d]
    end
    for s=2:S
        calc_next_k!(ks, ytmp, y, s, fn, t, dt, dof, btab)
        for d=1:dof
            ytrial[d] += btab.b[1,s]*ks[s][d]
            yerr[d]   += btab.b[2,s]*ks[s][d]
        end
    end
    for d=1:dof
        yerr[d]   = dt * (ytrial[d]-yerr[d])
        ytrial[d] = y[d] + dt * ytrial[d]
    end
end

function stepsize_hw92!(dt::T, tdir, x0, xtrial, xerr, order,
                    timeout, dof, abstol, reltol, maxstep, norm) where T
    # Estimates the error and a new step size following Hairer &
    # Wanner 1992, p167 (with some modifications)
    #
    # If timeout>0 no step size increase is allowed, timeout is
    # decremented in here.
    #
    # Returns the error, newdt and the number of timeout-steps
    #
    # The work-array xerr is modified in-place
    #
    # TODO:
    # - allow component-wise reltol and abstol?
    # - allow other norms

    # Needed interface:
    # On components: isoutofdomain
    # On y0 container: norm, get/setindex

    timout_after_nan = 5
    fac = T[0.8, 0.9, 0.25^(1/(order+1)), 0.38^(1/(order+1))][1]
    facmax = T(5) # maximal step size increase. 1.5-5
    facmin = T(1)/T(facmax)  # maximal step size decrease. ?
    maxstep = T(maxstep)
    abstol = T(abstol)
    reltol = T(reltol)

    # in-place calculate xerr./tol
    for d=1:dof
        # if outside of domain (usually NaN) then make step size smaller by maximum
        isoutofdomain(xtrial[d]) && return T(10), dt*facmin, timout_after_nan
        xerr[d] = xerr[d]/(abstol + max(norm(x0[d]), norm(xtrial[d]))*reltol) # Eq 4.10
    end
    err = norm(xerr, 2) # Eq. 4.11
    newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(T(1/(order+1)))) ) # Eq 4.13 modified
    if timeout>0
        newdt = min(newdt, dt)
        timeout -= 1
    end
    return err, tdir*newdt, timeout
end

function calc_next_k!(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) where Ty
    # Calculates the next ks and puts it into ks[s]
    # - ks and ytmp are modified inside this function.

    # Needed interface:
    # On components: +, *
    # On y0 container: setindex!, getindex, fn

    ytmp[:] = y
    for ss=1:s-1, d=1:dof
        ytmp[d] += dt * ks[ss][d] * btab.a[s,ss]
    end
    ks[s] = fn(t + btab.c[s]*dt, ytmp)::Ty
    nothing
end

# Helper functions:
function allocate!(vec::Vector{T}, y0, dof) where T
    # Allocates all vectors inside a Vector{Vector} using the same
    # kind of container as y0 has and element type eltype(eltype(vec)).
    for s=1:length(vec)
        vec[s] = similar(y0, eltype(T), dof)
    end
end
function index_or_push!(vec, i, val)
    # Fills in the vector until there is no space, then uses push!
    # instead.
    if length(vec)>=i
        vec[i] = val
    else
        push!(vec, val)
    end
end
vcat_nosplat(y) =  eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting

function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1)
    # For dense output see Hairer & Wanner p.190 using Hermite
    # interpolation. Updates y in-place.
    #
    # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 )
    # this is O(3). TODO for higher order.

    theta = (tquery-t)/dt
    for i=1:length(y0)
        y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) *
                ((1-2*theta)*(y1[i]-y0[i]) + (theta-1)*dt*f0[i] + theta*dt*f1[i]) )
    end
    nothing
end
function hermite_interp(tquery,t,dt,y0,y1,f0,f1)
    # Returns the y instead of in-place
    y = similar(y0)
    hermite_interp!(y,tquery,t,dt,y0,y1,f0,f1)
    return y
end


text:
function ParameterDriver(arg0::JString, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return ParameterDriver((JString, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function add_observer(obj::ParameterDriver, arg0::ParameterObserver)
    return jcall(obj, "addObserver", void, (ParameterObserver,), arg0)
end

function get_max_value(obj::ParameterDriver)
    return jcall(obj, "getMaxValue", jdouble, ())
end

function get_min_value(obj::ParameterDriver)
    return jcall(obj, "getMinValue", jdouble, ())
end

function get_name(obj::ParameterDriver)
    return jcall(obj, "getName", JString, ())
end

function get_normalized_value(obj::ParameterDriver)
    return jcall(obj, "getNormalizedValue", jdouble, ())
end

function get_observers(obj::ParameterDriver)
    return jcall(obj, "getObservers", List, ())
end

function get_reference_date(obj::ParameterDriver)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function get_reference_value(obj::ParameterDriver)
    return jcall(obj, "getReferenceValue", jdouble, ())
end

function get_scale(obj::ParameterDriver)
    return jcall(obj, "getScale", jdouble, ())
end

function get_value(obj::ParameterDriver)
    return jcall(obj, "getValue", jdouble, ())
end

function get_value(obj::ParameterDriver, arg0::jint, arg1::Map)
    return jcall(obj, "getValue", Gradient, (jint, Map), arg0, arg1)
end

function is_selected(obj::ParameterDriver)
    return jcall(obj, "isSelected", jboolean, ())
end

function remove_observer(obj::ParameterDriver, arg0::ParameterObserver)
    return jcall(obj, "removeObserver", void, (ParameterObserver,), arg0)
end

function replace_observer(obj::ParameterDriver, arg0::ParameterObserver, arg1::ParameterObserver)
    return jcall(obj, "replaceObserver", void, (ParameterObserver, ParameterObserver), arg0, arg1)
end

function set_max_value(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setMaxValue", void, (jdouble,), arg0)
end

function set_min_value(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setMinValue", void, (jdouble,), arg0)
end

function set_name(obj::ParameterDriver, arg0::JString)
    return jcall(obj, "setName", void, (JString,), arg0)
end

function set_normalized_value(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setNormalizedValue", void, (jdouble,), arg0)
end

function set_reference_date(obj::ParameterDriver, arg0::AbsoluteDate)
    return jcall(obj, "setReferenceDate", void, (AbsoluteDate,), arg0)
end

function set_reference_value(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setReferenceValue", void, (jdouble,), arg0)
end

function set_scale(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setScale", void, (jdouble,), arg0)
end

function set_selected(obj::ParameterDriver, arg0::jboolean)
    return jcall(obj, "setSelected", void, (jboolean,), arg0)
end

function set_value(obj::ParameterDriver, arg0::jdouble)
    return jcall(obj, "setValue", void, (jdouble,), arg0)
end

function to_string(obj::ParameterDriver)
    return jcall(obj, "toString", JString, ())
end

ext: jl


text:
vec_proj(vec_in::SVector{3, TD}, n̂::SVector{3,TD}) where {TD} = n̂ * dot(n̂, vec_in)
function vec_sub_vec_proj(v::SVector{3,T}, n̂::SVector{3,T}) where T
	t = -dot(v, n̂)
	return SVector{3,T}(muladd(t, n̂[1], v[1]),
						muladd(t, n̂[2], v[2]),
						muladd(t, n̂[3], v[3]))
end

function a_dot_one_pad_b(a::Union{SMatrix{1,4,T1,4},SVector{4,T1}}, b::SVector{3,T2}) where {T1,T2}
	d =    muladd(a[1], b[1], a[4])
	d =    muladd(a[2], b[2], d)
	return muladd(a[3], b[3], d)
end


text:
"""
Un object de classe Solution encapsule les éléments d'une solution

Args:
    inst (Instance): instance du problème
    planes (Vector{Plane}): Avions
    x (Vector{Int}): Dates d'atterrissage des avions
    costs (Vector{Float64}): Coût de chaque avion
    cost (Float64): Coût total de la solution
    timing_algo_solver
"""
mutable struct Solution
    inst::Instance
    planes::Vector{Plane}
    x::Vector{Int}
    costs::Vector{Float64}
    cost::Float64
    timing_algo_solver     # :earliest, :lp, ...
    # lpTimingSolver       # solver lp si c'est :lp qui est choisi
    # earliestTimingSolver # solver dp si c'est :dp qui est choisi

    # l'attribut contiendra l'objet solver dont le type dépend du symbole
    # contenu dans l'attribut timing_algo_solver précédent
    solver # e.g. objet de type LpTimingSolver ou EarliestimingSolver

    # Le constructeur peut prendre
    # - soit une instance qui sert à initialiser la solution
    # - soit un autre objet solution pour en faire une copie
    #
    function Solution(inst::Instance; update::Bool=true, algo=:ARGS)
        this = new()
        this.inst = inst
        this.planes = copy(inst.planes)
        # @show inst.nb_planes
        this.x = zeros(Int, inst.nb_planes)
        this.costs = zeros(inst.nb_planes)
        this.cost = typemax(Float64) # valeur Inf pour un Float64)

        # @show algo
        # @show Args.t_algos()
        if algo == :ARGS
            this.timing_algo_solver = Args.get("timing_algo_solver")
        else
            if String(algo) in Args.t_algos()
                this.timing_algo_solver = algo
            else
                error("timing_algo_solver \"$(algo)\" inconnu !")
            end
        end
        # this.solver = :nothing # sera initialisé lors du premier solve!()
        init_solver(this, this.timing_algo_solver)

        ln3("Solution: Construction d'une solution avec $(this.timing_algo_solver)")
        if update
            solve!(this)
        end
        return this
    end

    # ATTENTION la référence au timingSolver est copiées superficiellement !
    function Solution(sol::Solution)
        this = new()
        # On recopie une solution existante
        this.inst = sol.inst
        this.planes = copy(sol.planes)
        this.x = copy(sol.x)
        this.costs = copy(sol.costs)
        this.cost = sol.cost
        this.timing_algo_solver = sol.timing_algo_solver
        this.solver = sol.solver
        # solve!(this)  # ceci est en principe inutile
        return this
    end

    function Solution(instance::Instance, planes::Vector{Plane}, x)
        solution = new(instance, planes, x, zeros(Int, instance.nb_planes), 0,
                       "none", "none")
        update_costs!(solution)
        return solution
    end
end

"""Retourne la solution """
function Solution(inst::Instance, filename::String; update::Bool=true, algo=:ARGS)
    this = Solution(inst, update=false, algo=:ARGS)
    readsol(this, filename)
end

"""Recopie dans la solution *sol* le contenu de la solution *other*."""
function copy!(sol::Solution, other::Solution)
    sol.inst = other.inst
    copyto!(sol.planes, other.planes)
    copyto!(sol.x,      other.x)
    copyto!(sol.costs,  other.costs)
    sol.cost = other.cost  # ceci maintient la cohérence du coût avec la solution
    sol.timing_algo_solver = other.timing_algo_solver
    sol.solver = other.solver # ATTENTION : COPIE SUPERFICIELLE
    return sol   # pratique pour enchainer les appels de méthode (en POO !)
end

"""
Crée et affecte l'attribut solver pour la résolution des dates d'atterrissage.

Le type du solver créée est mémorisé dans l'attribut symbole
``sol.timing_algo_solver``
(par exemple :ealiest, :lp, ...)
"""
function init_solver(sol::Solution, algo::Symbol)
    if algo == :lp
        sol.solver = LpTimingSolver(sol.inst)
    elseif algo == :earliest
        # sol.solver = nothing # pas de TimingSolver (inutile) pour :earliest
        sol.solver = EarliestTimingSolver(sol.inst)
    else
        error("init_solver:  algo invalide : $(algo)")
    end
    sol.timing_algo_solver = algo
end

"""
    current_algo_symbol(sol::Solution)

Return le symbol de l'algo à partir du TimingSolver mémorisé dans la
variable sol.solver
Si l'objet n'est pas un solver, le symbol :earliest est retourné
"""
function current_algo_symbol(sol::Solution)
    if sol.solver !== nothing
        return symbol(sol.solver)
    else
        return :earliest
    end
end

"""
    hash(sol::Solution)

Calcul du hachage d'une solution.
Deux solutions identiques doivent donner le même hachage.
Deux solutions sont considérées comme identique ssi elles ont les avions
dans le même ordre, et avec les mêmes dates d'atterrissage.
HYP: les deux solutions sont à jour (par solve!)
"""
function hash(sol::Solution)
    # on construit une chaine de caractères représentative de la solution
    # (ce qui revient plus ou moins à sérialiser l'objet sol)
    io = IOBuffer()
    for p in sol.planes
        print(io, p.name, ":", round(Int, sol.x[p.id]), "-")
    end
    return Base.hash(String(take!(io)))
end

"""Retourne le tableaux des noms d'avions dans l'ordre de la solution."""
function get_names(sol::Solution)
    # RUBY: @planes.map{| plane |  plane.name}
    map(p::Plane->p.name, sol.planes)
end

"""
Affecte la solution en fonction de l'ordre des noms d'avion définit
en paramètre

Example:
   ``set_from_names!(mysol, [3, 4, 5, 6, 8, 9, 7, 1, 10, 2])``

Attention:
    Couteux O(n^2), mais pratique pour affichage et tests.
"""
function set_from_names!(sol::Solution, names)
    if length(names) != length(sol.planes)
        error("\nERREUR : taille de names $(length(names)) incorrecte " *
              "attendue : $(length(sol.planes))")
    end
    # sol.planes .= get_plane_from_name.(sol.inst, names)
    for i in 1:length(names)
        sol.planes[i] = get_plane_from_name(sol.inst, names[i])
    end
    return sol
end

"""
Trie les avions de l'objet sol selon leur date d'atterrissage croissante.

La solution n'est pas recalculée, mais les vecteurs *x* et *planes* sont
triés de manière cohérente.
"""
function sort!(sol::Solution)
    # Pour cela, on récupère l'ordre correspondant au tri des dates d'atterrissage
    # (le résultat est la permutation correspondant au tri à effectuer)
    permu = sortperm(sol.x)

    # On applique cette permutation aux deux tableaux de notre solution pour les
    # trier dans le même ordre
    sol.x = sol.x[permu]
    sol.planes = sol.planes[permu]
    sol.costs = sol.costs[permu]
end

function Base.show(io::IO, sol::Solution)
    Base.write(io, to_s(sol))
end

"""
Convertit une Solution en chaine

Example:
    ``"cost=123.456  :[p2,p5,p1,...,p10]"``

"""
function to_s(sol::Solution)
    # prc = Args.args[:cost_precision] + 2 # cmt le 30/10/2018 car 7 par défaut
    prc = Args.args[:cost_precision]
    cost = round(sol.cost, digits=prc)

    txt = "cost=" * rpad(cost, 8) * " :" *
    # "[" * sol.planes.map {|plane| plane.name }.join(",") << "]"
    "[" * join(map( p->p.name, sol.planes), ",") * "]"
    return txt
end

function to_s_long(sol::Solution)
    io = IOBuffer()
    prc = Args.args[:cost_precision]
    cost = round(sol.cost, digits=prc)
    viol_cost = get_viol_penality(sol)

    println(io, "# ALP solution version 1.0")
    println(io, "name $(sol.inst.name)  # nb_planes=$(length(sol.planes))")
    if viol_cost != 0.0
        println(io, "# Solution **invalide** ")
        println(io, "viols_cost=$(viol_cost)")
    else
        println(io, "# Solution valide")
    end
    println(io, "timestamp $(Dates.now())")
    println(io, "cost $(cost)")
    println(io, "order [", join([p.name for p in sol.planes], ","), "]")
    println(io)
    println(io, "#       name     t   dt  cost        # comments")
    for i in 1:length(sol.planes)
        p = sol.planes[i]
        print(io, @sprintf("landing %4s ", p.name) )
        print(io, @sprintf(" %4s ", sol.x[i]) )
        print(io, @sprintf("%4s  ", sol.x[i] - p.target) )

        # Ceci est horrible, mais la chaine de formatage de @sprintf ne peut
        # pas contenir de variable !
        #   print(io, @sprintf("% 7.$(prc)f ", sol.costs[i]) ) # CECI PLANTE
        if prc <= 5
            print(io, @sprintf("%9.5f ", sol.costs[i]) )
        else
            print(io, @sprintf("%11.7f ", sol.costs[i]) )
        end

        print(io, "  # ")
        print(io, @sprintf("E=%-4d ", p.lb) )
        print(io, @sprintf("T=%-4d ", p.target) )
        print(io, @sprintf("L=%-4d ", p.ub))
        print(io, "  ")
        print(io, @sprintf(" ep=%3.2f  ", p.ep) )
        print(io, @sprintf(" tp=%3.2f  ", p.tp) )

        # On affiche les temps d'écart et temps attendu entre les trois derniers
        # prédessesseurs pour vérifier le respect des temps d'écart.
        print(io, " sep ")
        isbad = false
        # sep_string = ""
        for j in i-1:-1:max(i-3,1)
            p_prev = sol.planes[j]
            sep_real = sol.x[i]-sol.x[j]
            sep_th = get_sep(sol.inst, p_prev, p)
            # sep_string *= " $sep_real($sep_th)"
            if sep_real < sep_th
                isbad = true
            end
            print(io, " $sep_real($sep_th)")
        end
        if isbad
            # print(io, to_sc(" VIOL "), sep_string)
            println(io, to_sc(" VIOL "))
        else
            # print(io, " -ok- ", sep_string)
            println(io, " -ok- ")
        end
    end
    println(io)
    String(take!(io))
end

"""
Vérifie la faisabilité de la solution par rapport aux dates d'atterrissage
réelle de la solution.

Retourne un coût de pénalité (nul si solution valide) si des contraintes sont
violées.
Ce coût prend en compte :

- l'amplitude des viols de la contrainte de séparation,
- l'amplitude des viols des bornes lb et ub de chaque avion.

Cette méthode suppose que l'information sol.x de la date d'atterrissage
réelle de chaque avion est à jour.

Les viols sont détectés même si l'inégalité triangulaire n'est pas respectée
dans l'instance (instance incorrecte par hypothèse)

Si level est suffisant : affiche message si solution non faisable
"""
function get_viol_penality(sol::Solution)
    total_penality = 0.0
    unit_penality = 1000.0
    SEP_MAX = sep_max(sol.inst)

    ### ASSERT : p1=sol.planes[1] est bien positionné sinon il

    # On commence par prendre en compte la pénalité de placement du premier avion
    # On calcule l'écart (viol sur les bornes) de la date d'atterissage x[1]
    # cet écart est positif sauf en cas de viol.
    p1 = sol.planes[1]
    dt = min(sol.x[1]-p1.lb, p1.ub-sol.x[1])
    if dt < 0
        total_penality += -dt*unit_penality
    end

    for i2 in 2:length(sol.planes)
        p2 = sol.planes[i2]

        # Étape 1 : calcul du viol des temps de séparation pour p2

        # parcours les prédecesseurs en partant de de p2
        for i1 in i2-1:-1:1
            if sol.x[i2] > sol.x[i1] + SEP_MAX
                # tous les temps d'écarts des prédécesseurs sont respectés
                break
            end
            p1 = sol.planes[i1]
            # s_val = sol.inst.sep_mat[p1.kind, p2.kind]
            s_val = get_sep(sol.inst, p1, p2)
            if sol.x[i2] < sol.x[i1] + s_val
                total_penality += unit_penality * (sol.x[i1] + s_val  - sol.x[i2])
                # il faudrait pousser j
                feasible = false
                if lg4()>=4
                    println("ERREUR écart entre type des avions $(p1.name)->$(p2.name)")
                    println("   x1=$(sol.x[i1]) pour $(p1)")
                    println("   x2=$(sol.x[i2]) pour $(p2)")
                    println("   écart réel=$(sol.x[i2]-sol.x[i1]) au lieu de $(s_val)")
                end
            end
        end

        # Étape 2 : calcul du viol de la Upper Bound
        if sol.x[i2] > p2.ub
            total_penality += unit_penality * (sol.x[i2] - p2.ub)
            lg4() && println("UB infaisable pour l'avion $(p2.name)")
        end

        # Étape 3 : calcul du viol de la Lower Bound
        if sol.x[i2] < p2.lb
            total_penality += unit_penality *  (p2.lb - sol.x[i2])
            lg4() && println("LB infaisable pour l'avion $(p2.name)")
        end
    end
    if lg4()>=4 && total_penality==0.0
        println("solution faisable.")
    end
    feasible = (total_penality==0.0)
    if !feasible && lg4()
        println("get_viol_penality: Infaisable : total_penality=", total_penality)
    end
    return total_penality
end

"""
Vérifie la faisabilité dates d'atterrissage x de la solution par rapport
aux contraintes de l'instance.

- Si pas de viol, retourne une chaine vide
- Si viol, retourne un texte avec une ligne par contrainte violée.

Les contraintes testées sont les suivante :

- la date d'atterissage doit entre dans l'intervale permi pour l'avion
- tous les temps de séparation doivent être respertés

Cette méthode suppose que l'information *sol.x* de la date d'atterrissage
réelle de chaque avion est à jour.

Les viols sont détectés même si l'inégalité triangulaire n'est pas respectée
dans l'instance (instance incorrecte par hypothèse).
"""
function get_viol_description(sol::Solution)
    io = IOBuffer()

    nbviols = 0
    for i in 1:length(sol.planes)
        p = sol.planes[i]
        x = sol.x[i] ; lb = p.lb ; ub = p.ub
        if sol.x[i] - p.lb < 0
            println(io, p.name, " atterrit trop tôt  x=$x < lb=$lb")
            nbviols += 1
        end
        if p.ub - sol.x[i] < 0
            println(io, p.name, " atterrit trop tard ub=$ub < x=$x")
            nbviols += 1
        end
    end

    SEP_MAX = sep_max(sol.inst)
    for i2 in 2:length(sol.planes)
        p2 = sol.planes[i2]
        x2 = sol.x[i2]
        # parcours les prédecesseurs en partant de de p2
        for i1 in i2-1:-1:1
            p1 = sol.planes[i1]
            x1 = sol.x[i1]
            # On pourrait ne pas explorer tout si l'on se contentait de détecter
            # une invalidité. Mais on veut afficher **tous** les viols
            #   if x2 > x1 + SEP_MAX
            #       # tous les temps d'écarts des prédécesseurs sont respectés
            #       # (s'il n'y a aucun e viol avant x1)
            #       break
            #   end

            sep = get_sep(sol.inst, p1, p2)
            if x2 < x1 + sep
                print(io,  "$(p1.name)->$(p2.name) : écart insuffisant ")
                print(io,  "x[$(p1.name)]=$x1->x[$(p2.name)]=$x2 => ")
                println(io, "sep($(p1.name),$(p2.name))=$(sol.x[i2]-sol.x[i1]) ",
                "au lieu de ", get_sep(sol.inst, p1, p2) )
                nbviols += 1
            end
        end
    end
    violtxt = String(take!(io))
    return (nbviols, violtxt)
end

"""
Calcul la faisabilité intrinsèque de la solution en fonction de l'ordre
imposé pour les avions.

Ne modifie pas la solution.

N'utilise pas l'attribut `x[i]` de la solution associée au placement
de chaque avion.

Note:
    Si on peut placer les avions au plus tôt en respectant les temps de
    séparation sans violer la borne sup du placement des avions suivants,
    alors la solution est faisable.
"""
function is_feasable(sol::Solution; param_bidon=nothing)
    # ASSUME_TRINEQ = Args.get(:assume_trineq)
    ASSUME_TRINEQ = false # L'inégalité triangulaire n'est pas présupposée

    n = length(sol.planes)
    # x : vecteur des placements au plus tôt des avions
    x = Vector{Int}(undef,n)

    i1 = 1
    p1 = sol.planes[i1]
    x[i1] = p1.lb
    # @assert (x[i1] <= p1.ub), "Instance infaisable pour avion $(p1) !"
    @assert x[i1] <= p1.ub

    for i2 in 2:length(sol.planes)
        p2 = sol.planes[i2]
        x[i2] = max(p2.lb, x[i1] + get_sep(sol.inst, p1, p2))
        if !ASSUME_TRINEQ
            # inégalité triangulaire non vérifiée : on doit vérifier
            # les prédécesseurs antérieurs à p1
            for i_prev in i1-1 : -1 : 1
                p_prev - sol.planes[i_prev]
                this_sep = x[i2] - x[i_prev]
                if this_sep >= MAX_SEP
                    # c'est bon : p2 est forcément suffisamment loin devant
                    # tous ces prédécesseurs
                    break
                end
                # on retarde p2 pour satisfaire cette séparation
                x[i2] = max(x[i2], x[i_prev] + get_sep(sol.inst, p_prev, p2))
            end
        end
        # maintenant que p2 est placé au plus tôt en respectant les contraintes
        # de séparation, on vérifie qu'il n'est pas hors borne supérieure
        if x[i2] > p2.ub
            # if lg1() println("viol x[$(i2)]=$(x[i2]) mais $(p2).ub=$(p2.ub)") end
            return false
        end
    end
    # Tous les avions ont pu être placés au plus tôt.
    # La solution est donc timing-faisable
    return true
end

"""

TODO:
    prévoir paramètre from_index=1 pour permettre un recalcul partiel.

do_update_cost (false par défaut) : si vrai, met à jours les coûts.

Notes:
    - l'option do_update_cost n'est utile que si ce n'est pas déjà fait
      par le TimingSolver.
    - Elle ne semble plus utilisée (sauf par solve_to_earliest!(sol)).
    - Il est conseillé d'appeler explicitement la méthode update_costs!(sol)
      plutot que d'utiliser cette option.
    - Si aucun TimingSolver n'a été spécifié lors de la création de la solution
      alors on utilise la résolution au plus tot.
"""
function solve!(sol::Solution; do_update_cost::Bool=false)
    if sol.solver === nothing
        init_solver(sol, sol.timing_algo_solver)
    end
    solve!(sol.solver, sol)
    if do_update_cost
        update_costs!(sol)
    end
    return sol
end

"""
Mélange la solution et mets à jour (par défaut) les timings

Args:
    do_update (Bool): si vrai, résoud le pb de timing

Note:
    L'option do_update ne semble plus utilisée.
"""
function shuffle!(sol::Solution; do_update::Bool=true)
    Random.shuffle!(sol.planes)
    if do_update
        solve!(sol) # NEW 03/04/2019
    end
    return sol
end

"""Remue légèrement la solution et mets à jour (par défaut) les timings"""
function disturb!(sol::Solution;
                  idx_first=1, idx_last=-1, shift_max=1,
                  nb_shift=10, do_update::Bool=true,
                  )
    # VÉRIFICATION DU DOMAINE DES PARAMÈTRES
    if idx_last == -1   idx_last = length(sol.planes)   end
    @assert 1 <  idx_last  <= length(sol.planes)
    @assert 1 <= idx_first <  idx_last

    # shift_max == -1 => on prend le maximum possible
    if shift_max == -1  shift_max = idx_last-idx_first  end
    @assert shift_max >= 1 "shift_max doit être stricement positif (ou -1 pour max)"

    # APPLICATION DES MOUVEMENTS (de type shift)
    for i in 1:nb_shift
        dist = rand(1:shift_max)
        idx1 = rand(idx_first:idx_last)
        idx2_min = max(idx1-dist, idx_first)
        idx2_max = min(idx1+dist, idx_last)
        idx2 = rand(idx2_min:idx2_max)
        shift!(sol, idx1, idx2, do_update=false)
    end

    # RECALCUL ÉVENTUELLE DES COÛTS DE LA SOLUTION
    if do_update
        # do_update_cost est utile ici par la solution obtenue peut-être infaisable
        # et il faut pouvoir évaluer le coût de cette solution
        solve!(sol, do_update_cost=true) # mise à jour après le dernier shift
    end
end

"""
Met à jour les coûts de chaque avion (sol.costs) et le coût global (sol.cost)
à partir des dates d'atterrissage supposées connues (sol.x).

Caution:
    Précondition: le Vector s.costs est dans le même ordre que s.planes
"""
function update_costs!(sol::Solution; add_viol_penality = true )

    # if !isnothing(x)
    #     copy!(sol.x, x)
    # end

    # Mise à jour des coûts individuels de chaque avion
    for i in 1:length(sol.planes)
        p = sol.planes[i]
        sol.costs[i] = get_cost(p, sol.x[i])
    end

    # Mise à jour du coût global de la solution
    #
    sol.cost = sum(sol.costs)
    # @show Args.args[:cost_precision]
    # @show sol.cost
    sol.cost = round(sol.cost, digits=Args.args[:cost_precision]) # 30/10/2018
    # @show sol.cost

    if lg4()>=4 println("get_viol_penality=$(get_viol_penality(sol))") end
    if add_viol_penality
        sol.cost += get_viol_penality(sol)
    end
    return sol
end

"""
Met à jour la date d'atterrissage de chaque avion au plus tôt, de façon à
respecter les contraintes de précédence au mieux.

Important:
  Les contraintes de séparation sont **toujours respectées**.
  Au besoin on viole la borne  d'atterrissage des avions successeurs (p.ub)
"""
function solve_to_earliest!(sol::Solution; do_update_cost=true)
    # println("="^70)
    # println("=== BEGIN solve_to_earliest! pour $(to_s(s))")
    # println("")
    SEP_MAX = sep_max(sol.inst)
    # Méthode : on parcourt tous les successeurs (plane2) et on s'assure qu'ils
    # sont suffisamment éloignés de **tous** leurs prédécesseurs (les plane1).
    for i2 in 1:length(sol.planes)
        plane2 = sol.planes[i2]
        # @show plane2
        sol.x[i2] = plane2.lb # date au plus tot
        # parcours les prédecesseurs de plane2 du plus proche au plus éloigné
        for i1 in i2-1:-1:1
            # @show sol.planes[i1]
            if (sol.x[i2] - sol.x[i1]) >= SEP_MAX
                # i2 est suffisamment en avance quels que soit le type de
                # l'avions i1 qui le précède
                break
            end
            plane1 = sol.planes[i1]
            if sol.x[i2] < sol.x[i1] + sol.inst.sep_mat[plane1.kind, plane2.kind]
                # il faut pousser plane2
                sol.x[i2] = sol.x[i1] + sol.inst.sep_mat[plane1.kind, plane2.kind]
            end
        end
        if sol.x[i2] > plane2.ub
            if lg4()>=4 println("pb infaisable pour l'avion $(plane2.name)" ) end
            # TODO: On pourrait ajouter un coût de pénalité pour les
            # contraintes de bornes non respectées
        end
    end
    if do_update_cost
        update_costs!(sol)
    end
    return sol
end

"""
Permute deux avions d'indices `idx1` et  `idx2` dans la solution.

Si un indice n'est pas défini, il est choisi aléatoirement.

Returns:
    indices des avions permutés
"""
function swap!(sol::Solution, idx1=-1, idx2=-1; do_update=true)
    # RUBY: idx1 = rand(@planes.size)  if not idx1
    if  idx1 == -1
        idx1 = rand(1:length(sol.planes))
    end
    if  idx2 == -1
        idx2 = rand(1:length(sol.planes))
    end
    # Il se peut que idx1 == idx2 mais tant pis
    sol.planes[idx1], sol.planes[idx2] = sol.planes[idx2], sol.planes[idx1]
    if idx1 != idx2 && do_update
        solve!(sol)
    end
    return (idx1, idx2)
end

"""
Déplace avions dans la solution

Args:
    idx1: indice de l'avion d'origine
    idx2: indice de l'avion cible

Si un indice n'est pas défini, il est choisi aléatoirement.

Returns:
    indices des avions permutés
"""
function shift!(sol::Solution, idx1=-1, idx2=-1; do_update=true)
    if  idx1 == -1
        idx1 = rand(1:length(sol.planes))
    end
    if  idx2 == -1
        idx2 = rand(1:length(sol.planes))
    end
    # Il se peut que idx1 == idx2 mais tant pis
    if idx1==idx2
        do_update && solve!(sol)
        return (idx1, idx1)
    elseif abs(idx1-idx2)==1
        return swap!(sol, idx1, idx2, do_update=do_update)
    end
    # println("AVANT SHIFT $idx1->$idx2 : ", to_s(sol))
    # 03/05/2019 : shift! est définie dans le fichier utils/array.jl
    shift!(sol.planes, idx1, idx2)

    if do_update
        solve!(sol)
    end
    # println("APRES SHIFT $idx1->$idx2 : ", to_s(sol))
    return (idx1, idx2)
end

"""
Permute des avions dans la solution.

Args:
    indices1: tableau des indices d'avions à permuter
    indices2: tableau des nouveaux indices des avions permutés
    do_update : recalcule le coût (true par défaut)

HYPOTHESE : indices2 est une permutation valide de indices1

Les avions de la solution ont pour indice 0:nb_planes-1
Par exemple, si indices1 vaut [0, 2, 7, 12]
et si indices2 vaut [12, 7, 0, 2],
L'opération consiste à permuter les avions positionnés en :
   ``0, 2, 7, 12``
de façon à ce que leur nouvelle position soit en :
  ``12, 7, 0, 2``

Exemple:
    On veut transformer la liste d'avion :
    ``
    @plane = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
                0   1*  2   3   4   5   6*  7*  8    9
    ``
    en :
    ``
    @plane = [10, 70, 30, 40, 50, 60, 80, 20, 90, 100]
                0   1*  2   3   4   5   6*  7*  8    9
    ``
    Pour cela l'appel doit être de la forme :
    ``mysol.permu! [1, 6, 7], [7, 1, 6]``

TODO: effectuer une mise à jour différencielle du coût
"""
function permu!(sol::Solution, indices1, indices2; do_update=true)
    sol.planes[indices1] = sol.planes[indices2]
    # sol.planes[[indices1...]] = sol.planes[[indices2...]] # avec des tuples iiii
    if do_update
        solve!(sol)
    end
    return indices2
end

"""
Effectue un tri de la solution courante selon le critère passer en paramètre.

Si aucun critère presort n'est imposé en argumant, celui spécifié par
l'option --presort est utilisé

Le timing de la solution est mis à jour.
"""
function initial_sort!(sol::Solution; presort=:ARGS, solve=true)
    if presort == :ARGS
        presort = Args.get("presort")
    end

    if presort == :target
        Base.sort!(sol.planes, by=p->p.target)
    elseif presort == :rtarget
        Base.sort!(sol.planes, by=p->p.target, rev=true)
    elseif presort == :lb
        Base.sort!(sol.planes, by=p->p.lb)
    elseif presort == :rlb
        Base.sort!(sol.planes, by=p->p.lb, rev=true)
    elseif presort == :ub
        Base.sort!(sol.planes, by=p->p.ub)
    elseif presort == :rub
        Base.sort!(sol.planes, by=p->p.ub, rev=true)
    elseif presort == :mean
        Base.sort!(sol.planes, by=p->(p.lb + p.ub))
    elseif presort == :bary
        Base.sort!(sol.planes, by=p->(p.lb + p.target + p.ub))
    elseif presort == :shuffle
        # Base.shuffle!(sol.planes)
        Random.shuffle!(sol.planes)
    elseif presort==:none
        # rien à faire
    else
        error("Valeur incorrecte pour presort : ", presort)
    end

    # Autre solution pour trier de type d'avions (A GARDER COMME EXEMPLES) :
    # On passe la fonction qui compare deux éléments :
    # Base.sort!(sol.planes, lt=(p1,p2)->p1.lb < p2.lb) # Marche aussi (comparator)
    # Base.sort!(sol.planes, rev=true,lt=(p1,p2)->p1.lb < p2.lb) # ok aussi THE WORST CASE
    solve!(sol, do_update_cost=true)
    return nothing
end

"""
Enregistre la solution dans un fichier.

Par défaut, le nom est construit automatiquement à partir des caractéristiques
de l'instance et du coût de la solution (arrondi à deux décimales)
"""
function write(sol::Solution, filename="")
    dir = Args.args[:outdir] # e.g "_tmp"
    if !isdir(dir)
        dir = "."
    end
    cost = round(sol.cost, digits=2) # cost intervient dans le nom du fichier
    if filename == ""
        # filename = "seqata_p$(sol.inst.nb_planes)=$(cost).sol"
        filename = guess_solname(sol)
    end
    # solpath = dir*"/"*filename
    solpath = "$dir/$filename"
    open(solpath, "w") do fh
        # Base.write(fh, to_s(sol))
        # Base.write(fh, "\n")
        Base.write(fh, to_s_long(sol))
    end
end

"""
Infère le nom de la solution.

Le nom est de la forme:
``alp_01_p10=700.0.sol``
"""
function  guess_solname(sol::Solution)
    prc = Args.args[:cost_precision]
    cost = round(sol.cost, digits=prc)
    return "$(sol.inst.name)=$(cost).sol"
end

"""Fonction d'impression de haut niveau."""
function print_sol(sol::Solution, msg::AbstractString= "")
    if msg!=""  println(msg) end
    println("\nmeilleure solution trouvée :")
    println("="^70)
    println(to_s(sol))
    println("="^70)
    println(to_s_long(sol))
    println("="^70)
end

# END TYPE Solution


text:
using ITensorLocalMPO
using Test

@testset "ITensorLocalMPO.jl" begin
    # Write your tests here.
end

max_stars_repo_name: mtfishman/ITensorLocalMPO

max_issues_repo_name: mtfishman/ITensorLocalMPO

max_forks_repo_name: mtfishman/ITensorLocalMPO


text:
module Interop

import ..CUDA: CuArray, CuSubArray, AbstractCuArray

immutable Array{T,N}
    ptr::Ptr{T}
    dims::N
    strides::N
    continuous::Bool
end

Array{T,N}(x::CuArray{T,N}) = Array(pointer(x),cint(size(x)),cint(strides(x)),true)
Array{T,N}(x::CuSubArray{T,N}) = Array(pointer(x),cint(size(x)),cint(strides(x)),false)

immutable Cint1
    i1::Cint
end
immutable Cint2
    i1::Cint
    i2::Cint
end
immutable Cint3
    i1::Cint
    i2::Cint
    i3::Cint
end
immutable Cint4
    i1::Cint
    i2::Cint
    i3::Cint
    i4::Cint
end

cint(t::NTuple{1,Int}) = Cint1(t[1])
cint(t::NTuple{2,Int}) = Cint2(t[1],t[2])
cint(t::NTuple{3,Int}) = Cint3(t[1],t[2],t[3])
cint(t::NTuple{4,Int}) = Cint4(t[1],t[2],t[3],t[4])

end


text:
using Flux, CUDA, Test
using Flux: pullback

@testset "CUDNN BatchNorm" begin
    @testset "4D Input" begin
        x = rand(Float32, 2, 2, 3, 4)
        m = BatchNorm(3)
        gx = gpu(x)
        gm = gpu(m)

        y, back = pullback((m, x) -> m(x), m, x)
        gy, gback = pullback((m, x) -> m(x), gm, gx)

        @test cpu(gy) ≈ y

        Δ = randn(Float32, size(y))
        dm, dx = back(Δ)
        gdm, gdx = gback(gpu(Δ))

        @test dm[].γ ≈ cpu(gdm[].γ)
        @test dm[].β ≈ cpu(gdm[].β)
        @test dx ≈ cpu(gdx)
    end

    @testset "2D Input" begin
        x = rand(Float32, 3, 4)
        m = BatchNorm(3)
        gx = gpu(x)
        gm = gpu(m)

        y, back = pullback((m, x) -> m(x), m, x)
        gy, gback = pullback((m, x) -> m(x), gm, gx)

        @test cpu(gy) ≈ y

        Δ = randn(Float32, size(y))
        dm, dx = back(Δ)
        gdm, gdx = gback(gpu(Δ))

        @test dm[].γ ≈ cpu(gdm[].γ)
        @test dm[].β ≈ cpu(gdm[].β)
        @test dx ≈ cpu(gdx)
    end
end


text:
mutable struct HTMLTemplater <: Rewriter
    assets::Assets
    template::XNode
    templatepath::AbstractPath
    inlineincludes::Bool
    insertpos::Position
end

Base.show(io::IO, templater::HTMLTemplater) = print(io, "HTMLTemplater($(templater.assets)), $(templater.templatepath)")

function HTMLTemplater(
        templatepath::AbstractPath,
        includes::Vector{<:AbstractPath}=Path[];
        assetdir = Path("."),
        inlineincludes=false,
        insertpos=FirstChild(SelectTag(:body)),
    )
    template = Pollen.parse(templatepath)
    assets = Assets(Dict(joinpath(p"template/", p) => absolute(joinpath(assetdir, p)) for p in includes))
    if inlineincludes
        template = inlineintemplate(
            template,
            [absolute(joinpath(assetdir, p)) for p in includes])
    else
        template = includeintemplate(template,  [joinpath(p"template/", p) for p in includes])
    end

    return HTMLTemplater(assets, template, templatepath, inlineincludes, insertpos)
end

function rewritedoc(templater::HTMLTemplater, p, doc)
    # Include the document in the template
    doc = insertfirst(templater.template, doc, templater.insertpos)
    return doc
end

function includeintemplate(template::XNode, includes)
    for p in includes
        ext = extension(p)
        if ext == "css"
            x = XNode(:link, Dict(:rel => "stylesheet", :href => "/" * string(p)))
        elseif ext == "js"
            x = XNode(:script, Dict(:src => "/" * string(p)))
        else
            continue
        end
        template = insertfirst(template, x, FirstChild(SelectTag(:head)))
    end

    return template
end

function inlineintemplate(template::XNode, includes)
    for p in includes
        ext = extension(p)
        if ext == "css"
            x = XNode(:style, [XLeaf(read(p, String))])
        elseif ext == "js"
            x = XNode(:script, [XLeaf(read(p, String))])
        end
        template = insertfirst(template, x, FirstChild(SelectTag(:head)))
    end
    return template
end

function getfilehandlers(templater::HTMLTemplater, project, dir, builder)
    handlers =  [
        # When template changes, reload it and rebuild every file
        (
            absolute(templater.templatepath),
            () -> onupdatetemplate(templater, project, builder)
        ),
    ]

    if templater.inlineincludes
        handlers = vcat(handlers, [
            (p, () -> onupdatetemplate(templater, project, builder)) for p in values(templater.assets.assets)
        ])
    else
        handlers = vcat(handlers, getfilehandlers(templater.assets, project, dir, builder))
    end
    return handlers
end

function onupdatetemplate(templater, project, builder)
    template = Pollen.parse(templater.templatepath)
    if templater.inlineincludes
        templater.template = inlineintemplate(template, collect(values(templater.assets.assets)))
    else
        templater.template = includeintemplate(template, collect(keys(templater.assets.assets)))
    end
    rebuild(builder, project)
end

function postbuild(templater::HTMLTemplater, project, builder)
    postbuild(templater.assets, project, builder)
end

reset!(templater::HTMLTemplater) = reset!(templater.assets)


text:
using StructTypes, URIs, HTTP.Messages
using CodeInfoTools: code_inferred

import StructTypes: StructType, NoStructType
import Base: |, ==

export Headers, Query, Body, Params, MissingHeaders, MissingCookies,
    GET, POST, PUT, DELETE, OPTIONS, CONNECT, TRACE, PATCH, ALL

include("dasherize.jl")

abstract type HttpParameter end

abstract type HttpMethod end
struct Get <: HttpMethod end
struct Post <: HttpMethod end
struct Put <: HttpMethod end
struct Delete <: HttpMethod end
struct Connect <: HttpMethod end
struct Options <: HttpMethod end
struct Trace <: HttpMethod end
struct Patch <: HttpMethod end
struct All <: HttpMethod end

Base.String(::HttpMethod) = "get"
Base.String(::Post) = "post"
Base.String(::Put) = "put"
Base.String(::Delete) = "delete"
Base.String(::Connect) = "connect"
Base.String(::Options) = "option"
Base.String(::Trace) = "trace"
Base.String(::Patch) = "patch"
Base.String(::All) = "*"

const GET = Get()
const POST = Post()
const PUT = Put()
const DELETE = Delete()
const OPTIONS = Options()
const CONNECT = Connect()
const TRACE = Trace()
const PATCH = Patch()
const ALL = All()

struct Cookies{T} <: HttpParameter
    t::Type{T}
end

struct Headers{T} <: HttpParameter
    t::Type{T}
    val::Union{T,Nothing}
end

function kw_constructor_data_type(T; kwargs...)
    k = []
    v = []
    for (x, y) in kwargs
        @assert y isa DataType "Argument is not a DataType"
        push!(k, x)
        push!(v, y)
    end
    t = NamedTuple{tuple(k...),Tuple{v...}}
    T(t)
end

function kw_constructor(T; kwargs...)
    k = []
    v = []
    has_datatype = false

    for (x, y) in kwargs

        if y isa DataType
            has_datatype = true
        end

        push!(k, x)
        push!(v, y)
    end

    if has_datatype
        @assert all(map(x -> x isa DataType, v))
        t = NamedTuple{tuple(k...),Tuple{v...}}
        return T(t, nothing)
    else
        nt = values(kwargs)
        T(typeof(nt), nt)
    end
end

Headers(;kwargs...) = kw_constructor(Headers; kwargs...)
Headers(t::DataType) =  Headers(t, nothing)

function fieldnames_to_header(T)
    fields = string.(fieldnames(T))
    l = []
    for i in fields
        push!(l, dasherize(i))
    end
    l
end

headerize(s) = dasherize(string(s))

struct Query{T} <: HttpParameter
    t::Type{T}
    val::Union{T,Nothing}
end

Query(t::DataType) = Query(t, nothing)
Query(; kwargs...) = kw_constructor_data_type(Query; kwargs...)

struct Body{T} <: HttpParameter
    t::Type{T}
    val::Union{T,Nothing}
end

Body(t::DataType) = Body(t, nothing)
Body(;kwargs...) = kw_constructor(Body; kwargs...)
Body(t) = Body(typeof(t), t)

function parameter_type(t::Type{<:HttpParameter})
    t.parameters[1]
end

struct Params{T} <: HttpParameter
    t::Type{T}
    val::Union{T,Nothing}
end

Params(; kwargs...) = kw_constructor_data_type(Params; kwargs...)
Params(t::DataType) = Params(t, nothing)

# https://www.juliabloggers.com/the-emergent-features-of-julialang-part-ii-traits/

function (c::Cookies{T})(stream) where {T}
    hs = headers(stream)
    cs = HTTP.Cookies.readcookies(hs, "")

    # TODO:
    # support fields other than strings!
    # how to handle additional cookie information such as
    # maxage, expires e.g
    cookies::Dict{String,Any} = Dict(c.name => c.value for c in cs)
    fields = fieldnames(T)
    d::Dict{Symbol,Any} = Dict(i =>
        get(cookies, string(i), missing)
                               for i in fields
    )

    try
        convert_numbers!(d, T)
        StructTypes.constructfrom(T, d)
    catch e
        rethrow(e)
    end
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
@testset "Scan" begin
    show(scanner)
    show(params)

    @test wave_x.version == v"1"
    @test wave_xy.version == v"1"
    @test wave_xyz.version == v"1"

    @test wave_x.axes == "x"
    @test wave_xy.axes == "xy"
    @test wave_xyz.axes == "xyz"

    @test isnothing(wave_x.metrics)
    @test isnothing(wave_xy.metrics)
    @test isnothing(wave_xyz.metrics)
end


text:
@testset "72.edit-distance.jl" begin
    @test min_distance("horse", "ros") == 3
    @test min_distance("intention", "execution") == 5
    @test min_distance("", "") == 0
end


text:
module PythonToJuliaHelper

export  PYTHON2JULIA_WORDS,
        python2julia

"""
List of `Pair`s to `replace` Python's words by Julia's words.
In the Regex, `\b` indicates a word boundary.
"""
const PYTHON2JULIA_WORDS = Pair{Regex, String}[
    # Regex: \b indicates a word boundary
    r"\bNone\b"    => "nothing",
    r"\bTrue\b"    => "true",
    r"\bFalse\b"   => "false",
    r"\btype\b\("  => "typeof(",
    r"\bid\b\("    => "object_id(",
    r"\blen\b\("   => "length(",
    r"\bprint\b\(" => "println(",
    r"\bdef\b"     => "function",
    r"\belif\b"    => "elseif",
    r">>>"         => "julia>"
    ]

"""
Return the indentation of the line.

```jldoctest
julia> using PythonToJuliaHelper

julia> PythonToJuliaHelper.get_indentation("        if value:")
"        "

julia> PythonToJuliaHelper.get_indentation("def fun(arg):")
""

```
"""
function get_indentation(line)
    m = match(r"^(\s+).", line)
    if m !== nothing
        return m[1]
    else
        return ""
    end
end

"It adds an `end` line at the end of an indented block."
function add_end_to_block!(lines, start_position)
    position = start_position
    look_for_indentation = false
    block_indentation = nothing
    end_indentation = get_indentation(lines[start_position])
    non_empty_lines = Int[]
    while position ≤ length(lines)
        line = lines[position]
        if length(strip(line)) == 0
            position += 1
            continue
        else
            push!(non_empty_lines, position)
        end
        if block_indentation === nothing && match(r"\:\s*(#.+)?$", line) !== nothing
            look_for_indentation = true
            position += 1
            continue
        end
        if look_for_indentation
            block_indentation = get_indentation(line)
            look_for_indentation = false
        end
        if block_indentation !== nothing
            line_indentation = get_indentation(line)
            if length(line_indentation) < length(block_indentation)
                if match(r"^\s*else\s*\:\s*(#.+)?$", line) === nothing &&
                   match(r"^\s*elif\s+.*\:\s*(#.+)?$", line) === nothing

                    insert!(lines, non_empty_lines[end-1]+1,
                            string(end_indentation, "end"))
                    block_indentation = nothing
                    break
                end
            end
        end
        position += 1
    end
    if block_indentation !== nothing
        insert!(lines, non_empty_lines[end]+1, string(end_indentation, "end"))
    end
    lines
end

"It deletes the colon at the end of the line."
function delete_colon!(lines, start_position)
    position = start_position
    while position ≤ length(lines)
        line = lines[position]
        if match(r"\:\s*(#.+)?$", line) !== nothing
            lines[position] = replace(line, r"\:(\s*(#.+)?)$" => s"\1")
            break
        end
        position += 1
    end
    lines
end

"""
Python | lambda x, y : y + x * 2
Julia  | (x,y) -> y + x * 2

```jldoctest
julia> using PythonToJuliaHelper

julia> PythonToJuliaHelper.anonymous_function("lambda x, y : y + x * 2")
"(x, y) -> y + x * 2"

julia> PythonToJuliaHelper.anonymous_function("lambda x : x * 2")
"x -> x * 2"

```
"""
function anonymous_function(line)
    m = match(r"lambda\s+(.+)\s*\:\s*", line)
    if m !== nothing
        python_args = rstrip(m.captures[1])
        julia_args = occursin(',', python_args) ? string("(", python_args, ")") : python_args
        line = replace(line, m.match => string(julia_args, " -> "))
    end
    line
end

"""
Python | range(stop)
Python | range(start, stop)
Python | range(start, stop, step)
Julia  | start:stop
Julia  | start:step:stop

Python's `range(3)` and `range(0, 3)` (i.e. 0,1,2) are going to be Julia's `0:3` (i.e. 0,1,2,3)

Python's `range` and `xrange` are translated to ranges unless `collect=true`.

```jldoctest
julia> using PythonToJuliaHelper

julia> PythonToJuliaHelper.ranges("range(stop)")
"0:stop"

julia> PythonToJuliaHelper.ranges("range(start, stop)")
"start:stop"

julia> PythonToJuliaHelper.ranges("range(start,stop,step)")
"start:step:stop"

julia> PythonToJuliaHelper.ranges("range(start,stop,step)", collect=true)
"collect(start:step:stop)"

julia> PythonToJuliaHelper.ranges("xrange(start,stop,step)", collect=true)
"start:step:stop"

```
"""
function ranges(line; collect::Bool=false)
    m = match(r"\bx?range\b\((.*)\)", line)
    if m !== nothing
        if collect && !occursin("xrange", m.match)
            prefix = "collect("
            suffix = ")"
        else
            prefix = ""
            suffix = ""
        end
        python_args = [ strip(arg) for arg in split(m.captures[1], ',') ]
        n_args = length(python_args)
        if n_args == 1
            line = replace(line, m.match => "$(prefix)0:$(python_args[1])$(suffix)")
        elseif n_args == 2
            line = replace(line, m.match => "$(prefix)$(python_args[1]):$(python_args[2])$(suffix)")
        elseif n_args == 3
            line = replace(line, m.match => "$(prefix)$(python_args[1]):$(python_args[3]):$(python_args[2])$(suffix)")
        end
    end
    line
end

"It moves the docstring before the function definition."
function function_docstrings(lines, start_position)
    position = start_position
    docstring_lines = Int[]
    started = false
    ended = false
    look_for_docstring = false
    while position ≤ length(lines)
        line = lines[position]
        if match(r"\:\s*(#.+)?$", line) !== nothing
            look_for_docstring = true
            position += 1
            continue
        end
        if look_for_docstring && !started
            if match(r"^\s*\"\"\"", line) !== nothing
                started = true
            elseif match(r"^\s*\"", line) !== nothing
                started = true
                if match(r"\"\s*(#.+)?$", line) !== nothing
                    ended = true
                else
                    break
                end
            elseif match(r"\S+", line) !== nothing
                break
            end
        end
        if started && match(r"\"\"\"\s*(#.+)?$", line) !== nothing
            ended = true
        end
        if started && !ended
            push!(docstring_lines, position)
        elseif started && ended
            push!(docstring_lines, position)
            docs = lines[docstring_lines]
            min_indentation = get_indentation(docs[1])
            for line in docs
                line_indentation = get_indentation(line)
                if length(strip(line)) > 0 && length(line_indentation) < length(min_indentation)
                    min_indentation = line_indentation
                end
            end
            docs = [ replace(line, min_indentation => "", count=1) for line in docs ]
            lines = vcat(
                lines[1:start_position-1],
                docs,
                lines[setdiff(start_position:length(lines), docstring_lines)]
                )
            break
        end
        position += 1
    end
    lines
end

"It moves all the docstrings before their function definitions."
function move_all_docstrings(lines)
    position = 1
    while position < length(lines)
        if match(r"^\s*def\s*", lines[position]) !== nothing
            lines = function_docstrings(lines, position)
        end
        position += 1
    end
    lines
end

"""
It takes Python code lines and transform it to something closer to Julia code.

1. It moves docstrings before their function definitions.
2. It adds an end line to for, if and function blocks.
3. It deletes colons at the end of for, if, def, else and elif statements.
4. It changes anonymous function syntax.
5. It changes (x)range syntax without modification of the values.
6. It replaces python words using the pairs in `PYTHON2JULIA_WORDS`.
"""
function python2julia(lines; collect::Bool=false)
    lines = move_all_docstrings(lines)
    position = 1
    while position < length(lines)
        if (match(r"^\s*if\s*", lines[position]) !== nothing) ||
                (match(r"^\s*for\s*", lines[position]) !== nothing) ||
                (match(r"^\s*def\s*", lines[position]) !== nothing)

            add_end_to_block!(lines, position)
            delete_colon!(lines, position)
        end
        if (match(r"^\s*else\s*\:\s*(#.+)?$", lines[position]) !== nothing) ||
                (match(r"^\s*elif\s+.*\:\s*(#.+)?$", lines[position]) !== nothing)
            delete_colon!(lines, position)
        end
        lines[position] = anonymous_function(lines[position])
        lines[position] = ranges(lines[position], collect=collect)
        for pair in PYTHON2JULIA_WORDS
            lines[position] = replace(lines[position], pair)
        end
        position += 1
    end
    lines
end

function python2julia(filename_input::String,
        filename_output::String = replace(filename_input, r"\.py$" => ".jl");
        kargs...)

    lines = python2julia(readlines(filename_input); kargs...)
    open(filename_output, "w") do fh
        for line in lines
            println(fh, line)
        end
    end
    lines
end

end # module


text:
#=
rotate_string:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-30
=#

module Rotation
    function rotate_left(str, num)
        str_rot = join([str[num+1:end], str[1:num]], "")
        return str_rot
    end

    function rotate_right(str, num)
        return rotate_left(str, length(str) - num)
    end

    function main()
        println(rotate_left("GeeksForGeeks", 2))
        println(rotate_right("GeeksForGeeks", 2))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Rotation
    Rotation.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using DrWatson

quickactivate(@__DIR__)

using PolarizationFramework
using DifferentialEquations

g = 5
attr = BinaryAttributes(g)

calc_heider_attr(45, 0.0, 1000.0, Heider7!, AutoTsit5(Rodas5(autodiff = false)), true, attr)

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
module DOT

using ParserCombinator
using AutoHashEquals
import Base: ==

export Statement, Statements, ID, StringID, NumericID, HtmlID, Attribute,
       Attributes, Graph, Port, NodeID, Node, EdgeNode, Edge, GraphAttributes,
       NodeAttributes, EdgeAttributes, SubGraph, parse_dot, nodes, edges

# i've gone with a very literal parsing, which returns a structure that is
# pretty much what is described in the grammar at
# http://www.graphviz.org/content/dot-language.  this is because it is not
# clear to me exactly how the syntax is connected to the semantics of graph
# layout.  for example, a simple graph like the one at
# http://www.graphviz.org/Gallery/gradient/colors.gv.txt has the node
# attributes specified before the node ID.  if you look at the syntax
# carefully that means that the attributes are not associated with that
# particular node (such attributes come after the node ID) and instead seem to
# be "global" node attributes in some sense.  so the format is more like a
# "program" (where order is important, with the "global" node attributes
# changing as the file is "processed") than a specification.

# see test/dot/examples.jl for examples accessing fields in this structure

abstract type Statement end

const Statements = Vector{Statement}

abstract type ID end

@auto_hash_equals struct StringID <: ID
    id::AbstractString
end

@auto_hash_equals struct NumericID <: ID
    id::AbstractString
end

@auto_hash_equals struct HtmlID <: ID
    id::AbstractString
end

@auto_hash_equals struct Attribute <: Statement
    name::ID
    value::ID
end

const Attributes = Vector{Attribute}

@auto_hash_equals struct Graph
    strict::Bool
    directed::Bool
    id::Union{ID, Nothing}
    stmts::Statements
    Graph(s::Bool, d::Bool, id::ID, st::Statements) = new(s, d, id, st)
    Graph(s::Bool, d::Bool, st::Statements) = new(s, d, nothing, st)
end

@auto_hash_equals struct SubGraph <: Statement
    id::Union{ID, Nothing}
    stmts::Statements
    SubGraph(id::ID, s::Statements) = new(id, s)
    SubGraph(s::Statements) = new(nothing, s)
end

@auto_hash_equals struct Port
    id::Union{ID, Nothing}
    point::Union{AbstractString, Nothing}
    Port(id::ID, p::Union{AbstractString, Nothing}=nothing) = new(id, p)
    Port(p::AbstractString) = new(nothing, p)
end

@auto_hash_equals struct NodeID
    id::ID
    port::Union{Port, Nothing}
    NodeID(id::ID, p::Union{Port, Nothing}=nothing) = new(id, p)
end

@auto_hash_equals struct Node <: Statement
    id::NodeID
    attrs::Attributes
    Node(id::NodeID, a::Attributes) = new(id, a)
    Node(id::NodeID) = new(id, Attribute[])
end

const EdgeNode = Union{NodeID, SubGraph}
const EdgeNodes = Vector{EdgeNode}

@auto_hash_equals struct Edge <: Statement
    nodes::EdgeNodes
    attrs::Attributes
    Edge(n::EdgeNodes, a::Attributes) = new(n, a)
    Edge(n::EdgeNodes) = new(n, Attribute[])
end

@auto_hash_equals struct GraphAttributes <: Statement
    attrs::Attributes
end

@auto_hash_equals struct NodeAttributes <: Statement
    attrs::Attributes
end

@auto_hash_equals struct EdgeAttributes <: Statement
    attrs::Attributes
end

@with_names begin

    # x can include alternatives that can also be discarded
    mkspc(x) = "([ \t]|\n#.*|\n|//.*|/\\*(.|\n)*?\\*/$x)*"

    spc = mkspc("")
    spc_init = ~Pattern(string("(#.*)?", spc))
    spc_star = ~Pattern(spc)

    wrd = p"[a-zA-Z\200-\377_][a-zA-Z\200-\377_0-9]*"

    # valid strings include:
    # 1 - "..."
    # 2 - "...\
    #      ..."
    # 3 = "..." + "..."
    # and all can contain quoted quotes
    unesc(s) = replace(s, "\\\"" => "\"")
    unesc_join(s...) = string(map(unesc, s)...)
    str_cont = Pattern("((?:[^\"\n]|\\\\\")*)\\\\\\\n", 1)
    str_end = p"([^\"\n]|\\\\\")*"
    str_one = Seq!(E"\"", Star!(str_cont), str_end, E"\"")
    str_many = Seq!(str_one, Star!(Seq!(spc_star, E"+", spc_star, str_one))) > unesc_join

    str_id = Alt!(wrd, str_many) > StringID

    num_id = p"-?(\.[0-9]+|[0-9]+(\.[0-9]*)?)" > NumericID

    # this is bracketed by <> (in addition to the xml).  not clear from
    # the grammar, but see test/dot/tictactoe.dot
    html_id = Pattern("<((:?[ \t\n]|<(:?[^<>]|\n)+>)*)>", 1) > HtmlID

    id = Alt!(str_id, num_id, html_id)

    cmp = p"(n|ne|e|se|s|sw|w|nw|c|_)"
    col = E":"

    # port grammar seeems to be ambiguous, since :ID could be :point
    # this is a best guess at what was meant
    port = Alt!(Seq!(col, spc_star, id, spc_star, col, spc_star, cmp),
                Seq!(col, spc_star, cmp),
                Seq!(col, spc_star, id)) > Port

    # this is a raw ID=ID, not the attr_stmt in the grammar
    attr = Seq!(id, spc_star, E"=", spc_star, id) > Attribute

    spc_attr = ~Pattern(mkspc("|;|,"))
    # for some reason 0.3 does not like |> Attributes here
    attr_list = PlusList!(Seq!(E"[", StarList!(attr, spc_attr), E"]"), spc_star) |> (x -> Attribute[x...])

    node_id = Seq!(id, spc_star, Opt!(port)) > NodeID
    node_stmt = Seq!(node_id, spc_star, Opt!(attr_list)) > Node

    NoCase(s) = Pattern(join(["[$(lowercase(c))$(uppercase(c))]" for c in s]))

    stmt = Delayed()

    # not a list because we can have a trailing ;
    # this eats trailing spaces, but i don't think it matters
    # the comma (",") below is not in the grammar at graphviz.org but is
    # needed to parse the examples in test/examples.jl
    # related, note that using a comma in this way seems to trigger bugs
    # in dot itself - see
    # https://github.com/JuliaGraphs/LightGraphs.jl/issues/107#issuecomment-131401430
    spc_stmt = ~Pattern(mkspc("|;|,"))
    # |> Statements but for 0.3
    stmt_list = Star!(Seq!(stmt, spc_stmt)) |> (x -> Statement[x...])
    stmt_brak = Seq!(E"{", spc_star, stmt_list, spc_star, E"}")

    sub_graph = Seq!(Opt!(~NoCase("subgraph")), spc_star, Opt!(id), spc_star, stmt_brak) > SubGraph

    # order important here, since we don't backtrack and "subgraph" could
    # be a node
    edge_node = Alt!(sub_graph, node_id)
    edge_sep = Seq!(spc_star, P"(--|->)", spc_star)
    # |> EdgeNodes but for 0.3
    edge_list = Seq!(edge_node, edge_sep, PlusList!(edge_node, edge_sep)) |> (x -> EdgeNode[x...])
    edge_stmt = Seq!(edge_list, spc_star, Opt!(attr_list)) > Edge

    attr_stmt = Alt!(Seq!(~NoCase("graph"), spc_star, attr_list) > GraphAttributes,
                     Seq!(~NoCase("node"), spc_star, attr_list) > NodeAttributes,
                     Seq!(~NoCase("edge"), spc_star, attr_list) > EdgeAttributes)

    # order important here as node_stmt can match almost anything
    stmt.matcher = Alt!(edge_stmt, attr_stmt, attr, sub_graph, node_stmt)

    strict = Alt!(Seq!(~NoCase("strict"), Insert(true)), Insert(false))
    direct = Alt!(Seq!(~NoCase("digraph"), Insert(true)),
                  Seq!(~NoCase("graph"), Insert(false)))
    graph = Seq!(strict, spc_star, direct, spc_star, Opt!(id), spc_star, stmt_brak) > Graph

    dot = Seq!(spc_init, Plus!(Seq!(graph, spc_star)), Eos())

end

# the file structured using the types above (returns an array of graphs)
function parse_dot(s; debug=false)
    try
        if debug
            parse_one_dbg(s, Trace(dot); debug=true)
        else
            parse_one(s, dot)
        end
    catch x
        if debug
            Base.show_backtrace(stdout, catch_backtrace())
        end
        rethrow()
    end
end

# set of all nodes
# expand both nodes and edges, then use sets to de-duplicate
nodes(g::Graph) = reduce(union, map(nodes, g.stmts))
nodes(s::Statement) = Set()
nodes(s::SubGraph) = reduce(union, map(nodes, s.stmts))
nodes(n::Node) = nodes(n.id)
nodes(n::NodeID) = Set([n.id.id])
nodes(e::Edge) = reduce(union, map(nodes, e.nodes))

# set of all node pairs that correspond to edges
fix_directed(g::Graph, e) = g.directed ? Set(e) : Set([Tuple(sort(collect(p))) for p in e])
edges(g::Graph) = fix_directed(g, reduce(vcat, map(edges_, g.stmts)))
edges_(s::Statement) = []
edges_(s::SubGraph) = reduce(vcat, map(edges_, s.stmts))
pair(n1::NodeID, n2::NodeID) = [(n1.id.id, n2.id.id)]
pair(n::NodeID, s::SubGraph) = vcat([(n.id.id, x) for x in nodes(s)], edges_(s))
pair(s::SubGraph, n::NodeID) = vcat([(x, n.id.id) for x in nodes(s)], edges_(s))
pair(s1::SubGraph, s2::SubGraph) =
    vcat(vec([(x, y) for x in nodes(s1), y in nodes(s2)]), edges_(s1), edges_(s2))
edges_(e::Edge) = reduce(vcat, [pair(a, b) for (a, b) in zip(e.nodes[1:end-1], e.nodes[2:end])])

end


text:
using WebDriver

capabilities = Capabilities("firefox")
wd = RemoteWebDriver(capabilities)
session = Session(wd)

function waitForPageLoaded(session)
    sleep(1)
    pageState = try
        script!(session, "return document.readyState")
    catch
        ""
    end
    for i = 1:30
        if pageState != "complete"
            sleep(1)
        end
    end
end

url = "https://proceedings.neurips.cc/"
navigate!(session, url)
waitForPageLoaded(session)
tableElem = Element(session, "css selector", "div.col-sm")
linkElems = Elements(tableElem, "css selector", "a")
links = [element_attr(linkElem, "href") for linkElem in linkElems]

file = open("DualAspectTopicRepresentation/data/NIPS.txt", "w")
for url in links
    textOfYear = split(url,"/")[end]
    navigate!(session, url)
    waitForPageLoaded(session)
    tableElem = Element(session, "css selector", "div.col")
    sublinkElems = Elements(tableElem, "css selector", "a")
    for sublinkElem in sublinkElems
        title = replace(replace(element_text(sublinkElem), "\n"=>" "), "￼"=>"")
        println(file, "$textOfYear. $title")
    end
end
close(file)

max_stars_repo_path: src/01 collect_data.jl

max_issues_repo_path: src/01 collect_data.jl

max_forks_repo_path: src/01 collect_data.jl


text: include("recalibration.jl")


text:
using StaticPermutations

using Test

@testset "StaticPermutations.jl" begin
    perm = Permutation(2, 3, 1)
    noperm = NoPermutation()

    @testset "Constructors" begin
        @test perm === Permutation((2, 3, 1))
        @inferred (() -> Permutation(2, 3, 1))()
        @inferred (() -> Permutation((2, 3, 1)))()
        @inferred NoPermutation()
    end

    @testset "getindex!" begin
        @test perm[2] == perm[Val(2)] == 3
        @test noperm[2] == noperm[Val(2)] == 2
        valgettwo(p) = Val(p[2])
        @inferred valgettwo(perm)
        @inferred valgettwo(noperm)
    end

    @testset "I/O" begin
        @test string(perm) == "Permutation(2, 3, 1)"
        @test string(noperm) == "NoPermutation()"
    end

    @testset "Base operations" begin
        @test Tuple(perm) === (2, 3, 1)
        @test Tuple(noperm) === nothing

        @test length(perm) === 3
        @test length(noperm) === nothing
    end

    @testset "Permutation checks" begin
        @test isperm(perm)
        @test isperm(noperm)
        @test !isperm(Permutation(2, 5, 3))

        # Check that result is fully inferred
        ispermval(p) = Val(isperm(p))
        @inferred ispermval(perm)
        @inferred ispermval(noperm)
        @inferred ispermval(Permutation(2, 5, 3))
    end

    @testset "Composition" begin
        p = Permutation(2, 3, 1)
        q = Permutation(1, 3, 2)
        @inferred p * q
        @inferred p * NoPermutation()
        @inferred NoPermutation() * p
        @test p * q === Permutation(3, 2, 1)
        @test q === p \ Permutation(3, 2, 1)
        @test q * p === Permutation(2, 1, 3)
        @test p === q \ Permutation(2, 1, 3)
        @test p * NoPermutation() === p
        @test NoPermutation() * p === p
        @test NoPermutation() * NoPermutation() === NoPermutation()
        @test p * inv(p) == NoPermutation()
        @test inv(p) * p == NoPermutation()
    end

    iperm = inv(perm)
    @testset "Inverse permutation" begin
        @test perm * iperm == NoPermutation()
        @test iperm * perm == NoPermutation()
        @test invperm(perm) === iperm
        @inferred inv(perm)
    end

    @testset "Identity permutation" begin
        @inferred identity_permutation(Val(4))
        id = identity_permutation(Val(4))
        @test id === Permutation(1, 2, 3, 4)
        @test id == NoPermutation()  # they're functionally equivalent
        idval(p) = Val(isidentity(p))
        @test (@inferred idval(perm)) === Val(false)
        @test (@inferred idval(noperm)) === Val(true)
        @test (@inferred idval(id)) === Val(true)
        @inferred idval(perm)
        @test !isidentity(perm)
        @test isidentity(NoPermutation())
        @test isidentity(id)
        @test isidentity(perm * iperm)
    end

    @testset "Left division operator" begin
        a = Permutation(2, 3, 1, 4)
        b = Permutation(3, 1, 4, 2)
        @inferred b / a
        r = b / a
        @test r * a == b
        np = NoPermutation()
        @test a / np === a
        @test np / a === inv(a)
        @test np / np === np
    end

    @testset "Comparisons" begin
        @test Permutation(1, 3, 2) == Permutation(1, 3, 2)
        @test Permutation(1, 3, 2) != Permutation(3, 1, 2)
        @test Permutation(2, 1, 3) != Permutation(2, 1)
        @test Permutation(1, 2, 3) == NoPermutation()
        @test NoPermutation() == Permutation(1, 2, 3)
        @test Permutation(1, 3, 2) != NoPermutation()
        @test NoPermutation() == NoPermutation()
    end

    @testset "Apply permutations" begin
        ind = (20, 30, 10)
        ind_perm = (30, 10, 20)
        permval(perm) = Val(perm * (:aaa, :bbb, :ccc))
        @inferred permval(perm)
        @inferred permval(iperm)
        @inferred permval(noperm)
        @test noperm * ind === ind
        @test ind === noperm \ ind
        @test perm * ind === ind_perm
        @test ind === perm \ ind_perm
        @test perm * CartesianIndex(ind) === CartesianIndex(ind_perm)
        @test noperm * CartesianIndex(ind) === CartesianIndex(ind)
    end

    @testset "Prepend / append" begin
        @inferred prepend(Permutation(2, 3, 1), Val(2))
        @inferred append(Permutation(2, 3, 1), Val(2))
        @test prepend(Permutation(2, 3, 1), Val(2)) === Permutation(1, 2, 4, 5, 3)
        @test prepend(NoPermutation(), Val(2)) === NoPermutation()
        @test append(Permutation(2, 3, 1), Val(2)) === Permutation(2, 3, 1, 4, 5)
        @test append(NoPermutation(), Val(2)) === NoPermutation()
    end

    @testset "PermutedDimsArray" begin
        x = rand(3, 5, 4)
        @inferred PermutedDimsArray(x, perm)
        @inferred PermutedDimsArray(x, noperm)

        # Compare new and original constructors
        @test PermutedDimsArray(x, perm) === PermutedDimsArray(x, Tuple(perm))
        @test PermutedDimsArray(x, noperm) === PermutedDimsArray(x, (1, 2, 3))
    end
end


text:
using StochasticBlockModel
using Test

@testset "Optimal assignments" begin

    A = [0 1 1 0 0 0 0 0;
         1 0 0 1 0 0 0 0;
         1 0 0 1 0 0 0 0;
         0 1 1 0 0 0 0 0;
         0 0 0 0 0 1 1 0;
         0 0 0 0 1 0 0 1;
         0 0 0 0 1 0 0 1;
         0 0 0 0 0 1 1 0]
    dataset = StochasticBlockModel.Dataset(A, 2)
    sbm = StochasticBlockModel.SBM([2.0 0.0; 0.0 2.0], "poisson")
    opt_results, x = StochasticBlockModel.optimalAssignments(dataset, sbm)
    @test x == [1 0; 1 0; 1 0; 1 0; 0 1; 0 1; 0 1; 0 1]
    @test opt_results.UB ≈ 2.4548225555204373 rtol = 1e-5
    @test opt_results.LB ≈ 2.4548225555204373 rtol = 1e-5
    @test opt_results.status == :Optimal

    estimator = StochasticBlockModel.SBMEstimator("ls1", 10.0, true, 1, false)
    sbm = StochasticBlockModel.SBM([1.3 0.2; 0.2 0.8], "poisson")
    opt_results1, x1 = StochasticBlockModel.optimalAssignments(estimator, dataset, sbm)
    opt_results2, x2 = StochasticBlockModel.optimalAssignments(dataset, sbm, verbose=true)
    @test x1 == x2
end


text:
using Test
using Comonicon.PATH
using Comonicon.Configurations
using Comonicon.Configurations: Install, Application, SysImg, Download, Precompile

module XYZ end

@test read_configs(XYZ) == Configurations.Comonicon(
    name="xyz",
    install = Install(
        path = "~/.julia",
        completion = true,
        quiet = false,
        compile = "yes",
        optimize = 2,
    ),
    sysimg = nothing,
    download = nothing,
    application = nothing,
)

@test read_configs(XYZ; name="zzz") == Configurations.Comonicon(
    name="zzz",
    install = Install(
        path = "~/.julia",
        completion = true,
        quiet = false,
        compile = "yes",
        optimize = 2,
    ),
    sysimg = nothing,
    download = nothing,
    application = nothing,
)

@test read_configs(XYZ; path="mypath") == Configurations.Comonicon(
    name="xyz",
    install = Install(
        path = "mypath",
        completion = true,
        quiet = false,
        compile = "yes",
        optimize = 2,
    ),
    sysimg = nothing,
    download = nothing,
    application = nothing,
)

@test_throws ArgumentError read_configs(XYZ; filter_stdlibs=true)
@test_throws ArgumentError Application(;path=pwd())
@test_throws ArgumentError SysImg(;path=pwd())

read_configs(XYZ; user="Roger-luo", repo="Foo") == Configurations.Comonicon(
    name="xyz",
    install = Install(
        path = "~/.julia",
        completion = true,
        quiet = false,
        compile = "yes",
        optimize = 2,
    ),
    sysimg = nothing,
    download = Download(
        user="Roger-luo",
        repo="Foo",
    ),
    application = nothing,
)

@test_throws ArgumentError read_configs(XYZ; abc=2)

@test read_configs(PATH.project("test", "Foo", "Comonicon.toml")) == Configurations.Comonicon(
    name = "foo",
    install = Install(
      path = "~/.julia",
      completion = true,
      quiet = false,
      compile = "min",
      optimize = 2,
    ),
    sysimg = SysImg(
      path = "deps",
      incremental = true,
      filter_stdlibs = false,
      cpu_target = "native",
      precompile = Precompile(
        execution_file = ["precopmile.jl"],
        statements_file = String[],
      ),
    ),
    download = Download(
      host = "github.com",
      user = "Roger-luo",
      repo = "Foo.jl",
    ),
    application = nothing,
  )

@test_throws ErrorException read_configs(PATH.project("test"))


text:
# test_GLS.jl

# Write tests corresponding to GLS here.

using PDESolver
using AdvectionEquationMod
using Input
using Base.Test
using ODLCommonTools
using ArrayViews

fname = "input_vals_GLS.jl"
opts = read_input(fname)
@testset "--- Check input arguments applied correctly ---" begin
  @test ( opts["use_GLS"] )== true
  @test ( opts["Q_transpose"] )== true
end

mesh, sbp, eqn, opts = solvePDE(fname)

@testset "--- Check functions in ../src/solver/advection/GLS.jl ---" begin

  @testset "Checking shape function derivatives" begin
    sbp2 = getTriSBPGamma{Float64}()
    shapefuncderiv = zeros(Tsol, sbp2.numnodes, sbp2.numnodes, Tdim)
    calcShapefuncDeriv(sbp2, shapefuncderiv)
    @test isapprox(shapefuncderiv[:,:,1] -  [-0.5  0.5  0.0
                                             -0.5  0.5  0.0
                                             -0.5  0.5  0.0]) atol=1e-13
    @test isapprox(shapefuncderiv[:,:,2] -  [-0.5  0.0  0.5
                                             -0.5  0.0  0.5
                                             -0.5  0.0  0.5]) atol=1e-13
  end

  @testset "Check calcAxidxi" begin
    Tdim = 2
    alpha_x = ones(Tsol, mesh.numNodesPerElement)
    alpha_y = ones(alpha_x)
    dxidx = sview(mesh.dxidx,:,:,:,1)
    shapefuncderiv = zeros(Tsol, sbp.numnodes, sbp.numnodes, Tdim)
    calcShapefuncDeriv(sbp, shapefuncderiv)
    AxiDxi = calcAxiDxi(mesh, dxidx, alpha_x, alpha_y, shapefuncderiv)
    @test isapprox(AxiDxi -[-0.5 0.25 0.25
                            -0.5 0.25 0.25
                            -0.5 0.25 0.25]) atol=1e-13
  end

  @testset "Check calcTau" begin
    fill!(eqn.params.alpha_x, 1.0)
    fill!(eqn.params.alpha_y, 1.0)
    tau = zeros(Tsol,mesh.numNodesPerElement, mesh.numEl)
    calcTau(mesh, sbp, eqn, tau)
    @test ( mesh.numEl )== 1
    @test isapprox( tau[:,1], [1.0,1.0,1.0])
  end

  @testset "Check GLS" begin
    fill!(eqn.params.alpha_x, 1.0)
    fill!(eqn.params.alpha_y, 1.0)
    fill!(eqn.res, 0.0)
    eqn.q[1,:,1] = eqn.q_vec
    GLS(mesh, sbp, eqn)
    @test isapprox( eqn.res[1,:,1], [8.0 -4.0 -4.0])

  end

end  # end facts


text:
"""
    rv(X, Y)
Compute the sample RV coefficient between matrices X and Y
* `X` : Matrix (n obs., p variables).
* `Y` : Matrix (n obs., q variables).

RV is bounded in [0, 1].

A dissimilarty measure between `X` and `Y` can be computed
by d = sqrt(2 * (1 - RV)).

## References

Escoufier, Y., 1973. Le Traitement des Variables Vectorielles. Biometrics 29, 751–760.
https://doi.org/10.2307/2529140

Josse, J., Pagès, J., Husson, F., 2008. Testing the significance of the RV coefficient.
Computational Statistics & Data Analysis 53, 82–91. https://doi.org/10.1016/j.csda.2008.06.012

Mayer, C.-D., Lorent, J., Horgan, G.W., 2011. Exploratory Analysis of Multiple Omics
Datasets Using the Adjusted RV Coefficient. Statistical Applications in Genetics and Molecular
Biology 10. https://doi.org/10.2202/1544-6115.1540

Smilde, A.K., Kiers, H.A.L., Bijlsma, S., Rubingh, C.M., van Erk, M.J., 2009.
Matrix correlations for high-dimensional data: the modified RV-coefficient.
Bioinformatics 25, 401–405. https://doi.org/10.1093/bioinformatics/btn634

Robert, P., Escoufier, Y., 1976. A Unifying Tool for Linear Multivariate Statistical Methods:
The RV-Coefficient. Journal of the Royal Statistical Society: Series C (Applied Statistics)
25, 257–265. https://doi.org/10.2307/2347233
"""

function rv(X, Y; centr = true)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    if centr
        X = center(X, mean(X, dims = 1))
        Y = center(Y, mean(Y, dims = 1))
    end
    YtX = Y' * X
    XtX = X' * X
    YtY = Y' * Y
    A = dot(YtX, YtX)
    B = dot(XtX, XtX)
    C = dot(YtY, YtY)
    rv = A / sqrt(B * C)
end


text:
import Pkg

include("templates.jl")

const LynxRepo = "https://github.com/jorge-brito/Lynx.jl"

function newapp(name::String, deps::Vector = [])
    path = joinpath(@__DIR__, "..", "src", name)
    cd(joinpath(@__DIR__, "..", "src"))
    Pkg.generate(name)
    @info "New app '$name' is located at $path"
    cd(path)

    write("README.md", README(name))
    write(joinpath("src", "$name.jl"), ENTRY_POINT(name))
    mkdir(".vscode")
    write(".vscode/tasks.json", VSCODE_TASKS(name))

    if !isempty(deps)
        Pkg.activate(".")
        for dep in deps
            @info "Adding dependency '$dep'"
            Pkg.add(dep)
        end
        Pkg.add(url=LynxRepo)
    end
end

if !isempty(ARGS)
    name = ARGS[1]
    deps = ARGS[2:end]
    if Meta.isidentifier(name)
        newapp(name, [deps..., "Luxor", "Colors"])
    else
        @error "The new app must have a valid name."
    end
else
    @error """No name argument specified

    Usage:

        julia ./newapp.jl <name> [dependencies]

    Example:

        julia ./newapp.jl MyApp
        julia ./newapp.jl AnotherApp Dep1 Dep2 Dep3

    """
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
module CVMData

using HTTP
using YearMonths

include("funds.jl")

end # module


text:
# ### m11.5t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

delim=';'
df = CSV.read(sr_datadir("UCBadmit.csv"), DataFrame; delim);

# Turing model

@model ppl11_5(admit, applications) = begin
    N=length(applications)
    θ ~ truncated(Exponential(1), 0, Inf)
    α ~ Normal(0,2)

    for i ∈ 1:N
        prob = logistic(α)

        # alpha and beta for the BetaBinomial must be provided.
        # The two parameterizations are related by
        # alpha = prob * theta, and beta = (1-prob) * theta.
        # See https://github.com/rmcelreath/rethinking/blob/master/man/dbetabinom.Rd
        alpha = prob * θ
        beta = (1 - prob) * θ

        admit[i] ~ BetaBinomial(applications[i], alpha, beta)
    end
end

# Sample

m11_5t = ppl11_5(df.admit, df.applications)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns11_5t = mapreduce(c -> sample(m11_5t, sampler, nsamples), chainscat, 1:nchains)

# Result rethinking

m11_5s = "
         mean   sd  5.5% 94.5% n_eff Rhat
θ        2.74 0.96  1.43  4.37  3583    1
a       -0.37 0.31 -0.87  0.12  3210    1
";

# End of `11/m811.5t.jl`


text:
# structure borrowed from Metalhead.jl's vgg19.jl
abstract type AutoEncoderModel end