text:
"""
    Gauging(level::Real, discharge::Real)

Construct an object of type Gauging
"""
struct Gauging
    level::Real
    discharge::Real
    Gauging(level::T,discharge::T) where {T<:Real} = new(level,discharge)
end

Gauging(level::Real, discharge::Real) = Gauging(promote(level, discharge)...)
Gauging(level::Integer, discharge::Integer,) = Gauging(float(level), float(discharge))
Gauging() = Gauging(0.0, 0.0)

"""
    bootstrap(G::Vector{Gauging})

Resample the gauging for bootstrap.

### Details

Sampling with replacement is performed to compose a gauging sample of the same size.
The smallest gauging is always selected so that the rating curve is valid over the entire gauged range.
"""
function bootstrap(G::Vector{Gauging})

    n = length(G)

    G_bootstrap = sample(G, n-1, replace=true)

    pushfirst!(G_bootstrap, minimum(G))

    return G_bootstrap

end

"""
    crcfit(G::Vector{Gauging})

Fit the compound rating curve corresponding to the gaugings `G`
"""
function crcfit(G::Vector{Gauging})

    h = level.(G)
    q = discharge.(G)

    hs = sort(unique(h))

    fobj(k) = sum( (log.(q) - logdischarge.(crcfit(G,k), h)).^2 )
    res = optimize(fobj, hs[3]+1e-2, hs[end-2]-1e-2)
    k = Optim.minimizer(res)

    crc = crcfit(G, k)

    return crc

end

"""
    crcfit(G::Vector{Gauging},k::Real)

Fit the compound rating curve corresponding to the gaugings `G` using the break-point `k`.

### Details

At least 3 gaugings are necessary by component of the compound rating curve.
"""
function crcfit(G::Vector{Gauging}, k::Real)

    h = level.(G)
    q = discharge.(G)

    ind = h.>k

    h2 = h[ind]
    q2 = q[ind]
    G₂ = G[ind]

    h1 = h[.!(ind)]
    q1 = q[.!(ind)]
    G₁ = G[.!(ind)]

    rc₂ = rcfit(G₂)

    # Continuity constraint at h = k
    rc₁ = rcfit(G₁, [k, discharge(rc₂, k)])

    crc = CompoundRatingCurve([k], [rc₁, rc₂] )

    return crc

end

"""
    discharge(G::Gauging)

Return the discharge of gauging `G`
"""
function discharge(G::Gauging)
    return G.discharge
end

"""
    Compute the rough initial values of the rating curve parameters.

### Details
These are rough estimates obtained by setting `c = 3/2`, `b = .9*min(h)`. The parameter `a` is computed with the median level.
"""
function getinitialvalues(G::Vector{Gauging})

    h = level.(G)
    q = discharge.(G)

    c₀ = 1.5
    b₀ = .9*minimum(h)

    n = length(h)
    p = sortperm(h)

    index = p[floor(Int64, n/2)]

    a₀ = q[index] / (h[index] - b₀)^c₀

    rc = RatingCurve(G, a₀, b₀, c₀)

    return rc

end

"""
    level(G::Gauging)

Return the level of gauging `G`
"""
function level(G::Gauging)
    return G.level
end

"""
    minimum(G::Vector{Gauging})

Return the gauging with the smallest level.
"""
function minimum(G::Vector{Gauging})

    ind = argmin(level.(G))

    return G[ind]
end

"""
    maximum(G::Vector{Gauging})

Return the gauging with the highest level.
"""
function maximum(G::Vector{Gauging})

    ind = argmax(level.(G))

    return G[ind]
end

"""
    rcfit(G::Vector{Gauging})

Fit the rating curve to the gaugings `G`.
"""
function rcfit(G::Vector{Gauging})

    h = level.(G)
    q = discharge.(G)

    rc₀ = RatingCurves.getinitialvalues(G)

    y = log.(q)

    fobj(b) = sum( (y - logdischarge.(rcfit(G,b[]), h)).^2 )

    res = optimize(fobj, [-Inf], [minimum(h)-1e-2], [rc₀.b])
    b = Optim.minimizer(res)[]

    rc = rcfit(G, b)

    return rc

end

"""
    rcfit(G::Vector{Gauging}, b::Real)

Fit the rating curve of parameter `b` to the gaugings `G`.
"""
function rcfit(G::Vector{Gauging}, b::Real)

    h = level.(G)
    q = discharge.(G)

    # @assert b < minimum(h)

    x = log.(h .- b)
    y = log.(q)

    X = hcat(ones(length(x)), x)

    β̂ = X\y

    c = β̂[2]

    a = exp(β̂[1])

    rc = RatingCurve(G, a, b, c)

    return rc

end

"""
    rcfit(G::Vector{Gauging}, constraint::AbstractVector{<:Real})

Fit the rating curve to the gaugings `G` passing through the point (h,q) specified in `constraint`.
"""
function rcfit(G::Vector{Gauging}, constraint::AbstractVector{<:Real})

    h = level.(G)
    q = discharge.(G)

    y = log.(q)

    rc₀ = RatingCurves.getinitialvalues(G)

    fobj(b) = sum( (y - logdischarge.(rcfit(G,b[],constraint), h)).^2 )

    res = optimize(fobj, [-Inf], [minimum(h)-1e-2], [rc₀.b])
    b = Optim.minimizer(res)[]

    rc = rcfit(G, b, constraint)

    return rc

end

"""
    rcfit(G::Vector{Gauging}, b::Real, constraint::AbstractVector{<:Real})

Fit the rating curve with parameter `b` to the gaugings `G` passing through the point (h,q) specified in `constraint`.
"""
function rcfit(G::Vector{Gauging}, b::Real, constraint::AbstractVector{<:Real})

    h = level.(G)
    q = discharge.(G)

    # @assert b <= minimum(h)

    h̃ = constraint[1]
    q̃ = constraint[2]

    x̃ = log(h̃ - b)
    ỹ = log(q̃)

    x = log.(h .- b)
    y = log.(q)

    z = x .- x̃

    c = z\(y .- ỹ)

    a = exp(ỹ - x̃*c)

    rc = RatingCurve(G, a, b, c)

    return rc

end

"""
    sort(G::Vector{Gauging}; rev::Bool=false)

Sort the gauging according to the level.

### Details

Sort in ascending order if `rev=false` (option by default) and in descending order if `rev=true`.
"""
function sort(G::Vector{Gauging}; rev::Bool=false)

    h = level.(G)

    Gs = G[sortperm(h, rev=rev)]

    return Gs

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Generic IO stubs -- all subtypes should implement these (if meaningful)

"""
    EOFError()

No more data was available to read from a file or stream.
"""
struct EOFError <: Exception end

"""
    SystemError(prefix::AbstractString, [errno::Int32])

A system call failed with an error code (in the `errno` global variable).
"""
struct SystemError <: Exception
    prefix::AbstractString
    errnum::Int32
    extrainfo
    SystemError(p::AbstractString, e::Integer, extrainfo) = new(p, e, extrainfo)
    SystemError(p::AbstractString, e::Integer) = new(p, e, nothing)
    SystemError(p::AbstractString) = new(p, Libc.errno())
end

lock(::IO) = nothing
unlock(::IO) = nothing
reseteof(x::IO) = nothing

const SZ_UNBUFFERED_IO = 65536
buffer_writes(x::IO, bufsize=SZ_UNBUFFERED_IO) = x

"""
    isopen(object) -> Bool

Determine whether an object - such as a stream or timer
-- is not yet closed. Once an object is closed, it will never produce a new event.
However, since a closed stream may still have data to read in its buffer,
use [`eof`](@ref) to check for the ability to read data.
Use the `FileWatching` package to be notified when a stream might be writable or readable.

# Examples
```jldoctest
julia> io = open("my_file.txt", "w+");

julia> isopen(io)
true

julia> close(io)

julia> isopen(io)
false
```
"""
function isopen end

"""
    close(stream)

Close an I/O stream. Performs a [`flush`](@ref) first.
"""
function close end
function flush end
function wait_connected end
function wait_readnb end
function wait_readbyte end
function wait_close end
function bytesavailable end

"""
    readavailable(stream)

Read all available data on the stream, blocking the task only if no data is available. The
result is a `Vector{UInt8,1}`.
"""
function readavailable end

"""
    isreadable(io) -> Bool

Return `true` if the specified IO object is readable (if that can be determined).

# Examples
```jldoctest
julia> open("myfile.txt", "w") do io
           print(io, "Hello world!");
           isreadable(io)
       end
false

julia> open("myfile.txt", "r") do io
           isreadable(io)
       end
true

julia> rm("myfile.txt")
```
"""
function isreadable end

"""
    iswritable(io) -> Bool

Return `true` if the specified IO object is writable (if that can be determined).

# Examples
```jldoctest
julia> open("myfile.txt", "w") do io
           print(io, "Hello world!");
           iswritable(io)
       end
true

julia> open("myfile.txt", "r") do io
           iswritable(io)
       end
false

julia> rm("myfile.txt")
```
"""
function iswritable end
function copy end
function eof end

"""
    read(io::IO, T)

Read a single value of type `T` from `io`, in canonical binary representation.

    read(io::IO, String)

Read the entirety of `io`, as a `String`.

# Examples
```jldoctest
julia> io = IOBuffer("JuliaLang is a GitHub organization");

julia> read(io, Char)
'J': ASCII/Unicode U+004a (category Lu: Letter, uppercase)

julia> io = IOBuffer("JuliaLang is a GitHub organization");

julia> read(io, String)
"JuliaLang is a GitHub organization"
```
"""
read(stream, t)

"""
    write(io::IO, x)
    write(filename::AbstractString, x)

Write the canonical binary representation of a value to the given I/O stream or file.
Return the number of bytes written into the stream.   See also [`print`](@ref) to
write a text representation (with an encoding that may depend upon `io`).

You can write multiple values with the same `write` call. i.e. the following are equivalent:

    write(io, x, y...)
    write(io, x) + write(io, y...)

# Examples
```jldoctest
julia> io = IOBuffer();

julia> write(io, "JuliaLang is a GitHub organization.", " It has many members.")
56

julia> String(take!(io))
"JuliaLang is a GitHub organization. It has many members."

julia> write(io, "Sometimes those members") + write(io, " write documentation.")
44

julia> String(take!(io))
"Sometimes those members write documentation."
```
"""
function write end

read(s::IO, ::Type{UInt8}) = error(typeof(s)," does not support byte I/O")
write(s::IO, x::UInt8) = error(typeof(s)," does not support byte I/O")

"""
    unsafe_write(io::IO, ref, nbytes::UInt)

Copy `nbytes` from `ref` (converted to a pointer) into the `IO` object.

It is recommended that subtypes `T<:IO` override the following method signature
to provide more efficient implementations:
`unsafe_write(s::T, p::Ptr{UInt8}, n::UInt)`
"""
function unsafe_write(s::IO, p::Ptr{UInt8}, n::UInt)
    written::Int = 0
    for i = 1:n
        written += write(s, unsafe_load(p, i))
    end
    return written
end

"""
    unsafe_read(io::IO, ref, nbytes::UInt)

Copy `nbytes` from the `IO` stream object into `ref` (converted to a pointer).

It is recommended that subtypes `T<:IO` override the following method signature
to provide more efficient implementations:
`unsafe_read(s::T, p::Ptr{UInt8}, n::UInt)`
"""
function unsafe_read(s::IO, p::Ptr{UInt8}, n::UInt)
    for i = 1:n
        unsafe_store!(p, read(s, UInt8)::UInt8, i)
    end
    nothing
end

# Generic wrappers around other IO objects
abstract type AbstractPipe <: IO end
function pipe_reader end
function pipe_writer end

write(io::AbstractPipe, byte::UInt8) = write(pipe_writer(io), byte)
unsafe_write(io::AbstractPipe, p::Ptr{UInt8}, nb::UInt) = unsafe_write(pipe_writer(io), p, nb)
buffer_writes(io::AbstractPipe, args...) = buffer_writes(pipe_writer(io), args...)
flush(io::AbstractPipe) = flush(pipe_writer(io))

