text:
"""
    has_constant_lags(integrator::DDEIntegrator)

Return if the DDE problem of the `integrator` contains constant delays.
"""
has_constant_lags(integrator::DDEIntegrator) = has_constant_lags(integrator.sol.prob)

"""
    has_dependent_lags(integrator::DDEIntegrator)

Return if the DDE problem of the `integrator` contains dependent delays.
"""
has_dependent_lags(integrator::DDEIntegrator) = has_dependent_lags(integrator.sol.prob)

"""
    has_constant_lags(prob::DDEProblem)

Return if the DDE problem `prob` contains constant delays.
"""
has_constant_lags(prob::DDEProblem) =
  prob.constant_lags !== nothing && !isempty(prob.constant_lags)

"""
    has_dependent_lags(prob::DDEProblem)

Return if the DDE problem `prob` contains dependent delays.
"""
has_dependent_lags(prob::DDEProblem) =
  prob.dependent_lags !== nothing && !isempty(prob.dependent_lags)

"""
    u_uprev(u0, alg; kwargs...)

Return state vectors `u` and `uprev` (possibly aliased) for solving the
differential equation problem for initial state `u0` with algorithm `alg`.
"""
function u_uprev(u0, alg;
                 alias_u0 = false,
                 adaptive = DiffEqBase.isadaptive(alg),
                 calck = false)
  if alias_u0
    u = u0
  else
    u = recursivecopy(u0)
  end

  # Some algorithms do not use `uprev` explicitly. In that case, we can save
  # some memory by aliasing `uprev = u`, e.g. for "2N" low storage methods.
  if OrdinaryDiffEq.uses_uprev(alg, adaptive) || calck
    uprev = recursivecopy(u)
  else
    uprev = u
  end

  u, uprev
end

"""
    u_uprev_uprev2(u0, alg; kwargs...)

Return state vectors `u`, `uprev`, and `uprev2` (possibly aliased) for solving the
differential equation problem for initial state `u0` with algorithm `alg`.
"""
function u_uprev_uprev2(u0, alg;
                        allow_extrapolation = alg_extrapolates(alg),
                        kwargs...)
  # compute u and uprev first
  u, uprev = u_uprev(u0, alg; kwargs...)

  if allow_extrapolation
    uprev2 = recursivecopy(u)
  else
    uprev2 = uprev
  end

  u, uprev, uprev2
end

