text:
function solve(model::Model; initial_values = nothing, output_writer = nothing, log = false)
    # TimerOutput handling
    reset_timer!()
    (log == false) && disable_timer!()

    # Reinitialize clock and output_writer # NO NEED, FORCE LOADING ALL INPUT FILE BEFORE SOLVING
    #reinit!(model.clock)
    #reinit!(output_writer,model)

    # iterate over time
    u = iterate(model,output_writer,initial_values)

    # print timings of code annotated with @timeit macro
    (log == true) && print_timer(title = "Analysis with $(getncells(model.grid)) elements", linechars = :ascii)
    println(model.material_state[1][1].temp_σ)
    return model, u
end

function iterate(model::Model{2,1,Nothing,Nothing,E,Nothing}, output_writer, initial_values = nothing) where {E<:Elasticity}
    @info "Rheology is purely elastic. Displacement field is solved once at the end of the requested time interval"
    c = model.clock
    c.current_time = c.tspan[2]
    # number of base functions per element
    nbasefuncs = getnbasefunctions.(model.cellvalues_tuple)

    # assembly of
    @timeit "assemble" doassemble!(model,nbasefuncs)

    #Apply Dirichlet boundary conditions
    update!(model.dirichlet_bc, c.current_time)
    apply!(model.K, model.RHS, model.dirichlet_bc)

    # Solve
    #u = Symmetric(model.K) \ model.RHS
    @timeit "linear_solve" u = linear_solve(model.K, model.RHS, model.solver)

    # check convergence :
    #fdofs = JuAFEM.free_dofs(dbc) no need for free dofs
    norm_res = norm(model.RHS - model.K*u)
    @info "residual L2 norm : $norm_res"
    # update material states :
    fill_state_from_u!(model,u) # fill temporary stress and strains
    update_material_state!(model) # set temporary values to converged type fields

    # output
    @timeit "export" write_output!(model, u, output_writer)

    return u
end
function iterate(model::Model{2,2,Nothing,Nothing,TE,Nothing}, output_writer, initial_values = nothing) where {TE<:Elasticity}
    @info "Rheology is purely elastic. Displacement field is solved once at the end of the requested time interval"
    c = model.clock
    c.current_time = c.tspan[2]
    # number of base functions per element
    nbasefuncs = getnbasefunctions.(model.cellvalues_tuple)

    # assembly of
    @timeit "assemble" doassemble!(model,nbasefuncs)

    #Apply Dirichlet boundary conditions
    update!(model.dirichlet_bc, c.current_time)
    apply!(model.K, model.RHS, model.dirichlet_bc)

    # Solve
    #u = Symmetric(model.K) \ model.RHS
    @timeit "linear_solve" u = linear_solve(model.K, model.RHS, model.solver)

    # update material states :
    fill_state_from_u!(model,u) # fill temporary stress and strains
    update_material_state!(model) # set temporary values to converged type fields

    # output
    @timeit "export" write_output!(model, u, output_writer)

    return u
end

function iterate(model::Model{2,2,Nothing,TV,TE,Nothing}, output_writer, initial_values = nothing) where {TV,TE<:Elasticity}
    @info "Rheology is purely elastic. Displacement field is solved once at the end of the requested time interval"
    clock = model.clock
    dbc, dh = model.dirichlet_bc, model.dofhandler
    # number of base functions per element
    nbasefuncs = getnbasefunctions.(model.cellvalues_tuple)

    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)
    u_converged = copy(u) # backup solution vector

    while clock.current_time <= clock.tspan[2]

        timestep!(clock) # update clock

        update!(dbc, clock.current_time) # evaluates the Dirichlet-bcs at time t
        apply!(u, dbc)  # set the prescribed dbcs values in the solution vector
        # assembly of
        @timeit "assemble" doassemble!(model, nbasefuncs, u, u_converged)

        #Apply Dirichlet boundary conditions
        apply!(model.K, model.RHS, dbc)

        # Solve
        #u = Symmetric(model.K) \ model.RHS
        @timeit "linear_solve" u .= linear_solve(model.K, model.RHS, model.solver)

        # update material states :
        fill_state_from_u!(model,u) # fill temporary stress and strains
        update_material_state!(model) # set temporary values to converged type fields

        # output
        @timeit "export" write_output!(model, u, output_writer)
        u_converged .= u
    end

    return u
end
#iterate(model::Model{2,1,Nothing,Nothing,E,Nothing}, output_writer, initial_values = nothing) where {E<:Elasticity} = iterate(model::Model{2,N,Nothing,Nothing,E,Nothing}, output_writer, initial_values) where {N,E<:Elasticity}
function iterate(model::Model{2,2,D,V,E,P}, output_writer, initial_values = nothing) where {N,D,V<:Viscosity,E,P<:ViscousDruckerPrager}

    # Unpack some model fields
    dh, dbc, cv, clock = model.dofhandler, model.dirichlet_bc, model.cellvalues_tuple, model.clock

    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)

    u_converged = copy(u) # backup solution vector
    δu = similar(u)

    while clock.current_time <= clock.tspan[2]
        @timeit "time iteration" begin

            timestep!(clock,Δt_max_damage) # update clock

            print("\n TIME ITERATION $(clock.iter)\n",
            " current simulation time = $(clock.current_time):\n",
            " timestep = $(clock.Δt)\n")

            restart_flag = false

            # Apply dirichlet bc and iteratively solve for u :
            # @timeit "nonlinear solve"
            restart_flag = nonlinear_solve!(u,u_converged,δu,model,restart_flag)

            if restart_flag == true
                u .= u_converged
                undo_timestep!(clock)
                clock.Δt *= clock.Δt_fact_down # decreased timestep
            else # converged
                update_material_state!(model) # update converged state values
                @timeit "export" write_output!(model, u, output_writer) # output

                clock.Δt *= clock.Δt_fact_up # increase timestep
                println(clock.Δt)
                u_converged .= u
            end

            clock.current_time == clock.tspan[2] && break # end time loop if requested end time is reached

        end
    end
    return u
end

function iterate(model::Model{2,2,D,V,E,P}, output_writer, initial_values = nothing) where {N,D<:Damage,V,E,P}

    # Unpack some model fields
    dh, dbc, cv, clock = model.dofhandler, model.dirichlet_bc, model.cellvalues_tuple, model.clock

    # Pre-allocate solution vectors, etc.
    dofs_D = get_field_dofs(:D,model)
    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)
    if initial_values == nothing
        u[dofs_D] .= log.(u[dofs_D].+1e-9) # ensure non zero starting Damage
    else
        set_initial_solution_vector!(u,dh,initial_values)
        u[dofs_D] .= log.(u[dofs_D])
    end

    u_converged = copy(u) # backup solution vector
    δu = similar(u)

    ###### TEST
    println("D initial extrema = ", extrema(exp.(u[dofs_D])))
    vtk_grid("TEST_u-logD_initial_fields", model.dofhandler) do vtkfile
        u2 = copy(u)
        u2[dofs_D] .= exp.(u2[dofs_D]) # coming back to damage from log(damage)
        vtk_point_data(vtkfile, model.dofhandler, u2)
    end
    ######

    Δt_max_damage = 1e9
    restart_flag = false

    while clock.current_time <= clock.tspan[2]
        @timeit "time iteration" begin

            # tdt = @elapsed Δt_max_damage = get_damage_constrained_Δt(model,u,0.3)
            # println("get_damage_constrained_Δt duration : ", tdt)
            # println("Δt max damage from last stress = ", Δt_max_damage)

            timestep!(clock,Δt_max_damage) # update clock

            print("\n TIME ITERATION $(clock.iter)\n",
            " current simulation time = $(clock.current_time):\n",
            " timestep = $(clock.Δt)\n")

            if restart_flag == true
                restart_flag = false
            end

            # Apply dirichlet bc and iteratively solve for u :
            # @timeit "nonlinear solve"
            restart_flag = nonlinear_solve!(u,u_converged,δu,model,restart_flag)
            println("restart_flag post nonlinear_solve : ", restart_flag)
            if restart_flag == true
                u .= u_converged
                undo_timestep!(clock)
                clock.Δt *= clock.Δt_fact_down # decreased timestep
            elseif restart_flag == :activate_damage
                undo_timestep!(clock)
                Δt_max_damage = get_damage_constrained_Δt(model,u,0.3)
                println("restarting with damage growth with constraint Δt of ", Δt_max_damage, "\n")
                u .= u_converged
            elseif restart_flag == false # converged
                update_material_state!(model) # update converged state values

                #test#
                display(model.material_state[1][1])
                ######

                @timeit "export" write_output!(model, u, output_writer) # output

                clock.Δt *= clock.Δt_fact_up # increase timestep
                println(clock.Δt)
                u_converged .= u

                #TEST#
                if clock.iter in (2,3)
                    #println("u cell 1 value after convergence and export : ")
                    #display(u[1:22])
                end
                ######
            end

            # ##### NLsolve VERSION :#####
            # perform_elastic_solve!(u,model)
            #
            # ###### TEST
            # # vtk_grid("TEST_u-D_after_elastic_solve_it_$(clock.iter)", model.dofhandler) do vtkfile
            # #     vtk_point_data(vtkfile, model.dofhandler, u)
            # # end
            # #write_output!(model, u, output_writer) # output
            # #println("temp_σyy = ", model.material_state[15][1].temp_σ[2,2])
            # #println("temp_ϵyy = ", model.material_state[15][1].temp_ϵ[2,2])
            # ######
            #
            # f!(res,u) = doassemble_res!(res, model, nbasefuncs, u, u_converged)
            # tt = @elapsed (sol = NLsolve.nlsolve(f!, u, method = :newton))
            # println("nlsolve time = ", tt)
            #
            # if NLsolve.converged(sol) == false
            #     println(sol)
            #     u .= u_converged
            #     undo_timestep!(clock)
            #     clock.Δt *= clock.Δt_fact_down # decreased timestep
            # else # converged
            #     u .= sol.zero # converged solution
            #     println(sol)
            #     update_material_state!(model) # update converged state values
            #     @timeit "export" write_output!(model, u, output_writer) # output
            #
            #     clock.Δt *= clock.Δt_fact_up # increase timestep
            #     println(clock.Δt)
            #     u_converged .= u
            # end

            ############################

            clock.current_time == clock.tspan[2] && break # end time loop if requested end time is reached

        end
    end
    return u
end

# TODO check the performances penalty magnitude of such an approach : one solution could
# be to build a container containing most EP_model fields and to reuse it for each subsequent elastic solve
# "The current strategy is to convert our damaged model into an equivalent undamaged model and to use the
# existing linear elastic solves implemented for undamaged rheologies"
# function perform_elastic_solve!(u,model)
#
#     tt = @elapsed begin
#     # model parameters
#     u_interp = typeof(model.dofhandler.field_interpolations[1]).parameters[3]
#     el_geom = typeof(model.grid.cells[1])
#     variables = PrimitiveVariables{1}((:u,), (u_interp,), el_geom)
#     dh = create_dofhandler(model.grid, variables)
#     dbc = modify_dirichlet_bc(dh, model.dirichlet_bc)
#     mp = convert_to_undamaged_material_properties(model.grid,model.material_properties)
#     K = create_sparsity_pattern(dh)
#     RHS = zeros(ndofs(dh))
#
#     EP_model = Model( model.grid,
#                       dh,
#                       dbc,
#                       model.neumann_bc,
#                       model.body_forces,
#                       (model.cellvalues_tuple[1],),
#                       model.facevalues,
#                       mp,
#                       model.material_state,
#                       K,
#                       RHS,
#                       model.clock,
#                       model.solver,
#                       model.multithreading )
#     end
#     println("EP_model building = ", tt)
#
#
#     dbc, K, res, solver, clock = EP_model.dirichlet_bc, EP_model.K, EP_model.RHS, EP_model.solver, EP_model.clock
#     nbasefuncs = getnbasefunctions(EP_model.cellvalues_tuple[1])
#     δu = similar(res)
#
#     update!(dbc, clock.current_time) #update newly formed dirichlet bc in EP_model
#     dofs_u = get_field_dofs(:u, model) # dofs associated with displacement field in u (u-D formulation)
#
#     tt = @elapsed doassemble!(EP_model, nbasefuncs, u[dofs_u] ; noplast = true)
#     println("assemble time : ", tt)
#     # compute residual norm
#     norm_res = norm(res[JuAFEM.free_dofs(dbc)])
#     print("First elastic iteration   \t residual: $(@sprintf("%.8f", norm_res))\n")
#
#     ### Linear Solve for δu ###
#     apply_zero!(K, res, dbc)
#     δu .= solver.linear_solver(K,-res,EP_model)
#     # displacement correction
#
#     @assert length(dofs_u) == length(δu)
#     u[dofs_u] .+= δu
#
#     # update initial model temporary state
#     rcell, rqp = rand(1:getncells(model.grid)), rand(1:getnquadpoints(model.cellvalues_tuple[1]))
#     copy_temp_state!(model,EP_model.material_state)
#     @assert EP_model.material_state[rcell][rqp].temp_σ == model.material_state[rcell][rqp].temp_σ
#
#     return nothing
# end