read(io::AbstractPipe, byte::Type{UInt8}) = read(pipe_reader(io), byte)
unsafe_read(io::AbstractPipe, p::Ptr{UInt8}, nb::UInt) = unsafe_read(pipe_reader(io), p, nb)
read(io::AbstractPipe) = read(pipe_reader(io))
readuntil(io::AbstractPipe, arg::UInt8; kw...) = readuntil(pipe_reader(io), arg; kw...)
readuntil(io::AbstractPipe, arg::AbstractChar; kw...) = readuntil(pipe_reader(io), arg; kw...)
readuntil(io::AbstractPipe, arg::AbstractString; kw...) = readuntil(pipe_reader(io), arg; kw...)
readuntil(io::AbstractPipe, arg::AbstractVector; kw...) = readuntil(pipe_reader(io), arg; kw...)
readuntil_vector!(io::AbstractPipe, target::AbstractVector, keep::Bool, out) = readuntil_vector!(pipe_reader(io), target, keep, out)

for f in (
        # peek/mark interface
        :peek, :mark, :unmark, :reset, :ismarked,
        # Simple reader functions
        :readavailable, :isreadable)
    @eval $(f)(io::AbstractPipe) = $(f)(pipe_reader(io))
end

iswritable(io::AbstractPipe) = iswritable(pipe_writer(io))
isopen(io::AbstractPipe) = isopen(pipe_writer(io)) || isopen(pipe_reader(io))
close(io::AbstractPipe) = (close(pipe_writer(io)); close(pipe_reader(io)))
wait_readnb(io::AbstractPipe, nb::Int) = wait_readnb(pipe_reader(io), nb)
wait_readbyte(io::AbstractPipe, byte::UInt8) = wait_readbyte(pipe_reader(io), byte)
wait_close(io::AbstractPipe) = (wait_close(pipe_writer(io)); wait_close(pipe_reader(io)))

"""
    bytesavailable(io)

Return the number of bytes available for reading before a read from this stream or buffer will block.

# Examples
```jldoctest
julia> io = IOBuffer("JuliaLang is a GitHub organization");

julia> bytesavailable(io)
34
```
"""
bytesavailable(io::AbstractPipe) = bytesavailable(pipe_reader(io))

"""
    eof(stream) -> Bool

Test whether an I/O stream is at end-of-file. If the stream is not yet exhausted, this
function will block to wait for more data if necessary, and then return `false`. Therefore
it is always safe to read one byte after seeing `eof` return `false`. `eof` will return
`false` as long as buffered data is still available, even if the remote end of a connection
is closed.
"""
eof(io::AbstractPipe) = eof(pipe_reader(io))
reseteof(io::AbstractPipe) = reseteof(pipe_reader(io))

# Exception-safe wrappers (io = open(); try f(io) finally close(io))

write(filename::AbstractString, a1, args...) = open(io->write(io, a1, args...), filename, "w")

"""
    read(filename::AbstractString, args...)

Open a file and read its contents. `args` is passed to `read`: this is equivalent to
`open(io->read(io, args...), filename)`.

    read(filename::AbstractString, String)

Read the entire contents of a file as a string.
"""
read(filename::AbstractString, args...) = open(io->read(io, args...), filename)

read(filename::AbstractString, ::Type{T}) where {T} = open(io->read(io, T), filename)

"""
    read!(stream::IO, array::Union{Array, BitArray})
    read!(filename::AbstractString, array::Union{Array, BitArray})

Read binary data from an I/O stream or file, filling in `array`.
"""
function read! end

read!(filename::AbstractString, a) = open(io->read!(io, a), filename)

"""
    readuntil(stream::IO, delim; keep::Bool = false)
    readuntil(filename::AbstractString, delim; keep::Bool = false)

Read a string from an I/O stream or a file, up to the given delimiter.
The delimiter can be a `UInt8`, `AbstractChar`, string, or vector.
Keyword argument `keep` controls whether the delimiter is included in the result.
The text is assumed to be encoded in UTF-8.

# Examples
```jldoctest
julia> open("my_file.txt", "w") do io
           write(io, "JuliaLang is a GitHub organization.\\nIt has many members.\\n");
       end
57

julia> readuntil("my_file.txt", 'L')
"Julia"

julia> readuntil("my_file.txt", '.', keep = true)
"JuliaLang is a GitHub organization."

julia> rm("my_file.txt")
```
"""
readuntil(filename::AbstractString, args...; kw...) = open(io->readuntil(io, args...; kw...), filename)

"""
    readline(io::IO=stdin; keep::Bool=false)
    readline(filename::AbstractString; keep::Bool=false)

Read a single line of text from the given I/O stream or file (defaults to `stdin`).
When reading from a file, the text is assumed to be encoded in UTF-8. Lines in the
input end with `'\\n'` or `"\\r\\n"` or the end of an input stream. When `keep` is
false (as it is by default), these trailing newline characters are removed from the
line before it is returned. When `keep` is true, they are returned as part of the
line.

# Examples
```jldoctest
julia> open("my_file.txt", "w") do io
           write(io, "JuliaLang is a GitHub organization.\\nIt has many members.\\n");
       end
57

julia> readline("my_file.txt")
"JuliaLang is a GitHub organization."

julia> readline("my_file.txt", keep=true)
"JuliaLang is a GitHub organization.\\n"

julia> rm("my_file.txt")
```
"""
function readline(filename::AbstractString; keep::Bool=false)
    open(filename) do f
        readline(f, keep=keep)
    end
end

function readline(s::IO=stdin; keep::Bool=false)
    line = readuntil(s, 0x0a, keep=true)
    i = length(line)
    if keep || i == 0 || line[i] != 0x0a
        return String(line)
    elseif i < 2 || line[i-1] != 0x0d
        return String(resize!(line,i-1))
    else
        return String(resize!(line,i-2))
    end
end

"""
    readlines(io::IO=stdin; keep::Bool=false)
    readlines(filename::AbstractString; keep::Bool=false)

Read all lines of an I/O stream or a file as a vector of strings. Behavior is
equivalent to saving the result of reading [`readline`](@ref) repeatedly with the same
arguments and saving the resulting lines as a vector of strings.

# Examples
```jldoctest
julia> open("my_file.txt", "w") do io
           write(io, "JuliaLang is a GitHub organization.\\nIt has many members.\\n");
       end
57

julia> readlines("my_file.txt")
2-element Array{String,1}:
 "JuliaLang is a GitHub organization."
 "It has many members."

julia> readlines("my_file.txt", keep=true)
2-element Array{String,1}:
 "JuliaLang is a GitHub organization.\\n"
 "It has many members.\\n"

julia> rm("my_file.txt")
```
"""
function readlines(filename::AbstractString; kw...)
    open(filename) do f
        readlines(f; kw...)
    end
end
readlines(s=stdin; kw...) = collect(eachline(s; kw...))

## byte-order mark, ntoh & hton ##

let a = UInt32[0x01020304]
    endian_bom = GC.@preserve a unsafe_load(convert(Ptr{UInt8}, pointer(a)))
    global ntoh, hton, ltoh, htol
    if endian_bom == 0x01
        ntoh(x) = x
        hton(x) = x
        ltoh(x) = bswap(x)
        htol(x) = bswap(x)
        const global ENDIAN_BOM = 0x01020304
    elseif endian_bom == 0x04
        ntoh(x) = bswap(x)
        hton(x) = bswap(x)
        ltoh(x) = x
        htol(x) = x
        const global ENDIAN_BOM = 0x04030201
    else
        error("seriously? what is this machine?")
    end
end

"""
    ENDIAN_BOM

The 32-bit byte-order-mark indicates the native byte order of the host machine.
Little-endian machines will contain the value `0x04030201`. Big-endian machines will contain
the value `0x01020304`.
"""
ENDIAN_BOM

"""
    ntoh(x)

Convert the endianness of a value from Network byte order (big-endian) to that used by the Host.
"""
ntoh(x)

"""
    hton(x)

Convert the endianness of a value from that used by the Host to Network byte order (big-endian).
"""
hton(x)

"""
    ltoh(x)

Convert the endianness of a value from Little-endian to that used by the Host.
"""
ltoh(x)

"""
    htol(x)

Convert the endianness of a value from that used by the Host to Little-endian.
"""
htol(x)

"""
    isreadonly(io) -> Bool

Determine whether a stream is read-only.

# Examples
```jldoctest
julia> io = IOBuffer("JuliaLang is a GitHub organization");

julia> isreadonly(io)
true

julia> io = IOBuffer();

julia> isreadonly(io)
false
```
"""
isreadonly(s) = isreadable(s) && !iswritable(s)

## binary I/O ##

write(io::IO, x) = throw(MethodError(write, (io, x)))
function write(io::IO, x1, xs...)
    written::Int = write(io, x1)
    for x in xs
        written += write(io, x)
    end
    return written
end

@noinline unsafe_write(s::IO, p::Ref{T}, n::Integer) where {T} =
    unsafe_write(s, unsafe_convert(Ref{T}, p)::Ptr, n) # mark noinline to ensure ref is gc-rooted somewhere (by the caller)
unsafe_write(s::IO, p::Ptr, n::Integer) = unsafe_write(s, convert(Ptr{UInt8}, p), convert(UInt, n))
write(s::IO, x::Ref{T}) where {T} = unsafe_write(s, x, Core.sizeof(T))
write(s::IO, x::Int8) = write(s, reinterpret(UInt8, x))
function write(s::IO, x::Union{Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,Float16,Float32,Float64})
    return write(s, Ref(x))
end

write(s::IO, x::Bool) = write(s, UInt8(x))
write(to::IO, p::Ptr) = write(to, convert(UInt, p))

function write(s::IO, A::AbstractArray)
    if !isbitstype(eltype(A))
        error("`write` is not supported on non-isbits arrays")
    end
    nb = 0
    for a in A
        nb += write(s, a)
    end
    return nb
end

function write(s::IO, a::Array)
    if isbitstype(eltype(a))
        return GC.@preserve a unsafe_write(s, pointer(a), sizeof(a))
    else
        error("`write` is not supported on non-isbits arrays")
    end
end

function write(s::IO, a::SubArray{T,N,<:Array}) where {T,N}
    if !isbitstype(T) || !isa(a, StridedArray)
        return invoke(write, Tuple{IO, AbstractArray}, s, a)
    end
    elsz = sizeof(T)
    colsz = size(a,1) * elsz
    GC.@preserve a if stride(a,1) != 1
        for idxs in CartesianIndices(size(a))
            unsafe_write(s, pointer(a, idxs.I), elsz)
        end
        return elsz * length(a)
    elseif N <= 1
        return unsafe_write(s, pointer(a, 1), colsz)
    else
        for idxs in CartesianIndices((1, size(a)[2:end]...))
            unsafe_write(s, pointer(a, idxs.I), colsz)
        end
        return colsz * trailingsize(a,2)
    end
end

function write(io::IO, c::Char)
    u = bswap(reinterpret(UInt32, c))
    n = 1
    while true
        write(io, u % UInt8)
        (u >>= 8) == 0 && return n
        n += 1
    end
end
# write(io, ::AbstractChar) is not defined: implementations
# must provide their own encoding-specific method.

function write(io::IO, s::Symbol)
    pname = unsafe_convert(Ptr{UInt8}, s)
    return unsafe_write(io, pname, Int(ccall(:strlen, Csize_t, (Cstring,), pname)))
end

function write(to::IO, from::IO)
    n = 0
    while !eof(from)
        n += write(to, readavailable(from))
    end
    return n
end

@noinline unsafe_read(s::IO, p::Ref{T}, n::Integer) where {T} = unsafe_read(s, unsafe_convert(Ref{T}, p)::Ptr, n) # mark noinline to ensure ref is gc-rooted somewhere (by the caller)
unsafe_read(s::IO, p::Ptr, n::Integer) = unsafe_read(s, convert(Ptr{UInt8}, p), convert(UInt, n))
read!(s::IO, x::Ref{T}) where {T} = (unsafe_read(s, x, Core.sizeof(T)); x)

