text:
module LogRoller

using Dates
using Logging
using CodecZlib
using Logging
using JSON
using JSON.Serializations: CommonSerialization, StandardSerialization
using JSON.Writer: StructuralContext
import JSON: show_json

import Logging: shouldlog, min_enabled_level, catch_exceptions, handle_message
import Base: write, close, rawhandle
export RollingLogger, RollingFileWriter, RollingFileWriterTee, postrotate

const BUFFSIZE = 1024*16  # try and read 16K pages when possible
const DEFAULT_MAX_LOG_ENTRY_SIZE = 256*1024

include("limitio.jl")
include("log_utils.jl")

"""
A file writer that implements the `IO` interface, but only provides `write` methods.

Constructor parameters:
- filename: name (including path) of file to log into
- sizelimit: size of file (in bytes) after which the file should be rotated
- nfiles: number of rotated files to maintain
"""
mutable struct RollingFileWriter <: IO
    filename::String
    sizelimit::Int
    nfiles::Int
    filesize::Int
    stream::IO
    lck::ReentrantLock
    procstream::Union{Nothing,Pipe}
    procstreamer::Union{Nothing,Task}
    procstreamteelogger::Union{Nothing,AbstractLogger}
    assumed_level::LogLevel
    postrotate::Union{Nothing,Function}

    function RollingFileWriter(filename::String, sizelimit::Int, nfiles::Int; rotateOnInit=false)
        stream = open(filename, "a")
        filesize = stat(stream).size
        io = new(filename, sizelimit, nfiles, filesize, stream, ReentrantLock(), nothing, nothing, nothing, Logging.Info, nothing)
        if rotateOnInit && filesize > 0
            rotate_file(io)
        end
        return io
    end
end

"""
Register a function to be called with the rotated file name just after the current log file is rotated.
The file name of the rotated file is passed as an argument. The function is blocking and so any lengthy
operation that needs to be done should be done asynchronously.
"""
function postrotate(fn::Function, io::RollingFileWriter)
    io.postrotate = fn
    nothing
end

"""
Tee all lines to the provided logger
"""
function tee(io::RollingFileWriter, logger::AbstractLogger, level::LogLevel)
    io.procstreamteelogger = logger
    io.assumed_level = level
    io
end

"""
Close any open file handle and streams.
A closed object must not be used again.
"""
function close(io::RollingFileWriter)
    if io.procstream !== nothing
        close(io.procstream)
        lock(io.lck) do
            io.procstream = nothing
            io.procstreamer = nothing
            io.procstreamteelogger = nothing
        end
    end
    close(io.stream)
end

"""
Write into the underlying stream, rolling over as and when necessary.
"""
write(io::RollingFileWriter, byte::UInt8) = _write(io, byte)
write(io::RollingFileWriter, str::Union{SubString{String}, String}) = _write(io, str)
write(io::RollingFileWriter, buff::Vector{UInt8}) = _write(io, buff)

function _write(io::RollingFileWriter, args...)
    lock(io.lck) do
        io.filesize += write(io.stream, args...)
        flush(io.stream)
        if io.filesize >= io.sizelimit
            with_logger(NullLogger()) do
                rotate_file(io)
            end
        end
    end
end

"""
Rotate files as below with increasing age:
    - <filename> : active file
    - <filename>_1.gz : last rotated file
    - <filename>_2.gz : previous <filename>_1.gz rotated to <filename>_2.gz
    - <filename>_3.gz : previous <filename>_2.gz rotated to <filename>_3.gz
    - ...
    - <filename>_n.gz : last rotated file is discarded when rotated
"""
function rotate_file(io::RollingFileWriter)
    for N in io.nfiles:-1:1
        nthlogfile = string(io.filename, "_", N, ".gz")
        if isfile(nthlogfile)
            if N == io.nfiles
                rm(nthlogfile) # discard
            else
                nextlogfile = string(io.filename, "_", N+1, ".gz")
                mv(nthlogfile, nextlogfile) # rotate
            end
        end
    end
    close(io.stream)
    nthlogfile = string(io.filename, "_", 1, ".gz")
    open(nthlogfile, "w") do fw
        cfw = GzipCompressorStream(fw)
        try
            open(io.filename, "r") do fr
                buff = Vector{UInt8}(undef, BUFFSIZE)
                while !eof(fr)
                    nbytes = readbytes!(fr, buff)
                    (nbytes > 0) && write(cfw, (nbytes < BUFFSIZE) ? view(buff, 1:nbytes) : buff)
                end
            end
        finally
            flush(cfw)
            close(cfw)
        end
    end
    io.stream = open(io.filename, "w")
    io.filesize = 0

    # call postrotate hook if one is registered
    (io.postrotate === nothing) || io.postrotate(nthlogfile)

    nothing
end

"""
Tees raw log entries made a RollingFileWriter on to a provided Julia AbstractLogger.

Each line of text is taken as a single log message.

All log entries are made with the same log level, which can be provided during construction. It leaves
further examination/parsing of log messages (to extract parameters, or detect exact log levels) to the
downstream logger.
"""
function RollingFileWriterTee(filename::String, sizelimit::Int, nfiles::Int, logger::AbstractLogger, assumed_level::LogLevel=Logging.Info)
    io = RollingFileWriter(filename, sizelimit, nfiles)
    RollingFileWriterTee(io, logger, assumed_level)
end

function RollingFileWriterTee(io::RollingFileWriter, logger::AbstractLogger, assumed_level::LogLevel=Logging.Info)
    tee(io, logger, assumed_level)
end

"""
RollingLogger(filename, sizelimit, nfiles, min_level=Info; timestamp_identifier::Symbol=:time, format::Symbol=:console)
Log into a log file. Rotate log file based on file size. Compress rotated logs.

Logs can be formatted as JSON by setting the optional keyword argument `format` to `:json`. A JSON formatted log entry
is a JSON object. It should have these keys (unless they are empty):
The message part can contain the following keys unless they are empty:
- `metadata`: event metadata e.g. timestamp, line, filename, ...
- `message`: the log message string
- `keywords`: any keywords provided
"""
mutable struct RollingLogger <: AbstractLogger
    stream::RollingFileWriter
    min_level::LogLevel
    message_limits::Dict{Any,Int}
    timestamp_identifier::Symbol
    format::Symbol
    entry_size_limit::Int
end
function RollingLogger(filename::String, sizelimit::Int, nfiles::Int, level=Logging.Info; timestamp_identifier::Symbol=:time, format::Symbol=:console, entry_size_limit::Int=DEFAULT_MAX_LOG_ENTRY_SIZE, rotateOnInit=false)
    stream = RollingFileWriter(filename, sizelimit, nfiles,rotateOnInit=rotateOnInit)
    RollingLogger(stream, level, Dict{Any,Int}(), timestamp_identifier, format, entry_size_limit)
end

"""
Register a function to be called with the rotated file name just after the current log file is rotated.
The file name of the rotated file is passed as an argument. The function is blocking and so any lengthy
operation that needs to be done should be done asynchronously.
"""
postrotate(fn::Function, io::RollingLogger) = postrotate(fn, io.stream)

"""
Close any open file handle and streams.
A closed object must not be used again.
"""
close(logger::RollingLogger) = close(logger.stream)

shouldlog(logger::RollingLogger, level, _module, group, id) = get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::RollingLogger) = logger.min_level

catch_exceptions(logger::RollingLogger) = false

function get_timestamp(logger::RollingLogger, kwargs)
    try
        for (key, val) in kwargs
            if key === logger.timestamp_identifier
                if isa(val, DateTime)
                    return (val, true)
                else
                    return (Dates.unix2datetime(Float64(val)), true)
                end
            end
        end
    catch
        # could not convert val to DateTime, fallback to now()
    end
    (now(), false)
end