"""
    get_abstol(u, tspan, alg; abstol = nothing)

Return the absolute tolerance for solving the differential equation problem with state
variable `u` and time span `tspan` with algorithm `alg`.
"""
function get_abstol(u, tspan, alg; abstol = nothing)
  if typeof(alg) <: FunctionMap
    _abstol = real.(zero.(u))
  elseif abstol === nothing
    uBottomEltype = recursive_bottom_eltype(u)
    uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

    if uBottomEltypeNoUnits == uBottomEltype
      _abstol = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1//10^6))
    else
      _abstol = real.(oneunit.(u).*1//10^6)
    end
  else
    _abstol = real.(abstol)
  end

  _abstol
end

"""
    get_reltol(u, tspan, alg; reltol = nothing)

Return the relative tolerance for solving the differential equation problem with state
variable `u` and time span `tspan` with algorithm `alg`.
"""
function get_reltol(u, tspan, alg; reltol = nothing)
  if typeof(alg) <: FunctionMap
    _reltol = real.(zero(first(u)/t))
  elseif reltol === nothing
    uBottomEltype = recursive_bottom_eltype(u)
    uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

    if uBottomEltypeNoUnits == uBottomEltype
      _reltol = real(convert(uBottomEltype, oneunit(uBottomEltype) * 1//10^3))
    else
      _reltol = real.(oneunit.(u).*1//10^3)
    end
  else
    _reltol = real.(reltol)
  end

  _reltol
end

"""
    callback_set_and_cache(prob, callback)

Return set of callbacks and its cache for the differential equation problem `prob` and the
user-provided `callback`.
"""
function callback_set_and_cache(prob, callback)
  callback_set = CallbackSet(callback)

  max_len_cb = DiffEqBase.max_vector_callback_length(callback_set)
  if max_len_cb isa VectorContinuousCallback
    uBottomEltype = recursive_bottom_eltype(prob.u0)
    callback_cache = DiffEqBase.CallbackCache(max_len_cb.len, uBottomEltype, uBottomEltype)
  else
    callback_cache = nothing
  end

  callback_set, callback_cache
end

"""
    rate_prototype_of(u0, tspan)

Return prototype of rates for a given differential equation problem with state `u` and
time span `tspan`.
"""
rate_prototype_of(u0, tspan) =
  DiffEqBase.@.. u0 * $(inv(oneunit(eltype(tspan))))

"""
    solution_arrays(u, tspan, rate_prototype; kwargs...)

Return arrays of saved time points, states, and rates, initialized with the solution at the
first time point if `save_start = true` (the default).
"""
function solution_arrays(u, tspan, rate_prototype;
                         timeseries_init = typeof(u)[],
                         ts_init = eltype(tspan)[],
                         ks_init = [],
                         save_idxs = nothing,
                         save_start = true)
  # determine types of time and state
  uType = typeof(u)
  tType = eltype(tspan)

  # initialize vector of saved time points
  ts = convert(Vector{tType}, ts_init)

  # initialize vector of saved states
  if save_idxs === nothing
    timeseries = convert(Vector{uType}, timeseries_init)
  else
    u_initial = u[save_idxs]
    timeseries = convert(Vector{typeof(u_initial)}, timeseries_init)
  end

  # initialize vector of saved rates
  if save_idxs === nothing
    ksEltype = Vector{typeof(rate_prototype)}
  else
    ks_prototype = rate_prototype[save_idxs]
    ksEltype = Vector{typeof(ks_prototype)}
  end
  ks = convert(Vector{ksEltype}, ks_init)

  # save solution at initial time point
  if save_start
    copyat_or_push!(ts, 1, first(tspan))
    if save_idxs === nothing
      copyat_or_push!(timeseries, 1, u)
      copyat_or_push!(ks, 1, [rate_prototype])
    else
      u_initial = u[save_idxs]
      copyat_or_push!(timeseries, 1, u_initial, Val{false})
      copyat_or_push!(ks, 1, [ks_prototype])
    end
  end

  ts, timeseries, ks
end

"""
    sizehint!(sol::DESolution, n)

Suggest that solution `sol` reserves capacity for at least `n` elements.
"""
function Base.sizehint!(sol::DESolution, n)
  sizehint!(sol.u, n)
  sizehint!(sol.t, n)
  sizehint!(sol.k, n)

  nothing
end

"""
    sizehint!(sol::DESolution, alg, tspan, tstops, saveat; kwargs...)

Suggest that solution `sol` reserves capacity for a number of elements that
depends on the parameter settings of the numerical solver.
"""
function Base.sizehint!(sol::DESolution, alg, tspan, tstops, saveat;
                        save_everystep = isempty(saveat),
                        adaptive = isadaptive(alg),
                        internalnorm = DiffEqBase.ODE_DEFAULT_NORM,
                        dt = zero(eltype(tspan)))
  # obtain integration time
  t0 = first(tspan)
  integrationtime = last(tspan) - t0

  if !adaptive && save_everystep && !isinf(integrationtime)
    # determine number of steps if known a priori
    if iszero(dt)
      steps = length(tstops)
    else
      steps = ceil(Int, internalnorm(integrationtime / dt, t0))
    end

    sizehint!(sol, steps + 1)
  elseif save_everystep
    sizehint!(sol, 50)
  elseif !isempty(saveat)
    sizehint!(sol, length(saveat) + 1)
  else
    sizehint!(sol, 2)
  end

  nothing
end

function build_history_function(prob, alg, rate_prototype, reltol;
                                dt = zero(eltype(prob.tspan)),
                                adaptive = DiffEqBase.isadaptive(alg.alg),
                                calck = false,
                                internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
  @unpack f, u0, tspan, p = prob

  t0 = first(tspan)
  tType = eltype(tspan)
  tTypeNoUnits = typeof(one(tType))
  tdir = sign(last(tspan) - t0)

  uEltypeNoUnits = recursive_unitless_eltype(u0)
  uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u0)

  # bootstrap an ODE integrator
  # - whose solution captures the dense history of the simulation
  # - that is used for extrapolation of the history for time points past the
  #   already fixed history
  # - that is used for interpolation of the history for time points in the
  #   current integration step (so the interpolation is fixed while updating the stages)
  # we wrap the user-provided history function such that function calls during the setup
  # of the integrator do not fail
  ode_f = ODEFunctionWrapper(f, prob.h)
  ode_prob = ODEProblem{isinplace(prob)}(ode_f, u0, tspan, p)

  # get states of ODE integrator (do not alias uprev)
  ode_u, ode_uprev = u_uprev(u0, alg; alias_u0 = false, calck = true)

  # initialize output arrays
  ode_k = typeof(rate_prototype)[]
  ode_ts, ode_timeseries, ode_ks = solution_arrays(ode_u, tspan, rate_prototype;
                                                   save_idxs = nothing,
                                                   save_start = true)

  # obtain cache (we alias uprev2 and uprev)
  ode_cache = OrdinaryDiffEq.alg_cache(alg.alg, ode_u, rate_prototype, uEltypeNoUnits,
                                       uBottomEltypeNoUnits, tTypeNoUnits, ode_uprev,
                                       ode_uprev, ode_f, t0, zero(tType), reltol, p, calck,
                                       Val(isinplace(prob)))

  # build dense interpolation of history
  if iscomposite(alg)
    ode_alg_choice = Int[]
    ode_id = OrdinaryDiffEq.CompositeInterpolationData(ode_f, ode_timeseries, ode_ts, ode_ks,
                                                       ode_alg_choice, true, ode_cache) # dense = true
    ode_sol = DiffEqBase.build_solution(ode_prob, alg.alg, ode_ts, ode_timeseries;
                                        dense = true, k = ode_ks, interp = ode_id,
                                        alg_choice = ode_alg_choice,
                                        calculate_error = false, destats = DiffEqBase.DEStats(0))
  else
    ode_id = OrdinaryDiffEq.InterpolationData(ode_f, ode_timeseries, ode_ts, ode_ks, true, ode_cache) # dense = true
    ode_sol = DiffEqBase.build_solution(ode_prob, alg.alg, ode_ts, ode_timeseries;
                                        dense = true, k = ode_ks, interp = ode_id,
                                        calculate_error = false, destats = DiffEqBase.DEStats(0))
  end

  # reserve capacity
  sizehint!(ode_sol, alg.alg, tspan, (), ();
            save_everystep = true, adaptive = adaptive, internalnorm = internalnorm, dt = tType(dt))

  # create simple integrator
  tdirType = typeof(sign(zero(tType)))
  ode_integrator = HistoryODEIntegrator{typeof(alg.alg),isinplace(prob),typeof(prob.u0),
                                        tType,tdirType,typeof(ode_k),
                                        typeof(ode_sol),typeof(ode_cache)}(
                                          ode_sol, ode_u, ode_k, t0, zero(tType), ode_uprev,
                                          t0, alg.alg, zero(tType), tdir, 1, 1, ode_cache)

  # combine the user-provided history function and the ODE integrator with dense solution
  # to a joint dense history of the DDE
  # we use this history information to create a problem function of the DDE with all
  # available history information that is of the form f(du,u,p,t) or f(u,p,t) such that
  # ODE algorithms can be applied
  HistoryFunction(prob.h, ode_integrator)
end

"""
    initialize_solution!(integrator::DDEIntegrator)

Initialize the solution of an integrator by adjusting the cache for composite algorithms.
"""
function initialize_solution!(integrator::DDEIntegrator)
  if iscomposite(integrator.alg)
    copyat_or_push!(integrator.integrator.sol.alg_choice, 1, integrator.cache.current)
    if integrator.opts.save_start
      copyat_or_push!(integrator.sol.alg_choice, 1, integrator.cache.current)
    end
  end

  nothing
end

function unwrap_alg(integrator::DDEIntegrator, is_stiff)
  alg = integrator.alg
  iscomp = typeof(alg) <: CompositeAlgorithm
  if !iscomp
    return alg
  elseif typeof(alg.choice_function) <: AutoSwitch
    num = is_stiff ? 2 : 1
    return alg.algs[num]
  else
    return alg.algs[integrator.cache.current]
  end
end

DiffEqBase.nlsolve_f(integrator::DDEIntegrator) =
  DiffEqBase.nlsolve_f(integrator.f, unwrap_alg(integrator, true))


text:
export update_cell_list!

"""
    update_cell_list!(cell_list)

Updates cell list with current positions of select particles stored under
`cell_list.particles`.
"""
function update_cell_list!(cell_list::CellList)
    fill!(cell_list.start_pid, -1)
    fill!(cell_list.next_pid, -1)

    @inbounds for (n, particle) in enumerate(cell_list.particles)
        i = trunc(Int64, particle.x / cell_list.cell_spacing_x) + 1
        j = trunc(Int64, particle.y / cell_list.cell_spacing_y) + 1

        if cell_list.start_pid[i, j] > 0
            cell_list.next_pid[n] = cell_list.start_pid[i, j]
        end
        cell_list.start_pid[i, j] = n
    end
end

ext: jl


text:
using DataDeps

register(DataDep(
  "home-mortgage",
  "http://www.stat.ufl.edu/~winner/data/myield.txt",
  "http://www.stat.ufl.edu/~winner/data/myield.dat",
  "b2a3205acc023b99e2e76db74973209cf9dbbf0aef3f418a48b8fc8d3d6661df",
  post_fetch_method=(path -> begin
    # IDs have spaces in them which is the delimiter
    lines = readlines(path)
    map!(x -> lstrip(x[21:end]), lines, lines)
    write(path, join(lines, '\n'))

    UCIData.process_dataset(path,
        target_index=7,
        feature_indices=1:6,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using HorizonSideRobots
include("../basics.jl")

md"""
    ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля

    РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью, следующий - весь,
                за исключением одной последней клетки на Востоке,следующий - за исключением двух последних клеток на Востоке, и т.д.
"""

"""Двигается в переданную сторону до границы и возвращает пройденное расстояние"""
function moveNcount_to_border!(r, side)
    cnt = 1
    while isborder(r, side) == 0
        move!(r, side)
        cnt += 1
    end
    return cnt
end
move_to_corner_and_get_start!(r) = (moveNcount_to_border!(r, West), moveNcount_to_border!(r, Sud))  # Двигает робота в начальный угол и возвращает координаты страта

"""Двигаться на {steps} шагов или до "упора" в сторону {side}"""
function move_for_x!(r, side::HorizonSide, steps::Int)
    for _ in 1:steps
        if isborder(r, side)
            return false
        end
        move!(r, side)
        putmarker!(r)
    end
    return true
end

"""Двигает робота на {steps} клеток вверх, возвращается назад и двигается направо на одну клетку, если может."""
function do_movements!(r, steps::Int)
    move_for_x!(r, Nord, steps)
    move_for_x!(r, Sud, steps)
    if isborder(r, Ost) == 0
        move!(r, Ost)
        return true
    else
        return false
    end
end

 """Решение задачи 4"""
 function solve_problem_4!(r)
    x, y = move_to_corner_and_get_start!(r)
    i = 0
    while do_movements!(r, i)
        i += 1
    end
    while ismarker(r)
        move!(r, West)
    end
    putmarker!(r)
    move_to!(r, x, y)
    show(r)
end


text:
using MultivariateStats
using Base.Test

srand(34568)

## PCA with zero mean

X = randn(5, 10)
Y = randn(3, 10)

P = qr(randn(5, 5))[1][:, 1:3]
pvars = [5., 4., 3.]
M = PCA(Float64[], P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == zeros(5)
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test_approx_eq transform(M, X[:,1]) P'X[:,1]
@test_approx_eq transform(M, X) P'X

@test_approx_eq reconstruct(M, Y[:,1]) P * Y[:,1]
@test_approx_eq reconstruct(M, Y) P * Y

## PCA with non-zero mean

mv = rand(5)
M = PCA(mv, P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == mv
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test_approx_eq transform(M, X[:,1]) P' * (X[:,1] .- mv)
@test_approx_eq transform(M, X) P' * (X .- mv)

@test_approx_eq reconstruct(M, Y[:,1]) P * Y[:,1] .+ mv
@test_approx_eq reconstruct(M, Y) P * Y .+ mv

## prepare training data

d = 5
n = 1000

R = qr(randn(d, d))[1]
@test_approx_eq R'R eye(5)
scale!(R, sqrt([0.5, 0.3, 0.1, 0.05, 0.05]))

X = R'randn(5, n) .+ randn(5)
mv = vec(mean(X, 2))
Z = X .- mv

C = cov(X; vardim=2)
pvs0 = sort(eigvals(C); rev=true)
tv = sum(diag(C))

## pcacov (default using cov when d < n)

M = fit(PCA, X)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mv
@test_approx_eq P'P eye(5)
@test_approx_eq C * P scale(P, pvs)
@test issorted(pvs; rev=true)
@test_approx_eq pvs pvs0
@test_approx_eq tvar(M) tv
@test_approx_eq sum(pvs) tvar(M)
@test_approx_eq reconstruct(M, transform(M, X)) X

M = fit(PCA, X; mean=mv)
@test_approx_eq projection(M) P

M = fit(PCA, Z; mean=0)
@test_approx_eq projection(M) P

M = fit(PCA, X; maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test_approx_eq P'P eye(3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test_approx_eq P'P eye(3)
@test issorted(pvs; rev=true)

## pcastd

M = fit(PCA, X; method=:svd)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mv
@test_approx_eq P'P eye(5)
@test_approx_eq_eps C * P scale(P, pvs) 1.0e-3
@test issorted(pvs; rev=true)
@test_approx_eq_eps pvs pvs0 1.0e-3
@test_approx_eq_eps tvar(M) tv 1.0e-3
@test_approx_eq sum(pvs) tvar(M)
@test_approx_eq reconstruct(M, transform(M, X)) X

M = fit(PCA, X; method=:svd, mean=mv)
@test_approx_eq projection(M) P

M = fit(PCA, Z; method=:svd, mean=0)
@test_approx_eq projection(M) P

M = fit(PCA, X; method=:svd, maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test_approx_eq P'P eye(3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; method=:svd, pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test_approx_eq P'P eye(3)
@test issorted(pvs; rev=true)


text:
"""
    simulate(bg::Number, kappa::Float64, kern::Function, maxT::Number)

Simulate a Hawkes process between 0 and `maxT` with parameters `bg`, `kappa`, `kern`.

# Arguments

* `bg`: The background rate of the Hawkes process. Constant or positive function.
* `kappa`: The kappa parameter of the Hawkes process.
* `kern` : The kernel function of the Hawkes process.
* `maxT` : Maximum time that the Hawkes process will be simulated for.

# Notes

* kappa must be between 0 and 1 for a stable Hawkes process.

# Examples

``julia
kern_f(x) = pdf.(Distributions.Exponential(1/0.5), x)

simevents = simulate(0.5, 0.5, kern_f, 100)
``

"""
function simulate(bg, kappa::Float64, kern::Function, maxT::Number)

    childfunction(x) = kappa * kern(x)
    lambda_max = childfunction(0)
    childsimulation(x) = poisson_simulation(childfunction, maxT - x, lambda_max) .+ x

    backgroundEvents = poisson_simulation(bg, maxT)

    eventbus = backgroundEvents
    allevents = backgroundEvents

    while !isempty(eventbus)
        newevents = map(childsimulation, eventbus)
        childevents = vcat(filter(!isempty, newevents)...)
        append!(allevents, childevents)
        eventbus = childevents
    end
    sort(allevents)
end

function simulate(bg, kappa::Float64, kernD::Distributions.Distribution, maxT::Number)
    kern, kern_cdf = generate_kernel_functions(kernD)
    return simulate(bg, kappa, kern, maxT)
end

function poisson_simulation(rate::Number, maxT::Number)
    nevents = rand(Distributions.Poisson(rate*maxT))
    event_times = sort(rand(Distributions.Uniform(0, maxT), nevents))
    return event_times
end

function poisson_simulation(rate::Function, maxT::Number)
    lambda_max = 1.1*maximum(rate(range(0, maxT, length=1000)))
    poisson_simulation(rate, maxT, lambda_max)
end

function poisson_simulation(rate::Function, maxT::Number, lambda_max::Number)
    tstar = poisson_simulation(lambda_max, maxT)
    accept_prob = rate(tstar) / lambda_max
    tstar[rand(length(accept_prob)) .< accept_prob]
end


text:
#Struct to hold information about wavelet-based multigrid preconditioner.
@with_kw struct Preconditioner{T <: Real, N <: Integer} <: AbstractPreconditioner{T,N}
    op::LinearMap{T}
    op_diag::Vector{T} = diag(sparse(op))
    nsmooth::N = 5
    sweep::Vector{N}
    sweep_order::Vector{N} = unique(sweep)
    smoother_omega::T = 1.0
    restrict::LinearMap{T}
    prolong::LinearMap{T}
    op_coarse::LinearMap{T} = restrict*op*prolong
    correction_coarse::Vector{T} = zeros(T,size(op_coarse,2))
    tol_coarse::T = 1e-7
    maxiter_coarse::N = 1000
end


text:
const DVISVGM_PATH = Ref{String}()
const LIBGS_PATH = Ref{String}()

function dvisvg()
    if !isassigned(DVISVGM_PATH)
        if Sys.isunix()
           DVISVGM_PATH[] = readchomp(`which dvisvgm`)
       end
    end
    return DVISVGM_PATH[]
end

function libgs()
    if !isassigned(LIBGS_PATH) && !haskey(ENV, "LIBGS")
        if Sys.isapple()
           LIBGS_PATH[] = joinpath(readchomp(`brew --cellar ghostscript`), "9.56.1", "lib", "libgs.dylib.9.56")
       end
   elseif !isassigned(LIBGS_PATH) && haskey(ENV, "LIBGS")
        LIBGS_PATH[] = ENV["LIBGS"]
    end
    return LIBGS_PATH[]
end

# The main compilation method - compiles arbitrary LaTeX documents
function compile_latex(
        document::AbstractString;
        tex_engine = `lualatex`,
        options = `-file-line-error -halt-on-error`,
        format = "dvi",
        read_format = format
    )

    use_tex_engine=tex_engine

    # First, we do some input checking.
    if !(format ∈ ("dvi", "pdf"))
        @error "Format must be either dvi or pdf; was $format"
    end
    formatcmd = ``
    if format == "dvi"
        if use_tex_engine==`lualatex`
            use_tex_engine=`dvilualatex`
        end
        # produce only dvi not pdf
        formatcmd = `-dvi -pdf-`
    else # format == "pdf"
        formatcmd = `-pdf`
    end

    # Unfortunately for us, Latexmk (which is required for any complex LaTeX doc)
    # does not like to compile straight to stdout, OR take in input from stdin;
    # it needs a file. We make a temporary directory for it to output to,
    # and create a file in there.
    return mktempdir() do dir
        cd(dir) do

            # First, create the tex file and write the document to it.
            touch("temp.tex")
            file = open("temp.tex", "w")
            print(file, document)
            close(file)

            # Now, we run the latexmk command in a pipeline so that we can redirect stdout and stderr to internal containers.
            # First we establish these pipelines:
            out = Pipe()
            err = Pipe()

            latex_cmd = `latexmk $options --shell-escape -latex=$use_tex_engine -cd -interaction=nonstopmode --output-format=$format $formatcmd temp.tex`

            latex_pipeline = pipeline(ignorestatus(latex_cmd), stdout=out, stderr=err)

            try
                latex = run(latex_pipeline)
                suc = success(latex)
                close(out.in)
                close(err.in)
                if !isfile("temp.$read_format")
                    println("Latex did not write temp.$(read_format)!")
                    println("Files in temp directory are:\n" * join(readdir(), ','))
                    printstyled("Stdout\n", bold=true)
                    println(read(out, String))
                    printstyled("Stderr\n", bold=true)
                    println(read(err, String))
                    return
                end
            finally
                if format == "pdf"

                    if pdf_num_pages("temp.pdf") > 1
                        @warn("The PDF has more than 1 page!  Choosing the first page.")
                    end

                    crop_engine = replace(string(tex_engine)[2:end-1], "la" => "")

                    pdfcrop = joinpath(@__DIR__, "pdfcrop.pl")
                    redirect_stderr(devnull) do
                        redirect_stdout(devnull) do
                            Ghostscript_jll.gs() do gs_exe
                                Perl_jll.perl() do perl_exe
                                    run(`$perl_exe $pdfcrop --margin '1 1 1 1' --$crop_engine --gscmd $gs_exe temp.pdf temp_cropped.pdf`)
                                end
                            end
                        end
                    end
                    return read("temp_cropped.pdf", String)
                else
                    return read("temp.$read_format", String)
                end
            end
        end
    end
end

compile_latex(document::TeXDocument; kwargs...) = compile_latex(convert(String, document); kwargs...)

latex2dvi(args...; kwargs...) = compile_latex(args...; format = "dvi", kwargs...)
latex2pdf(args...; kwargs...) = compile_latex(args...; format = "pdf", kwargs...)

# DVI pipeline - LaTeX → DVI → SVG → Rsvg → Cairo
function dvi2svg(
        dvi::String;
        bbox = .2, # minimal bounding box
        options = ``
    )
    # dvisvgm will allow us to convert the DVI file into an SVG which
    # can be rendered by Rsvg.  In this case, we are able to provide
    # dvisvgm a DVI file from stdin, and receive a SVG string from
    # stdout.  This greatly simplifies the pipeline, and anyone with
    # a working TeX installation should have these utilities available.

    dvisvgm_cmd = Cmd(`$(dvisvg()) --bbox=$bbox --no-fonts --stdin --stdout $options`, env = ("LIBGS" => libgs(),))

    # We create an IO buffer to redirect stderr
    err = Pipe()

    dvisvgm = open(dvisvgm_cmd, "r+")

    redirect_stderr(err) do
        write(dvisvgm, dvi)
        close(dvisvgm.in)
    end

    return read(dvisvgm.out, String) # read the SVG in as a String

end

# PDF pipeline - LaTeX → PDF → SVG → Rsvg → Cairo
# Better in general
function pdf2svg(pdf::Vector{UInt8}; kwargs...)
    pdftocairo = Poppler_jll.pdftocairo() do exe
        open(`$exe -f 1 -l 1 -svg - -`, "r+")
    end

    write(pdftocairo, pdf)

    close(pdftocairo.in)

    return read(pdftocairo.out, String)
end

pdf2svg(pdf::String) = pdf2svg(Vector{UInt8}(pdf))

# SVG/RSVG functions
# Real simple stuff

function svg2rsvg(svg::String, dpi = 72.0)
    handle = Rsvg.handle_new_from_data(svg)
    Rsvg.handle_set_dpi(handle, Float64(dpi))
    return handle
end

function rsvg2recordsurf(handle::Rsvg.RsvgHandle)
    surf = Cairo.CairoRecordingSurface()
    ctx  = Cairo.CairoContext(surf)
    Rsvg.handle_render_cairo(ctx, handle)
    return (surf, ctx)
end

function rsvg2img(handle::Rsvg.RsvgHandle, dpi = 72.0)
    Rsvg.handle_set_dpi(handle, Float64(dpi))

    # We can find the final dimensions (in pixel units) of the Rsvg image.
    # Then, it's possible to store the image in a native Julia array,
    # which simplifies the process of rendering.
    d = Rsvg.handle_get_dimensions(handle)

    # Cairo does not draw "empty" pixels, so we need to fill here
    w, h = d.width, d.height
    img = fill(Colors.ARGB32(1,1,1,0), w, h)

    # Cairo allows you to use a Matrix of ARGB32, which simplifies rendering.
    cs = Cairo.CairoImageSurface(img)
    c = Cairo.CairoContext(cs)

    # Render the parsed SVG to a Cairo context
    Rsvg.handle_render_cairo(c, handle)

    # The image is rendered transposed, so we need to flip it.
    return rotr90(permutedims(img))
end

function svg2img(svg::String, dpi = 72.0)

    # First, we instantiate an Rsvg handle, which holds a parsed representation of
    # the SVG.  Then, we set its DPI to the provided DPI (usually, 300 is good).
    handle = Rsvg.handle_new_from_data(svg)

    return rsvg2img(handle, dpi)
end

function render_surface(ctx::CairoContext, surf)
    Cairo.save(ctx)

    Cairo.set_source(ctx, surf, 0.0, 0.0)

    Cairo.paint(ctx)

    Cairo.restore(ctx)
    return
end

# Utility functions

function pdf_num_pages(filename)
    metadata = Poppler_jll.pdfinfo() do exe
        read(`$exe $filename`, String)
    end

    infos = split(metadata, '\n')

    ind = findfirst(x -> contains(x, "Pages"), infos)

    pageinfo = infos[ind]

    return parse(Int, split(pageinfo, ' ')[end])
end


text:
#==========================================================================================#

# CROSS PRODUCTS

crossprod(x::AbstractMatrix)                          = x' * x
crossprod(x::AbstractMatrix, w::UnitWeights)          = x' * x
crossprod(x::AbstractMatrix, w::AbstractVector)       = x' * (Diagonal(w) * x)
crossprod(x::AbstractMatrix, w::AbstractSparseMatrix) = x' * Matrix(w * x)
crossprod(x::AbstractMatrix, w::AbstractMatrix)       = x' * (w * x)

#==========================================================================================#

# FIRST AND SECOND STAGES

first_stage(obj::TwoStageModel)  = obj.first_stage
second_stage(obj::TwoStageModel) = obj.second_stage

switch_stage(obj::OneStageModel) = obj
switch_stage(obj::TwoStageModel) = obj.second_stage
switch_stage(obj::ParEstimate)   = obj

#==========================================================================================#

# RETRIEVAL

getvector(obj::ParModel, x::Symbol) = getvector(switch_stage(obj).sample, x)
getmatrix(obj::ParModel, args...)   = getmatrix(switch_stage(obj).sample, args...)
getnames(obj::ParModel, args...)    = getnames(switch_stage(obj).sample, args...)
getcorr(obj::ParModel)              = getcorr(switch_stage(obj).sample)
getnonmissing(obj::AnyModel)        = getnonmissing(switch_stage(obj).sample)
getweights(obj::AnyModel)           = getweights(switch_stage(obj).sample)

#==========================================================================================#

# REWEIGHTING

reweight(w::UnitWeights, v::ProbabilityWeights)     = v
reweight(w::AbstractWeights, v::ProbabilityWeights) = pweights(w .* v)

#==========================================================================================#

# FORMATTER

function frmtr(X::Array{<: Real}, d::Int)

    table = format.("{:.$(d)f}", X)

    for c = 1:size(table, 2)
        l = maximum(length.(table[:, c]))
        table[:, c] .= format.("{:>$(l)s}", table[:, c])
    end

    return table
end


text:
struct BrentMin{T}
    division::T
    detect_flatness::Bool
end
BrentMin(; division = (3 - sqrt(5)) / 2, detect_flatness = false) =
    BrentMin(division, detect_flatness)

function solve(
    problem::OptimizationProblem,
    approach::BrentMin,
    options::OptimizationOptions,
)
    _solve(problem, approach, options)
end

function _solve(prob, bm::BrentMin, options)
    time0 = time()
    a, b = bounds(prob)
    T = typeof(a)
    t = 1e-8
    c = bm.division
    v = w = x = a + c * (b - a)

    e = d = 0 * x

    f0 = fv = fw = fx = value(prob, x)
    p = q = r = 0 * x
    brent_iter = 0
    while brent_iter < options.maxiter
        brent_iter += 1
        m = (a + b) / 2
        tol = eps(T) * abs(x) + t
        if abs(x - m) > 2 * tol - (b - a) / 2 # stopping crit
            # fit parabola
            if abs(e) > tol
                r = (x - w) * (fx - fv)
                q = (x - v) * (fx - fw)
                p = (x - v) * q - (x - w) * r
                q = 2 * (q - r)
                if q > 0
                    p = -p
                else
                    q = -q
                end
                r = e
                e = d
            end

            if abs(p) < abs(q * r / 2) && p < q * (a - x) && p < q * (b - x)
                # then do the parapolic interpolation
                d = p / q
                u = x + d
                if u - a < 2 * tol || b - u < 2 * tol
                    d = sign(m-x)*tol
                end
            else
                # do golden section
                if x < m
                    e = b - x
                else
                    e = a - x
                end
                d = c * e
            end

            if abs(d) >= tol
                u = x + d
            else
                u = x + sign(d) * tol
            end

            fu = value(prob, u)
            if fu <= fx
                if u < x
                    b = x
                else
                    a = x
                end
                v = w
                fv = fw
                w = x
                fw = fx
                x = u
                fx = fu
            else
                if u < x
                    a = u
                else
                    b = u
                end
                if fu <= fw || w == x
                    v = w
                    fv = fw
                    w = u
                    fw = fu
                elseif fu <= fv || v == x || v == w
                    v = u
                    fv = fu
                end
            end
        else
            break
        end
    end

    return ConvergenceInfo(bm, (;x, f0, minimum=fx, time=time()-time0, iter=brent_iter), options)
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
module Davidson
using LinearAlgebra
export eigdav
function eigdav(A, eigs, k, kmax, tol)
    n = size(A, 1)
    V = zeros((n, n))
    theta = 0
    w = 0
    theta_old = 0
    I = Diagonal(ones(n, n))
    t = Diagonal(randn(n, n))
    for m = k:k:kmax
        if m <= k
            for j = 1:1:k
                V[:, j] = t[:, j]#/norm(t[:,j])
            end
            theta_old = ones(eigs)
        else
            theta_old = theta[1:eigs]
        end
        F = qr(V)
        V = Array{Float64}(F.Q)
        @views T = transpose(V[:, 1:(m+1)]) * A * V[:, 1:(m+1)]
        theta, s = eigen(T)

        #THETA = eigvals(T)
        #S = eigvecs(T)
        #idx = sortperm(THETA)
        #theta = THETA[idx]
        #s = S[:,idx]
        #println(theta)
        for j = 1:1:k
            @views w = (A - theta[j] * I) * V[:, 1:(m+1)] * s[:, j]
            println("DIFF ", theta[j] - A[j, j])
            q = w / (theta[j] - A[j, j])
            V[:, m + j + 1] = q
        end
        normm = norm(theta[1:eigs] - theta_old)
        println("NORM ", normm)
        if normm < tol
            return theta[1:eigs]
        end

    end
    return theta[1:eigs]
end
end #module Davidson


text:
ENV["CUDA_VISIBLE_DEVICES"] = "2,3"
using Revise
using ADSeismic
using ADCME
using MAT
using PyPlot
using DelimitedFiles
using LineSearches
using Random
using Mamba
using ProgressMeter

# matplotlib.use("Agg")
close("all")
use_gpu(0)

data_dir = "../data/acoustic"
if !ispath(data_dir)
  mkpath(data_dir)
end
figure_dir = "figure/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/"
if !ispath(result_dir)
  mkpath(result_dir)
end
# if isfile(joinpath(result_dir, "loss.txt"))
#   rm(joinpath(result_dir, "loss.txt"))
# end

################### Inversion using Automatic Differentiation #####################
model_name = "../models/marmousi2-model-smooth.mat"

## load model setting
params = load_params(model_name)

results = readdlm("Data/UQ_VP.txt")
@show size(results)
M = reshape(mean(results[1:end,:], dims = 1), params.NX+2, params.NY+2)
V = reshape(std(results[1:end,:], dims = 1), params.NX+2, params.NY+2)

close("all")
visualize_model(M, params)
savefig("VP_mean.png")

close("all")
visualize_model(V, params)
savefig("VP_std.png")


text:
"""
    handle_fault(...)

An internal function to handle a worker dying or being killed by the OS.
Attempts to determine which `Thunk`s were running on (or had their results
cached on) the dead worker, and stores them in a "deadlist". It uses this
deadlist to correct the scheduler's internal `ComputeState` struct to recover
from the fault.

Note: The logic for this functionality is not currently perfectly robust to
all failure modes, and is only really intended as a last-ditch attempt to
repair and continue executing. While it should never cause incorrect execution
of DAGs, it *may* cause a `KeyError` or other failures in the scheduler due to
the complexity of getting the internal state back to a consistent and proper
state.
"""
function handle_fault(ctx, state, deadproc)
    @assert !isempty(procs(ctx)) "No workers left for fault handling!"

    deadlist = Thunk[]

    # Evict cache entries that were stored on the worker
    for t in keys(state.cache)
        v = state.cache[t]
        if v isa Chunk && v.handle isa DRef && v.handle.owner == deadproc.pid
            push!(deadlist, t)
            pop!(state.cache, t)
        end
    end
    # Remove thunks that were running on the worker
    for t in collect(keys(state.running_on))
        pid = state.running_on[t].pid
        if pid == deadproc.pid
            push!(deadlist, t)
            delete!(state.running_on, t)
            pop!(state.running, t)
        end
    end
    # Clear thunk.cache_ref
    for t in deadlist
        t.cache_ref = nothing
    end

    # Remove thunks from state.ready that have inputs on the deadlist
    for idx in length(state.ready):-1:1
        rt = state.ready[idx]
        if any((input in deadlist) for input in rt.inputs)
            deleteat!(state.ready, idx)
        end
    end

    # Reschedule inputs from deadlist
    seen = Set{Thunk}()
    for t in deadlist
        reschedule_inputs!(state, t, seen)
    end
end


text:
export Curve
export Bernstein
export arclength

abstract type Curve{D, T} end

include("curves/bernstein.jl")


text:
module Source

import ..Points
import ..Blobs
import ..Elements: circulation, flux, kind, seed_zeros, seed_position, seed_strength,
                    seed_position!, seed_strength!

import ..Utils: dualize, ComplexGradientConfig, seed!, Partials

#== Wrapper for a point source ==#

"""
    Source.Point(z::ComplexF64, S::Float64)

A point source located at `z` with strength `S`.

A new point source can be created from an existing one by treating the
existing source as a function and passing in the parameter you
want to change as keyword arguments.
For example,
```jldoctest
julia> p = Source.Point(1.0, 1.0)
Source.Point(1.0 + 0.0im, 1.0)

julia> p()
Source.Point(1.0 + 0.0im, 1.0)

julia> p(S = 2.0)
Source.Point(1.0 + 0.0im, 2.0)
```
"""
const Point = Points.Point{T,R} where {T<:Complex, R <: Real}
Point(z::Complex{R},S::T) where {T<:Real,R<:Real} = Points.Point{Complex{T}}(z,S)
Point(z::Real,S::T) where {T} = Points.Point{Complex{T}}(complex(z),S)

(p::Point)(; z = p.z, S = imag(p.S)) = Point(z, S)

function Base.show(io::IO, s::Point)
    if iszero(real(s.S))
        print(io, "Source.Point($(s.z), $(imag(s.S)))")
    else
        print(io, "Points.Point($(s.z), $(s.S))")
    end
end
flux(p::Point) = imag(p.S)
circulation(::Point) = 0.0

#== Wrapper for a blob source ==#

"""
    Source.Blob(z::ComplexF64, S::Float64, δ::Float64)

A regularized point source located at `z` with strength `S` and blob radius `δ`.

A new blob source can be created from an existing one by treating the
existing blob as a function and passing in the parameter you want to
change as keyword arguments.
For example,
```jldoctest
julia> b = Source.Blob(1.0, 1.0, 0.1)
Source.Blob(1.0 + 0.0im, 1.0, 0.1)

julia> b()
Source.Blob(1.0 + 0.0im, 1.0, 0.1)

julia> b(S = 2.0, δ = 0.01)
Source.Blob(1.0 + 0.0im, 2.0, 0.01)
```
"""
const Blob = Blobs.Blob{T,R} where {T<:Complex, R <: Real}
Blob(z::Complex{R},S::T,δ::Float64) where {T<:Real,R<:Real} = Blobs.Blob{Complex{T}}(z,S,δ)
Blob(z::Real,S::T,δ) where {T<:Real} = Blobs.Blob{Complex{T}}(complex(z),S,δ)

(b::Blob)(; z = b.z, S = imag(b.S), δ = b.δ) = Blob(z, S, δ)

function Base.show(io::IO, s::Blob)
    if iszero(real(s.S))
        print(io, "Source.Blob($(s.z), $(imag(s.S)), $(s.δ))")
    else
        print(io, "Blobs.Blob($(s.z), $(imag(s.S)), $(s.δ))")
    end
end
circulation(::Blob) = 0.0
flux(b::Blob) = imag(b.S)

include("source/diff.jl")

end