read(s::IO, ::Type{Int8}) = reinterpret(Int8, read(s, UInt8))
function read(s::IO, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    return read!(s, Ref{T}(0))[]::T
end

read(s::IO, ::Type{Bool}) = (read(s, UInt8) != 0)
read(s::IO, ::Type{Ptr{T}}) where {T} = convert(Ptr{T}, read(s, UInt))

function read!(s::IO, a::Array{UInt8})
    GC.@preserve a unsafe_read(s, pointer(a), sizeof(a))
    return a
end

function read!(s::IO, a::Array{T}) where T
    if isbitstype(T)
        GC.@preserve a unsafe_read(s, pointer(a), sizeof(a))
    else
        for i in eachindex(a)
            a[i] = read(s, T)
        end
    end
    return a
end

function read(io::IO, ::Type{Char})
    b0 = read(io, UInt8)
    l = 8(4-leading_ones(b0))
    c = UInt32(b0) << 24
    if l < 24
        s = 16
        while s ≥ l && !eof(io)
            peek(io) & 0xc0 == 0x80 || break
            b = read(io, UInt8)
            c |= UInt32(b) << s
            s -= 8
        end
    end
    return reinterpret(Char, c)
end
# read(io, T) is not defined for other AbstractChar: implementations
# must provide their own encoding-specific method.

# readuntil_string is useful below since it has
# an optimized method for s::IOStream
readuntil_string(s::IO, delim::UInt8, keep::Bool) = String(readuntil(s, delim, keep=keep))

function readuntil(s::IO, delim::AbstractChar; keep::Bool=false)
    if delim ≤ '\x7f'
        return readuntil_string(s, delim % UInt8, keep)
    end
    out = IOBuffer()
    while !eof(s)
        c = read(s, Char)
        if c == delim
            keep && write(out, c)
            break
        end
        write(out, c)
    end
    return String(take!(out))
end

function readuntil(s::IO, delim::T; keep::Bool=false) where T
    out = (T === UInt8 ? StringVector(0) : Vector{T}())
    while !eof(s)
        c = read(s, T)
        if c == delim
            keep && push!(out, c)
            break
        end
        push!(out, c)
    end
    return out
end

# requires that indices for target are the integer unit range from firstindex to lastindex
# returns whether the delimiter was matched
# uses the Knuth–Morris–Pratt_algorithm, with the first and second cache entries unrolled
# For longer targets, the cache improves the big-O efficiency of scanning of sequences
# with repeated patterns
# Each cache entry tells us which index we should start the search at.
# We assume this is unlikely, so we only lazy-initialize as much of the cache as we need to use
# When we allocate the cache, we initialize it to 0 (and offset by the first index afterwards).
# Suppose target is:
#    Index:  1245689
#    Value: "aδcaδcx"
# We would set the cache to
#    0 0 0 1 2 3 4 0
# So after if we mismatch after the second aδc sequence,
# we can immediately jump back to index 5 (4 + 1).
function readuntil_vector!(io::IO, target::AbstractVector{T}, keep::Bool, out) where {T}
    first = firstindex(target)
    last = lastindex(target)
    len = last - first + 1
    if len < 1
        return true
    end
    pos = 0 # array-offset
    max_pos = 1 # array-offset in cache
    local cache # will be lazy initialized when needed
    output! = (isa(out, IO) ? write : push!)
    while !eof(io)
        c = read(io, T)
        # Backtrack until the next target character matches what was found
        while true
            c1 = target[pos + first]
            if c == c1
                pos += 1
                break
            elseif pos == 0
                break
            elseif pos == 1
                if !keep
                    output!(out, target[first])
                end
                pos = 0
            else
                # grow cache to contain up to `pos`
                if !@isdefined(cache)
                    cache = zeros(Int, len)
                end
                while max_pos < pos
                    ci = target[max_pos + first]
                    b = max_pos
                    max_pos += 1
                    while b != 0
                        b = cache[b]
                        cb = target[b + first]
                        if ci == cb
                            cache[max_pos] = b + 1
                            break
                        end
                    end
                end
                # read new position from cache
                pos1 = cache[pos]
                if !keep
                    # and add the removed prefix from the target to the output
                    # if not always keeping the match
                    for b in 1:(pos - pos1)
                        output!(out, target[b - 1 + first])
                    end
                end
                pos = pos1
            end
        end
        if keep || pos == 0
            output!(out, c)
        end
        pos == len && return true
    end
    if !keep
        # failed early without finishing the match,
        # add the partial match to the output
        # if not always keeping the match
        for b in 1:pos
            output!(out, target[b - 1 + first])
        end
    end
    return false
end

function readuntil(io::IO, target::AbstractString; keep::Bool=false)
    # small-string target optimizations
    isempty(target) && return ""
    c, rest = Iterators.peel(target)
    if isempty(rest) && c <= '\x7f'
        return readuntil_string(io, c % UInt8, keep)
    end
    # convert String to a utf8-byte-iterator
    if !(target isa String) && !(target isa SubString{String})
        target = String(target)
    end
    target = codeunits(target)::AbstractVector
    return String(readuntil(io, target, keep=keep))
end

function readuntil(io::IO, target::AbstractVector{T}; keep::Bool=false) where T
    out = (T === UInt8 ? StringVector(0) : Vector{T}())
    readuntil_vector!(io, target, keep, out)
    return out
end

"""
    readchomp(x)

Read the entirety of `x` as a string and remove a single trailing newline
if there is one. Equivalent to `chomp(read(x, String))`.

# Examples
```jldoctest
julia> open("my_file.txt", "w") do io
           write(io, "JuliaLang is a GitHub organization.\\nIt has many members.\\n");
       end;

julia> readchomp("my_file.txt")
"JuliaLang is a GitHub organization.\\nIt has many members."

julia> rm("my_file.txt");
```
"""
readchomp(x) = chomp(read(x, String))

# read up to nb bytes into nb, returning # bytes read

"""
    readbytes!(stream::IO, b::AbstractVector{UInt8}, nb=length(b))

Read at most `nb` bytes from `stream` into `b`, returning the number of bytes read.
The size of `b` will be increased if needed (i.e. if `nb` is greater than `length(b)`
and enough bytes could be read), but it will never be decreased.
"""
function readbytes!(s::IO, b::AbstractArray{UInt8}, nb=length(b))
    require_one_based_indexing(b)
    olb = lb = length(b)
    nr = 0
    while nr < nb && !eof(s)
        a = read(s, UInt8)
        nr += 1
        if nr > lb
            lb = nr * 2
            resize!(b, lb)
        end
        b[nr] = a
    end
    if lb > olb
        resize!(b, nr) # shrink to just contain input data if was resized
    end
    return nr
end

"""
    read(s::IO, nb=typemax(Int))

Read at most `nb` bytes from `s`, returning a `Vector{UInt8}` of the bytes read.
"""
function read(s::IO, nb::Integer = typemax(Int))
    # Let readbytes! grow the array progressively by default
    # instead of taking of risk of over-allocating
    b = Vector{UInt8}(undef, nb == typemax(Int) ? 1024 : nb)
    nr = readbytes!(s, b, nb)
    return resize!(b, nr)
end

read(s::IO, ::Type{String}) = String(read(s))
read(s::IO, T::Type) = error("The IO stream does not support reading objects of type $T.")

## high-level iterator interfaces ##

struct EachLine{IOT <: IO}
    stream::IOT
    ondone::Function
    keep::Bool
    EachLine(stream::IO=stdin; ondone::Function=()->nothing, keep::Bool=false) =
        new{typeof(stream)}(stream, ondone, keep)
end

"""
    eachline(io::IO=stdin; keep::Bool=false)
    eachline(filename::AbstractString; keep::Bool=false)

Create an iterable `EachLine` object that will yield each line from an I/O stream
or a file. Iteration calls [`readline`](@ref) on the stream argument repeatedly with
`keep` passed through, determining whether trailing end-of-line characters are
retained. When called with a file name, the file is opened once at the beginning of
iteration and closed at the end. If iteration is interrupted, the file will be
closed when the `EachLine` object is garbage collected.

# Examples
```jldoctest
julia> open("my_file.txt", "w") do io
           write(io, "JuliaLang is a GitHub organization.\\n It has many members.\\n");
       end;

julia> for line in eachline("my_file.txt")
           print(line)
       end
JuliaLang is a GitHub organization. It has many members.

julia> rm("my_file.txt");
```
"""
function eachline(stream::IO=stdin; keep::Bool=false)
    EachLine(stream, keep=keep)::EachLine
end

function eachline(filename::AbstractString; keep::Bool=false)
    s = open(filename)
    EachLine(s, ondone=()->close(s), keep=keep)::EachLine
end

function iterate(itr::EachLine, state=nothing)
    eof(itr.stream) && return (itr.ondone(); nothing)
    (readline(itr.stream, keep=itr.keep), nothing)
end

eltype(::Type{<:EachLine}) = String

IteratorSize(::Type{<:EachLine}) = SizeUnknown()

# IOStream Marking
# Note that these functions expect that io.mark exists for
# the concrete IO type. This may not be true for IO types
# not in base.

"""
    mark(s)

Add a mark at the current position of stream `s`. Return the marked position.

See also [`unmark`](@ref), [`reset`](@ref), [`ismarked`](@ref).
"""
function mark(io::IO)
    io.mark = position(io)
end

"""
    unmark(s)

Remove a mark from stream `s`. Return `true` if the stream was marked, `false` otherwise.

See also [`mark`](@ref), [`reset`](@ref), [`ismarked`](@ref).
"""
function unmark(io::IO)
    !ismarked(io) && return false
    io.mark = -1
    return true
end

"""
    reset(s)

Reset a stream `s` to a previously marked position, and remove the mark. Return the
previously marked position. Throw an error if the stream is not marked.

See also [`mark`](@ref), [`unmark`](@ref), [`ismarked`](@ref).
"""
function reset(io::T) where T<:IO
    ismarked(io) || throw(ArgumentError("$(T) not marked"))
    m = io.mark
    seek(io, m)
    io.mark = -1 # must be after seek, or seek may fail
    return m
end

"""
    ismarked(s)

Return `true` if stream `s` is marked.

See also [`mark`](@ref), [`unmark`](@ref), [`reset`](@ref).
"""
ismarked(io::IO) = io.mark >= 0

# Make sure all IO streams support flush, even if only as a no-op,
# to make it easier to write generic I/O code.

"""
    flush(stream)

Commit all currently buffered writes to the given stream.
"""
flush(io::IO) = nothing

"""
    skipchars(predicate, io::IO; linecomment=nothing)

Advance the stream `io` such that the next-read character will be the first remaining for
which `predicate` returns `false`. If the keyword argument `linecomment` is specified, all
characters from that character until the start of the next line are ignored.

# Examples
```jldoctest
julia> buf = IOBuffer("    text")
IOBuffer(data=UInt8[...], readable=true, writable=false, seekable=true, append=false, size=8, maxsize=Inf, ptr=1, mark=-1)

julia> skipchars(isspace, buf)
IOBuffer(data=UInt8[...], readable=true, writable=false, seekable=true, append=false, size=8, maxsize=Inf, ptr=5, mark=-1)

julia> String(readavailable(buf))
"text"
```
"""
function skipchars(predicate, io::IO; linecomment=nothing)
    while !eof(io)
        c = read(io, Char)
        if c === linecomment
            readline(io)
        elseif !predicate(c)
            skip(io, -ncodeunits(c))
            break
        end
    end
    return io
end

"""
    countlines(io::IO; eol::AbstractChar = '\\n')

Read `io` until the end of the stream/file and count the number of lines. To specify a file
pass the filename as the first argument. EOL markers other than `'\\n'` are supported by
passing them as the second argument.  The last non-empty line of `io` is counted even if it does not
end with the EOL, matching the length returned by [`eachline`](@ref) and [`readlines`](@ref).

# Examples
```jldoctest
julia> io = IOBuffer("JuliaLang is a GitHub organization.\\n");

julia> countlines(io)
1

julia> io = IOBuffer("JuliaLang is a GitHub organization.");

julia> countlines(io)
1

julia> countlines(io, eol = '.')
0
```
"""
function countlines(io::IO; eol::AbstractChar='\n')
    isascii(eol) || throw(ArgumentError("only ASCII line terminators are supported"))
    aeol = UInt8(eol)
    a = Vector{UInt8}(undef, 8192)
    nl = nb = 0
    while !eof(io)
        nb = readbytes!(io, a)
        @simd for i=1:nb
            @inbounds nl += a[i] == aeol
        end
    end
    if nb > 0 && a[nb] != aeol
        nl += 1 # final line is not terminated with eol
    end
    nl
end

countlines(f::AbstractString; eol::AbstractChar = '\n') = open(io->countlines(io, eol = eol), f)::Int

ext: jl

max_issues_repo_name: josecarlosmp/julia

max_forks_repo_name: josecarlosmp/julia


text:
# using Statistics, StatsBase
# using PyCall
# xr = pyimport("xarray")
# using PyPlot
#
# path = "/Users/milan/cams"
# allvars = ["no2","go3","so2","aermr04","aermr05","aermr06","ch4","co"]
#
# M = Array{Float32,3}(undef,length(allvars),6,137)
#
# for (i,vari) in enumerate(allvars)
#     filelist = filter(x->endswith(x,vari*".grib"),readdir(path))
#     X = xr.open_dataset(joinpath(path,filelist[1]),engine="cfgrib")[vari].data
#
#     M[i,1,:] = mean(X,dims=2)[:,1]
#     M[i,2,:] = median(X,dims=2)[:,1]
#     M[i,3,:] = minimum(X,dims=2)[:,1]
#     M[i,4,:] = maximum(X,dims=2)[:,1]
#     M[i,5,:] = [percentile(vec(X[i,:]),10) for i in level]
#     M[i,6,:] = [percentile(vec(X[i,:]),90) for i in level]
# end

# plotting
level = 1:137
alphabet = ["a","b","c","d","e","f","g","h","i","j"]

fig,axs = subplots(2,4,figsize=(10,8),sharey=true)
axs[1,1].invert_yaxis()

for (i,ax) in enumerate(axs)
    ax.semilogx(M[i,1,:],level,"C0",lw=3,label="mean")
    ax.plot(M[i,2,:],level,"C0",lw=1.5,label="median")
    ax.fill_betweenx(level,M[i,3,:],M[i,4,:],color="C0",alpha=0.2,label="Min-max range")
    ax.fill_betweenx(level,M[i,5,:],M[i,6,:],color="C0",alpha=0.5,label="10-90% range")
    ax.set_title(allvars[i],loc="left")
    ax.set_title(alphabet[i],loc="right",fontweight="bold")
end

axs[2,1].legend(loc=2)

axs[1,1].set_ylim(137,1)
axs[1,1].set_ylabel("model level")
axs[2,1].set_ylabel("model level")
axs[2,1].set_xlabel("values")
axs[2,2].set_xlabel("values")
axs[2,3].set_xlabel("values")
axs[2,4].set_xlabel("values")

tight_layout()


text:
using Pkg, Libdl

function shared_lib_filename()
    if Sys.iswindows()
        return "blpapi3_64.dll"
    elseif Sys.islinux() || Sys.isapple()
        return "libblpapi3_64.so"
    end
    error("Unsupported system.")
end

libblpapi3_path = "" # global var with path to the shared library
libblpapi3_handle = C_NULL # global var to be set at __init__() -> check_deps()
const libblpapi3 = shared_lib_filename()

function get_artifact_dir()
    if VERSION < v"1.6"
        #artifact_dir = abspath(artifact"blpapi") # didn't work
        artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
        @assert isfile(artifacts_toml) "Couldn't find $artifacts_toml"
        artifact_dict = Pkg.Artifacts.load_artifacts_toml(artifacts_toml)
        return Pkg.Artifacts.do_artifact_str("blpapi", artifact_dict, artifacts_toml, @__MODULE__)
    else
        return Pkg.Artifacts.@artifact_str("blpapi")
    end
end

get_libblpapi_artifact_filepath() = joinpath(get_artifact_dir(), libblpapi3)

# called by __init__()
function check_deps()

    global libblpapi3_path
    global libblpapi3_handle

    if Sys.islinux()
        try
            libblpapi3_handle = dlopen(libblpapi3)
        catch err
            if isa(err, ErrorException) && endswith(err.msg, "$libblpapi3: cannot open shared object file: No such file or directory")
                @warn("### WARNING ### \nCouldn't find $libblpapi3 shared lib.\nCopy the file $(get_libblpapi_artifact_filepath()) to your `LD_LIBRARY_PATH` and restart Julia.")
                return
            else
                rethrow(err)
            end
        end

        if libblpapi3_handle == C_NULL
            error("$libblpapi3 cannot be opened.")
        end
    else
        libblpapi3_path = get_libblpapi_artifact_filepath()
        if !isfile(libblpapi3_path)
            error("Couldn't find shared library artifact at $libblpapi3_path.")
        end

        libblpapi3_handle = dlopen(libblpapi3_path)
        if libblpapi3_handle == C_NULL
            error("$libblpapi3_path cannot be opened.")
        end
    end
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
#################################
### Convert LLA Bounds to ENU ###
#################################

# there's not an unambiguous conversion, but for now,
# returning the minimum bounds that contain all points contained
# by the input bounds
function ENU(bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, lla_ref::OpenStreetMapX.LLA = OpenStreetMapX.center(bounds), datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)

    max_x = max_y = -Inf
    min_x = min_y = Inf

    xs = [bounds.min_x, bounds.max_x]
    ys = [bounds.min_y, bounds.max_y]
    if bounds.min_y < 0.0 < bounds.max_y
        push!(ys, 0.0)
    end
    ref_x = OpenStreetMapX.getX(lla_ref)
    if bounds.min_x < ref_x < bounds.max_x ||
       (bounds.min_x > bounds.max_x && !(bounds.min_x >= ref_x >= bounds.max_x))
        push!(xs, ref_x)
    end

    for x_lla in xs, y_lla in ys
        pt = OpenStreetMapX.ENU(OpenStreetMapX.LLA(y_lla, x_lla), lla_ref, datum)
        x, y = OpenStreetMapX.getX(pt), OpenStreetMapX.getY(pt)

        min_x, max_x = min(x, min_x), max(x, max_x)
        min_y, max_y = min(y, min_y), max(y, max_y)
    end

    return OpenStreetMapX.Bounds{OpenStreetMapX.ENU}(min_y, max_y, min_x, max_x)
end

#########################################
### Get Center Point of Bounds Region ###
#########################################

function center(bounds::OpenStreetMapX.Bounds{OpenStreetMapX.ENU})
    x_mid = (bounds.min_x + bounds.max_x) / 2
    y_mid = (bounds.min_y + bounds.max_y) / 2

    return OpenStreetMapX.ENU(x_mid, y_mid)
end

function center(bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA})
    x_mid = (bounds.min_x + bounds.max_x) / 2
    y_mid = (bounds.min_y + bounds.max_y) / 2

    if bounds.min_x > bounds.max_x
        x_mid = x_mid > 0 ? x_mid - 180 : x_mid + 180
    end

    return OpenStreetMapX.LLA(y_mid, x_mid)