function iterate(model::Model{2,1,D,V,E,P}, output_writer, initial_values = nothing) where {D,V,E,P}

    # Unpack some model fields
    dh, dbc, cv, clock = model.dofhandler, model.dirichlet_bc, model.cellvalues_tuple, model.clock

    # Pre-allocate solution vectors, etc.
    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)  # solution vector
    u_converged = copy(u) # backup solution vector
    #u_grav_clean = copy(u)
    δu = zeros(n_dofs)  # displacement correction

    while clock.current_time <= clock.tspan[2]
        @timeit "time iteration" begin

            timestep!(clock) # update clock

            print("\n TIME ITERATION $(clock.iter)\n",
            " current simulation time = $(clock.current_time):\n",
            " timestep = $(clock.Δt)\n")

            restart_flag = false

            # Apply dirichlet bc and iteratively solve for u :
            # @timeit "nonlinear solve"
            restart_flag = nonlinear_solve!(u,u_converged,δu,model,restart_flag)

            if restart_flag == true
                u .= u_converged
                undo_timestep!(clock)
                clock.Δt *= clock.Δt_fact_down # decreased timestep
            else # converged

                #### TODO clear gravity displacement field nicer
                # if clock.current_time == 0
                #     u_grav_clean .= u
                # end
                update_material_state!(model) # update converged state values
                @timeit "export" write_output!(model, u, output_writer) # output

                clock.Δt *= clock.Δt_fact_up # increase timestep
                println(clock.Δt)
                u_converged .= u
            end

            clock.current_time == clock.tspan[2] && break # end time loop if requested end time is reached

        end
    end
    return u
end
# iterate(model::Model{2,1,Nothing,Nothing,E,P}, output_writer) where {E<:Elasticity,P<:Plasticity} = iterate(model::Model{2,1,Nothing,V,E,P}, output_writer) where {V,E,P}
# iterate(model::Model{2,1,D,Nothing,E,Nothing}, output_writer) where {D<:Damage,E<:Elasticity} = iterate(model::Model{2,1,Nothing,V,E,P}, output_writer) where {V,E,P}
# iterate(model::Model{2,1,D,Nothing,E,P}, output_writer) where {D<:Damage,E<:Elasticity,P<:Plasticity} = iterate(model::Model{2,1,Nothing,V,E,P}, output_writer) where {V,E,P}
# iterate(model::Model{2,1,Nothing,V,E,Nothing}, output_writer) where {V<:Viscosity,E<:Elasticity} = iterate(model::Model{2,1,Nothing,V,E,P}, output_writer) where {V,E,P}

function fill_state_from_u!(model,u)
    @inbounds for (i, cell) in enumerate(CellIterator(model.dofhandler))
        cell_states, r = model.material_state[i], model.material_properties[i]
        eldofs = celldofs(cell)
        nu = getnbasefunctions(model.cellvalues_tuple[1])
        ue = u[eldofs][1:nu]
        fill_cell_state!(model,cell,cell_states,r,ue)
    end
end

function fill_cell_state!(model,cell,cell_states,r,ue)
    cvu = model.cellvalues_tuple[1]
    n_basefuncs = getnbasefunctions(cvu)
    reinit!(cvu, cell)
    Dᵉ = r.elasticity.Dᵉ
    @inbounds for q_point in 1:getnquadpoints(cvu)
        state = cell_states[q_point]
        # For each integration point, compute stress and material stiffness
        ϵ2D = function_symmetric_gradient(cvu, q_point, ue)
        # Total strain recomputed each time because of the newton correction
        state.temp_ϵ = SymmetricTensor{2,3}((i,j)->get_3D_func(i,j,ϵ2D))
        state.temp_σ = Dᵉ ⊡ state.temp_ϵ
    end
end

export_vtk(model::Model{2,Nothing,V,E,P},u,filename) where {V,E,P} = nothing
function export_vtk(model::Model{2,Nothing,V,E,Nothing},u,filename) where {V,E}
    p_values = zeros(getncells(model.grid))
    σxx_values = zeros(getncells(model.grid))
    σyy_values = zeros(getncells(model.grid))
    σzz_values = zeros(getncells(model.grid))
    σxy_values = zeros(getncells(model.grid))
    ϵxx_values = zeros(getncells(model.grid))
    ϵyy_values = zeros(getncells(model.grid))
    ϵzz_values = zeros(getncells(model.grid))
    ϵxy_values = zeros(getncells(model.grid))
    ϵkk_values = zeros(getncells(model.grid))
    for (el, cell_states) in enumerate(model.material_state)
        for state in cell_states
            p_values[el] += - 1/3 * tr(state.σ)
            σxx_values[el] += state.σ[1,1]
            σyy_values[el] += state.σ[2,2]
            σzz_values[el] += state.σ[3,3]
            σxy_values[el] += state.σ[1,2]
            ϵxx_values[el] += state.ϵ[1,1]
            ϵyy_values[el] += state.ϵ[2,2]
            ϵzz_values[el] += state.ϵ[3,3]
            ϵxy_values[el] += state.ϵ[1,2]
            ϵkk_values[el] += tr(state.ϵ)
        end
        nqp = length(cell_states)
        p_values[el] /= nqp
        σxx_values[el] /= nqp
        σyy_values[el] /= nqp
        σzz_values[el] /= nqp
        σxy_values[el] /= nqp
        ϵxx_values[el] /= nqp
        ϵyy_values[el] /= nqp
        ϵzz_values[el] /= nqp
        ϵxy_values[el] /= nqp
        ϵkk_values[el] /= nqp
    end
    vtk_grid(filename, model.dofhandler) do vtkfile
        vtk_point_data(vtkfile, model.dofhandler, u)
        vtk_point_data(vtkfile, model.dirichlet_bc)
        vtk_cell_data(vtkfile, p_values, "pressure")
        vtk_cell_data(vtkfile, σxx_values, "sigma xx")
        vtk_cell_data(vtkfile, σyy_values, "sigma yy")
        vtk_cell_data(vtkfile, σzz_values, "sigma zz")
        vtk_cell_data(vtkfile, σxy_values, "sigma xy")
        vtk_cell_data(vtkfile, ϵxx_values, "epsilon xx")
        vtk_cell_data(vtkfile, ϵyy_values, "epsilon yy")
        vtk_cell_data(vtkfile, ϵzz_values, "epsilon zz")
        vtk_cell_data(vtkfile, ϵxy_values, "epsilon xy")
        vtk_cell_data(vtkfile, ϵkk_values, "volumetric strain")
    end
end


text:
"""
Records solution information from a mixed-integer linear programming (MILP)
problem.

    MILP(;
        termination_status::MOI.TerminationStatusCode,
        primal_status::MOI.ResultStatusCode,
        solve_time_sec::Float64;
        x::Union{Nothing,Vector{Float64}} = nothing,
        primal_bound::Float64 = NaN,
        dual_bound::Float64 = NaN,
        node_count::Int = -1,
        simplex_iters::Int = -1
    )

Records why the solver terminated (`termination_status`), if it had a feasible
primal solution (`primal_status`), what that solution was if available (`x`),
what the best primal cost was (`primal_bound`), the best available dual bound
(`dual_bound`), elapsed solve time in seconds (`solve_time_sec`), and
optionally the number of enumerated nodes (`node_count`) and total simplex
iterations (`simplex_iters`).

Some pieces of the [functional interface](@ref functional_interface) are
implemented:

    record_result(table::CSVRecord, result::MILPResult)

Additionally, there is a helper function for creating a `CSVRecord` with the
correct headers:

    CSVRecord(filename::String, ::Type{MILPResult})
"""
Base.@kwdef mutable struct MILPResult <: AbstractResult
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    solve_time_sec::Float64
    x::Union{Nothing,Vector{Float64}} = nothing
    primal_bound::Float64 = NaN
    dual_bound::Float64 = NaN
    node_count::Int = -1
    simplex_iters::Int = -1
end

# Drop :x as we don't want to record that in a CSVRecord
const _MILP_CSV_RESULT_FIELDS = [
    :termination_status,
    :primal_status,
    :primal_bound,
    :dual_bound,
    :solve_time_sec,
    :node_count,
    :simplex_iters,
]

function CSVRecord(filename::String, ::Type{MILPResult})
    # Don't write feasible solution to CSV
    return CSVRecord(
        filename,
        vcat(
            ["instance", "approach"],
            [string(field) for field in _MILP_CSV_RESULT_FIELDS],
        ),
    )
end

function record_result!(
    table::CSVRecord,
    result::MILPResult,
    instance_name::String,
    approach_name::String,
)
    println(
        table.fp,
        join(
            vcat(
                [instance_name, approach_name],
                [getfield(result, field) for field in _MILP_CSV_RESULT_FIELDS],
            ),
            ",",
        ),
    )
    return flush(table.fp)
end

ext: jl


text:
module RobertoMD

using Formatting
using Random, Distributions
using RoundingIntegers
using ProgressMeter
using LinearAlgebra
using Logging
using JSON
using MPI
using StaticArrays

include("core/bc/type_pbc.jl")
include("potential/bond/type_bond.jl")
include("thermostat/type_thermostat.jl")
include("core/velocity/type_velocity.jl")
include("core/integrator/type_integrator.jl")
include("core/dump/type_dump.jl")
include("core/log/type_log.jl")
include("core/configuration.jl")
include("core/restart/type_restart.jl")
include("core/MPI/balancingMPI.jl")
include("potential/field/type_field.jl")
include("core/MPI/initialize.jl")
include("core/atoms.jl")
include("core/system.jl")
include("util/array.jl")
include("util/random.jl")
include("core/force.jl")
include("core/dump/dump.jl")
include("core/restart/restart.jl")
include("potential/bond/force_bond.jl")
include("potential/field/force_field.jl")
include("core/integrator/integrate.jl")
include("thermostat/thermostat.jl")
include("core/thermo/thermo.jl")
include("core/log/log.jl")
include("core/simulation.jl")

end


text:
# template log-on and fetch graph

using Caesar
using IncrementalInference

include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
addrdict["session"] = "SESS??"
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

fg = Caesar.initfg(cloudgraph=cloudGraph, sessionname=addrdict["session"])
fullLocalGraphCopy!(fg)
savejld(fg, file="somefg.jld")

plotKDE(fg, :x1, dims=[1;2], api=localapi)
fncvar = getfnctype(getVert(fg,fg.fIDs[:x1710l200050]))
ret = whosNear2D(cloudGraph, "SESS??", x=21.7, y=-54.0)
plotLocalProduct(fg, :l200050, dims=[1;2], api=localapi)

#


text:
using PRM

simp = SimParam(Tm=10, epsilon=1e-6, x0=0, h0=0.2, t0=0, tf=20, ts=1);

run_sim(simp)

max_stars_repo_name: VeroAsghar/PRM

max_issues_repo_name: VeroAsghar/PRM

max_forks_repo_name: VeroAsghar/PRM


text:
module StringMatch

export lcsmatch

include("sais.jl")
include("suffixarray.jl")

isvalid(i::Int, j::Int, k::Int, l::Int) = i > 0 && i <= j && k > 0 && k <= l

function lcsmatch(seq1::Vector{Int}, seq2::Vector{Int})
    pairs = Tuple{Int,Int}[]
    buffer = [(1,length(seq1),1,length(seq2),true)]
    while !isempty(buffer)
        i,j,k,l,lmost = pop!(buffer)
        ii,jj,kk,ll = lcsubstring(seq1[i:j],seq2[k:l],lmost)
        isvalid(ii,jj,kk,ll) || continue
        ii += i - 1
        jj += i - 1
        kk += k - 1
        ll += k - 1
        append!(pairs, zip(ii:jj,kk:ll))

        left = i, ii-1, k, kk-1
        right = jj+1, j, ll+1, l
        isvalid(left...) && push!(buffer,(left...,false))
        isvalid(right...) && push!(buffer,(right...,true))
    end
    sort!(pairs)