function handle_message(logger::RollingLogger, level, message, _module, group, id, filepath, line; maxlog=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    if logger.format === :json
        timestamp, kwarg_timestamp = get_timestamp(logger, kwargs)
        log = (level=level, message=message, _module=_module, group=group, id=id, filepath=filepath, line=line, kwargs=kwargs)
        log = merge(log, [logger.timestamp_identifier=>timestamp])
        logentry = IndexedLogEntry(log, Symbol[])
        write(logger.stream, message_string(logentry, logger.entry_size_limit, true))
    else # if logger.format === :console
        buf = IOBuffer()
        lim = LimitIO(buf, logger.entry_size_limit)
        limited = false
        try
            iob = IOContext(lim, logger.stream)
            levelstr = level == Logging.Warn ? "Warning" : string(level)
            timestamp, kwarg_timestamp = get_timestamp(logger, kwargs)
            msglines = split(chomp(string(message)), '\n')
            println(iob, "┌ ", levelstr, ": ", timestamp, ": ", msglines[1])
            for i in 2:length(msglines)
                println(iob, "│ ", msglines[i])
            end
            for (key, val) in kwargs
                kwarg_timestamp && (key === logger.timestamp_identifier) && continue
                println(iob, "│   ", key, " = ", val)
            end
            println(iob, "└ @ ", something(_module, "nothing"), " ", something(filepath, "nothing"), ":", something(line, "nothing"))
        catch ex
            isa(ex, LimitIOException) || rethrow(ex)
            limited = true
        end
        write(logger.stream, take!(buf))
        limited && write(logger.stream, UInt8[0x0a])
    end

    nothing
end

function stream_process_logs(writer::RollingFileWriter)
    try
        while true
            logline = readline(writer.procstream; keep=true)
            if !isempty(logline)
                write(writer, logline)
                if writer.procstreamteelogger !== nothing
                    @logmsg(writer.assumed_level, strip(logline))
                end
            end
            eof(writer.procstream) && break
        end
    finally
        close(writer.procstream)
        lock(writer.lck) do
            writer.procstream = nothing
            writer.procstreamer = nothing
        end
    end
end

function rawhandle(writer::RollingFileWriter)
    lock(writer.lck) do
        if (writer.procstream === nothing) || !isopen(Base.pipe_writer(writer.procstream))
            writer.procstream = Pipe()
            Base.link_pipe!(writer.procstream)
            writer.procstreamer = @async begin
                if writer.procstreamteelogger !== nothing
                    with_logger(writer.procstreamteelogger) do
                        stream_process_logs(writer)
                    end
                else
                    stream_process_logs(writer)
                end
            end
        end
        return rawhandle(Base.pipe_writer(writer.procstream))
    end
end

end # module


text:
function w(x = Float16)
    return x*x*x - 6*x*x + 3*x - 0.149
end
function w1(x = Float32)
    return x*x*x - 6*x*x + 3*x - 0.149
end
function w2(x = Float64)
    return x*x*x - 6*x*x + 3*x - 0.149
end

function z(x = FLoat16)
    return ((x-6)*x+3)*x-0.149
end
function z1(x = FLoat32)
    return ((x-6)*x+3)*x-0.149
end
function z2(x = FLoat64)
    return ((x-6)*x+3)*x-0.149
end

wy = Float16(-14.636489)
wy1 = Float32(-14.636489)
wy2 = Float64(-14.636489)

x = Float16(4.71)
x1 = Float32(4.71)
x2 = Float64(4.71)

println((w(x) - wy2) / wy2)
println((w1(x1) - wy2) / wy2)
println((w2(x2) - wy2) / wy2)

println((z(x) - wy2) / wy2)
println((z1(x1) - wy2) / wy2)
println((z2(x2) - wy2) / wy2)

ext: jl

max_stars_repo_path: II rok/I semestr/ANM/Cwiczenia 01/zad3.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
module EfficientGlobalOptimization

using KrigingModel

export KrigingTuner
export sampling, LHS, GRID
export NoFeasibleInfill
export EGO, optimum, acquire, history, lowerbounds, upperbounds
export Acquire, EGOAquisition, ProbabilityOfImprovement, ExpectedImprovement
export constrains, ConstrainedAquisition
#export ego_save
export ParetoSet

@enum EGOSense Min=1 Max=-1
const _sym2sense = Dict(Symbol(i)=>i for i in instances(EGOSense))
function sym2sense(s::Symbol)
    try
        return _sym2sense[s]
    catch
        error("invalid sense, must be :Min or :Max")
    end
end

include("sampling/sampling.jl")
include("aquisition/aquisition.jl")
include("acquire/acquire.jl")
include("krgtuner.jl")
include("EGO.jl")
include("testfunction/testfunction.jl")
include("paretoset.jl")
end # module


text:
module comparisons
"""
    Code to facilitate the comparison of models with data.
"""

using Interpolations
import DataFrames: DataFrame
import MyterialColors: salmon, blue, green, purple, teal, indigo_dark
import Statistics: mean, std, median

import jcontrol: Track, trim, interpolate_wrt_to, closest_point_idx, int, FULLTRACK

export ComparisonPoint, ComparisonPoints, TrackSegment, track_segments

# ---------------------------------------------------------------------------- #
#                               COMPARISON POINTS                              #
# ---------------------------------------------------------------------------- #

"""
Store parametersof normal distribution
"""
struct Dist
    name::Union{String, Symbol}
    μ
    med
    σ
    values::Vector{Float64}
end
Base.show(io::IO, dist::Dist) = print(io, "Dist: $(dist.name)(μ: $(round(dist.μ; digits=2)), σ: $(round(dist.σ; digits=2)))")
Dist(name, data::Vector) = Dist(name, mean(data), median(data), std(data), data)

"""
Get the number of standard deviations away from the mean
"""
σ(x::Number, dist::Dist; use=:μ) = (x - getfield(dist, use))/dist.σ

"""
Store values of kinematics variables from tracking data
"""
struct KinematicsValues
    x::Dist
    y::Dist
    θ::Dist
    ω::Dist
    u::Dist
    v::Dist
end

"""
Represents a point at which the comparison is made.

s: s value of CP
x,y: position of corresponding track point
θ: angle of corrisponding track point + π/2
w: width, used for plotting.
"""
struct ComparisonPoint
    s::Float64
    x::Float64
    y::Float64
    θ::Float64
    η::Float64
    w::Float64

    kinematics::Union{KinematicsValues, Nothing}
end

function Base.show(io::IO, p::ComparisonPoint)
    return print(
        io,
        "ComparisonPoint s=$(round(p.s; digits=1)) @ xy=($(round(p.x; digits=1)), $(round(p.y; digits=1)))",
    )
end

function ComparisonPoint(s, x, y, θ, η, w; trials::Union{Nothing, Vector}=nothing)
    isnothing(trials) && return ComparisonPoint(s, x, y, θ, η, w, nothing)

    # get kinematic variables values at comparison point
    data = Dict(
        :x => [],
        :y => [],
        :θ => [],
        :ω => [],
        :u => [],
        :v => [],
    )

    for trial in trials
        idx = closest_point_idx(trial.x, x, trial.y, y)

        for var in keys(data)
            push!(data[var], eval(:($trial.$var[$idx])))
        end
    end
    # Store mean and standard deviation
    kv = KinematicsValues(
        Dist(:x, data[:x]), Dist(:y, data[:y]), Dist(:θ, data[:θ]), Dist(:ω, data[:ω]), Dist(:u, data[:u]),  Dist(:v, data[:v])
    )
    return ComparisonPoint(s, x, y, θ, η, w, kv)
end

"""
Store a bunch of comparison points
"""
struct ComparisonPoints
    s::Vector{Float64}
    points::Vector{ComparisonPoint}
end

"""
Find values of s ∈ [0:S] at which to do comparisons.

If `mode=::equallyspaced` a bunch of s values are selected
along the entire track length, with spacing `δs`.
Otherwise hand-defined values are used.
"""
function ComparisonPoints(
    track::Track; δs=10.04, s₀=nothing, s₁=nothing, kwargs...
)::ComparisonPoints
    # get values of s for comparisons
    s₀ = isnothing(s₀) ? δs : (s₀-δs)
    s₁ = isnothing(s₁) ? track.S_f : s₁
    ŝ = s₀:δs:s₁

    ŝ = [s₀, ŝ..., s₁]

    return ComparisonPoints(track, collect(ŝ); kwargs...)
end

function ComparisonPoints(track::Track, svalues::Vector{Float64}; trials::Union{Nothing, Vector}=nothing)
    # create interpolation objects for each varianle
    x = interpolate_wrt_to(track.S, track.X)
    y = interpolate_wrt_to(track.S, track.Y)
    θ = interpolate_wrt_to(track.S, track.θ)
    η = interpolate_wrt_to(track.S, track.θ .+ π / 2)
    width = 2 # hardcoded

    # get points and return ComparisonPoints
    pts = map((_s) -> ComparisonPoint(_s, x(_s), y(_s), θ(_s), η(_s), width; trials=trials), svalues)

    return ComparisonPoints(svalues, pts)
end

# ---------------------------------------------------------------------------- #
#                                 TRACK SEGMENT                                #
# ---------------------------------------------------------------------------- #
struct TrackSegment
    track::Track  # stores a Track object with the segment track.
    s₀::Float64  # s-value of start
    s₁::Float64  # s-value of end
    checkpoints::Vector{ComparisonPoint}
    color::String
end

function TrackSegment(s₀::Float64, s₁::Float64, color::String; δs=5)
    # theres 262 waypoints in the track's npy file, the track's S_f is 259.
    # get track
    track = trim(FULLTRACK, s₀ * 259, (s₁ - s₀) * 259)
    # rack(; start_waypoint=wp₀, keep_n_waypoints=nwp)

    return TrackSegment(
        track,
        s₀,
        s₁,
        ComparisonPoints(FULLTRACK; δs=δs, s₀=(s₀+.01)*FULLTRACK.S_f, s₁=(s₁-.01)*FULLTRACK.S_f).points,
        color
    )
end

track_segments = [
    TrackSegment(0.15, .30, salmon),
    TrackSegment(0.24, .43, blue),
    TrackSegment(0.40, .61, green),
    TrackSegment(0.55, .72, purple),
    TrackSegment(0.71, .86, indigo_dark),
    TrackSegment(0.82, .975, teal),
]

end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
"""
```
    LinearStretching <: AbstractHistogramAdjustmentAlgorithm
    LinearStretching(; [src_minval], [src_maxval],
                       dst_minval=0, dst_maxval=1,
                       no_clamp=false)

    LinearStretching((src_minval, src_maxval) => (dst_minval, dst_maxval))
    LinearStretching((src_minval, src_maxval) => nothing)
    LinearStretching(nothing => (dst_minval, dst_maxval))

    adjust_histogram([T,] img, f::LinearStretching)
    adjust_histogram!([out,] img, f::LinearStretching)
```

Returns an image where the range of the intensities spans the interval [`dst_minval`, `dst_maxval`].

# Details

Linear stretching (also called *normalization*) is a contrast enhancing
transformation that is used to modify the dynamic range of the image. In
particular, suppose that the input image has gray values in the range [A,B] and
one wishes to change the dynamic range to [a,b] using a linear mapping, then the
necessary transformation is given by the relation

```math
f(x) = (x-A) \\frac{b-a}{B-A} + a.
```

# Options

Various options for the parameters of the `adjust_histogram` and
`LinearStretching` type  are described in more detail below.

## Choices for `img`

The function can handle a variety of input types. The returned
image depends on the input type.

For colored images, the input is converted to the
[YIQ](https://en.wikipedia.org/wiki/YIQ)  type and the intensities of the Y
channel are stretched to the specified range. The modified Y channel is then
combined with the I and Q channels and the resulting image converted to the same
type as the input.

## Choices for `dst_minval` and `dst_maxval`

If destination value range `dst_minval` and `dst_maxval` are specified then intensities are
mapped to the range [`dst_minval`, `dst_maxval`]. The default values are 0 and 1.

## Choices for `src_minval` and `src_maxval`

The source value range `src_minval` and `src_maxval` specifies the intensity range of input
image. By default, the values are `extrema(img)` (finite). If custom values are provided,
the output intensity value will be clamped to range `[dst_minval, dst_maxval]` if it exceeds that.

## `no_clamp`

Setting `no_clamp=true` to disable the automatic clamp even if the output intensity value
exceeds the range `[dst_minval, dst_maxval]`. Note that a clamp is still applied for types
that has limited value range, for example, if the input eltype is `N0f8`, then the output will
be clamped to `[0.0N0f8, 1.0N0f8]` even if `no_clamp==true`.

# Example

```julia
using ImageContrastAdjustment, TestImages

img = testimage("mandril_gray")
# Stretches the contrast in `img` so that it spans the unit interval.
imgo = adjust_histogram(img, LinearStretching(dst_minval = 0, dst_maxval = 1))
```

For convenience, Constructing a `LinearStretching` object using `Pair` is also supported

```julia
# these two constructors are equivalent
LinearStretching(src_minval=0.1, src_maxval=0.9, dst_minval=0.05, dst_maxval=0.95)
LinearStretching((0.1, 0.9) => (0.05, 0.95))

# replace the part with `nothing` to use default values, e.g.,
# specify only destination value range
LinearStretching(nothing => (0.05, 0.95))
# specify only source value range and use default destination value range, i.e., (0, 1)
LinearStretching((0.1, 0.9) => nothing)
```

# References
1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)

"""
@with_kw struct LinearStretching{T} <: AbstractHistogramAdjustmentAlgorithm
    src_minval::T = nothing
    src_maxval::T = nothing
    dst_minval::T = 0.0f0
    dst_maxval::T = 1.0f0
    minval::T = nothing
    maxval::T = nothing
    no_clamp::Bool = false
    function LinearStretching(src_minval::T1,
                              src_maxval::T2,
                              dst_minval::T3,
                              dst_maxval::T4,
                              minval::T5=nothing,
                              maxval::T6=nothing,
                              no_clamp::Bool=false) where {T1 <: Union{Nothing,Real,AbstractGray},
                                                 T2 <: Union{Nothing,Real,AbstractGray},
                                                 T3 <: Union{Nothing,Real,AbstractGray},
                                                 T4 <: Union{Nothing,Real,AbstractGray},
                                                 T5 <: Union{Nothing,Real,AbstractGray},
                                                 T6 <: Union{Nothing,Real,AbstractGray}}
        # in order to deprecate old fields we have to introduce new fields if we still want to use @with_kw
        # https://github.com/JuliaImages/ImageContrastAdjustment.jl/pull/28#discussion_r395751301
        if !isnothing(minval)
            dst_minval = minval
            Base.depwarn("deprecated: use `dst_minval` for keyword `minval`", :LinearStretching)
        end
        if !isnothing(maxval)
            dst_maxval = maxval
            Base.depwarn("deprecated: use `dst_maxval` for keyword `maxval`", :LinearStretching)
        end

        dst_minval <= dst_maxval || throw(ArgumentError("dst_minval $dst_minval should be less than dst_maxval $dst_maxval"))
        if !(isnothing(src_minval) || isnothing(src_maxval))
            src_minval <= src_maxval || throw(ArgumentError("src_minval $src_minval should be less than src_maxval $src_maxval"))
        end
        T = promote_type(T1, T2, T3, T4, T5, T6)
        new{T}(convert(T, src_minval), convert(T, src_maxval),
               convert(T, dst_minval), convert(T, dst_maxval),
               convert(T, dst_minval), convert(T, dst_maxval),
               no_clamp)
    end
end
function LinearStretching(rangemap::Pair{Tuple{T1, T2}, Tuple{T3, T4}}; no_clamp=false) where {T1, T2, T3, T4}
    LinearStretching(rangemap.first..., rangemap.second..., nothing, nothing, no_clamp)
end
function LinearStretching(rangemap::Pair{Nothing, Tuple{T3, T4}}; no_clamp=false) where {T3, T4}
    LinearStretching(nothing, nothing, rangemap.second..., nothing, nothing, no_clamp)
end
function LinearStretching(rangemap::Pair{Tuple{T1, T2}, Nothing}; no_clamp=false) where {T1, T2}
    LinearStretching(src_minval=rangemap.first[1], src_maxval=rangemap.first[2], no_clamp=no_clamp)
end

function (f::LinearStretching)(out::GenericGrayImage, img::GenericGrayImage)
    T = eltype(out)
    FT = eltype(floattype(T))
    img_min, img_max = minfinite(img), maxfinite(img)
    # explicit annotation is needed because the ?: line mixes three value types:
    # Nothing, T, and typeof(f.src_minval)
    src_minval::FT = isnothing(f.src_minval) ? img_min : f.src_minval
    src_maxval::FT = isnothing(f.src_maxval) ? img_max : f.src_maxval
    dst_minval::FT = f.dst_minval
    dst_maxval::FT = f.dst_maxval

    # the kernel operation `r * x - o` is equivalent to `(x-A) * ((b-a)/(B-A)) + a`
    # precalculate these and make inner loop contains only multiplication and addition
    # to get better performance
    r = (dst_maxval - dst_minval) / (src_maxval - src_minval)
    o = (src_minval*dst_maxval - src_maxval*dst_minval) / (src_maxval - src_minval)

    if 1 ≈ r && 0 ≈ o
        # when image intensity is already adjusted, there's no need to do it again
        # it's a trivial but common case in practice
        out === img || (out .= img)
        return out
    end

    # In most cases, we don't need to clamp the output
    # this is only used when user specifies custom parameters
    out_minval = r * img_min - o
    out_maxval = r * img_max - o
    do_clamp = f.no_clamp ? false : (out_minval < dst_minval) || (out_maxval > dst_maxval)
    # re-adjust clamp option to avoid ArgumentError
    if !do_clamp && (out_minval < typemin(T) || out_maxval > typemax(T))
        do_clamp = true
        dst_minval = convert(typeof(dst_minval), typemin(eltype(T)))
        dst_maxval = convert(typeof(dst_maxval), typemax(eltype(T)))
    else
        do_clamp = do_clamp || out_minval < typemin(T) || out_maxval > typemax(T)
        dst_minval = max(dst_minval, convert(typeof(dst_minval), typemin(eltype(T))))
        dst_maxval = min(dst_maxval, convert(typeof(dst_maxval), typemax(eltype(T))))
    end

    # although slightly faster, clamp before linear stretching has some rounding issue, e.g,
    # it's likely to get results like -9.313226f-10, which can't be directly converted to N0f8
    # thus the line below is commented out and not used
    # do_clamp && (img = clamp.(img, src_minval, src_maxval))

    # tweak the performance of FixedPoint by fusing operations into one broadcast
    # for Float32 the fallback implementation is faster
    if eltype(T) <: FixedPoint
        # ?: is faster than if-else
        @. out = do_clamp ? clamp(r * img - o, dst_minval, dst_maxval) : r * img - o
        return out
    end

    # fallback implementation
    @inbounds @simd for p in eachindex(img)
        val = img[p]
        if isnan(val)
            out[p] = val
        else
            newval = r * val - o
            do_clamp && (newval = clamp(newval, dst_minval, dst_maxval))
            out[p] = T <: Integer ? round(Int, newval) : newval
        end
    end
    out
end

function (f::LinearStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(out)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)
    adjust_histogram!(view(yiq_view,1,:,:), f)
    out .= convert.(T, yiq)
end

(f::LinearStretching)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
"""
    fitsprintr(X, y)

Input: data X, outcome y
Return m, model fit by Sprint.r
typeof(m): OrderedCollections.OrderedDict{Symbol, Any}
"""
function fitsprintr(X, y;
    #nfold = 5, square = FALSE, nlam = 100, lam_min_ratio = 0.01
    )
    @rput y X # put data from Julia to R.
    R"""
    library(sprintr);
    m <- cv.sprinter(x = X, y = y,
        nfold = 5, square = FALSE, nlam = 100, lam_min_ratio = 0.01)
    """
    @rget m
    delete!(m, :call)
    delete!(m[:fit], :call)
    m
end

"""
    predictsprintr(m, XH)

Input: fit object m, new data XH
Return pr, vector of predictions
"""
function predictsprintr(m, XH)
    @rput m XH # put data from Julia to R.
    R"""
    library(sprintr);
    class(m) <- 'cv.sprinter'
    pred <- predict(m, newdata = XH)
    """
    return @rget pred
end


text:
# ---
# title: 1239. Maximum Length of a Concatenated String with Unique Characters
# id: problem1239
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking, Bit Manipulation
# link: <https://leetcode.com/problems/maximum-length-of-a-concatenated-string-with-unique-characters/description/>
# hidden: true
# ---
#
# Given an array of strings `arr`. String `s` is a concatenation of a sub-
# sequence of `arr` which have **unique characters**.
#
# Return _the maximum possible length_ of `s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = ["un","iq","ue"]
#     Output: 4
#     Explanation: All possible concatenations are "","un","iq","ue","uniq" and "ique".
#     Maximum length is 4.
#
#
# **Example 2:**
#
#
#
#     Input: arr = ["cha","r","act","ers"]
#     Output: 6
#     Explanation: Possible solutions are "chaers" and "acters".
#
#
# **Example 3:**
#
#
#
#     Input: arr = ["abcdefghijklmnopqrstuvwxyz"]
#     Output: 26
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 16`
#   * `1 <= arr[i].length <= 26`
#   * `arr[i]` contains only lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "add_Gb" begin
    Gb_h = 0.0347
    @test add_Gb(Gb_h) == NamedTuple()
    #
    Gb = add_Gb(Gb_h, :N20 => 2, :CH4 => 4)
    @test keys(Gb) == (:Gb_N20, :Gb_CH4)
    @test ≈(Gb.Gb_N20, 0.0173, rtol = 1/100)
    @test ≈(Gb.Gb_CH4, 0.0109, rtol = 1/100)
    #
    Gbm = add_Gb(missing, :N20 => 2, :CH4 => 4)
    @test keys(Gbm) == (:Gb_N20, :Gb_CH4)
    @test all(ismissing.(values(Gbm)))
    #
    # DataFrame variant
    dfo = DataFrame(Gb_h=[Gb_h, missing, 0.055])
    df = copy(dfo)
    df2 = @inferred add_Gb!(df)
    @test isequal(df2, dfo)
    df2 = @inferred add_Gb!(df, :N20 => 2, :CH4 => 4)
    @test df2 === df # mutating
    @test propertynames(df2)[end-1:end] == [:Gb_N20, :Gb_CH4]
    @test df2.Gb_N20[1] == Gb.Gb_N20
end

@testset "compute_Gb Gb_constant_kB1 and compute_Gb_quantities" begin
    kB_h = 1.18
    ustar = 0.1
    Gb = @inferred Gb_constant_kB1(ustar, kB_h)
    @test Gb ≈ 0.0347 rtol = 1e-2
    Gbq = compute_Gb_quantities(Gb, ustar)
    @test keys(Gbq) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
    @test all(isapprox.(values(Gbq), values(
        (Rb_h = 28.8, Gb_h = 0.0347, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
    Gb = @inferred Gb_constant_kB1(missing, kB_h)
    @test ismissing(Gb)
    Gbq = compute_Gb_quantities(Gb, ustar)
    @test keys(Gbq) == (:Rb_h, :Gb_h, :kB_h, :Gb_CO2)
    @test all(ismissing.(values(getindex.(Ref(Gbq), SA[:Rb_h, :Gb_h, :Gb_CO2]))))
    #
    # DataFrame variant
    #dfo = DataFrame(ustar = SA[0.1,missing,0.3]) # not use SVector in DataFrame
    dfo = DataFrame(ustar = [0.1,missing,0.3])
    df = copy(dfo)
    @inferred compute_Gb!(df, ConstantKB1(); kB_h)
    @test propertynames(df) == [:ustar, :Gb_h]
    compute_Gb_quantities!(df)
    @test propertynames(df) == [:ustar, :Gb_h, :Rb_h, :kB_h, :Gb_CO2]
    @test all(isapprox.(values(df[1,2:end]), values(
        (Gb_h = 0.0347, Rb_h = 28.8, kB_h = 1.18, Gb_CO2 = 0.0264)), rtol = 1e-2))
end

@testset "compute_Gb Gb_Thom" begin
    ustar = 0.1
    Gb1 = @inferred Gb_Thom(ustar)
    @test Gb1 ≈ 0.0347 rtol = 1e-2
    Gbm = Gb_Thom(missing)
    @test ismissing(Gbm)
    #
    # DataFrame variant
    dfo = DataFrame(ustar = [ustar,missing,0.3])
    df = copy(dfo)
    @inferred compute_Gb!(df, Thom1972())
    @test propertynames(df) == [:ustar, :Gb_h]
    @test df.Gb_h[1] == Gb1
end

@testset "compute_Gb Gb_Choudhury and Gb_Su" begin
    zh, zr, LAI = thal.zh, thal.zr, thal.LAI
    leafwidth=0.1
    wind_zh = 2.1662688
    Gb_Choud = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh)
    @test Gb_Choud ≈ 0.157 rtol=1e-2 # from R
    Gbm = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh=missing)
    @test ismissing(Gbm)
    #
    Tair, pressure, ustar = values(tha48[1, SA[:Tair, :pressure, :ustar]])
    Dl=0.01; fc = (1-exp(-LAI/2))
    Gb_S = @inferred Gb_Su(Tair, pressure, ustar; wind_zh, Dl, fc)
    @test Gb_S ≈ 0.185 rtol=1e-2 # from R
    Gbm = @inferred Gb_Su(missing, pressure, ustar; wind_zh, Dl, fc)
    @test ismissing(Gbm)
    #
    # DataFrame variant
    df = tha48[:,Not(:Gb_h)]
    # sum(skipmissing(...)) not inferrable in Julia 1.6
    # @descend_code_warntype roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)
    # @code_warntype roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)
    # z0m = (@inferred roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)).z0m
    z0m = roughness_parameters(
        Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
        zh, zr).z0m
    # function f1(zh, ustar, z0m, Tair, pressure, H)
    #     wind_zh = wind_profile.(zh, ustar, 0.7*zh, z0m, Tair, pressure, H)
    #     wind_zh * 2
    # end
    # @code_warntype f(zh, df.ustar, z0m,df.Tair, df.pressure, df.H)
    wind_zh = wind_profile.(zh, df.ustar, 0.7*zh, z0m, df.Tair, df.pressure, df.H)
    @inferred compute_Gb!(df, Choudhury1988(); leafwidth, LAI, wind_zh)
    @test last(propertynames(df)) == :Gb_h
    @test df.Gb_h[1] ≈ Gb_Choud rtol=1e-6
    #
    df = tha48[:,Not(:Gb_h)]
    @inferred compute_Gb!(df, Su2001(); wind_zh, Dl, LAI)
    @test last(propertynames(df)) == :Gb_h
    @test df.Gb_h[1] ≈ Gb_S rtol=1e-6
end

@testset "compute_Gb Gb_Choudhury and Gb_Su Float32" begin
    zh, zr, LAI = thal.zh, thal.zr, thal.LAI
    leafwidth=0.1
    wind_zh = 2.1662688f0
    Gb_Choud = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh)
    @test typeof(Gb_Choud) == Float32
    @test Gb_Choud ≈ 0.157 rtol=1e-2 # from R
    Gbm = @inferred Gb_Choudhury(; leafwidth, LAI, wind_zh=missing)
    @test ismissing(Gbm)
    #
    Tair, pressure, ustar = Float32.(values(tha48[1, SA[:Tair, :pressure, :ustar]]))
    Dl=0.01; fc = (1-exp(-LAI/2))
    Gb_S = @inferred Gb_Su(Tair, pressure, ustar; wind_zh, Dl, fc)
    @test typeof(Gb_S) == Float32
    @test Gb_S ≈ 0.185 rtol=1e-2 # from R
    Gbm = @inferred Gb_Su(missing, pressure, ustar; wind_zh, Dl, fc)
    @test ismissing(Gbm)
    Gbm = @inferred Gb_Su(Tair, pressure, missing; wind_zh, Dl, fc)
    @test ismissing(Gbm)
    #
    # DataFrame variant
    df0 = tha48[:,Not(:Gb_h)]
    float_cols = names(df0, Not(:datetime))
    df0 = transform(df0, float_cols .=> ByRow(passmissing(Float32)) .=> float_cols)
    df = copy(df0)
    # sum(skipmissing(...)) not inferrable in Julia 1.6
    # @descend_code_warntype roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)
    # @code_warntype roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)
    # z0m = (@inferred roughness_parameters(
    #     Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
    #     zh, zr)).z0m
    z0m = roughness_parameters(
        Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H;
        zh, zr).z0m
    @test typeof(z0m) == Float32
    # function f1(zh, ustar, z0m, Tair, pressure, H)
    #     wind_zh = wind_profile.(zh, ustar, 0.7*zh, z0m, Tair, pressure, H)
    #     wind_zh * 2
    # end
    # @code_warntype f(zh, df.ustar, z0m,df.Tair, df.pressure, df.H)
    wind_zh = wind_profile.(zh, df.ustar, 0.7*zh, z0m, df.Tair, df.pressure, df.H)
    @test eltype(wind_zh) == Float32
    @inferred compute_Gb!(df, Choudhury1988(); leafwidth, LAI, wind_zh)
    @test last(propertynames(df)) == :Gb_h
    @test eltype(df.Gb_h) == Float32
    @test df.Gb_h[1] ≈ Gb_Choud rtol=1e-6
    #
    df = copy(df0)
    @inferred compute_Gb!(df, Su2001(); wind_zh, Dl, LAI)
    @test last(propertynames(df)) == :Gb_h
    @test df.Gb_h[1] ≈ Gb_S rtol=1e-6
    @test eltype(df.Gb_h) == Float32
end


text:
module BLS12381FieldProfile

using Profile

using TowerFields.BLS12381Field: powx, Q
using ..TestUtils: randfq12

const PROFILE_TIME = 0.5
const MINCOUNT = parse(Int, get(ENV, "MINCOUNT", "10"))
const FORMAT = Symbol(get(ENV, "FORMAT", "tree"))

# dynamically load the tested functions
function runprofiler(func::Function, args...)
    func(args...) # compile
    local stoptime = time() + PROFILE_TIME
    Profile.clear()
    while time() < stoptime
        @profile func(args...)
    end
    Profile.print(format=FORMAT, sortedby=:count, mincount=MINCOUNT)
end

const A2 = (Q - 1, Q - 2)
const B2 = (Q - 2, Q - 3)
const C2 = (Q - 3, Q - 4)
const A6 = (A2, B2, C2)
const B6 = (C2, A2, B2)
const A = (A6, B6)

# Current bottleneck is in the Final Exponention (as expected).
# Search for "requires 5 exponentiations by z" in:
# https://eprint.iacr.org/2012/232.pdf
runprofiler(powx, A)

end


text:
using DiffEqBase
using Test

@test DiffEqBase.undefined_exports(DiffEqBase) == []


text:
using Test
using LanguageServer

@testset "URI conversion" begin

if Sys.iswindows()
    @test LanguageServer.uri2filepath("file://SERVER/foo/bar") == "\\\\SERVER\\foo\\bar"
    @test LanguageServer.uri2filepath("file://wsl%24/foo/bar") == "\\\\wsl\$\\foo\\bar"
    @test LanguageServer.uri2filepath("file:///D:/foo/bar") == "d:\\foo\\bar"
    @test LanguageServer.uri2filepath("file:///foo/bar") == "\\foo\\bar"

    @test LanguageServer.filepath2uri("\\\\SERVER\\foo\\bar") == "file://SERVER/foo/bar"
    @test LanguageServer.filepath2uri("\\\\wsl\$\\foo\\bar") == "file://wsl%24/foo/bar"
    @test LanguageServer.filepath2uri("d:\\foo\\bar") == "file:///d%3A/foo/bar"
else
    @test LanguageServer.uri2filepath("file:///foo/bar") == "/foo/bar"

    @test LanguageServer.filepath2uri("/foo/bar") == "file:///foo/bar"
end

end

@testset "URI comparison" begin

if Sys.iswindows()
    @test LanguageServer.escape_uri("file:///c:/foo/bar") == "file:///c%3A/foo/bar"
    @test LanguageServer.escape_uri("file://wsl%24/foo/bar") == "file://wsl%24/foo/bar"
    @test LanguageServer.escape_uri("file:///D:/FOO/bar") == "file:///D%3A/FOO/bar"

    @test hash(LanguageServer.URI2("file:///D:/FOO/bar")) == hash(LanguageServer.URI2("file:///d%3A/FOO/bar"))
    @test hash(LanguageServer.URI2("file:///C:/foo/space bar")) == hash(LanguageServer.URI2("file:///c%3A/foo/space%20bar"))
    @test hash(LanguageServer.URI2("file://wsl\$/foo/bar")) == hash(LanguageServer.URI2("file://wsl%24/foo/bar"))

    @test LanguageServer.URI2("file:///D:/FOO/bar") == LanguageServer.URI2("file:///d%3A/FOO/bar")
    @test LanguageServer.URI2("file:///C:/foo/space bar") == LanguageServer.URI2("file:///c%3A/foo/space%20bar")
    @test LanguageServer.URI2("file://wsl\$/foo/bar") == LanguageServer.URI2("file://wsl%24/foo/bar")
else
    @test LanguageServer.escape_uri("file:///foo/bar") == "file:///foo/bar"

    @test hash(LanguageServer.URI2("file:///foo/bar")) == hash(LanguageServer.URI2("file:///foo/bar"))
    @test hash(LanguageServer.URI2("file:///foo/space bar")) == hash(LanguageServer.URI2("file:///foo/space%20bar"))

    @test LanguageServer.URI2("file:///foo/bar") == LanguageServer.URI2("file:///foo/bar")
    @test LanguageServer.URI2("file:///foo/space bar") == LanguageServer.URI2("file:///foo/space%20bar")
end

end


text:
# Julia wrapper for header: /usr/include/expat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function XML_SetElementDeclHandler(parser::XML_Parser,eldecl::XML_ElementDeclHandler)
    ccall((:XML_SetElementDeclHandler,libexpat),Void,(XML_Parser,XML_ElementDeclHandler),parser,eldecl)
end

function XML_SetAttlistDeclHandler(parser::XML_Parser,attdecl::XML_AttlistDeclHandler)
    ccall((:XML_SetAttlistDeclHandler,libexpat),Void,(XML_Parser,XML_AttlistDeclHandler),parser,attdecl)
end

function XML_SetXmlDeclHandler(parser::XML_Parser,xmldecl::XML_XmlDeclHandler)
    ccall((:XML_SetXmlDeclHandler,libexpat),Void,(XML_Parser,XML_XmlDeclHandler),parser,xmldecl)
end

function XML_ParserCreate(encoding)
    ccall((:XML_ParserCreate,libexpat),XML_Parser,(Ptr{XML_Char},),encoding)
end

function XML_ParserCreateNS(encoding,namespaceSeparator::XML_Char)
    ccall((:XML_ParserCreateNS,libexpat),XML_Parser,(Ptr{XML_Char},XML_Char),encoding,namespaceSeparator)
end

function XML_ParserCreate_MM(encoding,memsuite,namespaceSeparator)
    ccall((:XML_ParserCreate_MM,libexpat),XML_Parser,(Ptr{XML_Char},Ptr{XML_Memory_Handling_Suite},Ptr{XML_Char}),encoding,memsuite,namespaceSeparator)
end

function XML_ParserReset(parser::XML_Parser,encoding)
    ccall((:XML_ParserReset,libexpat),XML_Bool,(XML_Parser,Ptr{XML_Char}),parser,encoding)
end

function XML_SetEntityDeclHandler(parser::XML_Parser,handler::XML_EntityDeclHandler)
    ccall((:XML_SetEntityDeclHandler,libexpat),Void,(XML_Parser,XML_EntityDeclHandler),parser,handler)
end

function XML_SetElementHandler(parser::XML_Parser,start::XML_StartElementHandler,_end::XML_EndElementHandler)
    ccall((:XML_SetElementHandler,libexpat),Void,(XML_Parser,XML_StartElementHandler,XML_EndElementHandler),parser,start,_end)
end

function XML_SetStartElementHandler(parser::XML_Parser,handler::XML_StartElementHandler)
    ccall((:XML_SetStartElementHandler,libexpat),Void,(XML_Parser,XML_StartElementHandler),parser,handler)
end

function XML_SetEndElementHandler(parser::XML_Parser,handler::XML_EndElementHandler)
    ccall((:XML_SetEndElementHandler,libexpat),Void,(XML_Parser,XML_EndElementHandler),parser,handler)
end

function XML_SetCharacterDataHandler(parser::XML_Parser,handler::XML_CharacterDataHandler)
    ccall((:XML_SetCharacterDataHandler,libexpat),Void,(XML_Parser,XML_CharacterDataHandler),parser,handler)
end

function XML_SetProcessingInstructionHandler(parser::XML_Parser,handler::XML_ProcessingInstructionHandler)
    ccall((:XML_SetProcessingInstructionHandler,libexpat),Void,(XML_Parser,XML_ProcessingInstructionHandler),parser,handler)
end

function XML_SetCommentHandler(parser::XML_Parser,handler::XML_CommentHandler)
    ccall((:XML_SetCommentHandler,libexpat),Void,(XML_Parser,XML_CommentHandler),parser,handler)
end

function XML_SetCdataSectionHandler(parser::XML_Parser,start::XML_StartCdataSectionHandler,_end::XML_EndCdataSectionHandler)
    ccall((:XML_SetCdataSectionHandler,libexpat),Void,(XML_Parser,XML_StartCdataSectionHandler,XML_EndCdataSectionHandler),parser,start,_end)
end

function XML_SetStartCdataSectionHandler(parser::XML_Parser,start::XML_StartCdataSectionHandler)
    ccall((:XML_SetStartCdataSectionHandler,libexpat),Void,(XML_Parser,XML_StartCdataSectionHandler),parser,start)
end

function XML_SetEndCdataSectionHandler(parser::XML_Parser,_end::XML_EndCdataSectionHandler)
    ccall((:XML_SetEndCdataSectionHandler,libexpat),Void,(XML_Parser,XML_EndCdataSectionHandler),parser,_end)
end

function XML_SetDefaultHandler(parser::XML_Parser,handler::XML_DefaultHandler)
    ccall((:XML_SetDefaultHandler,libexpat),Void,(XML_Parser,XML_DefaultHandler),parser,handler)
end

function XML_SetDefaultHandlerExpand(parser::XML_Parser,handler::XML_DefaultHandler)
    ccall((:XML_SetDefaultHandlerExpand,libexpat),Void,(XML_Parser,XML_DefaultHandler),parser,handler)
end

function XML_SetDoctypeDeclHandler(parser::XML_Parser,start::XML_StartDoctypeDeclHandler,_end::XML_EndDoctypeDeclHandler)
    ccall((:XML_SetDoctypeDeclHandler,libexpat),Void,(XML_Parser,XML_StartDoctypeDeclHandler,XML_EndDoctypeDeclHandler),parser,start,_end)
end

function XML_SetStartDoctypeDeclHandler(parser::XML_Parser,start::XML_StartDoctypeDeclHandler)
    ccall((:XML_SetStartDoctypeDeclHandler,libexpat),Void,(XML_Parser,XML_StartDoctypeDeclHandler),parser,start)
end

function XML_SetEndDoctypeDeclHandler(parser::XML_Parser,_end::XML_EndDoctypeDeclHandler)
    ccall((:XML_SetEndDoctypeDeclHandler,libexpat),Void,(XML_Parser,XML_EndDoctypeDeclHandler),parser,_end)
end

function XML_SetUnparsedEntityDeclHandler(parser::XML_Parser,handler::XML_UnparsedEntityDeclHandler)
    ccall((:XML_SetUnparsedEntityDeclHandler,libexpat),Void,(XML_Parser,XML_UnparsedEntityDeclHandler),parser,handler)
end

function XML_SetNotationDeclHandler(parser::XML_Parser,handler::XML_NotationDeclHandler)
    ccall((:XML_SetNotationDeclHandler,libexpat),Void,(XML_Parser,XML_NotationDeclHandler),parser,handler)
end

function XML_SetNamespaceDeclHandler(parser::XML_Parser,start::XML_StartNamespaceDeclHandler,_end::XML_EndNamespaceDeclHandler)
    ccall((:XML_SetNamespaceDeclHandler,libexpat),Void,(XML_Parser,XML_StartNamespaceDeclHandler,XML_EndNamespaceDeclHandler),parser,start,_end)
end

function XML_SetStartNamespaceDeclHandler(parser::XML_Parser,start::XML_StartNamespaceDeclHandler)
    ccall((:XML_SetStartNamespaceDeclHandler,libexpat),Void,(XML_Parser,XML_StartNamespaceDeclHandler),parser,start)
end

function XML_SetEndNamespaceDeclHandler(parser::XML_Parser,_end::XML_EndNamespaceDeclHandler)
    ccall((:XML_SetEndNamespaceDeclHandler,libexpat),Void,(XML_Parser,XML_EndNamespaceDeclHandler),parser,_end)
end

function XML_SetNotStandaloneHandler(parser::XML_Parser,handler::XML_NotStandaloneHandler)
    ccall((:XML_SetNotStandaloneHandler,libexpat),Void,(XML_Parser,XML_NotStandaloneHandler),parser,handler)
end

function XML_SetExternalEntityRefHandler(parser::XML_Parser,handler::XML_ExternalEntityRefHandler)
    ccall((:XML_SetExternalEntityRefHandler,libexpat),Void,(XML_Parser,XML_ExternalEntityRefHandler),parser,handler)
end

function XML_SetExternalEntityRefHandlerArg(parser::XML_Parser,arg)
    ccall((:XML_SetExternalEntityRefHandlerArg,libexpat),Void,(XML_Parser,Ptr{Void}),parser,arg)
end

function XML_SetSkippedEntityHandler(parser::XML_Parser,handler::XML_SkippedEntityHandler)
    ccall((:XML_SetSkippedEntityHandler,libexpat),Void,(XML_Parser,XML_SkippedEntityHandler),parser,handler)
end

function XML_SetUnknownEncodingHandler(parser::XML_Parser,handler::XML_UnknownEncodingHandler,encodingHandlerData)
    ccall((:XML_SetUnknownEncodingHandler,libexpat),Void,(XML_Parser,XML_UnknownEncodingHandler,Ptr{Void}),parser,handler,encodingHandlerData)
end

function XML_DefaultCurrent(parser::XML_Parser)
    ccall((:XML_DefaultCurrent,libexpat),Void,(XML_Parser,),parser)
end

function XML_SetReturnNSTriplet(parser::XML_Parser,do_nst::Cint)
    ccall((:XML_SetReturnNSTriplet,libexpat),Void,(XML_Parser,Cint),parser,do_nst)
end

function XML_SetUserData(parser::XML_Parser,userData)
    ccall((:XML_SetUserData,libexpat),Void,(XML_Parser,Ptr{Void}),parser,userData)
end

function XML_SetEncoding(parser::XML_Parser,encoding)
    ccall((:XML_SetEncoding,libexpat),Cint,(XML_Parser,Ptr{XML_Char}),parser,encoding)
end

function XML_UseParserAsHandlerArg(parser::XML_Parser)
    ccall((:XML_UseParserAsHandlerArg,libexpat),Void,(XML_Parser,),parser)
end

function XML_UseForeignDTD(parser::XML_Parser,useDTD::XML_Bool)
    ccall((:XML_UseForeignDTD,libexpat),Cint,(XML_Parser,XML_Bool),parser,useDTD)
end

function XML_SetBase(parser::XML_Parser,base)
    ccall((:XML_SetBase,libexpat),Cint,(XML_Parser,Ptr{XML_Char}),parser,base)
end

function XML_GetBase(parser::XML_Parser)
    ccall((:XML_GetBase,libexpat),Ptr{XML_Char},(XML_Parser,),parser)
end

function XML_GetSpecifiedAttributeCount(parser::XML_Parser)
    ccall((:XML_GetSpecifiedAttributeCount,libexpat),Cint,(XML_Parser,),parser)
end

function XML_GetIdAttributeIndex(parser::XML_Parser)
    ccall((:XML_GetIdAttributeIndex,libexpat),Cint,(XML_Parser,),parser)
end

function XML_Parse(parser::XML_Parser,s,len::Cint,isFinal::Cint)
    ccall((:XML_Parse,libexpat),Cint,(XML_Parser,Cstring,Cint,Cint),parser,s,len,isFinal)
end

function XML_GetBuffer(parser::XML_Parser,len::Cint)
    ccall((:XML_GetBuffer,libexpat),Ptr{Void},(XML_Parser,Cint),parser,len)
end

function XML_ParseBuffer(parser::XML_Parser,len::Cint,isFinal::Cint)
    ccall((:XML_ParseBuffer,libexpat),Cint,(XML_Parser,Cint,Cint),parser,len,isFinal)
end

function XML_StopParser(parser::XML_Parser,resumable::XML_Bool)
    ccall((:XML_StopParser,libexpat),Cint,(XML_Parser,XML_Bool),parser,resumable)
end

function XML_ResumeParser(parser::XML_Parser)
    ccall((:XML_ResumeParser,libexpat),Cint,(XML_Parser,),parser)
end

function XML_GetParsingStatus(parser::XML_Parser,status)
    ccall((:XML_GetParsingStatus,libexpat),Void,(XML_Parser,Ptr{XML_ParsingStatus}),parser,status)
end

function XML_ExternalEntityParserCreate(parser::XML_Parser,context,encoding)
    ccall((:XML_ExternalEntityParserCreate,libexpat),XML_Parser,(XML_Parser,Ptr{XML_Char},Ptr{XML_Char}),parser,context,encoding)
end

function XML_SetParamEntityParsing(parser::XML_Parser,parsing::XML_ParamEntityParsing)
    ccall((:XML_SetParamEntityParsing,libexpat),Cint,(XML_Parser,XML_ParamEntityParsing),parser,parsing)
end

function XML_SetHashSalt(parser::XML_Parser,hash_salt::Culong)
    ccall((:XML_SetHashSalt,libexpat),Cint,(XML_Parser,Culong),parser,hash_salt)
end

function XML_GetErrorCode(parser::XML_Parser)
    ccall((:XML_GetErrorCode,libexpat),Cint,(XML_Parser,),parser)
end

function XML_GetCurrentLineNumber(parser::XML_Parser)
    ccall((:XML_GetCurrentLineNumber,libexpat),XML_Size,(XML_Parser,),parser)
end

function XML_GetCurrentColumnNumber(parser::XML_Parser)
    ccall((:XML_GetCurrentColumnNumber,libexpat),XML_Size,(XML_Parser,),parser)
end

function XML_GetCurrentByteIndex(parser::XML_Parser)
    ccall((:XML_GetCurrentByteIndex,libexpat),XML_Index,(XML_Parser,),parser)
end

function XML_GetCurrentByteCount(parser::XML_Parser)
    ccall((:XML_GetCurrentByteCount,libexpat),Cint,(XML_Parser,),parser)
end

function XML_GetInputContext(parser::XML_Parser,offset,size)
    ccall((:XML_GetInputContext,libexpat),Cstring,(XML_Parser,Ptr{Cint},Ptr{Cint}),parser,offset,size)
end

function XML_FreeContentModel(parser::XML_Parser,model)
    ccall((:XML_FreeContentModel,libexpat),Void,(XML_Parser,Ptr{XML_Content}),parser,model)
end

function XML_MemMalloc(parser::XML_Parser,size::Cint)
    ccall((:XML_MemMalloc,libexpat),Ptr{Void},(XML_Parser,Cint),parser,size)
end

function XML_MemRealloc(parser::XML_Parser,ptr,size::Cint)
    ccall((:XML_MemRealloc,libexpat),Ptr{Void},(XML_Parser,Ptr{Void},Cint),parser,ptr,size)
end

function XML_MemFree(parser::XML_Parser,ptr)
    ccall((:XML_MemFree,libexpat),Void,(XML_Parser,Ptr{Void}),parser,ptr)
end

function XML_ParserFree(parser::XML_Parser)
    ccall((:XML_ParserFree,libexpat),Void,(XML_Parser,),parser)
end

function XML_ErrorString(code::XML_Error)
    ccall((:XML_ErrorString,libexpat),Ptr{XML_LChar},(XML_Error,),code)
end

function XML_ExpatVersion()
    ccall((:XML_ExpatVersion,libexpat),Ptr{XML_LChar},())
end

function XML_ExpatVersionInfo()
    ccall((:XML_ExpatVersionInfo,libexpat),XML_Expat_Version,())
end

function XML_GetFeatureList()
    ccall((:XML_GetFeatureList,libexpat),Ptr{XML_Feature},())
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    schursolve(A, x₀, howmany, which, algorithm)

Compute a partial Schur decomposition containing `howmany` eigenvalues from the linear map
encoded in the matrix or function `A`. Return the reduced Schur matrix, the basis of Schur
vectors, the extracted eigenvalues and a `ConvergenceInfo` structure.

See also [`eigsolve`](@eigsolve) to obtain the eigenvectors instead. For real symmetric or
complex hermitian problems, the (partial) Schur decomposition is identical to the (partial)
eigenvalue decomposition, and `eigsolve` should always be used.

### Arguments:
The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or callable
object, that acts on vector like objects similar to `x₀`, which is the starting guess from
which a Krylov subspace will be built. `howmany` specifies how many Schur vectors should be
converged before the algorithm terminates; `which` specifies which eigenvalues should be targetted.
Valid specifications of `which` are
*   `LM`: eigenvalues of largest magnitude
*   `LR`: eigenvalues with largest (most positive) real part
*   `SR`: eigenvalues with smallest (most negative) real part
*   `LI`: eigenvalues with largest (most positive) imaginary part, only if `T <: Complex`
*   `SI`: eigenvalues with smallest (most negative) imaginary part, only if `T <: Complex`
*   [`ClosestTo(λ)`](@ref): eigenvalues closest to some number `λ`
!!! note "Note about selecting `which` eigenvalues"
    Krylov methods work well for extremal eigenvalues, i.e. eigenvalues on the periphery of
    the spectrum of the linear map. Even with `ClosestTo`, no shift and invert is performed.
    This is useful if, e.g., you know the spectrum to be within the unit circle in the complex
    plane, and want to target the eigenvalues closest to the value `λ = 1`.

The final argument `algorithm` can currently only be an instance of [`Arnoldi`](@ref), but
should nevertheless be specified. Since `schursolve` is less commonly used as `eigsolve`, no
convenient keyword syntax is currently available.

### Return values:
The return value is always of the form `T, vecs, vals, info = schursolve(...)` with
*   `T`: a `Matrix` containing the partial Schur decomposition of the linear map, i.e. it's
    elements are given by `T[i,j] = dot(vecs[i], f(vecs[j]))`. It is of Schur form, i.e.
    upper triangular in case of complex arithmetic, and block upper triangular (with at most
    2x2 blocks) in case of real arithmetic.
*   `vecs`: a `Vector` of corresponding Schur vectors, of the same length as `vals`. Note
    that Schur vectors are not returned as a matrix, as the linear map could act on any
    custom  Julia type with vector like behavior, i.e. the elements of the list `vecs` are
    objects that are typically similar to the starting guess `x₀`, up to a possibly
    different `eltype`. When the linear map is a simple `AbstractMatrix`, `vecs` will be
    `Vector{Vector{<:Number}}`. Schur vectors are by definition orthogonal, i.e.
    `dot(vecs[i],vecs[j]) = I[i,j]`. Note that Schur vectors are real if the problem (i.e.
    the linear map and the initial guess) are real.
*   `vals`: a `Vector` of eigenvalues, i.e. the diagonal elements of `T` in case of complex
    arithmetic, or extracted from the diagonal blocks in case of real arithmetic. Note that
    `vals` will always be complex, independent of the underlying arithmetic.
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: indicates how many eigenvalues and Schur vectors were
        actually converged to the specified tolerance (see below under keyword arguments)
    -   `info.residuals::Vector`: a list of the same length as `vals` containing the actual
        residuals
        ```julia
          info.residuals[i] = f(vecs[i]) - sum(vecs[j]*T[j,i] for j = 1:i+1)
        ```
        where `T[i+1,i]` is definitely zero in case of complex arithmetic and possibly zero
        in case of real arithmetic
    -   `info.normres::Vector{<:Real}`: list of the same length as `vals` containing the
        norm of the residual for every Schur vector, i.e.
        `info.normes[i] = norm(info.residual[i])`
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
    No warning is printed if not all requested eigenvalues were converged, so always check
    if `info.converged >= howmany`.

### Algorithm
The actual algorithm is an implementation of the Krylov-Schur algorithm, where the
[`Arnoldi`](@ref) algorithm is used to generate the Krylov subspace. During the algorith,
the Krylov subspace is dynamically grown and shrunk, i.e. the restarts are so-called thick
restarts where a part of the current Krylov subspace is kept.
"""
function schursolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
    T, U, fact, converged, numiter, numops = _schursolve(A, x₀, howmany, which, alg)
    if eltype(T) <: Real && howmany < length(fact) && T[howmany+1,howmany] != 0
        howmany += 1
    end
    if converged > howmany
        howmany = converged
    end
    TT = view(T,1:howmany,1:howmany)
    values = schur2eigvals(TT)
    vectors = let B = basis(fact)
        [B*u for u in cols(U, 1:howmany)]
    end
    residuals = let r = residual(fact)
        [mul!(similar(r), r, last(u)) for u in cols(U, 1:howmany)]
    end
    normresiduals = [normres(fact)*abs(last(u)) for u in cols(U, 1:howmany)]

    if alg.verbosity > 0
        if converged < howmany
            @warn """Arnoldi schursolve finished without convergence after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))"""
        else
            @info """Arnoldi schursolve finished after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))"""
        end
    end
    return TT, vectors, values, ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end

function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
    T, U, fact, converged, numiter, numops = _schursolve(A, x₀, howmany, which, alg)
    if eltype(T) <: Real && howmany < length(fact) && T[howmany+1,howmany] != 0
        howmany += 1
    end
    if converged > howmany
        howmany = converged
    end
    TT = view(T,1:howmany,1:howmany)
    values = schur2eigvals(TT)

    # Compute eigenvectors
    V = view(U, :, 1:howmany)*schur2eigvecs(TT)
    vectors = let B = basis(fact)
        [B*v for v in cols(V)]
    end
    residuals = let r = residual(fact)
        [mul!(similar(r, Base.promote_type(eltype(V), eltype(r))), r, last(v)) for v in cols(V)]
    end
    normresiduals = [normres(fact)*abs(last(v)) for v in cols(V)]

    if alg.verbosity > 0
        if converged < howmany
            @warn """Arnoldi eigsolve finished without convergence after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        else
            @info """Arnoldi eigsolve finished after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        end
    end
    return values, vectors, ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end

function _schursolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
    krylovdim = alg.krylovdim
    maxiter = alg.maxiter
    howmany > krylovdim && error("krylov dimension $(krylovdim) too small to compute $howmany eigenvalues")

    ## FIRST ITERATION: setting up
    numiter = 1
    # Compute arnoldi factorization
    iter = ArnoldiIterator(A, x₀, alg.orth)
    fact = initialize(iter; verbosity = alg.verbosity - 2)
    numops = 1
    sizehint!(fact, krylovdim)
    β = normres(fact)
    tol::eltype(β) = alg.tol
    if normres(fact) > tol || howmany > 1
        while length(fact) < krylovdim
            fact = expand!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
            normres(fact) < tol && length(fact) >= howmany && break
        end
    end

    # Process
    # allocate storage
    HH = fill(zero(eltype(fact)), krylovdim+1, krylovdim)
    UU = fill(zero(eltype(fact)), krylovdim, krylovdim)

    # initialize
    β = normres(fact)
    m = length(fact)
    H = view(HH, 1:m, 1:m)
    U = view(UU, 1:m, 1:m)
    f = view(HH, m+1, 1:m)
    copyto!(U, I)
    copyto!(H, rayleighquotient(fact))

    # compute dense schur factorization
    T, U, values = hschur!(H, U)
    by, rev = eigsort(which)
    p = sortperm(values, by = by, rev = rev)
    T, U = permuteschur!(T, U, p)
    f = mul!(f, view(U, m, :), β)
    converged = 0
    while converged < length(fact) && abs(f[converged+1]) < tol
        converged += 1
    end
    if eltype(T) <: Real && 0< converged < length(fact) && T[converged+1,converged] != 0
        converged -= 1
    end

    if alg.verbosity > 1
        msg = "Arnoldi schursolve in iter $numiter: "
        msg *= "$converged values converged, normres = ("
        msg *= @sprintf("%.2e", abs(f[1]))
        for i = 2:howmany
            msg *= ", "
            msg *= @sprintf("%.2e", abs(f[i]))
        end
        msg *= ")"
        @info msg
    end
    ## OTHER ITERATIONS: recycle
    while numiter < maxiter && converged < howmany
        numiter += 1

        # Determine how many to keep
        keep = div(3*m + 2*converged, 5) # strictly smaller than krylovdim since converged < howmany <= krylovdim, at least equal to converged
        if eltype(H) <: Real && H[keep+1,keep] != 0 # we are in the middle of a 2x2 block
            keep += 1 # conservative choice
            keep >= krylovdim && error("krylov dimension $(krylovdim) too small to compute $howmany eigenvalues")
        end

        # Restore Arnoldi form in the first keep columns
        @inbounds for j = 1:keep
            H[keep+1,j] = f[j]
        end
        @inbounds for j = keep:-1:1
            h, ν = householder(H, j+1, 1:j, j)
            H[j+1,j] = ν
            H[j+1,1:j-1] .= 0
            lmul!(h, H)
            rmul!(view(H, 1:j,:), h')
            rmul!(U, h')
        end
        copyto!(rayleighquotient(fact), H) # copy back into fact

        # Update B by applying U
        B = basis(fact)
        basistransform!(B, view(U, :, 1:keep))
        # for j = 1:m
        #     h, ν = householder(U, j:m, j)
        #     lmul!(h, view(U, :, j+1:krylovdim))
        #     rmul!(B, h')
        # end
        r = residual(fact)
        B[keep+1] = rmul!(r, 1/normres(fact))

        # Shrink Arnoldi factorization
        fact = shrink!(fact, keep)

        # Arnoldi factorization: recylce fact
        while length(fact) < krylovdim
            fact = expand!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
            normres(fact) < tol && length(fact) >= howmany && break
        end

        # post process
        β = normres(fact)
        m = length(fact)
        H = view(HH, 1:m, 1:m)
        U = view(UU, 1:m, 1:m)
        f = view(HH, m+1, 1:m)
        copyto!(U, I)
        copyto!(H, rayleighquotient(fact))

        # compute dense schur factorization
        T, U, values = hschur!(H, U)
        by, rev = eigsort(which)
        p = sortperm(values, by = by, rev = rev)
        T, U = permuteschur!(T, U, p)
        mul!(f, view(U,m,:), β)
        converged = 0
        while converged < length(fact) && abs(f[converged+1]) < tol
            converged += 1
        end
        if eltype(T) <: Real && 0 < converged < length(fact) && T[converged+1,converged] != 0
            converged -= 1
        end
        if alg.verbosity > 1
            msg = "Arnoldi schursolve in iter $numiter: "
            msg *= "$converged values converged, normres = ("
            msg *= @sprintf("%.2e", abs(f[1]))
            for i = 2:howmany
                msg *= ", "
                msg *= @sprintf("%.2e", abs(f[i]))
            end
            msg *= ")"
            @info msg
        end
    end
    return T, U, fact, converged, numiter, numops
end


text:
module Scoring

include("../Types.jl")

export run

function scorePeak(data, window_size)
        midpoint = Int(size(data)[1] / 2)
        diffLeft = 0.0
        diffRight = 0.0

        for i in range(1, midpoint; step=1)
            diffLeft += data[midpoint].magnitude - data[i].magnitude
        end
        for j in range(midpoint + 1, size(data)[1]; step=1)
            diffRight += data[midpoint].magnitude - data[j].magnitude
        end

        (diffRight + diffLeft) / window_size
end

function run(input:: Channel, window_size)
    window = []
    new_dp_index = Int(window_size / 2)
    return function (output:: Channel)
        for dp in input
            push!(window, dp)
            if size(window)[1] == window_size
                score = scorePeak(window, window_size)
                new_dp = Types.Datapoint(window[new_dp_index].time, score)
                put!(output, new_dp)
                popfirst!(window)
            end
        end
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
using KrylovMethods
using LinearOperators
using MatrixDepot
using Base.Test

println("=== Testing LSQR ===")

# small test
A   = [1 1; 0 -1; 1 0.]
r1  = A[:,1]   # compatible
r2  = randn(3) # incompatible
r3  = [1;1;-1.] # in the nullpace of A'

x1 = lsqr(A,r1,out=2,atol=1e-8,doBidiag=true,storeInterm=true)
@test norm(x1[1][:,end]-[1;0])<1e-10
@test x1[2]==1
x2 = lsqr(A,r2,out=1,atol=1e-8)
@test x2[2]==2
x3 = lsqr(A,r3,out=2)
@test all(x3[1].==0)

# test with sparse random matrix
A   = sprandn(100,10,.2)
rhs = randn(100)
xgt = full(A)\rhs
Af = (flag,x,a=0.0,v=0.0) -> ((flag=='F') ? A*x +a*v : A'*x+a*v)

x1  = lsqr(LinearOperator(A),rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10,out=1)
x2  = lsqr(A,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)
x3  = lsqr(Af,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)

@test norm(xgt-x1[1])/norm(xgt) < 1e-6
@test norm(xgt-x2[1])/norm(xgt) < 1e-6
@test norm(xgt-x3[1])/norm(xgt) < 1e-6
@test norm(x1[1]-x2[1])/norm(x1[1]) < 1e-12
@test norm(x2[1]-x3[1])/norm(x3[1]) < 1e-12

# test with dense hilbert matrix
A   = matrixdepot("hilb",3,3)
rhs = [1;0;0.]
xgt = A\rhs

x1  = lsqr(LinearOperator(A),rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)
x2  = lsqr(A,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)

@test norm(xgt-x1[1])/norm(xgt) < 1e-6
@test norm(xgt-x2[1])/norm(xgt) < 1e-6
@test norm(x1[1]-x2[1])/norm(x1[1]) < 1e-12
conA = norm(vec(A))*norm(vec(pinv(A)))
@test norm(conA-x1[3][end,end])/conA < 1e-12
@test norm(conA-x2[3][end,end])/conA < 1e-12

println("=== LSQR: All tests passed. ===")


text:
# standardized tables given in RFC6330

const RQ_DEGREE_TABLE = [
    0, 5243, 529531, 704294, 791675, 844104, 879057, 904023, 922747, 937311,
    948962, 958494, 966438, 973160, 978921, 983914, 988283, 992138, 995565,
    998631, 1001391, 1003887, 1006157, 1008229, 1010129, 1011876, 1013490,
    1014983, 1016370, 1017662, 1048576
]

"""table V2 of rfc6330"""
const V2 = [
    1629829892, 282540176, 2794583710, 496504798, 2990494426,
    3070701851, 2575963183, 4094823972, 2775723650, 4079480416,
    176028725, 2246241423, 3732217647, 2196843075, 1306949278,
    4170992780, 4039345809, 3209664269, 3387499533, 293063229,
    3660290503, 2648440860, 2531406539, 3537879412, 773374739,
    4184691853, 1804207821, 3347126643, 3479377103, 3970515774,
    1891731298, 2368003842, 3537588307, 2969158410, 4230745262,
    831906319, 2935838131, 264029468, 120852739, 3200326460,
    355445271, 2296305141, 1566296040, 1760127056, 20073893,
    3427103620, 2866979760, 2359075957, 2025314291, 1725696734,
    3346087406, 2690756527, 99815156, 4248519977, 2253762642,
    3274144518, 598024568, 3299672435, 556579346, 4121041856,
    2896948975, 3620123492, 918453629, 3249461198, 2231414958,
    3803272287, 3657597946, 2588911389, 242262274, 1725007475,
    2026427718, 46776484, 2873281403, 2919275846, 3177933051,
    1918859160, 2517854537, 1857818511, 3234262050, 479353687,
    200201308, 2801945841, 1621715769, 483977159, 423502325,
    3689396064, 1850168397, 3359959416, 3459831930, 841488699,
    3570506095, 930267420, 1564520841, 2505122797, 593824107,
    1116572080, 819179184, 3139123629, 1414339336, 1076360795,
    512403845, 177759256, 1701060666, 2239736419, 515179302,
    2935012727, 3821357612, 1376520851, 2700745271, 966853647,
    1041862223, 715860553, 171592961, 1607044257, 1227236688,
    3647136358, 1417559141, 4087067551, 2241705880, 4194136288,
    1439041934, 20464430, 119668151, 2021257232, 2551262694,
    1381539058, 4082839035, 498179069, 311508499, 3580908637,
    2889149671, 142719814, 1232184754, 3356662582, 2973775623,
    1469897084, 1728205304, 1415793613, 50111003, 3133413359,
    4074115275, 2710540611, 2700083070, 2457757663, 2612845330,
    3775943755, 2469309260, 2560142753, 3020996369, 1691667711,
    4219602776, 1687672168, 1017921622, 2307642321, 368711460,
    3282925988, 213208029, 4150757489, 3443211944, 2846101972,
    4106826684, 4272438675, 2199416468, 3710621281, 497564971,
    285138276, 765042313, 916220877, 3402623607, 2768784621,
    1722849097, 3386397442, 487920061, 3569027007, 3424544196,
    217781973, 2356938519, 3252429414, 145109750, 2692588106,
    2454747135, 1299493354, 4120241887, 2088917094, 932304329,
    1442609203, 952586974, 3509186750, 753369054, 854421006,
    1954046388, 2708927882, 4047539230, 3048925996, 1667505809,
    805166441, 1182069088, 4265546268, 4215029527, 3374748959,
    373532666, 2454243090, 2371530493, 3651087521, 2619878153,
    1651809518, 1553646893, 1227452842, 703887512, 3696674163,
    2552507603, 2635912901, 895130484, 3287782244, 3098973502,
    990078774, 3780326506, 2290845203, 41729428, 1949580860,
    2283959805, 1036946170, 1694887523, 4880696, 466000198,
    2765355283, 3318686998, 1266458025, 3919578154, 3545413527,
    2627009988, 3744680394, 1696890173, 3250684705, 4142417708,
    915739411, 3308488877, 1289361460, 2942552331, 1169105979,
    3342228712, 698560958, 1356041230, 2401944293, 107705232,
    3701895363, 903928723, 3646581385, 844950914, 1944371367,
    3863894844, 2946773319, 1972431613, 1706989237, 29917467,
    3497665928
]

"""table V3 of rfc6330"""
const V3 = [
    1191369816, 744902811, 2539772235, 3213192037, 3286061266,
    1200571165, 2463281260, 754888894, 714651270, 1968220972,
    3628497775, 1277626456, 1493398934, 364289757, 2055487592,
    3913468088, 2930259465, 902504567, 3967050355, 2056499403,
    692132390, 186386657, 832834706, 859795816, 1283120926,
    2253183716, 3003475205, 1755803552, 2239315142, 4271056352,
    2184848469, 769228092, 1249230754, 1193269205, 2660094102,
    642979613, 1687087994, 2726106182, 446402913, 4122186606,
    3771347282, 37667136, 192775425, 3578702187, 1952659096,
    3989584400, 3069013882, 2900516158, 4045316336, 3057163251,
    1702104819, 4116613420, 3575472384, 2674023117, 1409126723,
    3215095429, 1430726429, 2544497368, 1029565676, 1855801827,
    4262184627, 1854326881, 2906728593, 3277836557, 2787697002,
    2787333385, 3105430738, 2477073192, 748038573, 1088396515,
    1611204853, 201964005, 3745818380, 3654683549, 3816120877,
    3915783622, 2563198722, 1181149055, 33158084, 3723047845,
    3790270906, 3832415204, 2959617497, 372900708, 1286738499,
    1932439099, 3677748309, 2454711182, 2757856469, 2134027055,
    2780052465, 3190347618, 3758510138, 3626329451, 1120743107,
    1623585693, 1389834102, 2719230375, 3038609003, 462617590,
    260254189, 3706349764, 2556762744, 2874272296, 2502399286,
    4216263978, 2683431180, 2168560535, 3561507175, 668095726,
    680412330, 3726693946, 4180630637, 3335170953, 942140968,
    2711851085, 2059233412, 4265696278, 3204373534, 232855056,
    881788313, 2258252172, 2043595984, 3758795150, 3615341325,
    2138837681, 1351208537, 2923692473, 3402482785, 2105383425,
    2346772751, 499245323, 3417846006, 2366116814, 2543090583,
    1828551634, 3148696244, 3853884867, 1364737681, 2200687771,
    2689775688, 232720625, 4071657318, 2671968983, 3531415031,
    1212852141, 867923311, 3740109711, 1923146533, 3237071777,
    3100729255, 3247856816, 906742566, 4047640575, 4007211572,
    3495700105, 1171285262, 2835682655, 1634301229, 3115169925,
    2289874706, 2252450179, 944880097, 371933491, 1649074501,
    2208617414, 2524305981, 2496569844, 2667037160, 1257550794,
    3399219045, 3194894295, 1643249887, 342911473, 891025733,
    3146861835, 3789181526, 938847812, 1854580183, 2112653794,
    2960702988, 1238603378, 2205280635, 1666784014, 2520274614,
    3355493726, 2310872278, 3153920489, 2745882591, 1200203158,
    3033612415, 2311650167, 1048129133, 4206710184, 4209176741,
    2640950279, 2096382177, 4116899089, 3631017851, 4104488173,
    1857650503, 3801102932, 445806934, 3055654640, 897898279,
    3234007399, 1325494930, 2982247189, 1619020475, 2720040856,
    885096170, 3485255499, 2983202469, 3891011124, 546522756,
    1524439205, 2644317889, 2170076800, 2969618716, 961183518,
    1081831074, 1037015347, 3289016286, 2331748669, 620887395,
    303042654, 3990027945, 1562756376, 3413341792, 2059647769,
    2823844432, 674595301, 2457639984, 4076754716, 2447737904,
    1583323324, 625627134, 3076006391, 345777990, 1684954145,
    879227329, 3436182180, 1522273219, 3802543817, 1456017040,
    1897819847, 2970081129, 1382576028, 3820044861, 1044428167,
    612252599, 3340478395, 2150613904, 3397625662, 3573635640,
    3432275192
]

"""table 5.6 of rfc6330 with column labels [Kp, J(Kp), S(Kp), H(Kp), W(Kp)]"""
const RQ_parameter_table = [
    10 254 7 10 17
    12 630 7 10 19
    18 682 11 10 29
    20 293 11 10 31
    26 80 11 10 37
    30 566 11 10 41
    32 860 11 10 43
    36 267 11 10 47
    42 822 11 10 53
    46 506 13 10 59
    48 589 13 10 61
    49 87 13 10 61
    55 520 13 10 67
    60 159 13 10 71
    62 235 13 10 73
    69 157 13 10 79
    75 502 17 10 89
    84 334 17 10 97
    88 583 17 10 101
    91 66 17 10 103
    95 352 17 10 107
    97 365 17 10 109
    101 562 17 10 113
    114 5 19 10 127
    119 603 19 10 131
    125 721 19 10 137
    127 28 19 10 139
    138 660 19 10 149
    140 829 19 10 151
    149 900 23 10 163
    153 930 23 10 167
    160 814 23 10 173
    166 661 23 10 179
    168 693 23 10 181
    179 780 23 10 191
    181 605 23 10 193
    185 551 23 10 197
    187 777 23 10 199
    200 491 23 10 211
    213 396 23 10 223
    217 764 29 10 233
    225 843 29 10 241
    236 646 29 10 251
    242 557 29 10 257
    248 608 29 10 263
    257 265 29 10 271
    263 505 29 10 277
    269 722 29 10 283
    280 263 29 10 293
    295 999 29 10 307
    301 874 29 10 313
    305 160 29 10 317
    324 575 31 10 337
    337 210 31 10 349
    341 513 31 10 353
    347 503 31 10 359
    355 558 31 10 367
    362 932 31 10 373
    368 404 31 10 379
    372 520 37 10 389
    380 846 37 10 397
    385 485 37 10 401
    393 728 37 10 409
    405 554 37 10 421
    418 471 37 10 433
    428 641 37 10 443
    434 732 37 10 449
    447 193 37 10 461
    453 934 37 10 467
    466 864 37 10 479
    478 790 37 10 491
    486 912 37 10 499
    491 617 37 10 503
    497 587 37 10 509
    511 800 37 10 523
    526 923 41 10 541
    532 998 41 10 547
    542 92 41 10 557
    549 497 41 10 563
    557 559 41 10 571
    563 667 41 10 577
    573 912 41 10 587
    580 262 41 10 593
    588 152 41 10 601
    594 526 41 10 607
    600 268 41 10 613
    606 212 41 10 619
    619 45 41 10 631
    633 898 43 10 647
    640 527 43 10 653
    648 558 43 10 661
    666 460 47 10 683
    675 5 47 10 691
    685 895 47 10 701
    693 996 47 10 709
    703 282 47 10 719
    718 513 47 10 733
    728 865 47 10 743
    736 870 47 10 751
    747 239 47 10 761
    759 452 47 10 773
    778 862 53 10 797
    792 852 53 10 811
    802 643 53 10 821
    811 543 53 10 829
    821 447 53 10 839
    835 321 53 10 853
    845 287 53 10 863
    860 12 53 10 877
    870 251 53 10 887
    891 30 53 10 907
    903 621 53 10 919
    913 555 53 10 929
    926 127 53 10 941
    938 400 53 10 953
    950 91 59 10 971
    963 916 59 10 983
    977 935 59 10 997
    989 691 59 10 1009
    1002 299 59 10 1021
    1020 282 59 10 1039
    1032 824 59 10 1051
    1050 536 59 11 1069
    1074 596 59 11 1093
    1085 28 59 11 1103
    1099 947 59 11 1117
    1111 162 59 11 1129
    1136 536 59 11 1153
    1152 1000 61 11 1171
    1169 251 61 11 1187
    1183 673 61 11 1201
    1205 559 61 11 1223
    1220 923 61 11 1237
    1236 81 67 11 1259
    1255 478 67 11 1277
    1269 198 67 11 1291
    1285 137 67 11 1307
    1306 75 67 11 1327
    1347 29 67 11 1367
    1361 231 67 11 1381
    1389 532 67 11 1409
    1404 58 67 11 1423
    1420 60 67 11 1439
    1436 964 71 11 1459
    1461 624 71 11 1483
    1477 502 71 11 1499
    1502 636 71 11 1523
    1522 986 71 11 1543
    1539 950 71 11 1559
    1561 735 73 11 1583
    1579 866 73 11 1601
    1600 203 73 11 1621
    1616 83 73 11 1637
    1649 14 73 11 1669
    1673 522 79 11 1699
    1698 226 79 11 1723
    1716 282 79 11 1741
    1734 88 79 11 1759
    1759 636 79 11 1783
    1777 860 79 11 1801
    1800 324 79 11 1823
    1824 424 79 11 1847
    1844 999 79 11 1867
    1863 682 83 11 1889
    1887 814 83 11 1913
    1906 979 83 11 1931
    1926 538 83 11 1951
    1954 278 83 11 1979
    1979 580 83 11 2003
    2005 773 83 11 2029
    2040 911 89 11 2069
    2070 506 89 11 2099
    2103 628 89 11 2131
    2125 282 89 11 2153
    2152 309 89 11 2179
    2195 858 89 11 2221
    2217 442 89 11 2243
    2247 654 89 11 2273
    2278 82 97 11 2311
    2315 428 97 11 2347
    2339 442 97 11 2371
    2367 283 97 11 2399
    2392 538 97 11 2423
    2416 189 97 11 2447
    2447 438 97 11 2477
    2473 912 97 11 2503
    2502 1 97 11 2531
    2528 167 97 11 2557
    2565 272 97 11 2593
    2601 209 101 11 2633
    2640 927 101 11 2671
    2668 386 101 11 2699
    2701 653 101 11 2731
    2737 669 101 11 2767
    2772 431 101 11 2801
    2802 793 103 11 2833
    2831 588 103 11 2861
    2875 777 107 11 2909
    2906 939 107 11 2939
    2938 864 107 11 2971
    2979 627 107 11 3011
    3015 265 109 11 3049
    3056 976 109 11 3089
    3101 988 113 11 3137
    3151 507 113 11 3187
    3186 640 113 11 3221
    3224 15 113 11 3259
    3265 667 113 11 3299
    3299 24 127 11 3347
    3344 877 127 11 3391
    3387 240 127 11 3433
    3423 720 127 11 3469
    3466 93 127 11 3511
    3502 919 127 11 3547
    3539 635 127 11 3583
    3579 174 127 11 3623
    3616 647 127 11 3659
    3658 820 127 11 3701
    3697 56 127 11 3739
    3751 485 127 11 3793
    3792 210 127 11 3833
    3840 124 127 11 3881
    3883 546 127 11 3923
    3924 954 131 11 3967
    3970 262 131 11 4013
    4015 927 131 11 4057
    4069 957 131 11 4111
    4112 726 137 11 4159
    4165 583 137 11 4211
    4207 782 137 11 4253
    4252 37 137 11 4297
    4318 758 137 11 4363
    4365 777 137 11 4409
    4418 104 139 11 4463
    4468 476 139 11 4513
    4513 113 149 11 4567
    4567 313 149 11 4621
    4626 102 149 11 4679
    4681 501 149 11 4733
    4731 332 149 11 4783
    4780 786 149 11 4831
    4838 99 149 11 4889
    4901 658 149 11 4951
    4954 794 149 11 5003
    5008 37 151 11 5059
    5063 471 151 11 5113
    5116 94 157 11 5171
    5172 873 157 11 5227
    5225 918 157 11 5279
    5279 945 157 11 5333
    5334 211 157 11 5387
    5391 341 157 11 5443
    5449 11 163 11 5507
    5506 578 163 11 5563
    5566 494 163 11 5623
    5637 694 163 11 5693
    5694 252 163 11 5749
    5763 451 167 11 5821
    5823 83 167 11 5881
    5896 689 167 11 5953
    5975 488 173 11 6037
    6039 214 173 11 6101
    6102 17 173 11 6163
    6169 469 173 11 6229
    6233 263 179 11 6299
    6296 309 179 11 6361
    6363 984 179 11 6427
    6427 123 179 11 6491
    6518 360 179 11 6581
    6589 863 181 11 6653
    6655 122 181 11 6719
    6730 522 191 11 6803
    6799 539 191 11 6871
    6878 181 191 11 6949
    6956 64 191 11 7027
    7033 387 191 11 7103
    7108 967 191 11 7177
    7185 843 191 11 7253
    7281 999 193 11 7351
    7360 76 197 11 7433
    7445 142 197 11 7517
    7520 599 197 11 7591
    7596 576 199 11 7669
    7675 176 211 11 7759
    7770 392 211 11 7853
    7855 332 211 11 7937
    7935 291 211 11 8017
    8030 913 211 11 8111
    8111 608 211 11 8191
    8194 212 211 11 8273
    8290 696 211 11 8369
    8377 931 223 11 8467
    8474 326 223 11 8563
    8559 228 223 11 8647
    8654 706 223 11 8741
    8744 144 223 11 8831
    8837 83 223 11 8923
    8928 743 223 11 9013
    9019 187 223 11 9103
    9111 654 227 11 9199
    9206 359 227 11 9293
    9303 493 229 11 9391
    9400 369 233 11 9491
    9497 981 233 11 9587
    9601 276 239 11 9697
    9708 647 239 11 9803
    9813 389 239 11 9907
    9916 80 239 11 10009
    10017 396 241 11 10111
    10120 580 251 11 10223
    10241 873 251 11 10343
    10351 15 251 11 10453
    10458 976 251 11 10559
    10567 584 251 11 10667
    10676 267 257 11 10781
    10787 876 257 11 10891
    10899 642 257 12 11003
    11015 794 257 12 11119
    11130 78 263 12 11239
    11245 736 263 12 11353
    11358 882 269 12 11471
    11475 251 269 12 11587
    11590 434 269 12 11701
    11711 204 269 12 11821
    11829 256 271 12 11941
    11956 106 277 12 12073
    12087 375 277 12 12203
    12208 148 277 12 12323
    12333 496 281 12 12451
    12460 88 281 12 12577
    12593 826 293 12 12721
    12726 71 293 12 12853
    12857 925 293 12 12983
    13002 760 293 12 13127
    13143 130 293 12 13267
    13284 641 307 12 13421
    13417 400 307 12 13553
    13558 480 307 12 13693
    13695 76 307 12 13829
    13833 665 307 12 13967
    13974 910 307 12 14107
    14115 467 311 12 14251
    14272 964 311 12 14407
    14415 625 313 12 14551
    14560 362 317 12 14699
    14713 759 317 12 14851
    14862 728 331 12 15013
    15011 343 331 12 15161
    15170 113 331 12 15319
    15325 137 331 12 15473
    15496 308 331 12 15643
    15651 800 337 12 15803
    15808 177 337 12 15959
    15977 961 337 12 16127
    16161 958 347 12 16319
    16336 72 347 12 16493
    16505 732 347 12 16661
    16674 145 349 12 16831
    16851 577 353 12 17011
    17024 305 353 12 17183
    17195 50 359 12 17359
    17376 351 359 12 17539
    17559 175 367 12 17729
    17742 727 367 12 17911
    17929 902 367 12 18097
    18116 409 373 12 18289
    18309 776 373 12 18481
    18503 586 379 12 18679
    18694 451 379 12 18869
    18909 287 383 12 19087
    19126 246 389 12 19309
    19325 222 389 12 19507
    19539 563 397 12 19727
    19740 839 397 12 19927
    19939 897 401 12 20129
    20152 409 401 12 20341
    20355 618 409 12 20551
    20564 439 409 12 20759
    20778 95 419 13 20983
    20988 448 419 13 21191
    21199 133 419 13 21401
    21412 938 419 13 21613
    21629 423 431 13 21841
    21852 90 431 13 22063
    22073 640 431 13 22283
    22301 922 433 13 22511
    22536 250 439 13 22751
    22779 367 439 13 22993
    23010 447 443 13 23227
    23252 559 449 13 23473
    23491 121 457 13 23719
    23730 623 457 13 23957
    23971 450 457 13 24197
    24215 253 461 13 24443
    24476 106 467 13 24709
    24721 863 467 13 24953
    24976 148 479 13 25219
    25230 427 479 13 25471
    25493 138 479 13 25733
    25756 794 487 13 26003
    26022 247 487 13 26267
    26291 562 491 13 26539
    26566 53 499 13 26821
    26838 135 499 13 27091
    27111 21 503 13 27367
    27392 201 509 13 27653
    27682 169 521 13 27953
    27959 70 521 13 28229
    28248 386 521 13 28517
    28548 226 523 13 28817
    28845 3 541 13 29131
    29138 769 541 13 29423
    29434 590 541 13 29717
    29731 672 541 13 30013
    30037 713 547 13 30323
    30346 967 547 13 30631
    30654 368 557 14 30949
    30974 348 557 14 31267
    31285 119 563 14 31583
    31605 503 569 14 31907
    31948 181 571 14 32251
    32272 394 577 14 32579
    32601 189 587 14 32917
    32932 210 587 14 33247
    33282 62 593 14 33601
    33623 273 593 14 33941
    33961 554 599 14 34283
    34302 936 607 14 34631
    34654 483 607 14 34981
    35031 397 613 14 35363
    35395 241 619 14 35731
    35750 500 631 14 36097
    36112 12 631 14 36457
    36479 958 641 14 36833
    36849 524 641 14 37201
    37227 8 643 14 37579
    37606 100 653 14 37967
    37992 339 653 14 38351
    38385 804 659 14 38749
    38787 510 673 14 39163
    39176 18 673 14 39551
    39576 412 677 14 39953
    39980 394 683 14 40361
    40398 830 691 15 40787
    40816 535 701 15 41213
    41226 199 701 15 41621
    41641 27 709 15 42043
    42067 298 709 15 42467
    42490 368 719 15 42899
    42916 755 727 15 43331
    43388 379 727 15 43801
    43840 73 733 15 44257
    44279 387 739 15 44701
    44729 457 751 15 45161
    45183 761 751 15 45613
    45638 855 757 15 46073
    46104 370 769 15 46549
    46574 261 769 15 47017
    47047 299 787 15 47507
    47523 920 787 15 47981
    48007 269 787 15 48463
    48489 862 797 15 48953
    48976 349 809 15 49451
    49470 103 809 15 49943
    49978 115 821 15 50461
    50511 93 821 16 50993
    51017 982 827 16 51503
    51530 432 839 16 52027
    52062 340 853 16 52571
    52586 173 853 16 53093
    53114 421 857 16 53623
    53650 330 863 16 54163
    54188 624 877 16 54713
    54735 233 877 16 55259
    55289 362 883 16 55817
    55843 963 907 16 56393
    56403 471 907 16 56951
]

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
#=
Package adapted from https://github.com/luca-aki/TraceEstimation.jl
which has an MIT license as follows.

MIT License

Copyright (c) 2021 Akshay Jain and Contributors

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

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

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

module TraceEstimators

include("common.jl")
include("hutchinson.jl")
include("diagapp.jl")
include("slq.jl")
include("chebyhutch.jl")
include("diagonalapprox.jl")

end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using PGFPlotsX, FoodCachingPlotting, DataFrames, DataFramesMeta
import FoodCachingExperiments: bload, CLAYTON0103_EXPERIMENTS, EXPERIMENTS, summarize

include(joinpath(@__DIR__, "paths.jl"))

###
### Appendix result figures
###
results = bload(joinpath(datapath, "results"))[:results]

sort!(results, :logp_hat)
best = combine(groupby(results, [:model, :experiment]),
               names(results) .=> x -> [last(x)], renamecols = false)

for row in eachrow(best)
    @show row.model row.experiment
    f = plot_compare(row.experiment, row.avg, row.best, backend = :pgf)
    if row.experiment == :Clayton0103
        for (ei, fi) in zip(CLAYTON0103_EXPERIMENTS, f)
            pgfsave(joinpath(figpath, "$(row.model)_$(ei).tikz"), fi)
        end
    else
        pgfsave(joinpath(figpath, "$(row.model)_$(row.experiment).tikz"), f)
    end
end

###
### Appendix document
###

length_without_missing(x) = sum(1 .- ismissing.(x))
process_testval(val) = match(r"[0-9.]", string(val[1])) !== nothing ? "=" * val : val
texstring(s) = reduce(replace, ["_" => " ", "&" => "\\&", "×" => "\$\\times\$"], init = string(s))
figstring(m, s) = "\\begin{center}\\textbf{$m}\\end{center}\\input{$(m)_$s.tikz}"
function process_pval(val)
    pval = string(round(val, sigdigits = 2))
    m = match(r"(.*)e(.*)", pval)
    if m !== nothing
        "$(m.captures[1])\\cdot 10^{$(m.captures[2])}"
    else
        pval
    end
end
function experiment_summary()
    result = DataFrame(key = Symbol[],
                       ntests = Int[],
                       summary = String[],
                       comments = String[],
                       nvals = Int[],
                       atests = Array{String}[])
    for (key, experiment) in EXPERIMENTS
        key == :Clayton0103 && continue
        ntests = sum(length(test.pvalues) for (_, test) in experiment.tests)
        atests = vcat([["$k ($(test.keys[i])): \$$(test.statistic[i])($(test.df[i])) $(process_testval(test.value[i])), p = $(process_pval(test.pvalues[i]))\$"
                        for i in eachindex(test.labels) if test.labels[i] == "A"]
                       for (k, test) in experiment.tests]...)
        data = experiment.data
        if key ∈ (:Cheke11_specsat, :Cheke11_planning)
            data = summarize(experiment, data)
        end
        nvals = length_without_missing(data.μ)
        hasproperty(data, :sem) && (nvals += length_without_missing(data.sem))
        hasproperty(data, :firstinspection) && (nvals += length_without_missing(data.firstinspection))
        push!(result, [key, ntests, experiment.major_finding,
                       experiment.comments, nvals, atests])
    end
    result
end
function section(row)
    key = row.key
    s = "\\subsection{$(texstring(key))}\n"
    s *= "Published in \\cite{$(split(string(key), '_')[1])}."
    s *= "\\subsubsection*{Major Finding}\n$(texstring(row.summary))\n\n"
    if row.comments != ""
        s *= "\\subsubsection*{Comments}\n$(texstring(row.comments))\n"
    end
#     s *= "\\subsubsection*{Most Important Tests}\n"
#     s *= join(texstring.(row.atests), "\\newline\n")
    s *= "\\subsubsection*{Total Number of Statistics}\n$(row.ntests) p-values, $(row.nvals) means or sems.\n"
    s *= "\\subsubsection*{Figures}\\begin{center}"
    for model in ["Baseline", "MotivationalControl", "EpisodicLikeMemory",
                  "PlasticCaching", "ReplayAndPlan"]
        s *= figstring(model, key)
    end
    s *= raw"\end{center}"
    s
end

summary = experiment_summary()
summary = @transform summary ntotal = :ntests .+ :nvals
summary = @transform summary nAtests = length.(:atests)

include(joinpath(@__DIR__, "order.jl"))
sort!(summary, :key, by = x -> findfirst(==(x), experiment_order))
summary = join([section(summary[i, :]) for i in 1:nrow(summary)], "\n\n");
open(joinpath(figpath, "experiment_summary.tex"), "w") do io
    write(io, summary)
end


text:
#=
    Contents: Create, read, update and delete functions for server-side sessions that use a LoggedDict as the database.

    Schema: ld[session_id] => Dict(...)
=#
using LoggedDicts

"""
Init session_id => session on the database and set the specified cookie to session_id.
Return: session_id
"""
function session_create!(sessions::LoggedDict, username::AbstractString, res::Response, cookiename::AbstractString)
    get_n_sessions(username) >= config[:session][:max_n_sessions] && return
    session_id = generate_session_id()
    set!(sessions, session_id, Dict{AbstractString, Any}("username" => username))
    haskey(config[:session], :timeout) && update_lastreq!(sessions, session_id)
    add_sessionid_to_user!(username, session_id)
    write_to_cookie!(res, cookiename, session_id)
    session_id
end

"Delete session from database and set the specified cookie to an invalid state."
function session_delete!(sessions::LoggedDict, session_id::AbstractString)
    if haskey(sessions, session_id)
	username = get(sessions, session_id, "username")
	delete!(sessions, session_id)
	remove_sessionid_from_user!(username, session_id)
    end
end

function session_is_valid(sessions::LoggedDict, session_id::AbstractString)
    result = haskey(sessions, session_id)
    if haskey(config[:session], :timeout)
	if haskey(sessions, session_id, "lastreq")
	    dt_str = get(sessions, session_id, "lastreq")
	    dt     = DateTime(dt_str)
	    if dt + Dates.Second(config[:session][:timeout]) < now()
		result = false    # Session has timed out
	    end
	else
	    result = false
	end
    end
    result
end

function update_lastreq!(sessions::LoggedDict, session_id::AbstractString)
    set!(sessions, session_id, "lastreq", string(now()))
end

function session_get(sessions::LoggedDict, session_id, keys...)
    result = nothing
    if haskey(sessions, session_id, keys...)
	result = get(sessions, session_id, keys...)
    end
    result
end

function session_set!(sessions::LoggedDict, session_id, keys_value...)
    set!(sessions, session_id, keys_value...)
end

# EOF


text:
export preparation!

function preparation!( ua         :: UA,
                       dt         :: Float64,
                       particles  :: Particles,
                       xt         :: Array{ComplexF64,3},
                       yt         :: Array{ComplexF64,3})

    ε      = ua.ε
    ntau   = ua.ntau
    nbpart = particles.nbpart
    tau    = ua.tau
    ltau   = ua.ltau

    for m = 1:nbpart

        x1 = particles.x[1,m]
        x2 = particles.x[2,m]

        particles.b[m] = 1 + 0.5 * sin(x1) * sin(x2)
        particles.t[m] = dt * particles.b[m]

        ua.pl[1,m] = particles.t[m]
        ua.ql[1,m] = particles.t[m]^2 / 2

        t = particles.t[m]
        b = particles.b[m]

        for n=2:ntau
            elt = exp(-1im*ua.ltau[n]*t/ε)
            ua.pl[n,m] = ε * 1im*(elt-1)/ua.ltau[n]
            ua.ql[n,m] = ε * (ε*(1-elt) -1im*ua.ltau[n]*t)/ua.ltau[n]^2
        end

        ex,  ey  = particles.e[1:2,m]
        vx,  vy  = particles.v[1:2,m]
        vxb, vyb = vx/b, vy/b

        for n = 1:ntau

            τ = tau[n]

            h1 = ε * (sin(τ) * vxb - cos(τ) * vyb)
            h2 = ε * (sin(τ) * vyb + cos(τ) * vxb)

            xt1 = x1 + h1 + ε * vyb
            xt2 = x2 + h2 - ε * vxb

            xt[n,1,m] = xt1
            xt[n,2,m] = xt2

            interv=(1+0.5*sin(xt1)*sin(xt2)-b)/ε

            exb = ((  cos(τ)*vy - sin(τ)*vx) * interv + ex)/b
            eyb = (( -cos(τ)*vx - sin(τ)*vy) * interv + ey)/b

            ua.r[n,1] = cos(τ)* exb - sin(τ) * eyb
            ua.r[n,2] = sin(τ)* exb + cos(τ) * eyb

        end

        mul!(ua.r̃, ua.rtau, ua.r)

        for n = 2:ntau
            ua.r̃[n,1] = -1im * ua.r̃[n,1]/ltau[n]
            ua.r̃[n,2] = -1im * ua.r̃[n,2]/ltau[n]
        end

        ldiv!(ua.r, ua.rtau, ua.r̃)

        for n = 1:ntau
            yt[n,1,m] = vx + (ua.r[n,1] - ua.r[1,1]) * ε
            yt[n,2,m] = vy + (ua.r[n,2] - ua.r[1,2]) * ε
        end

    end

end

export update_particles_e!

function update_particles_e!( particles :: Particles,
                              et        :: Array{Float64,3},
                              fields    :: MeshFields,
                              ua        :: UA,
                              xt        :: Array{ComplexF64,3})

    interpol_eb_m6!( et, fields, xt, particles.nbpart, ua.ntau)

end

export update_particles_x!

function update_particles_x!( particles :: Particles,
                              fields    :: MeshFields,
                              ua        :: UA,
                              xt        :: Array{ComplexF64,3})

    compute_rho_m6!( fields, particles, xt, ua )

end

export compute_f!

function compute_f!( fx        :: Array{ComplexF64,3},
                     fy        :: Array{ComplexF64,3},
                     ua        :: UA,
                     particles :: Particles,
                     xt        :: Array{ComplexF64,3},
                     yt        :: Array{ComplexF64,3},
                     et        :: Array{Float64,3} )

    for m=1:particles.nbpart

        b = particles.b[m]

        for n=1:ua.ntau

            xt1 = real(xt[n,1,m])
            xt2 = real(xt[n,2,m])

            yt1 = yt[n,1,m]
            yt2 = yt[n,2,m]

            τ = ua.tau[n]

            fx[n,1,m] = ( cos(τ) * yt1 + sin(τ) * yt2)/b
            fx[n,2,m] = (-sin(τ) * yt1 + cos(τ) * yt2)/b

            interv = (1 + 0.5*sin(xt1)*sin(xt2)-b)/ua.ε

            tmp1 = et[n,1,m]+(  cos(τ)*yt2 - sin(τ)*yt1)*interv
            tmp2 = et[n,2,m]+(- cos(τ)*yt1 - sin(τ)*yt2)*interv

            fy[n,1,m] = (cos(τ)*tmp1-sin(τ)*tmp2)/b
            fy[n,2,m] = (sin(τ)*tmp1+cos(τ)*tmp2)/b

        end

    end

    fft!(fx,1)
    fft!(fy,1)

end

export ua_step!

function ua_step!( xt        :: Array{ComplexF64, 3},
                   x̃t        :: Array{ComplexF64, 3},
                   ua        :: UA,
                   particles :: Particles,
                   fx        :: Array{ComplexF64, 3} )

    for m=1:particles.nbpart

        t = particles.t[m]

        for n=1:ua.ntau

            elt = exp(-1im*ua.ltau[n]*t/ua.ε)
            xt[n,1,m] = elt * x̃t[n,1,m] + ua.pl[n,m] * fx[n,1,m]
            xt[n,2,m] = elt * x̃t[n,2,m] + ua.pl[n,m] * fx[n,2,m]

        end

    end

end

function ua_step!( xt        :: Array{ComplexF64, 3},
                   x̃t        :: Array{ComplexF64, 3},
                   ua        :: UA,
                   particles :: Particles,
                   fx        :: Array{ComplexF64, 3},
                   gx        :: Array{ComplexF64, 3} )

     for m=1:particles.nbpart

         t = particles.t[m]

         for n=1:ua.ntau

             elt = exp(-1im*ua.ltau[n]*t/ua.ε)
             fx1 = fx[n,1,m]
             fx2 = fx[n,2,m]
             xt1 = elt * x̃t[n,1,m] + ua.pl[n,m] * fx1
             xt2 = elt * x̃t[n,2,m] + ua.pl[n,m] * fx2
             xt1 += ua.ql[n,m] * (gx[n,1,m] - fx1) / t
             xt2 += ua.ql[n,m] * (gx[n,2,m] - fx2) / t

             xt[n,1,m] = xt1
             xt[n,2,m] = xt2

         end

    end

end

export compute_v!

function compute_v!( yt       :: Array{ComplexF64, 3},
                    particles :: Particles,
                    ua        :: UA )

        for m=1:particles.nbpart

            t = particles.t[m]

            px, py = 0.0im, 0.0im
            for n = 1:ua.ntau
                elt = exp(1im*ua.ltau[n]*t/ua.ε)
                px += yt[n,1,m]/ua.ntau * elt
                py += yt[n,2,m]/ua.ntau * elt
            end

            particles.v[1,m] = real(cos(t/ua.ε)*px+sin(t/ua.ε)*py)
            particles.v[2,m] = real(cos(t/ua.ε)*py-sin(t/ua.ε)*px)

        end

end


text:
# Autogenerated wrapper script for JasPer_jll for armv7l-linux-gnueabihf
export imgcmp, imginfo, jasper, libjasper

using JpegTurbo_jll
JLLWrappers.@generate_wrapper_header("JasPer")
JLLWrappers.@declare_executable_product(imgcmp)
JLLWrappers.@declare_executable_product(imginfo)
JLLWrappers.@declare_executable_product(jasper)
JLLWrappers.@declare_library_product(libjasper, "libjasper.so.4")
function __init__()
    JLLWrappers.@generate_init_header(JpegTurbo_jll)
    JLLWrappers.@init_executable_product(
        imgcmp,
        "bin/imgcmp",
    )

    JLLWrappers.@init_executable_product(
        imginfo,
        "bin/imginfo",
    )

    JLLWrappers.@init_executable_product(
        jasper,
        "bin/jasper",
    )

    JLLWrappers.@init_library_product(
        libjasper,
        "lib/libjasper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
InverseFunctions.inverse(::typeof(log1pexp)) = logexpm1
InverseFunctions.inverse(::typeof(logexpm1)) = log1pexp

InverseFunctions.inverse(::typeof(log1mexp)) = log1mexp

InverseFunctions.inverse(::typeof(log2mexp)) = log2mexp

InverseFunctions.inverse(::typeof(logit)) = logistic
InverseFunctions.inverse(::typeof(logistic)) = logit