end

#################################################
### Check Whether a Location is Within Bounds ###
#################################################

function inbounds(loc::OpenStreetMapX.ENU, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.ENU})
    x, y = OpenStreetMapX.getX(loc), OpenStreetMapX.getY(loc)
    bounds.min_x <= x <= bounds.max_x &&
    bounds.min_y <= y <= bounds.max_y
end

function inbounds(loc::OpenStreetMapX.LLA, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA})
    x, y = OpenStreetMapX.getX(loc), OpenStreetMapX.getY(loc)
    min_x, max_x = bounds.min_x, bounds.max_x
    (min_x > max_x ? !(max_x < x < min_x) : min_x <= x <= max_x) &&
    bounds.min_y <= y <= bounds.max_y
end

# only for points that have passed the inbounds test
function onbounds(loc::T, bounds::OpenStreetMapX.Bounds{T}) where T<:Union{OpenStreetMapX.LLA,OpenStreetMapX.ENU}
    x, y = OpenStreetMapX.getX(loc), OpenStreetMapX.getY(loc)
    x == bounds.min_x || x == bounds.max_x ||
    y == bounds.min_y || y == bounds.max_y
end

#############################################
### Find the Closest Point  Within Bounds ###
#############################################

# only for points where inbounds(p1) != inbounds(p2)
function boundary_point(p1::T, p2::T, bounds::OpenStreetMapX.Bounds{T}) where T<:Union{OpenStreetMapX.LLA,OpenStreetMapX.ENU}
    x1, y1 = OpenStreetMapX.getX(p1), OpenStreetMapX.getY(p1)
    x2, y2 = OpenStreetMapX.getX(p2), OpenStreetMapX.getY(p2)

    x, y = Inf, Inf

    if bounds.min_x >  bounds.max_x && x1*x2 < 0

        if x1 < bounds.min_x && x2 < bounds.max_x || x2 < bounds.min_x && x1 < bounds.max_x
            x = bounds.min_x
            y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
        elseif x1 > bounds.max_x && x2 > bounds.min_x || x2 > bounds.max_x && x1 > bounds.min_x
            x = bounds.max_x
            y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
        end

        p3 = T(OpenStreetMapX.XY(x, y))
        OpenStreetMapX.inbounds(p3, bounds) && return p3
    end

    # Move x to x bound if segment crosses boundary
    if x1 < bounds.min_x < x2 || x1 > bounds.min_x > x2
        x = bounds.min_x
        y = y1 + (y2 - y1) * (bounds.min_x - x1) / (x2 - x1)
    elseif x1 < bounds.max_x < x2 || x1 > bounds.max_x > x2
        x = bounds.max_x
        y = y1 + (y2 - y1) * (bounds.max_x - x1) / (x2 - x1)
    end

    p3 = T(OpenStreetMapX.XY(x, y))
    OpenStreetMapX.inbounds(p3, bounds) && return p3

    # Move y to y bound if segment crosses boundary
    if y1 < bounds.min_y < y2 || y1 > bounds.min_y > y2
        x = x1 + (x2 - x1) * (bounds.min_y - y1) / (y2 - y1)
        y = bounds.min_y
    elseif y1 < bounds.max_y < y2 || y1 > bounds.max_y > y2
        x = x1 + (x2 - x1) * (bounds.max_y - y1) / (y2 - y1)
        y = bounds.max_y
    end

    p3 = T(OpenStreetMapX.XY(x, y))
    OpenStreetMapX.inbounds(p3, bounds) && return p3

    error("Failed to find boundary point.")
end


text:
@testset "Renewable DCPLossLess FullDispatch" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFullDispatch, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFullDispatch, PM.DCPlosslessForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

#@testset "Renewable Testing" begin
#ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
#PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFullDispatch, PM.StandardACPForm, sys5b, time_range);
#end

@testset "Renewable DCPLossLess ConstantPowerFactor" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableConstantPowerFactor, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableConstantPowerFactor, PM.DCPlosslessForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

@testset "Renewable ACP ConstantPowerFactor" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableConstantPowerFactor, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 72
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableConstantPowerFactor, PM.StandardACPForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 72
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 72
end

@testset "Renewable DCPLossLess FixedOutput" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFixed, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFixed, PM.DCPlosslessForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

@testset "Renewable ACP FixedOutput" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFixed, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.RenewableGen, PSI.RenewableFixed, PM.StandardACPForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end


text:
module TestPipelines2

using MLJBase
using Test
using ..Models
using StableRNGs
using Tables

rng = StableRNG(698790187)

@testset "helpers" begin
    @test MLJBase.individuate([:x, :y, :x, :z, :y, :x]) ==
        [:x, :y, :x2, :z, :y2, :x3]
end

# # DUMMY MODELS

mutable struct MyTransformer2 <: Static
    ftr::Symbol
end

MLJBase.transform(transf::MyTransformer2, verbosity, X) =
    fill(:st, nrows(X))

mutable struct MyDeterministic <: Deterministic
    x::Symbol
end

MLJBase.fit(::MyDeterministic, args...) = nothing, nothing, nothing
MLJBase.transform(m::MyDeterministic, ::Any, Xnew) = fill(:dt, nrows(Xnew))
MLJBase.predict(m::MyDeterministic, ::Any, Xnew) = fill(:dp, nrows(Xnew))
MLJBase.supports_training_losses(::Type{<:MyDeterministic}) = true
MLJBase.iteration_parameter(::Type{<:MyDeterministic}) = :x
MLJBase.training_losses(::MyDeterministic, report) = ones(3)