end

function lcsubstring(seq1::Vector{Int}, seq2::Vector{Int}, leftmost::Bool)
    n1 = length(seq1)
    seq = copy(seq1)
    maxi = max(maximum(seq1), maximum(seq2))
    push!(seq, maxi+1)
    append!(seq, seq2)
    sa = SuffixArray(seq)
    lcparray = sa.lcparray
    maxlcp = 0
    for i = 2:length(lcparray)
        (sa[i-1] <= n1) == (sa[i] <= n1) && continue
        lcp = lcparray[i]
        lcp > maxlcp && (maxlcp = lcp)
    end
    maxlcp == 0 && return (0,0,0,0)

    bestp1, bestp2 = 0, 0
    for i = 1:length(lcparray)
        lcp = lcparray[i]
        lcp == maxlcp || continue
        i < length(lcparray) && lcparray[i+1] == maxlcp && continue

        s = findprev(x -> x != maxlcp, lcparray, i)
        posits = map(j -> sa[j], s:i)
        all(p -> (p <= n1) == (posits[1] <= n1), posits) && continue

        p1 = minimum(posits)
        p2 = minimum(filter(p -> p > n1, posits)) - (n1 + 1)
        if bestp1+bestp2 == 0
            (bestp1,bestp2) = (p1,p2)
        elseif leftmost && p1+p2 < bestp1+bestp2
            (bestp1,bestp2) = (p1,p2)
        elseif !leftmost && p1+p2 > bestp1+bestp2
            (bestp1,bestp2) = (p1,p2)
        end
    end
    (bestp1, bestp1+maxlcp-1, bestp2, bestp2+maxlcp-1)
end

end


text:
"""
Gère l'accès et la manipulation des données
"""

const DATAPATH = "../data";

# retourne un Array{Station,1} d'objets Station
function getStations(fileName::String)
    f = openFile(fileName) # ouverture du fichier
    nbStations = parse(Int,readline(f)) # nombre de stations
    stations = Vector{Station}(undef,nbStations)
    for indexStation in 1:nbStations
        line = readline(f)
        splitLine = split(line," ")
        stations[indexStation] = Station(indexStation, splitLine[1],parse(Float64,splitLine[3]),parse(Float64,splitLine[2])) # construction de la station
    end
    return stations
end

# retourne le fichier ouvert correspondant au nom donné
function openFile(fileName::String)
    actual_path = pwd()
    cd(joinpath(DATAPATH)) # déplacement dans le repertoire où sont stockés les fichiers source
    f = open(fileName)
    cd(actual_path)
    return f
end

# parse un fichier simulation et retourne la liste des requetes
function parseRequests(fileName::String)
    f = openFile(fileName)
    nbShuttles = parse(Int,readline(f))
    nbRequest = parse(Int,readline(f))
    requests = Vector{Request}(undef,nbRequest)
    for indexRequest in 1:nbRequest
        splittedLine = split(readline(f),' ')
        t = parse(Int,splittedLine[1])
        departureStation = string(splittedLine[2])
        arrivalStation = string(splittedLine[3])
        rq = Request(t,departureStation,arrivalStation)
        requests[indexRequest] = rq
    end
    return requests, nbShuttles
end


text:
module ExpressionUtils

using Base.Meta
using ..LintCompat

export split_comparison, simplify_literal, ispairexpr, isliteral,
       lexicaltypeof, lexicalfirst, lexicallast, lexicalvalue,
       withincurly, expand_trivial_calls, expand_assignment, COMPARISON_OPS,
       understand_import, dots, kind, path

"""
    withincurly(ex)

Get just the function part of a function declaration, or just the type head of
a parameterized type name.

```jldoctest
julia> using Lint.ExpressionUtils

julia> withincurly(:(Vector{T}))
:Vector

julia> withincurly(:((::Type{T}){T}))
:(::Type{T})
```
"""
withincurly(ex) = isexpr(ex, :curly) ? ex.args[1] : ex

"""
    split_comparison(::Expr)

Split a :comparison expression into an equivalent :&& expression, given that
each component of the comparison is pure.
"""
function split_comparison(ex)
    if !isexpr(ex, :comparison)
        throw(ArgumentError("expected a comparison expression, got a $(ex.head)"))
    end
    left = ex.args[1]
    op = ex.args[2]
    right = ex.args[3]
    remainder = ex.args[3:end]
    if length(remainder) == 1
        :($op($left, $right))
    else
        :($op($left, $right) &&
          $(split_comparison(Expr(:comparison, remainder...))))
    end
end

"""
    simplify_literal(::Expr)

Simplify certain macros into the literals they produce.

Simplifications performed include:

 - v"x.y.z" literals are simplified into `VersionNumber` objects.
"""
function simplify_literal(ex)
    if isexpr(ex, :macrocall) && ex.args[1] == Symbol("@v_str") &&
       isa(ex.args[2], AbstractString)
        VersionNumber(ex.args[2])
    else
        ex
    end
end

"""
Return `true` if `x` is an expression of the form `:(k => v)`.  Note that on
v0.5, this is `Expr(:(=>), ...)`, whereas on v0.6 it is `Expr(:call, :(=>),
...)`.
"""
ispairexpr(x) = isexpr(x, :(=>)) || isexpr(x, :call) && x.args[1] == :(=>)

"""
Return the first entry of the given pair expression, determined lexically.
Note that on v0.5, this is the first argument of the expression, whereas on
v0.6, this is the second argument of the expression.
"""
lexicalfirst(x) = VERSION < v"0.6.0-dev.2613" ? x.args[1] : x.args[2] # changed by julia PR #20327

"""
Return the last entry of the given pair expression, determined lexically.  Note
that on v0.5, this is the second argument of the expression, whereas on v0.6,
this is the third argument of the expression.
"""
lexicallast(x) = VERSION < v"0.6.0-dev.2613" ? x.args[2] : x.args[3] # changed by julia PR #20327

"""
Return `true` if the value represented by expression `x` is exactly `x` itself;
that is, `x` is not `Expr`, `QuoteNode`, or `Symbol`.
"""
isliteral(x) = !isa(x, Expr) && !isa(x, QuoteNode) && !isa(x, Symbol)

"""
    lexicalvalue(x) :: Nullable{Any}

If `x` is a literal, or a quoted literal, return that literal wrapped in a
`Nullable`. Otherwise, return `Nullable{Any}()`.
"""
function lexicalvalue(x)
    if isliteral(x)
        Nullable{Any}(x)
    elseif isexpr(x, :quote)
        if isexpr(x.args[1], :($))
            lexicalvalue(x.args[1].args[1])
        else
            Nullable{Any}(x.args[1])
        end
    elseif isa(x, QuoteNode)
        Nullable{Any}(x.value)
    else
        Nullable{Any}()
    end
end

"""
Return the most specific known lexical type of the given expression. Lexical
type is defined as

- If the expression is a literal, the type of that literal;
- If the expression is a `QuoteNode`, or an `Expr(:quote, ...)`, the type of
  whatever is quoted;
- Otherwise, `Any`.

That is, the maximal amount of information detectable from the lexical context
alone.
"""
lexicaltypeof(x) = get(BROADCAST(typeof, lexicalvalue(x)), Any)

"""
    expand_trivial_calls(x)

Expand the outer layer of trivial calls. Trivial calls are defined as
expression nodes that almost always lower to calls but are not represented as
such. The special case lowering of `A*B'` is neglected.
"""
function expand_trivial_calls(ex)
    if isexpr(ex, :(:))
        Expr(:call, :colon, ex.args...)
    elseif isexpr(ex, Symbol("'"))
        Expr(:call, :ctranspose, ex.args...)
    elseif isexpr(ex, Symbol(".'"))
        Expr(:call, :transpose, ex.args...)
    elseif isexpr(ex, :(=>))
        Expr(:call, :(=>), ex.args...)
    elseif isexpr(ex, :vect)
        Expr(:call, :(Base.vect), ex.args...)
    else
        ex
    end
end

# TODO: deal with dot-calls
"""
    expand_assignment(x)

Return a tuple `(LHS, RHS)` by expanding the expression as if it represents a
single assignment. For example, `x += y` is expanded to `x = x + y`, which is
returned as `(x, x + y)`. Return `Nullable()` if the argument could not be
interpreted as an assignment.
"""
function expand_assignment(expr::Expr)::Nullable
    op = expr.head
    if op in COMPARISON_OPS
        Nullable()
    elseif op == :(=)
        @assert length(expr.args) == 2
        Nullable((expr.args[1], expr.args[2]))
    else
        str = string(op)
        if str[end] == '='
            fop = Symbol(str[1:end-1])
            @assert length(expr.args) == 2
            Nullable((expr.args[1], Expr(:call, fop, expr.args[1],
                                         expr.args[2])))
        else
            Nullable()
        end
    end
end
expand_assignment(_) = Nullable()

const COMPARISON_OPS = [:(==), :(<), :(>), :(<=), :(>=), :(!=)]

immutable Import
    """
    The number of dots preceding the import. If `dots` is `0`, then this is a
    toplevel import (i.e., from Main, but additionally requiring the module if
    it is not already loaded). A value of `1` indicates importing something
    from this module (typically only useful for `using`). From there, any
    increase of `1` to `dots` will move up one level to the parent module,
    until it is `Main`, after which adding more dots does not change the
    result.
    """
    dots :: Int

    """
    After any dots, the series of modules that must be loaded to obtain the
    object to be imported. All but the last of these should represent a module.
    """
    path :: Vector{Symbol}

    """
    If this is `:import`, then the given symbol is simply imported into the
    current namespace. If it is `:importall`, then that is done, and if the
    loaded object is a module, all its exports are also imported. If it is
    `:using`, then the effect is like `:importall`, except the imported symbols
    are not available for method extension.
    """
    kind :: Symbol
end
dots(imp::Import) = imp.dots
path(imp::Import) = imp.path
kind(imp::Import) = imp.kind

"""
    understand_import(ex)::Nullable{Import}

Return an `Import` from extracting the important semantics from the given
expression `ex`, or `Nullable()` otherwise.
"""
function understand_import(ex)::Nullable{Import}
    if !isexpr(ex, [:using, :import, :importall])
        return Nullable()
    end

    kind = ex.head
    dots = 0
    for x in ex.args
        if x === :.
            dots += 1
        else
            break
        end
    end
    path = ex.args[dots+1:end]

    Import(dots, path, kind)
end

end


text:
using InteractiveDynamics, Observables, Test
using DynamicalSystems

ds = Systems.henonheiles()

potential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
energy(x,y,px,py) = 0.5(px^2 + py^2) + potential(x,y)
E = energy(get_state(ds)...)
function complete(y, py, x)
    V = potential(x, y)
    Ky = 0.5*(py^2)
    Ky + V ≥ E && error("Point has more energy!")
    px = sqrt(2(E - V - Ky))
    ic = [x, y, px, py]
    return ic
end

chaotic = get_state(ds)
stable = [0., 0.1, 0.5, 0.]
plane = (1, 0.0)
tf = 1000.0
idxs = SVector{2, Int}(2, 4)
u0 = get_state(ds); D = 4
Ttr = 200.0
rootkw = (xrtol = 1e-6, atol = 1e-6)

integ = integrator(ds, u0)
planecrossing = PlaneCrossing(plane, true)
f = (t) -> planecrossing(integ(t))
i = SVector{2, Int}(idxs)
data = poincaresos(integ, planecrossing, tf, Ttr, i, rootkw)
data0 = copy(data)

# Emulate clicking on the psos
click = Observable([0.0, 0.0])
Data = Observable(data)

on(click) do pos
    x, y = pos; z = plane[2]
    newstate = try
       complete(x, y, z)
    catch err
       @error "Could not get state, got error:" exception=err
       return
    end
    @assert length(newstate) == D

    reinit!(integ, newstate)

    Data[] = poincaresos(integ, planecrossing, tf, Ttr, i, rootkw)
    return data
end

# we click on the psos now
click[] = [0.2, 0.2]

@test Data[] != data0


text:
if nprocs() > 1
    rmprocs(workers())
end
sleep(1.0)