mutable struct MyProbabilistic <: Probabilistic
    x::Symbol
end

mutable struct MyUnsupervised <: Unsupervised
    x::Symbol
end

MLJBase.fit(::MyUnsupervised, args...) = nothing, nothing, nothing
MLJBase.transform(m::MyUnsupervised, ::Any, Xnew) = fill(:ut, nrows(Xnew))
MLJBase.predict(m::MyUnsupervised, ::Any, Xnew) = fill(:up, nrows(Xnew))

mutable struct MyInterval <: Interval
    x::Symbol
end

d = MyDeterministic(:d)
p = MyProbabilistic(:p)
u = MyUnsupervised(:u)
s = MyTransformer2(:s) # Static
i = MyInterval(:i)
m = MLJBase.matrix
t = MLJBase.table

@testset "pipe_named_tuple" begin
    @test_throws MLJBase.ERR_EMPTY_PIPELINE MLJBase.pipe_named_tuple((),())
    @test_throws(MLJBase.ERR_TOO_MANY_SUPERVISED,
                 MLJBase.pipe_named_tuple((:foo, :foo, :foo), (d, u, p)))
    _names      = (:trf, :fun, :fun, :trf, :clf)
    components = (u,    m,    t,    u,    d)
    @test MLJBase.pipe_named_tuple(_names, components) ==
        NamedTuple{(:trf, :fun, :fun2, :trf2, :clf),
                   Tuple{Unsupervised,
                         Any,
                         Any,
                         Unsupervised,
                         Deterministic}}(components)
end

@testset "public constructor" begin
    # un-named components:
    @test Pipeline(m, t, u) isa UnsupervisedPipeline
    @test Pipeline(m, t, u, p) isa ProbabilisticPipeline
    @test Pipeline(m, t, u, p, operation=predict_mean) isa DeterministicPipeline
    @test Pipeline(u, p, u, operation=predict_mean) isa DeterministicPipeline
    @test Pipeline(m, t) isa StaticPipeline
    @test Pipeline(m, t, s) isa StaticPipeline
    @test Pipeline(m, t, s, d) isa DeterministicPipeline
    @test Pipeline(m, t, i) isa IntervalPipeline
    @test_logs((:info, MLJBase.INFO_TREATING_AS_DETERMINISTIC),
               @test Pipeline(m, t, u, p, u) isa DeterministicPipeline)
    @test_logs((:info, MLJBase.INFO_TREATING_AS_DETERMINISTIC),
               @test Pipeline(m, t, u, i, u) isa DeterministicPipeline)

    # if "hidden" supervised model is already deterministic,
    # no need for warning:
    @test_logs @test Pipeline(m, t, u, d, u) isa DeterministicPipeline

    # named components:
    @test_throws MLJBase.ERR_USING_TARGET_KWARG Pipeline(target=u)
    @test Pipeline(c1=m, c2=t, c3=u) isa UnsupervisedPipeline
    @test Pipeline(c1=m, c2=t, c3=u, c5=p) isa ProbabilisticPipeline
    @test Pipeline(c1=m, c2=t) isa StaticPipeline
    @test Pipeline(c1=m, c2=t, c6=s) isa StaticPipeline
    @test Pipeline(c1=m, c2=t, c6=s, c7=d) isa DeterministicPipeline
    @test Pipeline(c1=m, c2=t, c8=i) isa IntervalPipeline
    @test_logs((:info, MLJBase.INFO_TREATING_AS_DETERMINISTIC),
               @test Pipeline(c1=m, c2=t, c3=u, c5=p, c4=u) isa
               DeterministicPipeline)
    @test(Pipeline(c1=m, c2=t, c3=u, c5=p, c4=u, prediction_type=:interval) isa
          IntervalPipeline)
    @test(Pipeline(c1=m, c2=t, c3=u, c5=p, c4=u,
                   prediction_type=:probabilistic) isa
          ProbabilisticPipeline)
    @test_logs((:info, MLJBase.INFO_TREATING_AS_DETERMINISTIC),
               @test Pipeline(c1=m, c2=t, c3=u, c8=i, c4=u) isa
               DeterministicPipeline)
    # if "hidden" supervised model is already deterministic,
    # no need for warning:
    @test_logs(@test Pipeline(c1=m, c2=t, c3=u, c7=d, c4=u) isa
               DeterministicPipeline)

    # errors and warnings:
    @test_throws MLJBase.ERR_MIXED_PIPELINE_SPEC Pipeline(m, mymodel=p)
    @test_throws(MLJBase.ERR_INVALID_OPERATION,
                 Pipeline(u, s, operation=cos))
    @test_throws(MLJBase.ERR_INVALID_PREDICTION_TYPE,
                 Pipeline(u=u, s=s, prediction_type=:ostrich))
    @test_logs((:warn, MLJBase.WARN_IGNORING_PREDICTION_TYPE),
               Pipeline(m, t, u, prediction_type=:deterministic))
    @test_throws(MLJBase.err_prediction_type_conflict(d, :probabilistic),
                 Pipeline(m, t, d, prediction_type=:probabilistic))

end

@testset "property access" begin
    pipe = Pipeline(m, u, u, s)

    # property names:
    @test propertynames(pipe) ===
        (:f, :my_unsupervised, :my_unsupervised2, :my_transformer2, :cache)

    # getproperty:
    @test pipe.my_unsupervised == u
    @test pipe.my_unsupervised2 == u
    @test pipe.cache == true

    # replacing a component with one whose abstract supertype is the same
    # or smaller:
    pipe.my_unsupervised = s
    @test pipe.my_unsupervised == s

    # attempting to replace a component with one whose abstract supertype
    # is bigger:
    @test_throws MethodError pipe.my_transformer2 = u

    # mutating the components themeselves:
    pipe.my_unsupervised.ftr = :z
    @test pipe.my_unsupervised.ftr == :z

    # or using MLJBase's recursive getproperty:
    MLJBase.recursive_setproperty!(pipe, :(my_unsupervised.ftr), :bonzai)
    @test pipe.my_unsupervised.ftr ==  :bonzai

    # more errors:
    @test_throws(MLJBase.err_pipeline_bad_property(pipe, :mount_fuji),
                 pipe.mount_fuji)
    @test_throws(MLJBase.err_pipeline_bad_property(pipe, :mount_fuji),
                 pipe.mount_fuji = 42)

end

@testset "show" begin
    io = IOBuffer()
    pipe = Pipeline(x-> x^2, m, t, p)
    show(io, MIME("text/plain"), pipe)
end

@testset "Front and extend" begin
    Xs = source(rand(3))
    ys = source(rand(3))
    front = MLJBase.Front(Xs, ys, false)
    @test front.predict() == Xs()
    @test front.transform() == ys()
    @test MLJBase.active(front)() == Xs()

    front = MLJBase.Front(Xs, Xs, true)
    front = MLJBase.extend(front, u, false, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:up, :up, :up]
    @test tnode() == [:ut, :ut, :ut]

    front = MLJBase.extend(front, s, true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:up, :up, :up]
    @test tnode() == [:st, :st, :st]

    front = MLJBase.extend(front, u, true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:up, :up, :up]
    @test tnode() == [:ut, :ut, :ut]

    front = MLJBase.extend(front, d, true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:dp, :dp, :dp]
    @test tnode() == [:dt, :dt, :dt]

    front = MLJBase.extend(front, x->string.(x), true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == string.([:dp, :dp, :dp])
    @test tnode() == [:dt, :dt, :dt]

    front = MLJBase.extend(front, u, true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:ut, :ut, :ut]
    @test tnode() == [:dt, :dt, :dt]

    front = MLJBase.extend(front, s, true, predict, ys)
    pnode, tnode = front.predict, front.transform
    fit!(pnode, verbosity=0)
    fit!(tnode, verbosity=0)
    @test pnode() == [:st, :st, :st]
    @test tnode() == [:dt, :dt, :dt]
end

NN = 7
X = MLJBase.table(rand(rng, NN, 3));
y = 2X.x1 - X.x2 + 0.05*rand(rng,NN);
Xs = source(X); ys = source(y)

broadcast_mode(v) = mode.(v)
doubler(y) = 2*y

@testset "pipeline_network_machine" begin
    t = MLJBase.table
    m = MLJBase.matrix
    f = FeatureSelector()
    h = OneHotEncoder()
    k = KNNRegressor()
    u = UnivariateStandardizer()
    c = ConstantClassifier()

    components = [f, k]
    mach = MLJBase.pipeline_network_machine(
        Deterministic, true, predict, components, Xs, ys)
    tree = mach.fitresult.predict |> MLJBase.tree
    @test mach.model isa DeterministicSurrogate
    @test tree.operation == predict
    @test tree.model == k
    @test tree.arg1.operation == transform
    @test tree.arg1.model == f
    @test tree.arg1.arg1.source == Xs
    @test tree.arg1.train_arg1.source == Xs
    @test tree.train_arg1 == tree.arg1
    @test tree.train_arg2.source == ys

    components = [f, h]
    mach = MLJBase.pipeline_network_machine(
        Unsupervised, true, predict, components, Xs)
    tree = mach.fitresult.transform |> MLJBase.tree
    @test mach.model isa UnsupervisedSurrogate
    @test tree.operation == transform
    @test tree.model == h
    @test tree.arg1.operation == transform
    @test tree.arg1.model == f
    @test tree.arg1.arg1.source == Xs
    @test tree.arg1.train_arg1.source == Xs
    @test tree.train_arg1 == tree.arg1

    components = [m, t]
    mach = MLJBase.pipeline_network_machine(
        Unsupervised, true, predict, components, Xs)
    tree = mach.fitresult.transform |> MLJBase.tree
    @test mach.model isa UnsupervisedSurrogate
    @test tree.operation == t
    @test tree.model == nothing
    @test tree.arg1.operation == m
    @test tree.arg1.model == nothing
    @test tree.arg1.arg1.source == Xs

    # check a probablistic case:
    components = [f, c]
    mach = MLJBase.pipeline_network_machine(
        Probabilistic, true, predict, components, Xs, ys)
    @test mach.model isa ProbabilisticSurrogate

    # check a static case:
    components = [m, t]
    mach = MLJBase.pipeline_network_machine(
        Static, true, predict, components, Xs, ys)
    @test mach.model isa StaticSurrogate

    # An integration test...

    # build a linear network for training:
    components = [f, k]
    mach = MLJBase.pipeline_network_machine(
        Deterministic, true, predict, components, Xs, ys)

    # build the same network by hand:
    fM = machine(f, Xs)
    Xt = transform(fM, Xs)
    uM = machine(u, ys)
    yt = transform(uM, ys)
    kM = machine(k, Xt, yt)
    zhat = predict(kM, Xt)
    N2 = inverse_transform(uM, zhat)

    # compare predictions
    fit!(mach, verbosity=0);
    fit!(N2, verbosity=0)
    yhat = predict(mach, X);
    @test yhat ≈ N2()
    k.K = 3; f.features = [:x3,]
    fit!(mach, verbosity=0);
    fit!(N2, verbosity=0)
    @test !(yhat ≈ predict(mach, X))
    @test predict(mach, X) ≈ N2()
    global hand_built = predict(mach, X);

end

@testset "training_losses" begin
    model = MyDeterministic(:bla)
    pipe = Standardizer() |> model

    # test helpers:
    @test MLJBase.supervised_component_name(pipe) == :my_deterministic
    @test MLJBase.supervised_component(pipe) == model

    @test supports_training_losses(pipe)
    _, _, rp = MLJBase.fit(pipe, 0, X, y)
    @test training_losses(pipe, rp) == ones(3)
    @test iteration_parameter(pipe) == :(my_deterministic.x)
end

# # INTEGRATION TESTS

@testset "integration 1" begin
    # check a simple pipeline prediction agrees with prediction of
    # hand-built learning network built earlier:
    p = Pipeline(FeatureSelector,
                 KNNRegressor,
                 prediction_type=:deterministic)
    p.knn_regressor.K = 3; p.feature_selector.features = [:x3,]
    mach = machine(p, X, y)
    fit!(mach, verbosity=0)
    @test MLJBase.tree(mach.fitresult.predict).model.K == 3
    MLJBase.tree(mach.fitresult.predict).arg1.model.features == [:x3, ]
    @test predict(mach, X) ≈ hand_built

    # Check target_scitype of a supervised pipeline is the same as the supervised component
    @test target_scitype(p) == target_scitype(KNNRegressor())

    # test cache is set correctly internally:
    @test all(fitted_params(mach).machines) do m
        MLJBase._cache_status(m)  == " caches data"
    end

    # test correct error thrown for inverse_transform:
    @test_throws(MLJBase.err_unsupported_operation(:inverse_transform),
                 inverse_transform(mach, 3))
end

@testset "integration 2" begin
    # a simple probabilistic classifier pipeline:
    X = MLJBase.table(rand(rng,7,3));
    y = categorical(collect("ffmmfmf"));
    Xs = source(X)
    ys = source(y)
    p = Pipeline(OneHotEncoder, ConstantClassifier, cache=false)
    mach = machine(p, X, y)
    fit!(mach, verbosity=0)
    @test p isa ProbabilisticComposite
    pdf(predict(mach, X)[1], 'f') ≈ 4/7

    # test cache is set correctly internally:
    @test all(fitted_params(mach).machines) do m
        MLJBase._cache_status(m) == " does not cache data"
    end

    # test invalid replacement of classifier with regressor throws
    # informative error message:
    p.constant_classifier = ConstantRegressor()
    @test_logs((:error, r"^Problem"),
           (:info, r"^Running type"),
           (:warn, r"The scitype of"),
           (:info, r"It seems"),
           (:error, r"Problem"),
               @test_throws Exception fit!(mach, verbosity=-1))
end

@testset "integration 3" begin
    # test a simple deterministic classifier pipeline:
    X = MLJBase.table(rand(rng,7,3))
    y = categorical(collect("ffmmfmf"))
    Xs = source(X)
    ys = source(y)
    p = Pipeline(OneHotEncoder, ConstantClassifier, broadcast_mode,
                 prediction_type=:probabilistic)
    mach = machine(p, X, y)
    fit!(mach, verbosity=0)
    @test predict(mach, X) == fill('f', 7)

    # Check target_scitype of a supervised pipeline is the same as the supervised component
    @test target_scitype(p) == target_scitype(ConstantClassifier())

    # test pipelines with weights:
    w = map(y) do η
        η == 'm' ? 100 : 1
    end
    mach = machine(p, X, y, w)
    fit!(mach, verbosity=0)
    @test predict(mach, X) == fill('m', 7)
end

age = [23, 45, 34, 25, 67]
X = (age = age,
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
height = [67.0, 81.5, 55.6, 90.0, 61.1]

mutable struct MyTransformer3 <: Static
    ftr::Symbol
end

MLJBase.transform(transf::MyTransformer3, verbosity, X) =
     selectcols(X, transf.ftr)

@testset "integration 4" begin
    #static transformers in pipelines
    p99 = Pipeline(X -> coerce(X, :age=>Continuous),
                   OneHotEncoder,
                   MyTransformer3(:age))
    mach = fit!(machine(p99, X), verbosity=0)
    @test transform(mach, X) == float.(X.age)

    # Check target_scitype of an unsupervised pipeline is Unknown
    @test target_scitype(p99) == Unknown
end

@testset "integration 5" begin
    # pure static pipeline:
    p = Pipeline(X -> coerce(X, :age=>Continuous),
                 MyTransformer3(:age))

    mach = fit!(machine(p), verbosity=0) # no training arguments!
    @test transform(mach, X) == X.age

    # and another:
    p = Pipeline(exp, log, x-> 2*x)
    mach = fit!(machine(p), verbosity=0)
    @test transform(mach, 20) ≈ 40
end

@testset "integration 6" begin
    # operation different from predict:
    p = Pipeline(OneHotEncoder,
                 ConstantRegressor,
                 operation=predict_mean)
    @test p isa Deterministic
    mach = fit!(machine(p, X, height), verbosity=0)
    @test scitype(predict(mach, X)) == AbstractVector{Continuous}
end

@testset "integration 7" begin
    # inverse transform:
    p = Pipeline(UnivariateBoxCoxTransformer,
                 UnivariateStandardizer)
    xtrain = rand(rng, 10)
    mach = machine(p, xtrain)
    fit!(mach, verbosity=0)
    x = rand(rng, 5)
    y = transform(mach, x)
    x̂ = inverse_transform(mach, y)
    @test isapprox(x, x̂)
    # Check target_scitype of an unsupervised pipeline is Unknown
    @test target_scitype(p) == Unknown
end

# A dummy clustering model:
mutable struct DummyClusterer <: Unsupervised
    n::Int
end
DummyClusterer(; n=3) = DummyClusterer(n)
function MLJBase.fit(model::DummyClusterer, verbosity::Int, X)
    Xmatrix = Tables.matrix(X)
    n = min(size(Xmatrix, 2), model.n)
    centres = Xmatrix[1:n, :]
    levels = categorical(1:n)
    report = (centres=centres,)
    fitresult = levels
    return fitresult, nothing, report
end
MLJBase.transform(model::DummyClusterer, fitresult, Xnew) =
    selectcols(Xnew, 1:length(fitresult))
MLJBase.predict(model::DummyClusterer, fitresult, Xnew) =
    [fill(fitresult[1], nrows(Xnew))...]

@testset "integration 8" begin
    # calling predict on unsupervised pipeline
    # https://github.com/JuliaAI/MLJClusteringInterface.jl/issues/10

    N = 20
    X = (a = rand(N), b = rand(N))

    p = Pipeline(PCA, DummyClusterer)
    mach = machine(p, X)
    fit!(mach, verbosity=0)
    y = predict(mach, X)
    @test y == fill(categorical(1:2)[1], N)
end

@testset "syntactic sugar" begin

    # recall u, s, p, m, are defined way above

    # unsupervised model |> static model:
    pipe1 = u |> s
    @test pipe1 == Pipeline(u, s)

    # unsupervised model |> supervised model:
    pipe2 = u |> p
    @test pipe2 == Pipeline(u, p)

    # pipe |> pipe:
    hose = pipe1 |> pipe2
    @test hose == Pipeline(u, s, u, p)

    # pipe |> model:
    @test Pipeline(u, s) |> p == Pipeline(u, s, p)

    # model |> pipe:
    @test u |> Pipeline(s, p) == Pipeline(u, s, p)

    # pipe |> function:
    @test Pipeline(u, s) |> m == Pipeline(u, s, m)

    # function |> pipe:
    @test m |> Pipeline(s, p) == Pipeline(m, s, p)

    # model |> function:
    @test u |> m == Pipeline(u, m)

    # function |> model:
    @test t |> u == Pipeline(t, u)

    @test_logs((:info, MLJBase.INFO_AMBIGUOUS_CACHE),
               Pipeline(u, cache=false) |> p)

    # with types
    @test PCA |> Standardizer() |> KNNRegressor ==
        Pipeline(PCA(), Standardizer(), KNNRegressor())
end

@testset "miscelleneous coverage" begin
    @test MLJBase.as_type(:unsupervised) == Unsupervised
end

end

true


text:
using WebIO
using JSExpr # you may need to install this package

function counter(start=0)
    scope = Scope()

    # updates to this update the UI
    count = Observable(scope, "count", start)

    onjs(count, # listen on JavaScript
         JSExpr.@js x->this.dom.querySelector("#count").textContent = x)

    on(count) do n # listen on Julia
        println(n > 0 ? "+"^n : "-"^abs(n))
    end

    btn(label, d) = dom"button"(
        label,
        events=Dict(
            "click" => JSExpr.@js () -> $count[] = $count[] + $d
        )
    )

    scope.dom = dom"div"(
        btn("increment", 1),
        btn("decrement", -1),
        dom"div#count"(string(count[])),
    )

    scope
end

# Display in whatever frontend is avalaible
function main()
    if @isdefined(IJulia) || @isdefined(Juno)
        return counter(1)
    elseif @isdefined(Blink)
        win = Window()
        body!(win, counter(1))
    elseif @isdefined(Mux)
        @sync webio_serve(page("/", req -> counter(1)), rand(8000:9000))
    else
        error("do one of using Mux, using Blink before running the
               example, or run it from within IJulia or Juno")
    end
end

main()


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Build a system image binary at sysimg_path.dlext.  By default, put the system image
# next to libjulia (except on Windows, where it goes in $JULIA_HOME\..\lib\julia)
# Allow insertion of a userimg via userimg_path.  If sysimg_path.dlext is currently loaded into memory,
# don't continue unless force is set to true.  Allow targeting of a CPU architecture via cpu_target
@unix_only function default_sysimg_path(debug=false)
    joinpath(dirname(Libdl.dlpath(debug ? "libjulia-debug" : "libjulia")),
             "julia", debug ? "sys-debug" : "sys")
end

@windows_only function default_sysimg_path(debug=false)
    joinpath(JULIA_HOME, "..", "lib", "julia", debug ? "sys-debug" : "sys")
end

function build_sysimg(sysimg_path=nothing, cpu_target="native", userimg_path=nothing; force=false, debug=false)
    if sysimg_path == nothing
        sysimg_path = default_sysimg_path(debug)
    end

    # Quit out if a sysimg is already loaded and is in the same spot as sysimg_path, unless forcing
    sysimg = Libdl.dlopen_e("sys")
    if sysimg != C_NULL
        if !force && Base.samefile(Libdl.dlpath(sysimg), "$(sysimg_path).$(Libdl.dlext)")
            info("System image already loaded at $(Libdl.dlpath(sysimg)), set force to override")
            return
        end
    end

    # Canonicalize userimg_path before we enter the base_dir
    if userimg_path != nothing
        userimg_path = abspath(userimg_path)
    end

    # Enter base/ and setup some useful paths
    base_dir = dirname(Base.find_source_file("sysimg.jl"))
    cd(base_dir) do
        julia = joinpath(JULIA_HOME, debug ? "julia-debug" : "julia")
        cc = find_system_compiler()

        # Ensure we have write-permissions to wherever we're trying to write to
        try
            touch("$sysimg_path.ji")
        catch
            err_msg =  "Unable to modify $sysimg_path.ji, ensure parent directory exists "
            err_msg *= "and is writable. Absolute paths work best.)"
            error( err_msg )
        end

        # Copy in userimg.jl if it exists...
        if userimg_path != nothing
            if !isfile(userimg_path)
                error("$userimg_path is not found, ensure it is an absolute path!")
            end
            if isfile("userimg.jl")
                error("$base_dir/userimg.jl already exists, delete manually to continue.")
            end
            cp(userimg_path, "userimg.jl")
        end
        try
            # Start by building inference0.{ji,o}
            inference0_path = joinpath(dirname(sysimg_path), "inference0")
            info("Building inference0.o...")
            println("$julia -C $cpu_target --output-ji $inference0_path.ji --output-o $inference0_path.o coreimg.jl")
            run(`$julia -C $cpu_target --output-ji $inference0_path.ji --output-o $inference0_path.o coreimg.jl`)

            # Bootstrap off off that to create inference.{ji,o}
            inference_path = joinpath(dirname(sysimg_path), "inference")
            info("Building inference.o...")
            println("$julia -C $cpu_target --output-ji $inference_path.ji --output-o $inference_path.o coreimg.jl")
            run(`$julia -C $cpu_target --output-ji $inference_path.ji --output-o $inference_path.o coreimg.jl`)

            # Bootstrap off off that to create sys.{ji,o}
            info("Building sys.o...")
            println("$julia -C $cpu_target --output-ji $sysimg_path.ji --output-o $sysimg_path.o -J $inference_path.ji --startup-file=no sysimg.jl")
            run(`$julia -C $cpu_target --output-ji $sysimg_path.ji --output-o $sysimg_path.o -J $inference_path.ji --startup-file=no sysimg.jl`)

            if cc != nothing
                link_sysimg(sysimg_path, cc, debug)
            else
                info("System image successfully built at $sysimg_path.ji")
            end

            if !Base.samefile("$default_sysimg_path.ji", "$sysimg_path.ji")
                if Base.isfile("$sysimg_path.$(Libdl.dlext)")
                    info("To run Julia with this image loaded, run: julia -J $sysimg_path.$(Libdl.dlext)")
                else
                    info("To run Julia with this image loaded, run: julia -J $sysimg_path.ji")
                end
            else
                info("Julia will automatically load this system image at next startup")
            end
        finally
            # Cleanup userimg.jl
            if userimg_path != nothing && isfile("userimg.jl")
                rm("userimg.jl")
            end
        end
    end
end

# Search for a compiler to link sys.o into sys.dl_ext.  Honor LD environment variable.
function find_system_compiler()
    if haskey( ENV, "CC" )
        if !success(`$(ENV["CC"]) -v`)
            warn("Using compiler override $(ENV["CC"]), but unable to run `$(ENV["CC"]) -v`")
        end
        return ENV["CC"]
    end

    # On Windows, check to see if WinRPM is installed, and if so, see if gcc is installed
    @windows_only try
        eval(Main, :(using WinRPM))
        winrpmgcc = joinpath(WinRPM.installdir,"usr","$(Sys.ARCH)-w64-mingw32",
            "sys-root","mingw","bin","gcc.exe")
        if success(`$winrpmgcc --version`)
            return winrpmgcc
        else
            throw()
        end
    catch
        warn("Install GCC via `Pkg.add(\"WinRPM\"); WinRPM.install(\"gcc\")` to generate sys.dll for faster startup times")
    end

    # See if `cc` exists
    try
        if success(`cc -v`)
            return "cc"
        end
    end

    warn( "No supported compiler found; startup times will be longer" )
end

# Link sys.o into sys.$(dlext)
function link_sysimg(sysimg_path=nothing, cc=find_system_compiler(), debug=false)
    if sysimg_path == nothing
        sysimg_path = default_sysimg_path(debug)
    end
    julia_libdir = dirname(Libdl.dlpath(debug ? "libjulia-debug" : "libjulia"))

    FLAGS = ["-L$julia_libdir"]

    push!(FLAGS, "-shared")
    push!(FLAGS, debug ? "-ljulia-debug" : "-ljulia")
    @windows_only push!(FLAGS, "-lssp")

    info("Linking sys.$(Libdl.dlext)")
    run(`$cc $FLAGS -o $sysimg_path.$(Libdl.dlext) $sysimg_path.o`)

    info("System image successfully built at $sysimg_path.$(Libdl.dlext)")
    @windows_only begin
        if convert(VersionNumber, Base.libllvm_version) < v"3.5.0"
            LLVM_msg = "Building sys.dll on Windows against LLVM < 3.5.0 can cause incorrect backtraces!"
            LLVM_msg *= " Delete generated sys.dll to avoid these problems"
            warn( LLVM_msg )
        end
    end
end

# When running this file as a script, try to do so with default values.  If arguments are passed
# in, use them as the arguments to build_sysimg above
if !isinteractive()
    if length(ARGS) > 5 || ("--help" in ARGS || "-h" in ARGS)
        println("Usage: build_sysimg.jl <sysimg_path> <cpu_target> <usrimg_path.jl> [--force] [--debug] [--help]")
        println("   <sysimg_path>    is an absolute, extensionless path to store the system image at")
        println("   <cpu_target>     is an LLVM cpu target to build the system image against")
        println("   <usrimg_path.jl> is the path to a user image to be baked into the system image")
        println("   --debug          Using julia-debug instead of julia to build the system image")
        println("   --force          Set if you wish to overwrite the default system image")
        println("   --help           Print out this help text and exit")
        println()
        println(" Example:")
        println("   build_sysimg.jl /usr/local/lib/julia/sys core2 ~/my_usrimg.jl --force")
        println()
        println(" Running this script with no arguments is equivalent to:")
        println("   build_sysimg.jl $(default_sysimg_path) native")
        return 0
    end

    debug_flag = "--debug" in ARGS
    filter!(x -> x != "--debug", ARGS)
    force_flag = "--force" in ARGS
    filter!(x -> x != "--force", ARGS)
    build_sysimg(ARGS...; force=force_flag, debug=debug_flag)
end


text:
function constructnetwork!(m::JuMP.AbstractModel, branch_models::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, netinjection::BalanceNamedTuple, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: CopperPlatePowerModel}

    copperplatebalance(m, netinjection, sys.time_periods)

end

function constructnetwork!(m::JuMP.AbstractModel, branch_models::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, netinjection::BalanceNamedTuple, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: AbstractFlowForm}

    for category in branch_models
        constructdevice!(m, netinjection, category.device, category.formulation, system_formulation, sys; args...)
    end

    nodalflowbalance(m, netinjection, system_formulation, sys)

end

function constructnetwork!(m::JuMP.AbstractModel, branch_models::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, netinjection::BalanceNamedTuple, system_formulation::Type{StandardPTDF}, sys::PowerSystems.PowerSystem; args...)
    if :PTDF in keys(args)
        PTDF = args[:PTDF]
    else
        PTDF = nothing
    end

    if !isa(PTDF,PTDFArray)
        @warn "no PTDF supplied"
        PTDF,  A = PowerSystems.buildptdf(sys.branches, sys.buses)
    end

    for category in branch_models
        constructdevice!(m, netinjection, category.device, category.formulation, system_formulation, sys; args..., PTDF=PTDF)
    end

    nodalflowbalance(m, netinjection, system_formulation, sys)

end

function constructnetwork!(m::JuMP.AbstractModel, branch_models::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, netinjection::BalanceNamedTuple, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: AbstractDCPowerModel}

    #= TODO: Needs to be generalized later for other branch models not covered by PM.
    for category in branch_models
        constructdevice!(m, netinjection, category.device, category.formulation, system_formulation, sys; args...)
    end
    =#

    nodalflowbalance(m, netinjection, system_formulation, sys)

    PM_F = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, system_formulation; kwargs...)

    PM_object = PS.build_nip_expr_model(m.ext[:PM_object], PM_F, jump_model=m);

    m = PM_object.model

end

function constructnetwork!(m::JuMP.AbstractModel, branch_models::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}, netinjection::BalanceNamedTuple, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: AbstractACPowerModel}

    #= TODO: Needs to be generalized later for other branch models not covered by PM.
    for category in branch_models
        constructdevice!(m, netinjection, category.device, category.formulation, system_formulation, sys; args...)
    end
    =#

    nodalflowbalance(m, netinjection, system_formulation, sys)

    PM_F = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, system_formulation; kwargs...)

    PM_object = PS.build_nip_expr_model(m.ext[:PM_object], PM_F, jump_model=m);

    m.ext[:PM_object] = PM_object