thisfolder = dirname(@__FILE__())
modulefolder = normpath(joinpath(thisfolder,".."))
sincosfolder = joinpath(modulefolder,"target","sincos")
springmassfolder = joinpath(modulefolder,"ode","springmass")
fitzhughnagumofolder = joinpath(modulefolder,"ode","fitzhughnagumo")
predatorpreyfolder = joinpath(modulefolder,"ode","predatorprey")

sincos1 = joinpath(sincosfolder,"scripts","SinCos1.jl")
sincos2 = joinpath(sincosfolder,"scripts","SinCos2.jl")
springmass1 = joinpath(springmassfolder,"scripts","SpringMass1.jl")
fitzhughnagumo1 = joinpath(fitzhughnagumofolder,"scripts","FitzHughNagumo1.jl")
predatorprey1 = joinpath(predatorpreyfolder,"scripts","PredatorPrey1.jl")

println()
println("============================================")
println("++++++++++++++++++++++++++++++++++++++++++++")
println("GMHExamples.jl tests running in main process")
println("++++++++++++++++++++++++++++++++++++++++++++")
println("============================================")
println()
include(sincos1)
include(sincos2)
include(springmass1)
include(fitzhughnagumo1)
include(predatorprey1)

addprocs(1)

println()
println("===========================================")
println("+++++++++++++++++++++++++++++++++++++++++++")
println("GMHExamples.jl tests running in one process")
println("+++++++++++++++++++++++++++++++++++++++++++")
println("===========================================")
println()
include(sincos1)
include(sincos2)
include(springmass1)
include(fitzhughnagumo1)
include(predatorprey1)

addprocs(1)

println()
println("=============================================")
println("+++++++++++++++++++++++++++++++++++++++++++++")
println("GMHExamples.jl tests running in two processes")
println("+++++++++++++++++++++++++++++++++++++++++++++")
println("=============================================")
println()
include(sincos1)
include(sincos2)
include(springmass1)
include(fitzhughnagumo1)
include(predatorprey1)

if nprocs() > 1
    rmprocs(workers())
end
sleep(1.0)


text:
immutable Beta <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function Beta(a::Real, b::Real)
        (a > zero(a) && b > zero(b)) || error("alpha and beta must be positive")
        new(float64(a), float64(b))
    end
end

Beta(a::Real) = Beta(a, a) # symmetric in [0, 1]
Beta() = Beta(1.0) # uniform

@_jl_dist_2p Beta beta

function entropy(d::Beta)
    o = lbeta(d.alpha, d.beta)
    o -= (d.alpha - 1.0) * digamma(d.alpha)
    o -= (d.beta - 1.0) * digamma(d.beta)
    o += (d.alpha + d.beta - 2.0) * digamma(d.alpha + d.beta)
    o
end

function kurtosis(d::Beta)
    α, β = d.alpha, d.beta
    num = 6.0 * ((α - β)^2 * (α + β + 1.0) - α * β * (α + β + 2.0))
    den = α * β * (α + β + 2.0) * (α + β + 3.0)
    num / den
end

mean(d::Beta) = d.alpha / (d.alpha + d.beta)

median(d::Beta) = quantile(d, 0.5)

function mode(d::Beta)
    α, β = d.alpha, d.beta
    if α >= 1.0
        if β > 1.0
            (α - 1.0) / (α + β - 2.0)
        elseif α == 1.0 && β == 1.0
            # Uniform[0,1]: what should be returned?
            0.5
        else
            1.0
        end
    else
        if β >= 1.0
            0.0
        else
            # not unique: return largest
            α > β ? 1.0 : 0.0
        end
    end
end

function modes(d::Beta)
    α, β = d.alpha, d.beta
    if α >= 1.0
        if β > 1.0
            [(α - 1.0) / (α + β - 2.0)]
        elseif α == 1.0 && β == 1.0
            # Uniform[0,1]: what should be returned?
            Float64[]
        else
            [1.0]
        end
    else
        if β >= 1.0
            [0.0]
        else
            [0.0,1.0]
        end
    end
end

function rand(d::Beta)
    u = randg(d.alpha)
    u / (u + randg(d.beta))
end

function rand!(d::Beta, A::Array{Float64})
    sa = GammaSampler(d.alpha)
    sb = GammaSampler(d.beta)

    for i = 1:length(A)
        u = rand(sa)
        v = rand(sb)
        @inbounds A[i] = u / (u + v)
    end

    return A
end

function skewness(d::Beta)
    num = 2.0 * (d.beta - d.alpha) * sqrt(d.alpha + d.beta + 1.0)
    den = (d.alpha + d.beta + 2.0) * sqrt(d.alpha * d.beta)
    num / den
end

function var(d::Beta)
    ab = d.alpha + d.beta
    d.alpha * d.beta / (ab * ab * (ab + 1.0))
end

function gradloglik(d::Beta, x::Float64)
  insupport(Beta, x) ? (d.alpha - 1.0) / x - (d.beta - 1.0) / (1 - x) : 0.0
end

### handling support

@continuous_distr_support Beta 0.0 1.0

## Fit model

# TODO: add MLE method (should be similar to Dirichlet)

# This is a moment-matching method (not MLE)
#
function fit(::Type{Beta}, x::Array)
    for xi in x
        insupport(Beta, xi) || error("Beta observations must be in [0,1]")
    end
    x_bar = mean(x)
    v_bar = varm(x, x_bar)
    α = x_bar * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    β = (1.0 - x_bar) * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    Beta(α, β)
end


text:
export PointLight, DistantLight

# ----- #
# Light #
# ----- #

"""
    Light

All objects emitting light should be a subtype of this. To add a custom light
source, two only two functions need to be defined.

* `get_direction(l::MyCustomLightSource, pt::Vec3)` - `pt` is the point receiving the
                                                      light from this source
* `get_intensity(l::MyCustomLightSource, pt::Vec3, dist::Array)` - `dist` is the array
                                                                   representing the distance
                                                                   of the point from the light
                                                                   source
"""
abstract type Light end

"""
    get_shading_info(l::Light, pt::Vec3)

Returns the `direction` of light incident from the light source onto that
Point and the intensity of that light ray. It computes these by internally
calling the `get_direction` and `get_intensity` functions and hence this
function never has to be modified for any `Light` subtype.

### Arguments:

* `l`  - The Object of the [`Light`](@ref) subtype
* `pt` - The Point in 3D Space for which the shading information is queried
"""
function get_shading_info(l::Light, pt::Vec3)
    dir = get_direction(l, pt)
    dist = sqrt.(l2norm(dir))
    intensity = get_intensity(l, pt, dist)
    return normalize(dir), intensity
end

"""
    get_direction(l::Light, pt::Vec3)

Returns the direction of light incident from the light source to `pt`.
"""
function get_direction(l::Light, pt::Vec3) end

"""
    get_intensity(l::Light, pt::Vec3, dist::Array)

Computed the intensity of the light ray incident at `pt`.
"""
function get_intensity(l::Light, pt::Vec3, dist::AbstractArray) end

# --------------- #
# - Point Light - #
# --------------- #

"""
    PointLight

A source of light which emits light rays from a point in 3D space. The
intensity for this light source diminishes as inverse square of the distance
from the point.

### Fields:

* `color`     - Color of the Light Rays emitted from the source.
* `intensity` - Intensity of the Light Source. This decreases as ``\\frac{1}{r^2}``
                where ``r`` is the distance of the point from the light source.
* `position`  - Location of the light source in 3D world space.

### Available Constructors:

* `PointLight(;color = Vec3(1.0f0), intensity::Real = 100.0f0, position = Vec3(0.0f0))`
* `PointLight(color, intensity::Real, position)`
"""
struct PointLight{T<:AbstractArray} <: Light
    color::Vec3{T}
    intensity::T
    position::Vec3{T}
end

PointLight(;color = Vec3(1.0f0), intensity::Real = 100.0f0, position = Vec3(0.0f0)) =
    PointLight(color, intensity, position)

PointLight(c, i::Real, p) = PointLight(clamp(c, 0.0f0, 1.0f0), [i], p)

show(io::IO, pl::PointLight) =
    print(io, "Point Light\n    Color - ", pl.color, "\n    Intensity - ",
          pl.intensity[], "\n    Position - ", pl.position)

@diffops PointLight

get_direction(p::PointLight, pt::Vec3) = p.position - pt

get_intensity(p::PointLight, pt::Vec3, dist::AbstractArray) =
    p.intensity[] * p.color / (4 .* (eltype(p.color.x))(π) .* (dist .^ 2))

# ----------------- #
# - Distant Light - #
# ----------------- #

"""
    DistantLight

A source of light that is present at infinity as such the rays from it are
in a constant direction. The intensity for this light source always remains
constant. This is extremely useful when we are trying to render a large scene
like a city.

### Fields:

* `color`     - Color of the Light Rays emitted from the Source.
* `intensity` - Intensity of the Light Source.
* `direction` - Direction of the Light Rays emitted by the Source.

### Available Constructors:

* `DistantLight(;color = Vec3(1.0f0), intensity::Real = 1.0f0,
                 direction = Vec3(0.0f0, 1.0f0, 0.0f0))`
* `DistantLight(color, intensity::Real, direction)`
"""
struct DistantLight{T<:AbstractArray} <: Light
    color::Vec3{T}
    intensity::T
    direction::Vec3{T}  # Must be normalized
end

DistantLight(;color = Vec3(1.0f0), intensity::Real = 1.0f0,
              direction = Vec3(0.0f0, 1.0f0, 0.0f0)) =
    DistantLight(color, intensity, direction)

DistantLight(c, i::Real, d) = DistantLight(clamp(c, 0.0f0, 1.0f0), [i], normalize(d))

show(io::IO, dl::DistantLight) =
    print(io, "Distant Light\n    Color - ", dl.color, "\n    Intensity - ",
          dl.intensity[], "\n    Direction - ", dl.direction)

@diffops DistantLight

get_direction(d::DistantLight, pt::Vec3) = d.direction

get_intensity(d::DistantLight, pt::Vec3, dist::AbstractArray) = d.intensity[]