end


text:
# Autogenerated wrapper script for Cairomm_jll for i686-w64-mingw32-cxx03
export libcariomm

using libsigcpp_jll
using Cairo_jll
JLLWrappers.@generate_wrapper_header("Cairomm")
JLLWrappers.@declare_library_product(libcariomm, "libcairomm-1.16-1.dll")
function __init__()
    JLLWrappers.@generate_init_header(libsigcpp_jll, Cairo_jll)
    JLLWrappers.@init_library_product(
        libcariomm,
        "bin\\libcairomm-1.16-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
        capacity(P::Drip;      # Drip structure
                 unit = "bit"  # optional: unit of capacity (bit or nat).
                 )
Returns the amount of information processes per unit of time in the steady state
of the DRIP `P`.
"""
function capacity(P::Drip;      # Drip structure
                  unit = "bit"  # optional: unit of capacity (bit or nat).
                  )
    if unit == "bit"
        κ = 0.5*log(det(P.ss.Σ_1)/det(P.ss.Σ_p))/log(2); #returns capacity in bits
    elseif unit == "nat"
        κ = 0.5*log(det(P.ss.Σ_1)/det(P.ss.Σ_p));        #returns capacity in nats
    else
        println("Invalid input for unit! Capacity is reported in bits.")
        κ = 0.5*log(det(P.ss.Σ_1)/det(P.ss.Σ_p))/log(2);
    end
    return(κ)
end

"""
    infinitesum(func; tol = 1e-6,maxit = 1000,start=0)
Returns the infinite sum `Σₓfunc(x)` starting from `x = start` up to tolderance
`tol` or max iteration `maxit`.
"""
function infinitesum(func; tol = 1e-6,maxit = 1000,start=0)
    diff  = 1.0
    infsum = func(start)
    it    = start + 1
    while (diff > tol) & (it < maxit)
        func_it = func(it)
        infsum += func_it
        diff = maximum(func_it)
        it += 1
    end
    return(infsum)
end


text:
# loading basic packages
using DelimitedFiles, Printf, LinearAlgebra, Random

# main function to convert a file with three dimensions coordinates in a pdb type file
"""
```
pdbCreate(inputFile::AbstractString, outputFile::AbstractString = "")

```
This is the main function in the program. To run this function, you need a file with a matrix (n, 3) of n real coordinates [x, y, z] of each atom.
    in inputFile parameter you must to inform the name of the file, in pwd() dir, that contains this matrix and, in outputFile, um can inform the name of generated pdb file.

## Example

```julia-repl
julia> pdbCreate("\\examples\\10.txt", "10.pdb")
```
if successful, a pdb file is provided in current directory, named as default "<name_of_inputFile>.pdb".
"""
function pdbCreate(inputFileName::AbstractString, outputFileName::AbstractString = "")
    inputFileName = string(pwd(),inputFileName)
    if (isdir(inputFileName))
        println("Invalid File Path")
        usage()
        return
    else
        realizations = readdlm(inputFileName)
        open(outputFileName == "" ? basename(splitext(inputFileName)[1]*".pdb") : outputFileName, "w") do f
            (m,n) = size(realizations)
            write(f, "HEADER    LAVOR-INSTANCES                         28-MAR-20   0000              \n")
            write(f, @sprintf "TITLE     THE THREE-DIMENSIONAL LAVOR INSTANCE OF LENGHT %9d             \n" m)
            write(f, "MODEL        1                                                                  \n")
            for i in 1:m
                    s = @sprintf "this is a %s %15.1f" "test" 34.567;
                    write(f, @sprintf "ATOM  %5d  H   MET A   1     %7.3f %7.3f %7.3f  1.00  0.00           H  \n" i realizations[i,1] realizations[i,2] realizations[i,3])
            end
            write(f, "ENDMDL                                                                          \n")
            write(f, @sprintf "MASTER        0    0    0    0    0    0    0 %5d    0    0    0              \n" m)
            write(f, "END                                                                             \n")
        end
        size(realizations)
    end
end

# function to create a file with a random instance
"""
```
randomInstanceCreate(n::Integer, outputFileName::AbstractString = ""; isLimited::Bool = true)

```
This function generate a proteic instance with `n` atoms randomly from a physical model that is close to reality, based on the paper

    LAVOR, C. . On generating instances for the molecular distance geometry
    problem. In: Leo Liberti; Nelson Maculan. (Org.). Global Optimization: from
    Theory to Implementation. Global Optimization: from Theory to
    Implementation. New York: Springer, 2006, v. 84, p. 405-414.

and that instance is saved in a delimited file named, optionally, just like the outputFileName parameter. By default, the file has the name "`n`.xyz", where `n` is the size of instance.

> Note: Since PDB files have a limited number of digits in their coordinates, it is of interest to limit the space of the molecule. For this, the `isLimited` parameter allows an internal loop to generate a limited protein.

## Example

```julia-repl
julia> randomInstanceCreate(10)
```
"""
function randomInstanceCreate(n::Integer, outputFileName::AbstractString = ""; isLimited::Bool = true)
    positions = zeros(n, 3)
    ω=[60;180;300]
    ωₑ=(-15:15)'

    θ=109.5*pi/180;
    r=1.526;

    B1=Array(1.0*I(4));
    B2=Array(1.0*I(4));
    B3=Array(1.0*I(4));

    B2[1,1]=-1;
    B2[3,3]=-1;
    B2[1,4]=-r;

    B3[1,1]=-cos(θ);
    B3[1,2]=-sin(θ);
    B3[1,4]=-r*cos(θ);
    B3[2,1]=sin(θ);
    B3[2,2]=-cos(θ);
    B3[2,4]=r*sin(θ);

    Orig=[0;0;0;1];

    B=B1*B2;
    v=B*Orig;
    positions[2,:]=(v[1:3])';
    B=B*B3;
    v=B*Orig;
    positions[3,:]=(v[1:3])';
    Omega=zeros(n-3,1);

    for i=4:n
        while true
            n1aux=randperm(3);
            n1=n1aux[1];
            n2aux=randperm(31);
            n2=n2aux[7];

            omega=(ω[n1]+ωₑ[n2])*pi/180;
            Omega[i-3]=omega;
            Bi=zeros(4,4);

            Bi[1,1]=-cos(θ);
            Bi[1,2]=-sin(θ);
            Bi[1,4]=-r*cos(θ);
            Bi[2,1]=sin(θ)*cos(omega);
            Bi[2,2]=-cos(θ)*cos(omega);
            Bi[2,3]=-sin(omega);
            Bi[2,4]=r*sin(θ)*cos(omega);
            Bi[3,1]=sin(θ)*sin(omega);
            Bi[3,2]=-cos(θ)*sin(omega);
            Bi[3,3]=cos(omega);
            Bi[3,4]=r*sin(θ)*sin(omega);
            Bi[4,4]=1;

            Btemp=B*Bi;
            v=Btemp*Orig;

            if (isLimited && (abs(v[1]) >= 99 || abs(v[2]) >= 99 || abs(v[3]) >= 99))
                continue
            end
            B = Btemp

            positions[i,:]=(v[1:3])';
            break
        end
    end

    open((outputFileName == "") ? string(n)*".xyz" : outputFileName, "w") do io
        writedlm(io, positions)
    end
end

# function for define the usage of pdbCreate main function
"""
```
usage()

```
This is a very simple function to inform the correct utilization of pdbCreate function

## Example

```julia-repl
julia> pdbCreate("")
usage: pdbCreate("<your_file_of_coordinates>", "<name_of_pdb_output_file>")
```
"""
function usage()
    print("usage: pdbCreate(\"<your_file_of_coordinates>\", \"<name_of_pdb_output_file>\")\n")
end

max_stars_repo_name: caomem/PdbFileCreator

max_issues_repo_name: caomem/PdbFileCreator

max_forks_repo_name: caomem/PdbFileCreator


text:
using Test, BitBasis

function swaprows!(v::AbstractVector, i::Int, j::Int)
    temp = v[i]
    v[i] = v[j]
    v[j] = temp
    v
end

function eye(b::Int)
    res = zeros(b, b)
    for i in 1:b
        res[i, i] = 1
    end
    res
end

mulrow!(v::AbstractVector, i::Int, f) = (v[i] *= f; v)

@testset "test group_shift! and lmove" begin
    @test group_shift!(5, [1, 2, 5]) == ([0, 15], [2, 1])
    @test group_shift!(5, [2, 3]) == ([1], [2])
    @test group_shift!(5, [1, 3, 5]) == ([0, 3, 15], [1, 1, 1])

    @test string(lmove(5, 1, 2), base = 2) == "10001"
end

@testset "test IterControl" begin
    v = randn(ComplexF64, 1 << 4)
    it = itercontrol(4, [3], [1])
    vec = Int[]
    it2 = itercontrol(4, [3, 4], [0, 0])
    for i in it2
        push!(vec, i)
    end
    @test vec == [0, 1, 2, 3]

    vec = Int[]
    it4 = itercontrol(4, [4, 2, 1], [1, 1, 1])
    for i in it4
        push!(vec, i)
    end
    @test vec == [11, 15]

    # test controlled operations
    @testset "test controlled operations" begin
        P0 = ComplexF64[1 0; 0 0]
        P1 = ComplexF64[0 0; 0 1]
        Z = ComplexF64[1 0; 0 -1]

        rrr = copy(v)
        controldo(x -> mulrow!(rrr, x + 1, -1.0), it4)
        M =
            kron(P1, eye(2), P1, Z) +
            kron(P0, eye(2), P0, eye(2)) +
            kron(P1, eye(2), P0, eye(2)) +
            kron(P0, eye(2), P1, eye(2))
        @test rrr ≈ M * v

        it = itercontrol(8, [3], [1])
        V = randn(ComplexF64, 1 << 8)
        res = kron(eye(1 << 5), ComplexF64[0 1; 1 0], eye(1 << 2)) * V
        rrr = copy(V)
        controldo(x -> swaprows!(rrr, x + 1, x - 3), it)
        @test rrr ≈ res

        rrr = copy(V)
        controldo(x -> mulrow!(rrr, x + 1, -1), itercontrol(8, [3, 7, 6], [1, 1, 1]))

        M =
            kron(eye(2), P1, Z, eye(4), P1, eye(4)) +
            kron(eye(2), P0, eye(8), P0, eye(4)) +
            kron(eye(2), P1, eye(8), P0, eye(4)) +
            kron(eye(2), P0, eye(8), P1, eye(4))

        @test rrr ≈ M * V
        ic = itercontrol(4, [2, 3], [0, 0])
        @test [ic...] == [ic[k] for k in 1:4]
    end
end


text:
using DataFrames

using GLM

using UnPack

using Feather
using RegressionDiscontinuity

using Statistics
using StatsBase
using StatsModels

using Test

lee08_path = joinpath(dirname(@__FILE__), "..", "data", "lee08.feather")
lee08 = Feather.read(lee08_path)

lee08_rdd = load_rdd_data(:lee08)

Z = lee08.margin ./ 100
Y = lee08.voteshare ./ 100

rdd_df = DataFrame(Z = Z, Y = Y)

N = length(Z)
left_idx = findall(Z .< 0)
right_idx = findall(Z .>= 0)
N_left = length(left_idx)
N_right = length(right_idx)

# Step 1: Density and conditional variance at 0
h₁ = 1.84 * std(Z) * N^(-1 / 5)

h₁_left_idx = findall(-h₁ .<= Z .< 0)
h₁_right_idx = findall(+h₁ .>= Z .>= 0)

N_h₁_left = length(h₁_left_idx)
N_h₁_right = length(h₁_right_idx)

Ȳ_h₁_left = mean(Y[h₁_left_idx])
sd_Y_h₁_left = std(Y[h₁_left_idx])
Ȳ_h₁_right = mean(Y[h₁_right_idx])
sd_Y_h₁_right = std(Y[h₁_right_idx])

f̂₀ = (N_h₁_left + N_h₁_right) / 2 / N / h₁

# Step 2: Estimation of second derivatives
global_cubic_lm = lm(@formula(Y ~ 1 + (Z >= 0) + Z + Z^2 + Z^3), rdd_df)

m̂₀_triple_prime = 6 * coef(global_cubic_lm)[5]

h₂_left = 3.56 * (sd_Y_h₁_left^2 / f̂₀ / m̂₀_triple_prime^2)^(1 / 7) * N_left^(-1 / 7)
h₂_right = 3.56 * (sd_Y_h₁_right^2 / f̂₀ / m̂₀_triple_prime^2)^(1 / 7) * N_right^(-1 / 7)

h₂_left_idx = findall(-h₂_left .<= Z .< 0)
h₂_right_idx = findall(+h₂_right .>= Z .>= 0)

N_h₂_left = length(h₂_left_idx)
N_h₂_right = length(h₂_right_idx)

quadratic_fit_left = lm(@formula(Y ~ 1 + Z + Z^2), rdd_df[h₂_left_idx, :])
m̂₀_double_prime_left = 2 * last(coef(quadratic_fit_left))

quadratic_fit_right = lm(@formula(Y ~ 1 + Z + Z^2), rdd_df[h₂_right_idx, :])
m̂₀_double_prime_right = 2 * last(coef(quadratic_fit_right))

# Step 3: Calculation of regularization terms

r̂_left = 2160 * sd_Y_h₁_left^2 / N_h₂_left / h₂_left^4
r̂_right = 2160 * sd_Y_h₁_right^2 / N_h₂_right / h₂_right^4

Ckernel = 3.4375

regularized_squared_diff =
    abs2(m̂₀_double_prime_right - m̂₀_double_prime_left) + r̂_left + r̂_right
ĥ_IK =
    Ckernel *
    ((sd_Y_h₁_left^2 + sd_Y_h₁_right^2) / f̂₀ / regularized_squared_diff)^(1 / 5) *
    N^(-1 / 5)

ik_triang_paper = 0.2939
point_est_triang_paper = 0.0799
se_triang_paper = 0.0083

bw_ik_triang = bandwidth(ImbensKalyanaraman(), SymTriangularDist(), lee08_rdd)
fit_triang = fit(
    NaiveLocalLinearRD(kernel = SymTriangularDist(), bandwidth = ImbensKalyanaraman()),
    lee08_rdd,
)

@test bw_ik_triang ≈ ik_triang_paper atol = 0.001
@test fit_triang.fitted_bandwidth == bw_ik_triang
@test fit_triang.tau_est ≈ point_est_triang_paper atol = 0.001
@test fit_triang.se_est ≈ se_triang_paper atol = 0.001

@test bw_ik_triang ≈ ĥ_IK atol = 0.001

ik_rect_paper = 0.4617
point_est_rect_paper = 0.0806
se_rect_paper = 0.0087

bw_ik_rect = bandwidth(ImbensKalyanaraman(), Rectangular(), lee08_rdd)
fit_rect = fit(
    NaiveLocalLinearRD(kernel = Rectangular(), bandwidth = ImbensKalyanaraman()),
    lee08_rdd,
)

@test bw_ik_rect ≈ ik_rect_paper atol = 0.001
@test fit_rect.fitted_bandwidth == bw_ik_rect
@test fit_rect.tau_est ≈ point_est_rect_paper atol = 0.001
@test fit_rect.se_est ≈ se_rect_paper atol = 0.001

Ys = lee08_rdd.Ys
ZsR = lee08_rdd.ZsR

bw_kernel = SymTriangularDist(0, bw_ik_triang)

support(bw_kernel)

#ZsR[support(bw_kernel)]

propertynames(lee08_rdd)
propertynames(lee08_rdd)

RegressionDiscontinuity.kernel_constant(ImbensKalyanaraman(), Rectangular())
RegressionDiscontinuity.kernel_constant(ImbensKalyanaraman(), SymTriangularDist())


text:
"""
    BasicHeuristic(selectValue::Function)

A BasicHeuristic is the type of ValueSelection that is not learning. It makes it possible to use
SeaPearl.jl as a classic CP Solver. This is useful in itself as it is a fully julia native CP Solver.
This is also useful to be able to compare the performances of the LearnedHeuristic to some handcrafted heuristics.

To create one, the user just as to give it a function which map an `AbstractIntVar` to a value of its domain.
"""
mutable struct BasicHeuristic <: ValueSelection
    selectValue::Function
    search_metrics::Union{Nothing, SearchMetrics}

end

"""
    BasicHeuristic()

Create the default `BasicHeuristic` that selects the maximum value of the domain
"""
BasicHeuristic() = BasicHeuristic((x; cpmodel=nothing) -> maximum(x.domain), nothing)
BasicHeuristic(selectValue::Function) = BasicHeuristic(selectValue, nothing)
"""
    (valueSelection::BasicHeuristic)(::LearningPhase, model, x, current_status)

Explains what the basicHeurstic should do at each step of the solving. This is useful to have a unified `search!` function working with both
BasicHeuristic and LearnedHeuristic. In the case of the BasicHeuristic, it is only called in the DecisionPhase where the selectValue function is used
to choose the value assigned.
"""
(valueSelection::BasicHeuristic)(::Type{InitializingPhase}, model::Union{Nothing, CPModel}=nothing) = (valueSelectionsearch_metrics = SearchMetrics(model))
(valueSelection::BasicHeuristic)(::Type{StepPhase}, model::Union{Nothing, CPModel}=nothing, current_status::Union{Nothing, Symbol}=nothing) = nothing
function (valueSelection::BasicHeuristic)(::Type{DecisionPhase}, model::Union{Nothing, CPModel}=nothing, x::Union{Nothing, AbstractIntVar}=nothing)
    model.statistics.lastVar = x
    return valueSelection.selectValue(x; cpmodel=model)
end
(valueSelection::BasicHeuristic)(::Type{EndingPhase}, model::Union{Nothing, CPModel}=nothing, current_status::Union{Nothing, Symbol}=nothing) = nothing

wears_mask(valueSelection::BasicHeuristic) = true