text:
function awss3_tests(config)
    bucket_name =
        "ocaws.jl.test." * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ"))

    @testset "Create Bucket" begin
        s3_create_bucket(config, bucket_name)
        @test bucket_name in s3_list_buckets(config)
        is_aws(config) && s3_enable_versioning(config, bucket_name)
        sleep(1)
    end

    @testset "Bucket Tagging" begin
        @test isempty(s3_get_tags(config, bucket_name))
        tags = Dict("A" => "1", "B" => "2", "C" => "3")
        s3_put_tags(config, bucket_name, tags)
        @test s3_get_tags(config, bucket_name) == tags
        s3_delete_tags(config, bucket_name)
        @test isempty(s3_get_tags(config, bucket_name))
    end

    @testset "Create Objects" begin
        s3_put(config, bucket_name, "key1", "data1.v1")
        s3_put(bucket_name, "key2", "data2.v1"; tags=Dict("Key" => "Value"))
        s3_put(config, bucket_name, "key3", "data3.v1")
        s3_put(config, bucket_name, "key3", "data3.v2")
        s3_put(config, bucket_name, "key3", "data3.v3"; metadata=Dict("foo" => "bar"))
        s3_put(config, bucket_name, "key4", "data3.v4"; acl="bucket-owner-full-control")
        s3_put_tags(config, bucket_name, "key3", Dict("Left" => "Right"))

        @test isempty(s3_get_tags(config, bucket_name, "key1"))
        @test s3_get_tags(config, bucket_name, "key2")["Key"] == "Value"
        @test s3_get_tags(config, bucket_name, "key3")["Left"] == "Right"
        s3_delete_tags(config, bucket_name, "key2")
        @test isempty(s3_get_tags(config, bucket_name, "key2"))

        @test s3_get(config, bucket_name, "key1") == b"data1.v1"
        @test s3_get(config, bucket_name, "key2") == b"data2.v1"
        @test s3_get(bucket_name, "key3") == b"data3.v3"
        @test s3_get(bucket_name, "key4") == b"data3.v4"
        @test s3_get_meta(bucket_name, "key3")["x-amz-meta-foo"] == "bar"
    end

    @testset "ASync Get" begin
        @sync begin
            for i in 1:2
                @async begin
                    @test s3_get(bucket_name, "key3") == b"data3.v3"
                end
            end
        end
    end

    @testset "Raw Return - XML" begin
        xml = "<?xml version='1.0'?><Doc><Text>Hello</Text></Doc>"
        s3_put(config, bucket_name, "file.xml", xml, "text/xml")
        @test String(s3_get(config, bucket_name, "file.xml"; raw=true)) == xml
        @test s3_get(config, bucket_name, "file.xml")["Text"] == "Hello"
    end

    @testset "Get byte range" begin
        teststr = "123456789"
        s3_put(config, bucket_name, "byte_range", teststr)
        range = 3:6
        @test String(s3_get(config, bucket_name, "byte_range"; byte_range=range)) ==
            teststr[range]
    end

    @testset "Object Copy" begin
        s3_copy(bucket_name, "key1"; to_bucket=bucket_name, to_path="key1.copy")
        @test s3_get(config, bucket_name, "key1.copy") == b"data1.v1"
    end

    @testset "Object exists" begin
        for key in ["key1", "key2", "key3", "key1.copy"]
            @test s3_exists(bucket_name, key)
        end
    end

    @testset "List Objects" begin
        for key in ["key1", "key2", "key3", "key1.copy"]
            @test key in [o["Key"] for o in s3_list_objects(config, bucket_name)]
        end
    end

    @testset "Object Delete" begin
        s3_delete(config, bucket_name, "key1.copy")
        @test !("key1.copy" in [o["Key"] for o in s3_list_objects(config, bucket_name)])
    end

    @testset "Check Metadata" begin
        meta = s3_get_meta(config, bucket_name, "key1")
        @test meta["ETag"] == "\"68bc8898af64159b72f349b391a7ae35\""
    end

    @testset "default Content-Type" begin
        # https://github.com/samoconnor/AWSS3.jl/issues/24
        ctype(key) = s3_get_meta(bucket_name, key)["Content-Type"]

        for k in ["file.foo", "file", "file_html", "file.d/html", "foobar.html/file.htm"]
            is_aws(config) && k == "file" && continue
            s3_put(config, bucket_name, k, "x")
            @test ctype(k) == "application/octet-stream"
        end

        for (k, t) in [
            ("foo/bar/file.html", "text/html"),
            ("x.y.z.js", "application/javascript"),
            ("downalods/foo.pdf", "application/pdf"),
            ("data/foo.csv", "text/csv"),
            ("this.is.a.file.txt", "text/plain"),
            ("my.log", "text/plain"),
            ("big.dat", "application/octet-stream"),
            ("some.tar.gz", "application/octet-stream"),
            ("data.bz2", "application/octet-stream"),
        ]
            s3_put(config, bucket_name, k, "x")
            @test ctype(k) == t
        end
    end

    @testset "Multi-Part Upload" begin
        MIN_S3_CHUNK_SIZE = 5 * 1024 * 1024 # 5 MB
        key_name = "multi-part-key"
        upload = s3_begin_multipart_upload(config, bucket_name, key_name)
        tags = Vector{String}()

        for part_number in 1:5
            push!(
                tags,
                s3_upload_part(config, upload, part_number, rand(UInt8, MIN_S3_CHUNK_SIZE)),
            )
        end

        s3_complete_multipart_upload(config, upload, tags)
        @test s3_exists(bucket_name, key_name)
    end

    # these tests are needed because lack of functionality of the underlying AWS API makes certain
    # seemingly inane tasks incredibly tricky: for example checking if an "object" (file or
    # directory) exists is very subtle
    @testset "path naming edge cases" begin
        # this seemingly arbitrary operation is needed because of the insanely tricky way we
        # need to check for directories
        s3_put(bucket_name, "testdir.", "") # create an empty file called `testdir.`
        s3_put(bucket_name, "testdir/", "") # create an empty file called `testdir/` which AWS will treat as an "empty directory"
        @test s3_exists(bucket_name, "testdir/")
        @test isdir(S3Path(bucket_name, "testdir/"))
        @test !isfile(S3Path(bucket_name, "testdir/"))
        @test s3_exists(bucket_name, "testdir.")
        @test isfile(S3Path(bucket_name, "testdir."))
        @test !isdir(S3Path(bucket_name, "testdir."))
        @test !s3_exists(bucket_name, "testdir")

        s3_put(bucket_name, "testdir/testfile.txt", "what up")
        @test s3_exists(bucket_name, "testdir/testfile.txt")
        @test isfile(S3Path(bucket_name, "testdir/testfile.txt"))
        # make sure the directory still "exists" even though there's a key in there now
        @test s3_exists(bucket_name, "testdir/")
        @test isdir(S3Path(bucket_name, "testdir/"))
        @test !isfile(S3Path(bucket_name, "testdir/"))

        # but it is still a directory and not an object
        @test !s3_exists(bucket_name, "testdir")
    end

    @testset "Version is empty" begin
        # Create the file to ensure we're only testing `version`
        k = "version_empty.txt"
        s3_put(config, bucket_name, k, "v1")
        s3_put(config, bucket_name, k, "v2")

        if is_aws(config)
            @test_throws AWSException s3_get(bucket_name, k; version="")
            @test_throws AWSException s3_get_meta(bucket_name, k; version="")
            @test_throws AWSException s3_exists(bucket_name, k; version="")
            @test_throws AWSException s3_delete(bucket_name, k; version="")
        else
            # Using an empty string as the version returns the latest version
            @test s3_get(bucket_name, k; version="") == "v2"
            @test s3_get_meta(bucket_name, k; version="") isa AbstractDict
            @test s3_exists(bucket_name, k; version="")
            @test s3_delete(bucket_name, k; version="") == UInt8[]
        end
    end

    @testset "Version is nothing" begin
        # Create the file to ensure we're only testing `version`
        k = "version_nothing.txt"
        s3_put(config, bucket_name, k, "v1")
        s3_put(config, bucket_name, k, "v2")

        # Using an empty string as the version returns the latest version
        @test s3_get(bucket_name, k; version=nothing) == "v2"
        @test s3_get_meta(bucket_name, k; version=nothing) isa AbstractDict
        @test s3_exists(bucket_name, k; version=nothing)
        @test s3_delete(bucket_name, k; version=nothing) == UInt8[]
    end

    is_aws(config) && @testset "Sign URL" begin
        for v in ["v2", "v4"]
            url = s3_sign_url(config, bucket_name, "key1"; signature_version=v)
            curl_output = ""

            @repeat 3 try
                curl_output = read(`curl -s -o - $url`, String)
            catch e
                @delay_retry if true
                end
            end

            @test curl_output == "data1.v1"

            fn = "/tmp/jl_qws_test_key1"
            if isfile(fn)
                rm(fn)
            end

            @repeat 3 try
                s3_get_file(config, bucket_name, "key1", fn)
            catch e
                sleep(1)
                @retry if true
                end
            end

            @test read(fn, String) == "data1.v1"
            rm(fn)
        end
    end

    is_aws(config) && @testset "Check Object Versions" begin
        versions = s3_list_versions(config, bucket_name, "key3")
        @test length(versions) == 3
        @test (
            s3_get(config, bucket_name, "key3"; version=versions[3]["VersionId"]) ==
            b"data3.v1"
        )
        @test (
            s3_get(config, bucket_name, "key3"; version=versions[2]["VersionId"]) ==
            b"data3.v2"
        )
        @test (
            s3_get(config, bucket_name, "key3"; version=versions[1]["VersionId"]) ==
            b"data3.v3"
        )

        tmp_file = joinpath(tempdir(), "jl_qws_test_key3")
        s3_get_file(config, bucket_name, "key3", tmp_file; version=versions[2]["VersionId"])
        @test read(tmp_file) == b"data3.v2"
    end

    is_aws(config) && @testset "Purge Versions" begin
        s3_purge_versions(config, bucket_name, "key3")
        versions = s3_list_versions(config, bucket_name, "key3")
        @test length(versions) == 1
        @test s3_get(config, bucket_name, "key3") == b"data3.v3"
    end

    if is_aws(config)
        @testset "Empty and Delete Bucket" begin
            AWSS3.s3_nuke_bucket(config, bucket_name)
            @test !in(bucket_name, s3_list_buckets(config))
        end

        @testset "Delete Non-Existant Bucket" begin
            @test_throws AWS.AWSException s3_delete_bucket(config, bucket_name)
        end
    end
end


text:
using Equations
consts=@equs(c=3e8, G=6.67e-11, massSun=2e30)
constsu=@equs(c=U(3e8,:m/:s), G=U(6.67e-11,:N*:m^2/:kg^2), massSun=U(2e30,:kg))

F1=-:G*:M*:m/:r^2
F2=:m*:v^2/:r
orbit=Equation(F2,F1)
v=sqrt(orbit/:m*:r) #orbital velocity

a=-:G*:M/:R^2

#Black holes
E=@equ m*v^2/2=G*M*m/R
R=E*:R/:m/:v^2*2
R_BH=R&@equs(v=c, M=massSun)&consts

R=U(6.67e-11,:m^3/:s^2/:kg)*U(2e30,:kg)/U(3e8,:m/:s)/U(3e8,:m/:s)*2
println(R)


text:
#=
push!(LOAD_PATH, "<path to directory WarpFlow>") # <--- use this if Julia can't find WarpFlow
include("scripts/example2.jl")
=#
using WarpFlow
using PyPlot

# set the data and the target
dim = 2
tmpx = shuffle([rand(50); randn(80)/10 .+ .8])
tmpy = tmpx + rand(size(tmpx)) .* 1.1
nPhi = length(tmpx)
X = Array{Float64,1}[[(tmpx[i]- minimum(tmpx))/maximum(tmpx), (tmpy[i]-minimum(tmpy))/maximum(tmpy)] for i=1:nPhi]
# target(x) = targetUnif2d(x; targSig = 0.1, limit = 0.4, center = 0.5)
target(x) = WarpFlow.targetNormal2d(x; targSig = 0.5, center = [0.5,0.5])

# generate Flow object
kappa     = Array{Float64,1}[X[Int(i)]  for i in 1:round(nPhi/2)]
# x_grd_kns, y_grd_kns =  meshgrid(linspace(-1.0,1.0, 10),linspace(-1.0,1.0, 10))
# append!(kappa, Array{Float64,1}[ [x_grd_kns[i], y_grd_kns[i]] for i in 1:length(x_grd_kns)])
nKap   = length(kappa)

y0 = Flow(kappa, WarpFlow.array1(dim, nKap), X, WarpFlow.array2eye(dim, nPhi), dim)

# function to save images
function plotrun(;fignum = 1, save = false)
	x_grd, y_grd =  meshgrid(linspace(-0.1, 1.1, 175),linspace(-0.1, 1.1, 175))
	phix_grd_0   = Array{Float64,1}[[x_grd[i], y_grd[i]] for i=1:length(x_grd)]
	Dphix_grd_0  = Array{Float64,2}[eye(2)               for i=1:length(x_grd)]
	yplt0        = WarpFlow.Flow(y0.kappa, y0.eta_coeff, phix_grd_0, Dphix_grd_0, dim)
	dydt(t,y)    = WarpFlow.d_forward_dt(y, sigma)
	(t1,yplt1)   = WarpFlow.ode23_abstract_end(dydt,[0,1], yplt0) # Flow y0 forward to time 1
	det_grd  = Float64[abs(det(yplt1.Dphix[i])) for i=1:length(x_grd)]
	den,     = target(yplt1.phix)
	est_den  = det_grd.*den

	fig = figure()
	scatter(Float64[point[1] for point in X], Float64[point[2] for point in X], c="b")
	contour(x_grd, y_grd, reshape(est_den,size(x_grd)), 30 )
	if save
		plt[:savefig]("scripts/example2_v$fignum.pdf",dpi=180)
		plt[:close](fig)
	end
end

#  gradient ascent on kappa and eta_coeff
lambda, sigma = 1.0, 0.1
for counter = 1:50
	@time z0 = WarpFlow.get_grad(y0, target, lambda, sigma)
	y0 = y0 + 0.005 * z0
end

plotrun()

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
# ---
# title: 543. Diameter of Binary Tree
# id: problem543
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/diameter-of-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, you need to compute the length of the diameter of the
# tree. The diameter of a binary tree is the length of the **longest** path
# between any two nodes in a tree. This path may or may not pass through the
# root.
#
# **Example:**
# Given a binary tree
#
#
#
#               1
#              / \
#             2   3
#            / \
#           4   5
#
#
# Return **3** , which is the length of the path [4,2,1,3] or [5,2,1,3].
#
# **Note:** The length of path between two nodes is represented by the number of
# edges between them.
#
#
## @lc code=start
using LeetCode

function diameter_of_binary_tree(root::TreeNode)::Int
    return max(
        diameter_of_binary_tree(root.left),
        diameter_of_binary_tree(root.right),
        _max_depth_to_leaf(root.left) + _max_depth_to_leaf(root.right) + 1,
    ) - 1
end

diameter_of_binary_tree(root::Nothing) = 0

function _max_depth_to_leaf(root::TreeNode)
    return max(1 + _max_depth_to_leaf(root.left), 1 + _max_depth_to_leaf(root.right))
end

_max_depth_to_leaf(::Nothing) = 0
## @lc code=end


text:
using StaticArrays

abstract type AbstractVocabulary{T} end

Base.eltype(::AbstractVocabulary{T}) where T = T

struct Vocab{T, A<:AbstractVector{T}} <: AbstractVocabulary{T}
    list::A
    unk::T
    unki::Int
end

"""
    Vocab(data::Vector{<:AbstractString}, unk::AbstractString="[UNK]")

Constructor for `Vocab`. `data` is the list of vocabulary word, can be nonunique.
 The actual list will be the unique version of `data` (i.e. `vocab.list = unique(data)`).
 `unk` is the indicator word for all unknown words. `unk` can be either in or not in `data`,
 depends on the use case.
"""
Vocab(data::AbstractVector, unk::AbstractString="[UNK]") = Vocab{String}(data, unk)

"""
    Vocab{T}(data::AbstractVector, unk) where T

construct Vocab with element type `T`. `unk` must be specified.
"""
function Vocab{T}(data::AbstractVector, unk) where T
    udata = Vector{T}(undef, length(data))
    unk = T(unk)
    unique!(map!(T, udata, data))
    list = SizedVector{length(udata)}(udata)
    i = findfirst(==(unk), list)
    unki = isnothing(i) ? 0 : i
    return Vocab(list, unk, unki)
end

Base.length(v::Vocab) = length(v.list)

function Base.show(io::IO, v::Vocab)
    print(io, "Vocab{", eltype(v), ", ", nameof(typeof(v.list)), '}')
    print(io, "(size = ", length(v))
    print(io, ", unk = ", v.unk)
    print(io, ", unki = ", v.unki, ')')
end

_lookup_index(list, unki, s) = (i = findfirst(==(s), list); isnothing(i) ? unki : i)
_lookup_word(list, unk, i) = 0 < i <= length(list) ? list[i] : unk

lookup(v::Vocab) = Base.Fix1(lookup,  v)
lookup(::Type{T}, v::Vocab) where T = lookup $ T $ v
lookup(::Type{I}, v::Vocab{T}, s::T) where {T, I<:Integer} = I(_lookup_index(v.list, v.unki, s))
lookup(::Type{I}, v::Vocab{<:AbstractString}, s::AbstractString) where I<:Integer = I(_lookup_index(v.list, v.unki, s))
lookup(::Type{T}, v::Vocab{T}, i::Integer) where T = _lookup_word(v.list, v.unk, i)
lookup(::Type{<:Integer}, v::Vocab{T}, i::Integer) where T = throw(DomainError(i, "Cannot lookup the value $i in the vocabulary: value should have the same type as Vocab's element type ($(eltype(v)))"))

lookup(v::Vocab{<:AbstractString}, s::AbstractString) = lookup(Int, v, s)
lookup(v::Vocab{T}, s::T) where T = lookup(Int, v, s)

lookup(v::Vocab, i::Integer) = lookup(eltype(v), v, i)
lookup(v::Vocab, i, j, k...) = (lookup(v, i), lookup(v, j), map(lookup(v), k)...)
lookup(v::Vocab, is::Union{AbstractArray, Tuple, NamedTuple}) = map(lookup(v), is)

lookup(::Type{T}, v::Vocab, i, j, k...) where T = lookup(T, v, (i, j, k...))
lookup(::Type{T}, v::Vocab, is::Union{AbstractArray, Tuple, NamedTuple}) where T = map(lookup(T, v), is)

# handle Vocab{Int}
# by default calling lookup(v, i) only do word lookup (i.e. v.list[i]), use lookup(Int, v, i) for index lookup
lookup(v::Vocab{T}, i::Integer) where T <: Integer = _lookup_word(v.list, v.unk, i)
lookup(::Type{I}, v::Vocab{<:Integer}, s::Integer) where I<:Integer = I(_lookup_index(v.list, v.unki, s))

# lookup(::Type{OneHot}, v::Vocab) = lookup $ OneHot $ v
lookup(::Type{OneHot}, v::Vocab, i) = lookup_onehot(v, i)
lookup(T::Type{OneHot}, v::Vocab, is::AbstractArray) = OneHotArray{length(v)}(lookup_onehot(v, is))
lookup(::Type{OneHot}, v::Vocab, is::Union{Tuple, NamedTuple}) = map(lookup(OneHot, v), is)

lookup_onehot(v::Vocab, i) = OneHot(length(v))(lookup(UInt32, v, i))
lookup_onehot(v::Vocab, is::AbstractArray) = map(lookup_onehot $ v, is)

function lookup(T::Type{OneHot}, v::Vocab, i, j, k...)
    c = 0
    si = lookup(T, v, i)
    c += si isa OneHot ? 1 : length(parent(si))
    sj = lookup(T, v, j)
    c += sj isa OneHot ? 1 : length(parent(sj))
    sk = map(lookup(T, v), k)
    c += sum(k->k isa OneHot ? 1 : length(parent(k)), sk; init=0)
    arr = Vector{OneHot(length(v))}(undef, c); empty!(arr)
    si isa OneHot ? push!(arr, si) : append!(arr, parent(si))
    sj isa OneHot ? push!(arr, sj) : append!(arr, parent(sj))
    for k in sk
        k isa OneHot ? push!(arr, k) : append!(arr, parent(k))
    end
    return OneHotArray(arr)
end

lookup(v::Vocab, i::OneHot) = lookup(v, Int(i))
lookup(v::Vocab, i::OneHotArray) = lookup(v, reinterpret(UInt32, i))
lookup(::Type{T}, v::Vocab{T}, i::OneHot) where T = lookup(v, i)
lookup(::Type{T}, v::Vocab{T}, i::OneHotArray) where T = lookup(v, i)

"""
    lookup(v::Vocab, x)

Lookup `x` in `v`. `lookup` words depends on the type of `x`. If `x` is an integer,
 return the `x`-th word on the vocabulary list (i.e. `v.list[x]`) and return the unknown word
 if `x` is out-of-bound (`v.unk`). If `x` is a string, return the indice of `x` in the vocabulary
 list (i.e `findfirst(==(x), v.list`) and return the unknown indice if `x` not found in the list.
 If the unknown word `v.unk` is in the list, the unknown indice is its indice, otherwise 0.

This function is bidirectional except for `Vocab{<:Integer}`. For integer vocabulary, this function
 only get the `x`-th word (`v.list[x]`). Use `lookup(Int, v, x)` for explicit indice lookup.

# Example
```julia
julia> vocab = Vocab(["a", "b", "c", "a", "b", "c"])
Vocab{String, StaticArrays.SizedVector{3, String, Vector{String}}}(size = 3, unk = [UNK], unki = 0)

julia> vocab_unk = Vocab(["a", "b", "xxx"], "xxx")
Vocab{String, StaticArrays.SizedVector{3, String, Vector{String}}}(size = 3, unk = xxx, unki = 3)

julia> lookup(vocab, "b")
2

julia> lookup(vocab, "d")
0

julia> lookup(vocab_unk, "d")
3

julia> lookup(vocab, 1)
"a"

julia> lookup(vocab, 10000)
"[UNK]"

julia> lookup(vocab_unk, 10000)
"xxx"

```
"""
function lookup end

@eval $((@macroexpand @doc """
    lookup(Int, v::Vocab, x)

The explicit version of `lookup(v, x)`. Lookup the indice of `x` in the vocabulary
 list. `x` should have the same type as Vocab's element type.

# Example
```julia
julia> vocab_unk = Vocab(["a", "b", "xxx"], "xxx")
Vocab{String, StaticArrays.SizedVector{3, String, Vector{String}}}(size = 3, unk = xxx, unki = 3)

julia> lookup(Int, vocab_unk, "b")
2

```
"""
function lookup(Int, v::Vocab, x) end
).args[2])

@eval $((@macroexpand @doc """
    lookup(::Type{T}, v::Vocab{T}, i::Integer) where T

The explicit version of `lookup(v, i)`. Lookup the word at index `i` on vocabulary
 list. `T` should be the same type as Vocab's element type. This method won't
 work on integer vocab, use `lookup(v, i)` directly.

# Example
```julia
julia> vocab_unk = Vocab(["a", "b", "xxx"], "xxx")
Vocab{String, StaticArrays.SizedVector{3, String, Vector{String}}}(size = 3, unk = xxx, unki = 3)

julia> lookup(String, vocab_unk, 1)
"a"

```
"""
lookup(::Type{T}, v::Vocab{T}, i::Integer) where T
).args[2])

@eval $((@macroexpand @doc """
    lookup(v::Vocab, is::AbstractArray)

recursively lookup value from `is`

# Example
```julia
julia> lookup(vocab, ["b", "c", "a", "A", "[UNK]"])
5-element Vector{Int64}:
 2
 3
 1
 0
 0

julia> lookup(vocab, [1, "a", 0, "A", "[UNK]"])
5-element Vector{Any}:
  "a"
 1
  "[UNK]"
 0
 0

```
"""
function lookup(v::Vocab, is::AbstractArray) end
).args[2])

@eval $((@macroexpand @doc """
    lookup(OneHot, v::Vocab, i)

lookup `i` and convert into one-hot representation.

# Example
```julia
julia> lookup(OneHot, vocab, "a")
3-element OneHot{3}:
 1
 0
 0

julia> lookup(OneHot, vocab, ["a" "b"; "c" "d"])
3x2x2 OneHotArray{3, 3, Matrix{OneHot{0x00000003}}}:
[:, :, 1] =
 1  0
 0  0
 0  1

[:, :, 2] =
 0  0
 1  0
 0  0

julia> lookup(OneHot, vocab, 3)
ERROR: DomainError with c:
cannot convert `lookup(::Vocab, 3)` = "c" into one-hot representation.
Stacktrace:
[...]

```
"""
function lookup(::Type{OneHot}, v::Vocab, i) end
).args[2])

@eval $((@macroexpand @doc """
    lookup(v::Vocab, i::OneHotArray)

convert the one-hot representation back into words.

# Example
```julia
julia> lookup(OneHot, vocab, ["a" "b"; "c" "d"])
3x2x2 OneHotArray{3, 3, Matrix{OneHot{0x00000003}}}:
[:, :, 1] =
 1  0
 0  0
 0  1

[:, :, 2] =
 0  0
 1  0
 0  0

julia> lookup(vocab, ans)
2×2 Matrix{String}:
 "a"  "b"
 "c"  "[UNK]"

```
"""
function lookup(v::Vocab, i::OneHotArray) end
).args[2])


text:
function add_pad_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(pad_gen_cb,b["pad_check_button"],"clicked",Void,(),false,(handles,))
    nothing
end

view_pad(b::Gtk.GtkBuilder)=getproperty(b["pad_check_button"],:active,Bool)

function pad_gen_cb(w::Ptr, user_data::Tuple{Tracker_Handles})

    han, = user_data

    if getproperty(han.b["pad_check_button"],:active,Bool)
        han.selection_mode = 10
    else
        han.selection_mode = 1
    end

    redraw_all(han)

    nothing
end


text:
module KernelCompiler

using Core.Compiler
using Core.Compiler: MethodInstance, NativeInterpreter, CodeInfo, CodeInstance, WorldView, OptimizationState
using Base.Meta

using LLVM
using LLVM.Interop

include("codecache.jl")
const CACHE = Dict{DataType, CodeCache}()
get_cache(ai::DataType) = CACHE[ai]

include("kernelinterpreter.jl")
include("optimize.jl")
include("codegen.jl")

# We have one global JIT and TM
const orc = Ref{LLVM.OrcJIT}()
const tm  = Ref{LLVM.TargetMachine}()

function __init__()
    CACHE[NativeInterpreter] = CodeCache(cpu_invalidate)

    opt_level = Base.JLOptions().opt_level
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end

    tm[] = LLVM.JITTargetMachine(optlevel=optlevel)
    LLVM.asm_verbosity!(tm[], true)

    orc[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
    atexit() do
        LLVM.dispose(orc[])
    end
end

include("pipeline.jl")
include("jit.jl")

end # module


text:
using SimilaritySearch

function main(n, m, dim)
    S = MatrixDatabase(rand(Float32, dim, n))
    Q = MatrixDatabase(rand(Float32, dim, m))

    seq = ExhaustiveSearch(SqL2Distance(), S)
    k = 10
    @info "running matrix knnresult"
    @time searchbatch(seq, Q, k)
    knnlist = [KnnResult(k) for _ in 1:m]
    @info "running vector knnresult"
    @time searchbatch(seq, Q, knnlist)
end

@info "warming"
@info size(main(100, 10, 3))
@info "large benchmark"
@info size(main(100_000, 1000, 3))

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
"""
```julia
    NegBinomialLikelihood(r::Int=10)
```

[Negative Binomial likelihood](https://en.wikipedia.org/wiki/Negative_binomial_distribution) with number of failures `r`
```math
    p(y|r,f) = binomial(y+r-1,y) (1-σ(f))ʳσ(f)ʸ
```
Where `σ` is the logistic function

"""
struct NegBinomialLikelihood{T<:Real} <: EventLikelihood{T}
    r::Int
    c::Vector{T}
    θ::Vector{T}
    function NegBinomialLikelihood{T}(r::Int) where {T<:Real}
        new{T}(r)
    end
    function NegBinomialLikelihood{T}(r,c,θ) where {T<:Real}
        new{T}(r,c,θ)
    end
end

function NegBinomialLikelihood(r::Int=10) where {T<:Real}
    NegBinomialLikelihood{Float64}(r)
end

function init_likelihood(likelihood::NegBinomialLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Int,nFeatures::Int) where T
    NegBinomialLikelihood{T}(
    likelihood.r,
    rand(T,nSamplesUsed),
    zeros(T,nSamplesUsed))
end

function pdf(l::NegBinomialLikelihood,y::Real,f::Real)
    pdf(NegativeBinomial(lr,get_p(l,f)),y)
end

function expec_count(l::NegBinomialLikelihood,f)
    broadcast((p,r)->p*r./(1.0.-p) ,get_p.(l,f),l.r)
end

function get_p(::NegBinomialLikelihood,f)
    logistic.(f)
end

function Base.show(io::IO,model::NegBinomialLikelihood{T}) where T
    print(io,"Negative Binomial Likelihood")
end

function compute_proba(l::NegBinomialLikelihood{T},μ::Vector{T},σ²::Vector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        x = pred_nodes.*sqrt(max(σ²[i],zero(T))).+μ[i]
        pred[i] = dot(pred_weights,get_p.(l,x))
    end
    return pred
end

## Local Updates ##

function local_updates!(l::NegBinomialLikelihood{T},y::AbstractVector,μ::AbstractVector,Σ::AbstractVector) where {T}
    l.c .= sqrt.(abs2.(μ) + Σ)
    l.θ .= (l.r.+y)./l.c.*tanh.(0.5*l.c)
end

## Global Updates ##

@inline ∇E_μ(l::NegBinomialLikelihood{T},::AOptimizer,y::AbstractVector) where {T} = (0.5*(y.-l.r),)
@inline ∇E_Σ(l::NegBinomialLikelihood{T},::AOptimizer,y::AbstractVector) where {T} = (0.5.*l.θ,)

## ELBO Section ##

AugmentedKL(l::NegBinomialLikelihood{T},y::AbstractVector) where {T} = PolyaGammaKL(l,y)

function logabsbinomial(n,k)
    log(binomial(n,k))
end

function expec_log_likelihood(l::NegBinomialLikelihood{T},i::AnalyticVI,y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    tot = sum(logabsbinomial.(y.+(l.r-1),y))-log(2.0)*sum((y.+l.r))
    tot += 0.5*dot(μ,(y.-l.r))-0.5*dot(l.θ,μ)-0.5*dot(l.θ,diag_cov)
    return tot
end

function PolyaGammaKL(l::NegBinomialLikelihood,y::AbstractVector) where {T}
    PolyaGammaKL(y.+l.r,l.c,l.θ)
end


text:
using ForestBiometrics

tl1 = [
Tree(10.0,65,"WO",1.0),
Tree(10.0,65,"WO",1.0),
Tree(10.0,65,"WO",1.0),
Tree(10.0,65,"WO",1.0)]

tl2 = [
Tree(10.0,65,"WO",1.0),
Tree(10.0,60,"WO",1.0),
Tree(10.0,65,"WO",1.0),
Tree(10.0,65,"WO",1.0)]

tl3 = [
Tree(10.0,65,"WO",1.0),
Tree(10.0,60,"WO",1.0),
Tree(10.0,60,"WO",1.0),
Tree(10.0,65,"WO",1.0),
Tree(10.0,70,"WO",1.0)]

@test length(collapse(tl1)) == 1
@test length(collapse(tl2)) == 2
@test length(collapse(tl3)) == 3

t1 = Tree(12.0,65,"WO",4.25)
t2 = Tree(12.0,65,"WO",10)

tl4 = [Tree(12.0,60,"WO",5),
Tree(10.0,65,"WO",5),
Tree(12.0,60,"WO",1),
Tree(10.0,65,"WO",1)]

@test length(expand(t1)) == 5
@test length(expand(t1, min_exp = 2)) == 3
@test length(expand(tl4)) == 12

@test length(t2 / 3.0) == 3

s1 = Stand(tl4)
expand!(s1)

@test length(s1.treelist) == 12


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

using Markdown
using InteractiveUtils

# ╔═╡ e7e7c0e4-47e1-4856-9e71-11dfd12dfcd3
using StatsBase

# ╔═╡ 8c0c6652-aa5d-4cb8-a4cf-018c80b380cf
student_names = "Triana Forment";

# ╔═╡ f3aa9d9c-75e1-4948-b7dd-7f0596486e83
md"""
# Project: Hungarian algorithm

**STMO**

2021-2022

project by $student_names

## Introduction

Optimal transportation leads with two probability distributions and there is a cost function for moving elements from one distribution to the other distribution. And you have to find a transportation scheme to map one distribution into the other one. The original version of this problem of optimal transport dates back to Gaspard Monge in 1781. The Monge's problem has a discrete n number iron mines and n factories and, the cost is the distance between the mine and the factory. The solution is to find which mine suplies which factory, finding the minimum average distance transported.

The Hungarian algorithm solves an assignment problem which can be considered as a special type of transportation problem in which the number of sources and sinks are equal. The capacity of each source as well as the requirement of each sink is taken as 1. In this project, the Hungarian algorithm is implemented using an adjacency matrix.
This algorithm can be used to find the map that minimizes the transport cost, given two distributions and a cost function. For two discrete probability vectors, $$\mathbf{a} ∈ Σ_n$$ and $$\mathbf{b} ∈ Σ_m$$, we have a $$n \times m$$ cost matrix $$C$$. It can be represented as a adjacency matrix, with elements in $$a$$ as rows and elements in $$b$$ as columns, and their weights as entries in the matrix.

An example of an assignment problem, can be, given $$n$$ agents and the money they ask for performing each task, finding which agent performs which task in order to minimize the cost of performing all $$m$$ tasks.
When there are the same number of agents than tasks, the the problem is called balanced assignment, otherwise is called unbalanced assignment. And if the total cost is the sum for each agent performing its task, the problem is called linear assignment, because the cost function to be optimized as well as all the constraints contain only linear terms. This implementation can be used to solve linear balanced assignment problems.
"""

# ╔═╡ bf771d76-099c-449b-bebd-5201ed72ebb6
md"""
# Implementation
"""

# ╔═╡ ed527906-da70-4186-a6db-da9906d72efa
md"""
## Outline of the approach

1. Subtract the smallest entry in each row from all the entries in the row, thus making the smallest entry in the row now equal to `0`. The same for each column.

2. Mark the rows and columns which have 0 entries, such that the fewest lines possible are drawn.

3. If there are $m$ marked rows and columns, an optimal assignment of zeros is possible and the algorithm is finished. The cost value and the optimal agent-task combination can be calculated. If the number of lines is less than $$m$$, then the optimal number of zeroes is not yet reached and the next step needs to be performed.

4. Find the smallest entry not marked by any marked row and by any column. Subtract this entry from each row which is not marked, and then add it to each column which is marked. Then, go back to step 2.
"""

# ╔═╡ 6be9cccf-160a-4d05-88bd-9f5e33ad3d05
md"""
**`Hungarian_algorithm`** is the funtion that uses as an input the cost matrix and gives as an output the cost value and the matrix solution. This function performs steps 1 to 4, by calling the previous functions. Once the marked rows and columns sum the number of tasks to perform, it calculates the cost value by adding each cost value of the different agent-task selected. The matrix solution has all zeros, except for each agent-task pair selected, where it has its original value.
"""

# ╔═╡ 23020cd7-333d-4d97-aabb-c13614d57c65
md"""
**`find_min_row`** is a function that gets the row with the fewest zeros. For that purpose, it needs as an input a boolean matrix, having `true` if the original value after subtracting the smallest entry in each row from all the other entries in the row was `0`, and a `false` if it was different from `0`. It also needs a vector `zero_list` that will store the coordinates of the `0` in the row with fewest zeros. After the first finding, it reajusts the boolean matrix so the whole row and column from the last coordinate stored is set to `false`. This function is used in the next function `mark_matrix`.
"""

# ╔═╡ d7c22bc5-088d-4a74-a00d-5833430bf219
"""
	find_min_row(matrix_zero, zero_list)

Inputs:
	- `matrix_zero`: boolean matrix
	- `zero_list`: array to store the coordinates of the zeros in the row with fewest zeros.

Outputs:
	- `min_row`: tuple with the first element the number of zeros in the row with fewest zeros and the index of that row as second element
	- `zero_list`: list of coordinates tuples of the zeros found
"""
function find_min_row(matrix_zero, zero_list)
	min_row = [Inf, -1]

	for row in 1:size(matrix_zero)[1]
		#if true in matrix_zero[row, :] # Needs a 0 in the row
			# If the number of zeros < than the last min stored
			if sum(matrix_zero[row,:]) > 0 && min_row[1] > sum(matrix_zero[row,:])
				#stores the number of zeros and the row index
				min_row = [sum(matrix_zero[row,:]), row]
			end
	end
	# Get the column index
	zero_index = findall(x->x==true, matrix_zero[min_row[2],:])[1]
	# Store tuple in zero_list the tuple of coordinates
	append!(zero_list, [(min_row[2], zero_index)])
	# Mark the specific row and column as false.
	matrix_zero[min_row[2], :] .= false
	matrix_zero[:, zero_index] .= false
	return min_row, zero_list
end

# ╔═╡ 4f82e659-944a-4ff2-aa77-db81187da001
md"""
**`mark_matrix`** needs as an input the matrix with the already subtracted the smallest entry in each row from all the other entries and it returns the `zero_list` and the marked row and column indexes. For that, it creates the boolean matrix and calls the function `find_min_row`. It performs steps 2 and 3.
"""

# ╔═╡ 8d04677a-5f97-4e76-937c-e79f69c14788
"""
	mark_matrix(mat)

Perform steps 2 and 3: It marks the rows and columns that have the 0 entries such that the fewest lines possible are drawn. If there are "m" marked rows and columns, an optimal assignment solution can be found.

Input:
	- `mat` : cost matrix already modified by step 1 or an adjusted matrix by `adjust_matrix`

Outputs:
	- `zero_list` : coordinates of the zeros found by `find_min_row`
	- `marked_rows`: array of marked rows
	- `marked_cols`: array of marked columns
"""
function mark_matrix(mat)

	#Transform the matrix to boolean matrix(0 = true, others = false)
	cur_mat = mat
	boolean_mat = iszero.(mat)
	boolean_mat_copy = copy(boolean_mat)
	#Recording possible answer positions by zero_list
    zero_list = []
    while true in boolean_mat_copy
        find_min_row(boolean_mat_copy, zero_list)
	end
    #Recording the row indexes
    zero_list_row = []
    for i in 1:length(zero_list)
        append!(zero_list_row, zero_list[i][1])
	end

	# Get non marked rows
    non_marked_row = (x->collect(x))(setdiff(Set(1:size(cur_mat)[1]), Set(zero_list_row)))

    marked_cols = []
    flag = 0
    while flag == 0 # Enter in the loop until there's no more unmarked cols and rows
        flag = 1
        for i in 1:length(non_marked_row)
            row_array = boolean_mat[non_marked_row[i], :]
            for j in 1:length(row_array)
				# Find unmarked 0 elements in the corresponding column
				if row_array[j] == true && j ∉ marked_cols
                    # Store column index in "marked_cols"
                    append!(marked_cols, j)
                    flag = 0
				end
			end
		end
		zero_rows = [t[1] for t in zero_list]
		zero_cols = [t[2] for t in zero_list]
		for i in 1:length(zero_list)
			# If in "zero_list" coordinates there was a 0 marked in cols,
			# but not in marked_rows, add its row index to "non_marked_row"
            if zero_rows[i] ∉ non_marked_row && zero_cols[i] in marked_cols
                append!(non_marked_row, zero_rows[i])
                flag = 0
			end
		end
	end
	# Add those the indexes not stored in "non_marked_row" to "marked_rows"
	marked_rows = (x->collect(x))(setdiff(Set(1:size(cur_mat)[1]), Set(non_marked_row)))
    return zero_list, marked_rows, marked_cols
end

# ╔═╡ 05e77871-d6d1-4d7b-b368-89aeeed878a4
md"""
**`adjust_matrix`** is a function needed when, after performing steps 1, 2 and 3, the number of marked rows plus the number of marked columns is not equal to the number of tasks/activities to perform by the agents. Then, it adjustes and returns the matrix as described in step 4.
"""

# ╔═╡ f94e0b59-80d0-4d26-bac3-d75c30e4d3ee
"""
	adjust_matrix(mat, marked_rows, marked_cols)

Perform step 4: Adjust the matrix by finding the smallest entry not marked by any marked row and by any column. Subtract this entry from each row that isn’t marked, and then add it to each column that is marked.

Inputs:
	- `mat` : boolean matrix created and modified by `mark_matrix`
	- `marked_rows`: array of marked rows produced by `mark_matrix`
	- `marked_cols`: array of marked columns produced by `mark_matrix`

Output:
	- `adjusted`: the adjusted matrix
"""
function adjust_matrix(mat, marked_rows, marked_cols)
	adjusted = mat
    min_values = []
    n = size(adjusted)[1]

    # Find the minimum value for that is not in marked_rows and not in marked_cols
    for row in 1:n
        if row ∉ marked_rows
            for i in 1:n
                if i ∉ marked_cols
                    append!(min_values, adjusted[row,i])
				end
			end
		end
	end
    min_value = minimum(min_values)
	# Substract that min value to all elements not in marked rows and columns
    for row in 1:n
        if row ∉ marked_rows
            for i in 1:n
                if i ∉ marked_cols
                    adjusted[row, i] = adjusted[row, i] - min_value
				end
			end
		end
	end
    # Add the min value to elements that are in both marked row and columns
    for row in marked_rows
        for col in marked_cols
            adjusted[row, col] = adjusted[row,col] + min_value
		end
	end
    return adjusted
end

# ╔═╡ d4435d9e-1c3a-4d0d-80b3-8c0061336975
"""
	Hungarian_algorithm(mat)

Solves a linear balanced assignment problem.

Input:
	- 'mat`: the cost matrix

Outputs:
	- `cost`: the resulting cost value
	- `matrix solution`: the matrix with the selected agent-task as its original cost value and the rest zeros.
"""
function Hungarian_algorithm(mat)
	num_rows = size(mat)[1]
	num_cols = size(mat)[1]
		 @assert num_rows == num_cols throw(DimensionMismatch("It is not balanced. Number of sources and sinks need to match"))

	matrix = copy(mat)
	#Step 1: subtract its internal minimum from every column and every row
	matrix .-= minimum(matrix, dims=2) # from rows
	matrix .-= minimum(matrix, dims=1) # from columns

	num_zeros = 0

	# Step 2: mark rows and columns that have the 0 entries
	pos = []
    while num_zeros < num_rows
        pos, rows, cols = mark_matrix(matrix)
		rows = reverse(rows)
        num_zeros = length(rows) + length(cols)
		# Step 3: if marked rows & columns is less than num_rows, go to step 4.
        if num_zeros < num_rows
			# Step 4: adjust the matrix
            matrix = adjust_matrix(matrix, rows, cols)
		end
	end
	# Calculate total cost and matrix solution
	cost = 0
	matrix_solution = zeros(Int64, size(mat)[1], size(mat)[2])
	for i in 1:num_rows
		cost += mat[pos[i][1], pos[i][2]]
		matrix_solution[pos[i][1], pos[i][2]] = mat[pos[i][1], pos[i][2]]
	end
	return cost, matrix_solution
end

# ╔═╡ d8dd3d3f-2ed6-434c-8878-c3af69abc45b
md"""
### Example and solution of assignment problem

Suppose we have 5 companies having different costs to realize 5 different tasks. You want to select a different company for each task in the way that we spend the minimum money. In order to find that combination, you just have to create a matrix with the cost each company asks for performing each activity. Here we have the following cost matrix with the companies as rows, the tasks as columns, and the values in thousand euros unit.
"""

# ╔═╡ 4cbb69df-9573-4b9f-a03a-4d1e3c53b8ff
cost_matrix = [7 6 2 9 2;
               1 2 5 3 9;
               5 3 9 6 5;
               9 2 5 8 7;
               2 5 3 6 1]

# ╔═╡ 8f59c095-4cab-4d2f-bfe5-85ab121070bd
md"""
After applying the function `Hungarian_algorithm`, you get the minimum cost (12000 €),
and the solution matrix, which shows the company-task combination to get that minimum cost. The first company will perform the third task for 2000 €, the second company the first task for 1000 €, etc.
"""

# ╔═╡ 545a38bf-4edc-4c16-9f1a-0c26ff2e5f33
Hungarian_algorithm(cost_matrix)

# ╔═╡ 011e0f2f-ac2f-4ca6-94a9-912dee4b21ae
md"""
In case you want to maximize the values (e.g. instead of cost values you have benefit values), you just have to create the benefit matrix with negative values. This would solve an assigment problem maximizing the solution.
"""

# ╔═╡ 0ad00f08-c712-4dde-892f-e48a8b7f02d4
begin
	benefit_matrix = copy(cost_matrix)
	benefit_matrix .= - cost_matrix
end

# ╔═╡ 50b784ac-64a1-41eb-ba01-fd1af65f6fe7
md"""
After applying `Hungarian_agorithm` to the benefit matrix, you get a maximum benefit of 41000 € by selecting the first company performing the fourth task, the second company the fifth, etc.
"""

# ╔═╡ b76c5a89-fb4c-4840-b5e7-ce816ba4b4fc
Hungarian_algorithm(benefit_matrix)

# ╔═╡ dfd65b50-b0a7-4ad1-83d8-7bfa0bedfed7
md"""
## References

- H. W. Kuhn. The Hungarian Method for the Assignment Problem. (2010), from http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Lec10-HungarianMethod-Kuhn.pdf
- M. Thorpe. Introduction to Optimal Transport. (2018), from https://www.math.cmu.edu/~mthorpe/OTNotes
- Chapter 6 - Selected Topic in Mathematical Optimization course by Michiel Stock
- https://brilliant.org/wiki/hungarian-matching/
- https://python.plainenglish.io/hungarian-algorithm-introduction-python-implementation-93e7c0890e15

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

[compat]
StatsBase = "~0.33.13"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.2"

[[ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2bb0cb32026a66037360606510fca5984ccc6b75"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.13"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e7e7c0e4-47e1-4856-9e71-11dfd12dfcd3
# ╟─f3aa9d9c-75e1-4948-b7dd-7f0596486e83
# ╟─8c0c6652-aa5d-4cb8-a4cf-018c80b380cf
# ╟─bf771d76-099c-449b-bebd-5201ed72ebb6
# ╟─ed527906-da70-4186-a6db-da9906d72efa
# ╟─6be9cccf-160a-4d05-88bd-9f5e33ad3d05
# ╠═d4435d9e-1c3a-4d0d-80b3-8c0061336975
# ╠═23020cd7-333d-4d97-aabb-c13614d57c65
# ╠═d7c22bc5-088d-4a74-a00d-5833430bf219
# ╠═4f82e659-944a-4ff2-aa77-db81187da001
# ╠═8d04677a-5f97-4e76-937c-e79f69c14788
# ╟─05e77871-d6d1-4d7b-b368-89aeeed878a4
# ╠═f94e0b59-80d0-4d26-bac3-d75c30e4d3ee
# ╟─d8dd3d3f-2ed6-434c-8878-c3af69abc45b
# ╠═4cbb69df-9573-4b9f-a03a-4d1e3c53b8ff
# ╟─8f59c095-4cab-4d2f-bfe5-85ab121070bd
# ╠═545a38bf-4edc-4c16-9f1a-0c26ff2e5f33
# ╟─011e0f2f-ac2f-4ca6-94a9-912dee4b21ae
# ╠═0ad00f08-c712-4dde-892f-e48a8b7f02d4
# ╟─50b784ac-64a1-41eb-ba01-fd1af65f6fe7
# ╠═b76c5a89-fb4c-4840-b5e7-ce816ba4b4fc
# ╟─dfd65b50-b0a7-4ad1-83d8-7bfa0bedfed7
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: triforment/STMOZOO

max_issues_repo_name: triforment/STMOZOO

max_forks_repo_name: triforment/STMOZOO


text:
# ---
# title: 1014. Best Sightseeing Pair
# id: problem1014
# author: Indigo
# date: 2021-01-30
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/best-sightseeing-pair/description/>
# hidden: true
# ---
#
# Given an array `A` of positive integers, `A[i]` represents the value of the
# `i`-th sightseeing spot, and two sightseeing spots `i` and `j` have distance
# `j - i` between them.
#
# The _score_  of a pair (`i < j`) of sightseeing spots is (`A[i] + A[j] + i -
# j)` : the sum of the values of the sightseeing spots, **minus** the distance
# between them.
#
# Return the maximum score of a pair of sightseeing spots.
#
#
#
# **Example 1:**
#
#
#
#     Input: [8,1,5,2,6]
#     Output: 11
#     Explanation: i = 0, j = 2, A[i] + A[j] + i - j = 8 + 5 + 0 - 2 = 11
#
#
#
#
# **Note:**
#
#   1. `2 <= A.length <= 50000`
#   2. `1 <= A[i] <= 1000`
#
#
## @lc code=start
using LeetCode

function max_score_sightseeing_pair(A::Vector{Int})
    res, mx = 0, A[1] + 0
    for j in 2:length(A)
        res = max(res, A[j] - j + 1 + mx)
        mx = max(mx, A[j] + j - 1)
    end
    res
end
## @lc code=end


text:
using Test, Nonconvex

@test_throws ArgumentError using NonconvexIpopt
Nonconvex.@load Ipopt
IpoptAlg()

@test_throws ArgumentError using NonconvexNLopt
Nonconvex.@load NLopt
NLoptAlg(:LD_MMA)

@test_throws ArgumentError using NonconvexJuniper
Nonconvex.@load Juniper
JuniperIpoptAlg()
IpoptAlg()

@test_throws ArgumentError using NonconvexMMA
Nonconvex.@load MMA
MMA87()
MMA02()

@test_throws ArgumentError using NonconvexPavito
Nonconvex.@load Pavito
PavitoIpoptCbcAlg()
IpoptAlg()

@test_throws ArgumentError using NonconvexPercival
Nonconvex.@load AugLag
AugLag()

@test_throws ArgumentError using NonconvexBayesian
Nonconvex.@load Bayesian
BayesOptAlg(IpoptAlg())

#@test_throws ArgumentError using NonconvexAugLagLab
#Nonconvex.@load AugLag2
#AugLag2()

@test_throws ArgumentError using NonconvexSemidefinite
Nonconvex.@load Semidefinite
SDPBarrierAlg(IpoptAlg())

@test_throws ArgumentError using NonconvexSearch
Nonconvex.@load Search
MTSAlg()
LS1Alg()

@test_throws ArgumentError using NonconvexMultistart
Nonconvex.@load Multistart
HyperoptAlg(IpoptAlg())
DeflatedAlg(IpoptAlg())