text:
"""
    allocate_input(p::PencilFFTPlan)          -> PencilArray
    allocate_input(p::PencilFFTPlan, dims...) -> Array{PencilArray}
    allocate_input(p::PencilFFTPlan, Val(N))  -> NTuple{N, PencilArray}

Allocate uninitialised
[`PencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.PencilArray)
that can hold input data for the given plan.

The second and third forms respectively allocate an array of `PencilArray`s of
size `dims`, and a tuple of `N` `PencilArray`s.

!!! note "In-place plans"

    If `p` is an in-place plan, a
    [`ManyPencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.ManyPencilArray)
    is allocated. This
    type holds `PencilArray` wrappers for the input and output transforms (as
    well as for intermediate transforms) which share the same space in memory.
    The input and output `PencilArray`s should be respectively accessed by
    calling [`first(::ManyPencilArray)`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#Base.first-Tuple{ManyPencilArray}) and
    [`last(::ManyPencilArray)`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#Base.last-Tuple{ManyPencilArray}).

    #### Example

    Suppose `p` is an in-place `PencilFFTPlan`. Then,
    ```julia
    @assert is_inplace(p)
    A = allocate_input(p) :: ManyPencilArray
    v_in = first(A)       :: PencilArray  # input data view
    v_out = last(A)       :: PencilArray  # output data view
    ```

    Also note that in-place plans must be performed directly on the returned
    `ManyPencilArray`, and not on the contained `PencilArray` views:
    ```julia
    p * A       # perform forward transform in-place
    p \\ A       # perform backward transform in-place
    # p * v_in  # not allowed!!
    ```
"""
function allocate_input end

# Out-of-place version
function allocate_input(p::PencilFFTPlan{T,N,false} where {T,N})
    T = eltype_input(p)
    pen = pencil_input(p)
    PencilArray{T}(undef, pen, p.extra_dims...)
end

# In-place version
function allocate_input(p::PencilFFTPlan{T,N,true} where {T,N})
    pencils = map(pp -> pp.pencil_in, p.plans)

    # Note that for each 1D plan, the input and output pencils are the same.
    # This is because the datatype stays the same for in-place transforms
    # (in-place real-to-complex transforms are not supported!).
    @assert pencils === map(pp -> pp.pencil_out, p.plans)

    T = eltype_input(p)
    ManyPencilArray{T}(undef, pencils...; extra_dims=p.extra_dims)
end

allocate_input(p::PencilFFTPlan, dims...) =
    _allocate_many(allocate_input, p, dims...)

"""
    allocate_output(p::PencilFFTPlan)          -> PencilArray
    allocate_output(p::PencilFFTPlan, dims...) -> Array{PencilArray}
    allocate_output(p::PencilFFTPlan, Val(N))  -> NTuple{N, PencilArray}

Allocate uninitialised [`PencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.PencilArray) that can hold output data for the
given plan.

If `p` is an in-place plan, a [`ManyPencilArray`](https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.ManyPencilArray) is allocated.

See [`allocate_input`](@ref) for details.
"""
function allocate_output end

# Out-of-place version.
function allocate_output(p::PencilFFTPlan{T,N,false} where {T,N})
    T = eltype_output(p)
    pen = pencil_output(p)
    PencilArray{T}(undef, pen, p.extra_dims...)
end

# For in-place plans, the output and input are the same ManyPencilArray.
allocate_output(p::PencilFFTPlan{T,N,true} where {T,N}) = allocate_input(p)

allocate_output(p::PencilFFTPlan, dims...) =
    _allocate_many(allocate_output, p, dims...)

_allocate_many(allocator::Function, p::PencilFFTPlan, dims::Vararg{Int}) =
    [allocator(p) for I in CartesianIndices(dims)]

_allocate_many(allocator::Function, p::PencilFFTPlan, ::Val{N}) where {N} =
    ntuple(n -> allocator(p), Val(N))


text:
include("SimpleNumField/Field.jl")
include("SimpleNumField/Subfields.jl")

ext: jl


text:
using CodecBGZF
import TranscodingStreams:
    TranscodingStream,
    test_roundtrip_read,
    test_roundtrip_write,
    test_roundtrip_lines,
    test_roundtrip_transcode
using Base.Test

@testset "BGZF Codec" begin
    # Test the EOF trailer (28 bytes).
    data = UInt8[]
    buffer = IOBuffer(data, true, true)
    stream = BGZFCompressorStream(buffer)
    write(stream, "foobar")
    data = buffer.data
    close(stream)
    @test sizeof(data) ≥ 28
    @test data[end-27:end] == [
        0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00,
        0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x02, 0x00,
        0x1b, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00,
        0x00, 0x00, 0x00, 0x00]
    @test read(BGZFDecompressorStream(IOBuffer(data))) == b"foobar"

    # test.bgz is created using the bgzip command
    bgzfile = joinpath(dirname(@__FILE__), "test.bgz")
    stream = BGZFDecompressorStream(open(bgzfile))
    @test read(stream) == b"abracadabra"
    close(stream)

    data = read(bgzfile)
    stream = BGZFDecompressorStream(IOBuffer(vcat(data, data, data)))
    seek(stream, UInt64(sizeof(data) << 16))
    @test read(stream, 5) == b"abrac"
    seek(stream, UInt64(sizeof(data) << 16) | UInt64(2))
    @test read(stream, 5) == b"racad"
    # Seeking a stream to an unaligned offset will fail.
    @test_throws Exception seek(stream, UInt64((sizeof(data) + 1) << 16) | UInt64(2))
    close(stream)

    test_roundtrip_read(BGZFCompressorStream, BGZFDecompressorStream)
    test_roundtrip_write(BGZFCompressorStream, BGZFDecompressorStream)
    test_roundtrip_lines(BGZFCompressorStream, BGZFDecompressorStream)
    test_roundtrip_transcode(BGZFCompressor, BGZFDecompressor)
end


text:
# EBI Proteins
# ======
#
# APIs for EBI Proteins (UniProt)
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
EBI Proteins module.

The `EBIProteins` module provides a programming interface to the Uniprot
databases at EMBL-EBI. Ten functions are exported from this module:
1. `proteins`: protein sequences
2. `features`: functional annotations
3. `variation`: variation data
4. `proteomics`: proteomics datasets
5. `antigen`: antigens
6. `proteomes`: proteome data
7. `genecentric`: proteome data grouped by under a gene
8. `taxonomy`: taxonomic relationships
9. `coordinates`: genomic coordinates
10. `uniparc` : non-redundant protein sequences

See "EMBL-EBI Proteins API"
(https://www.ebi.ac.uk/proteins/api/doc/) for more details. The implemented
APIs are based on the document as of March 24, 2021.

Note: size = -1 allows all results to be streamed; this is not currently
implemented.
"""
module EBIProteins

export
    ebiproteins,
    ebifeatures,
    ebivariation,
    ebiproteomics,
    ebiantigen,
    ebiproteomes,
    ebigenecentric,
    ebitaxonomy,
    ebicoordinates,
    ebiuniparc

import XMLDict
import JSON
import HTTP
import BioServices.request

const baseURL = "https://www.ebi.ac.uk/proteins/api/"

# APIs of UniProt Proteins
# -------------------

"""
    ebiproteins(params...)

Retrieve protein data.

Calls the endpoint specified by "<operation>/<accession>/<subset>". If any argument is
unset it will be excluded.

If `dbtype` and `dbid` are set, calls the endpoint "<dbtype>:<dbid>".

If any query parameters are provided, calls the base endpoint.

Key Parameters: contenttype, operation, accession, subset, dbtype, dbid

`contenttype` may be "text/x-fasta", "text/x-flatfile", "application/json", "application/xml"

Query Parameters: offset, size, reviewed, isoform, goterms, keywords, ec, gene, exact_gene,
protein, organism, taxid, pubmed, seqLength, md5

"""
function ebiproteins(; contenttype = nothing, accession = nothing, operation = nothing, subset = nothing,
                  dbtype = nothing, dbid = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "proteins"
    if !isnothing(dbtype) || !isnothing(dbid)
        all(!isnothing, [dbid, dbtype]) || throw(ArgumentError("Both `dbtype` and `dbid` must be set."))
        path *= "/" * HTTP.escapeuri(dbtype) * ":" * HTTP.escapeuri(dbid)
    elseif isempty(params)
        any(!isnothing, [operation, accession, subset]) || throw(ArgumentError("At least one parameter must be set."))
        if !isnothing(operation)
            path *= "/" * HTTP.escapeuri(operation)
        end
        if !isnothing(accession)
            path *= "/" * HTTP.escapeuri(accession)
        end
        if !isnothing(subset)
            path *= "/" * HTTP.escapeuri(subset)
        end
    elseif !isnothing(accession)
        params[:accession] = string(accession)
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebifeatures(params...)

Retrieve protein data.

Calls the endpoint at "<accession>".

If `type` is set, calls the endpoint "type/<type>".

If any query parameters other than `categories` or `types` are provided,
or `accession` is not set, calls the base endpoint.

Key Parameters: contenttype, type, accession

`contenttype` may be "text/x-gff", "application/json", "application/xml"

Query Parameters: offset, size, reviewed, gene, exact_gene, protein,
organism, taxid, categories, types, terms
"""
function ebifeatures(; contenttype = nothing, accession = nothing, type = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "features"
    !isnothing(type) && !isnothing(accession) && throw(ArgumentError("Cannot set both `type` and `accession`"))
    if !isnothing(type)
        haskey(params, :terms) && params[:terms] != "" || throw(ArgumentError("At least one term is required."))
        path *= "/type/" * HTTP.escapeuri(type)
    elseif !isnothing(accession) && isempty(filter(kv -> first(kv) ∉ [:categories, :types], params))
        path *= "/" * HTTP.escapeuri(accession)
    elseif isnothing(accession) && isempty(params)
        throw(ArgumentError("At least one parameter must be set."))
    elseif !isnothing(accession)
        params[:accession] = string(accession)
    else
        any(haskey.(Ref(params), [:gene, :protein, :exact_gene, :organism, :taxid])) ||
            throw(ArgumentError("Must set at least one of `gene`, `exact_gene`, `protein`, `organism`, or `taxid`."))
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebivariation(params...)

Retrieve variant data.

Calls the endpoint at "<accession>".

If `dbid` is set, calls the endpoint "dbsnp/<dbid>".

If `hgvs` is set, calls the endpoint "hgvs/<hgvs>".

If any query parameters other than `sourcetype`, `consequencetype`,
`wildtype`, `alternativesequence`, or `location` are provided,
or `accession` is not set, calls the base endpoint.

Key Parameters: contenttype, accession, dbid, hgvs

`contenttype` may be "text/x-gff", "text/x-peff", "application/json", "application/xml"

Query Parameters: offset, size, sourcetype, consequencetype, wildtype,
alternativesequence, location, disease, omim, evidence, taxid, dbtype, dbid
"""
function ebivariation(; contenttype = nothing, accession = nothing, dbid = nothing, hgvs = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "variation"
    !isnothing(dbid) && !isnothing(hgvs) && throw(ArgumentError("Cannot set both `dbid` and `hgvs`"))
    if !isnothing(dbid)
        path *= "/dbsnp/" * HTTP.escapeuri(dbid)
    elseif !isnothing(hgvs)
        path *= "/hgvs/" * HTTP.escapeuri(hgvs)
    elseif !isnothing(accession) && isempty(filter(kv -> first(kv) ∉ [:sourcetype, :consequencetype, :wildtype, :alternativesequence, :location], params))
        path *= "/" * accession
    elseif isnothing(accession) && isempty(params)
        throw(ArgumentError("At least one parameter must be set."))
    elseif !isnothing(accession)
        params[:accession] = string(accession)
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebiproteomics(params...)

Retrieve proteomics data.

Calls the endpoint at "<accession>".

If any query parameters are set or `accession` is not set, calls the base endpoint.

Key Parameters: contenttype, accession

`contenttype` may be "text/x-gff", "application/json", "application/xml"

Query Parameters: offset, size, taxid, upid, datasource, peptide, unique
"""
function ebiproteomics(; contenttype = nothing, accession = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "proteomics"
    if isempty(params)
        isnothing(accession) && throw(ArgumentError("At least one term is required."))
        path *= "/" * HTTP.escapeuri(accession)
    elseif !isnothing(accession)
        params[:accession] = string(accession)
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebiantigen(params...)

Retrieve proteomics data.

Calls the endpoint at "<accession>"

If any query parameters are set or `accession` is not set, calls the base endpoint.

Key Parameters: contenttype, accession

`contenttype` may be "text/x-gff", "application/json", "application/xml"

Query Parameters: offset, size, antigen_sequence, antigen_id, ensembl_ids, match_score
"""
function ebiantigen(; contenttype = nothing, accession = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "antigen"
    if isempty(params)
        isnothing(accession) && throw(ArgumentError("At least one term is required."))
        path *= "/" * HTTP.escapeuri(accession)
    elseif !isnothing(accession)
        params[:accession] = string(accession)
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebiproteomes(params...)

Retrieve proteome data.

Calls the endpoint at "<operation>/<upid>". If any argument is unset it will be excluded.

If any other query parameters are set or `upid` is not set, calls the base endpoint.

Key Parameters: contenttype, upid, reviewed

`contenttype` may be "application/json", "application/xml"

Query Parameters: offset, size, upid, name, taxid, keyword, xref, genome_acc,
is_ref_proteome, is_redundant, reviewed
"""
function ebiproteomes(; contenttype = nothing, operation = nothing, upid = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "proteomes"
    if isempty(filter(kv -> first(kv) != :reviewed, params))
        !isnothing(upid) || throw(ArgumentError("Must set `upid`"))
        if !isnothing(operation)
            path *= "/" * HTTP.escapeuri(operation)
        end
        path *= "/" * HTTP.escapeuri(upid)
    elseif !isnothing(upid)
        params[:upid] = HTTP.escapeuri(string(upid))
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebigenecentric(params...)

Retrieve proteomics data.

Calls the endpoint at "<accession>"

If any query parameters are set or `accession` is not set, calls the base endpoint.

Key Parameters: contenttype, accession

`contenttype` may be "application/json", "application/xml"

Query Parameters: offset, size, upid, gene
"""
function ebigenecentric(; contenttype = nothing, accession = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "genecentric"
    if isempty(params)
        isnothing(accession) && throw(ArgumentError("At least one term is required."))
        path *= "/" * HTTP.escapeuri(accession)
    elseif !isnothing(accession)
        params[:accession] = HTTP.escapeuri(string(accession))
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebitaxonomy(params...)

Retrieve taxonomy data.

Calls the endpoint at "<operation>/<ids>/<subset>". If `getnode` is true, calls the endpoint
at "<operation>/<ids>/<subset>/node". If any argument is unset it will be excluded.

If `operation` is "path" or "relationship", calls those respective endpoints: "path", "path/nodes",
"relationship".

Key Parameters: contenttype, operation, ids, subset

`contenttype` may be "application/json", "application/xml"

Query Parameters: pageNumber, pageSize, searchType, fieldName, depth, direction, from, to
"""
function ebitaxonomy(; contenttype = nothing, operation = nothing, ids = nothing, subset = nothing,
                     getnode = false, params...)
    params = process_parameters(params)
    !isnothing(operation) || throw(ArgumentError("Must set `operation`."))
    path = baseURL * "taxonomy/" * HTTP.escapeuri(operation)
    if operation == "path"
        haskey(params, :direction) || throw(ArgumentError("Must set `direction`."))
        !isnothing(ids) || throw(ArgumentError("Must set `ids`."))
        params[:id] = HTTP.escapeuri(string(ids))
        if getnode
            path *= "/nodes"
        end
    elseif operation == "relationship"
        haskey(params, :from) && haskey(params, :to) || throw(ArgumentError("Must set `from` and `to`."))
    else
        !isnothing(ids) || throw(ArgumentError("Must set `ids`."))
        path *= "/" * HTTP.escapeuri(string(ids))
        if !isnothing(subset)
            path *= "/" * HTTP.escapeuri(subset)
        end
        if getnode
            path *= "/node"
        end
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebicoordinates(params...)

Retrieve coordinates data.

If `accession` is set, calls the endpoint at "<accession>". If a position is also provided,
calls the endpoint at either "location/<accession>:<pPosition>" or "location/<accession>:<pStart>-<pEnd>".

If `dbtype` and `dbid` are provided, calls the endpoint at "<dbtype>:<dbid>".

If `taxonomy` and `locations` are provided, calls the endpoint at "<taxonomy>/<locations>". If `getfeature`
is true, access the endpoint at "<taxonomy>/<locations>/feature".

If any query parameters not approptiate to the specific endpoint are set, calls the base endpoint.

Key Parameters: contenttype, accession, dbtype, dbid, pPosition, pEnd, taxonomy, locations, getfeature

`contenttype` may be "application/json", "application/xml"

Query Parameters: offset, size, in_range
"""
function ebicoordinates(; accession = nothing, dbtype = nothing, dbid = nothing, pPosition = nothing,
                        pStart = nothing, pEnd = nothing, locations = nothing, taxonomy = nothing,
                        getfeature = false, contenttype = nothing, params...)
    params = process_parameters(params)
    path = baseURL * "coordinates"
    if !isnothing(accession) && isempty(params)
        isnothing(dbtype) && isnothing(dbid) && isnothing(taxonomy) && isnothing(locations) ||
            throw(ArgumentError("Cannot set `dbtype`, `dbid`, `taxonomy`, or `locations` with `accession`."))
        path *= "/"
        if !isnothing(pPosition) || !isnothing(pStart)
            isnothing(pStart) && isnothing(pEnd) || (pPosition = pStart)
            path *= "location/" * HTTP.escapeuri(accession) * ":" * HTTP.escapeuri(pPosition)
            if !isnothing(pEnd)
                path *= "-" * HTTP.escapeuri(pEnd)
            end
        else
            path *= HTTP.escapeuri(accession)
        end
    elseif !isnothing(dbtype) || !isnothing(dbid)
        isnothing(accession) && isnothing(taxonomy) && isnothing(locations) ||
            throw(ArgumentError("Cannot set `accession`, `taxonomy`, or `locations` with `dbtype` and `dbid`."))
        all(!isnothing, [dbid, dbtype]) || throw(ArgumentError("Both `dbtype` and `dbid` must be set."))
        path *= "/" * HTTP.escapeuri(dbtype) * ":" * HTTP.escapeuri(dbid)
    elseif !isnothing(taxonomy) || !isnothing(locations)
        isnothing(accession) ||
            throw(ArgumentError("Cannot set `accession` with `taxonomy` and `locations`."))
        all(!isnothing, [taxonomy, locations]) || throw(ArgumentError("Both `taxonomy` and `locations` must be set."))
        path *= "/" * HTTP.escapeuri(taxonomy) * "/" * HTTP.escapeuri(locations)
        getfeature && (path *= "/feature")
    elseif !isnothing(accession)
        params[:accession] = HTTP.escapeuri(string(accession))
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

"""
    ebiuniparc(params...)

Retrieve proteomics data.

Calls the endpoint at "accession/<accession>", "dbreference/<dbid>", "uniparc/upi/<upi>",
or "sequence", depending on provided values.

If `bestguess` is true, calls the endpoint at "bestguess".

If any excess query parameters are set for the specific endpoint, calls the base endpoint.

Key Parameters: contenttype, accession, upid, upi, dbid, sequence, sequencecontenttype,
bestguess

`contenttype` may be "text/x-fasta", "application/json", "application/xml"

`sequencecontenttype` may be "text/plain" (default), "application/json", "application/xml", if set

Query Parameters: offset, size, dbtype, gene, protein, taxid, organism, sequencechecksum,
ipr, signaturetype, signatureid, seqLength, rdDdtype, rfDbid, rfActive, rfTaxId
"""
function ebiuniparc(; contenttype = nothing, accession = nothing, dbid = nothing, upid = nothing,
                    upi = nothing, sequence = nothing, sequencecontenttype = "text/plain", bestguess = false,
                    params...)
    params = process_parameters(params)
    path = baseURL * "uniparc"
    if !isnothing(sequence)
        path *= "/sequence"
        header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype, "Content-Type" => sequencecontenttype]
        body = HTTP.escapeuri(sequence)
        return request("POST", path, body = body, header, query = params, retry_non_idempotent = true)
    elseif !isnothing(accession) && isempty(filter(kv -> first(kv) ∉ [:rfDdtype, :rfDbid, :rfActive, :rfTaxId], params))
        path *= "/accession/" * HTTP.escapeuri(accession)
    elseif !isnothing(upid) && isempty(filter(kv -> first(kv) ∉ [:offset, :size, :rfDdtype, :rfDbid, :rfActive, :rfTaxId], params))
        path *= "/proteome/" * HTTP.escapeuri(upid)
    elseif !isnothing(dbid) && isempty(filter(kv -> first(kv) ∉ [:offset, :size, :rfDdtype, :rfDbid, :rfActive, :rfTaxId], params))
        path *= "/dbreference/" * HTTP.escapeuri(dbid)
    elseif !isnothing(upi) && isempty(filter(kv -> first(kv) ∉ [:rfDdtype, :rfDbid, :rfActive, :rfTaxId], params))
        path *= "/upi/" * HTTP.escapeuri(upi)
    elseif bestguess
        path *= "/bestguess"
    else
        isnothing(accession) || (params[:accession] = HTTP.escapeuri(accession))
        isnothing(upid) || (params[:upid] = HTTP.escapeuri(upid))
        isnothing(dbid) || (params[:dbid] = HTTP.escapeuri(dbid))
        isnothing(upi) || (params[:upi] = HTTP.escapeuri(upi))
    end
    header = isnothing(contenttype) ? Pair{String, String}[] : ["Accept" => contenttype]
    return request("GET", path, header, query = params)
end

# Process query parameters.
function process_parameters(params)
    # stringify all values
    params = Dict{Symbol, Any}(params)
    for (key, val) in params
        params[key] = string(val)
    end

    return params
end

if VERSION < v"1.1"
isnothing(x) = x === nothing
end

end

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
"""
    normalize!(pcloud::PointCloud)

Normalize the PointCloud `pcloud` with mean centered at origin and unit standard deviation
and overwrite the pcloud with normalized PointCloud.

See also: [`normalize`](@ref)

### Examples:

```julia
julia> p = PointCloud(rand(1024,3))
julia> normalize!(p)
```
"""
function normalize!(pcloud::PointCloud)
    centroid = mean(pcloud.points, dims = 2)
    pcloud.points =
        (pcloud.points .- centroid) ./
        (std(pcloud.points, mean = centroid, dims = 2) .+ EPS)
    return pcloud
end

"""
    normalize(pcloud::PointCloud)

Normalize the PointCloud `pcloud` with mean centered at origin and unit standard deviation

See also: [`normalize!`](@ref)

### Examples:

```julia
julia> p = PointCloud(rand(1024,3))
julia> p = normalize(p)
```
"""
function normalize(pcloud::PointCloud)
    p = deepcopy(pcloud)
    normalize!(p)
    return p
end

"""
    scale!(pcloud::PointCloud, factor::Number)

Scale the PointCloud `pcloud` by scaling factor `factor`
and overwrite `pcloud` with scaled PointCloud.

Scaling factor `factor` should be strictly greater than `0.0`.

See also: [`scale`](@ref)

### Examples:
```julia
julia> p = PointCloud(rand(1024,3))
julia> scale!(p, 1.0)
```
"""
function scale!(pcloud::PointCloud, factor::Float32)
    (factor > 0.0) || error("factor must be greater than 0.0")
    lmul!(factor, pcloud.points)
    return pcloud
end

scale!(pcloud::PointCloud, factor::Number) = scale!(pcloud, Float32(factor))

"""
    scale(pcloud::PointCloud, factor::Number)

Scale the PointCloud `pcloud` by scaling factor `factor`.

Scaling factor `factor` should be strictly greater than `0.0`.

See also: [`scale!`](@ref)

### Examples:
```julia
julia> p = PointCloud(rand(1024,3))
julia> p = scale(1.0)
```
"""
function scale(pcloud::PointCloud, factor::Float32)
    p = deepcopy(pcloud)
    scale!(p, factor)
    return p
end

scale(pcloud::PointCloud, factor::Number) = scale(pcloud, Float32(factor))

"""
    rotate!(pcloud::PointCloud, rotmat::AbstractArray{Number,2})
    rotate!(pcloud::PointCloud, rotmat::AbstractArray{Number,3})

Rotate the PointCloud `pcloud` by rotation matrix `rotmat`
and overwrite `pcloud` with rotated PointCloud.

Rotation matrix `rotmat` should be of size `(3,3)` or `(3,3,B)`
where B is the batch size of PointCloud.

See also: [`rotate`](@ref)

### Examples:
```julia
julia> p = PointCloud(rand(1024,3))
julia> rotmat = rand(3,3)
julia> rotate!(p, rotmat)
```
"""
function rotate!(pcloud::PointCloud, rotmat::AbstractArray{Float32,2})
    size(rotmat) == (3, 3) ||
        error("rotmat must be (3, 3) array, but instead got $(size(rotmat)) array")
    size(pcloud.points, 1) == 3 || error("dimension of points in PointCloud must be 3")
    points_packed = transpose(rotmat) * reshape(pcloud.points, 3, :)
    pcloud.points = reshape(points_packed, size(pcloud.points)...)
    return pcloud
end

function rotate!(pcloud::PointCloud, rotmat::AbstractArray{Float32,3})
    _B = size(pcloud.points, 3)
    size(rotmat) == (3, 3, _B) ||
        error("rotmat must be (3, 3, $(_B)) array, but instead got $(size(rotmat)) array")
    size(pcloud.points, 1) == 3 || error("dimension of points in PointCloud must be 3")
    pcloud.points = Flux.batched_mul(Flux.batched_transpose(rotmat), pcloud.points)
    return pcloud
end

rotate!(pcloud::PointCloud, rotmat::AbstractArray{<:Number}) =
    rotate!(pcloud, Float32.(rotmat))

"""
    rotate(pcloud::PointCloud, rotmat::Array{<:Number})

Rotate the PointCloud `pcloud` by rotation matrix `rotmat`.

Rotation matrix `rotmat` should be of size `(3,3)`

See also: [`rotate!`](@ref)

### Examples:
```julia
julia> p = PointCloud(rand(1024,3))
julia> rotmat = rand(3,3)
julia> p = rotate(p, rotmat)
```
"""
function rotate(pcloud::PointCloud, rotmat::AbstractArray{<:Number})
    p = deepcopy(pcloud)
    rotate!(p, rotmat)
    return p
end

"""
    realign!(src::PointCloud, tgt::PointCloud)
    realign!(src::PointCloud, tgt_min::AbstractArray{<:Number,2}, tgt_max::AbstractArray{<:Number,2})

Re-Align the PointCloud `src` with the axis aligned bounding box of PointCloud `tgt`
and overwrite `pcloud` with re-aligned PointCloud.

PointCloud `src` and `tgt` should be of same dimension.

See also: [`realign`](@ref)

### Examples:
```julia
julia> src = PointCloud(rand(1024,3))
julia> tgt = PointCloud(rand(1024,3))
julia> realign!(src, tgt)
```
"""
function realign!(
    src::PointCloud,
    tgt_min::AbstractArray{Float32,2},
    tgt_max::AbstractArray{Float32,2},
)
    size(src.points, 1) == size(tgt_max, 1) ||
        error("source and target pointcloud dimension mismatch")
    src_min = minimum(src.points, dims = 2)
    src_max = maximum(src.points, dims = 2)
    src.points =
        ((src.points .- src_min) ./ (src_max - src_min .+ EPS)) .* (tgt_max - tgt_min) .+
        tgt_min
    return src
end

realign!(
    src::PointCloud,
    tgt_min::AbstractArray{<:Number,2},
    tgt_max::AbstractArray{<:Number,2},
) = realign!(src, Float32.(tgt_min), Float32.(tgt_max))

function realign!(src::PointCloud, tgt::PointCloud, index::Number = 1)
    points = tgt[index]
    tgt_min = minimum(points, dims = 2)
    tgt_max = maximum(points, dims = 2)
    realign!(src, tgt_min, tgt_max)
    return src
end

"""
    realign(src::PointCloud, tgt::PointCloud)
    realign(src::PointCloud, tgt_min::AbstractArray{<:Number,2}, tgt_max::AbstractArray{<:Number,2})

Re-Align the PointCloud `src` with the axis aligned bounding box of PointCloud `tgt`.

PointCloud `src` and `tgt` should be of same dimension.

See also: [`realign!`](@ref)

### Examples:
```julia
julia> src = PointCloud(rand(1024,3))
julia> tgt = PointCloud(rand(1024,3))
julia> src = realign!(src, tgt)
```
"""
function realign(
    src::PointCloud,
    tgt_min::AbstractArray{<:Number,2},
    tgt_max::AbstractArray{<:Number,2},
)
    p = deepcopy(src)
    realign!(p, tgt_min, tgt_max)
    return p
end

function realign(src::PointCloud, tgt::PointCloud, index::Number = 1)
    p = deepcopy(src)
    realign!(p, tgt, index)
    return p
end


text:
struct FlattenRNNHidden <: SlothLayer; end
(l::FlattenRNNHidden)(x) = ndims(x) == 1 ? x : reshape(x, size(x,1), :)

function reset!(r::RNN; h=0.0f0, c=0.f0)
    r.h, r.c = h, c
    nothing
end

function RNN(; input_size::Int, hidden_size::Int, kwargs...)
    RNN(input_size, hidden_size; kwargs...)
end

ReluRNN(; kwargs...) = RNN(; rnnType=:relu, kwargs...)
TanhRNN(; kwargs...) = RNN(; rnnType=:tanh, kwargs...)
LSTM(; kwargs...) = RNN(; kwargs...)
GRU(; kwargs...) = RNN(; rnnType=:gru, kwargs...)
BiLSTM(; kwargs...) = RNN(; bidirectional=true, kwargs...)


text:
# Number of vector pointers in Product struct
function Base.length(p::Product)::Int
    length(p.xss)
end

# Number of combinations in a cartesian inner product
function Base.size(p::Product)::Int
    l = 1
    for i = 1:length(p)
        l *= length(p.xss[i])
    end
    return l
end

# Length * Width
area(v)::Int = prod(size(v))

# Maps a function f over a vector, returning NaN if vector is empty
function safemap(f::Function, v)
    function f_safe(x)
        try
            f(x)
        catch
            NaN
        end
    end
    map(f_safe, v)
end

# Converts vector into a vector of tuples
function to_tuple(v)
    map(x -> (x,), v)
end

# Compares arrays against a tuple
function opt_compare(varr, p::T where T <: Tuple)
    if length(varr) !== length(p)
        throw("length of vectors does not equal length of keys")
    end

    bv = BitArray{1}(length(varr[1]))
    for i=1:length(p)
        if i===1
            bv .= varr[i] .== p[i]
        else
            bv .= bv .& varr[i] .== p[i]
        end
    end

    return bv
end

# Returns indices where arrays match key combinations
function find_matches(varr, k)
    buffer = Array{Array{Int, 1}, 1}(area(k))
    bvecs = Array{BitArray{1}, 1}(length(k))
    bvec = BitArray{1}(length(varr[1]))
    i = 0
    for p in k
        i += 1
        #bvecs .= [varr[j] .== ifelse(length(p)==1, p, p[j]) for j=1:length(p)]
        #bvec .= sum(bvecs) .== length(p)
        bvec .= opt_compare(varr, p)
        buffer[i] = find(bvec)
    end
    return(buffer)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SshConfiguration <: SwaggerModel
    publicKeys::Any # spec type: Union{ Nothing, Vector{SshPublicKey} } # spec name: publicKeys

    function SshConfiguration(;publicKeys=nothing)
        o = new()
        validate_property(SshConfiguration, Symbol("publicKeys"), publicKeys)
        setfield!(o, Symbol("publicKeys"), publicKeys)
        o
    end
end # type SshConfiguration

const _property_map_SshConfiguration = Dict{Symbol,Symbol}(Symbol("publicKeys")=>Symbol("publicKeys"))
const _property_types_SshConfiguration = Dict{Symbol,String}(Symbol("publicKeys")=>"Vector{SshPublicKey}")
Base.propertynames(::Type{ SshConfiguration }) = collect(keys(_property_map_SshConfiguration))
Swagger.property_type(::Type{ SshConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SshConfiguration[name]))}
Swagger.field_name(::Type{ SshConfiguration }, property_name::Symbol) =  _property_map_SshConfiguration[property_name]

function check_required(o::SshConfiguration)
    true
end

function validate_property(::Type{ SshConfiguration }, name::Symbol, val)
end


text:
"""
    ncon(tensorlist, indexlist, [conjlist, sym]; order = ..., output = ...)

Contract the tensors in `tensorlist` (of type `Vector` or `Tuple`) according to the network
as specified by `indexlist`. Here, `indexlist` is a list (i.e. a `Vector` or `Tuple`) with
the same length as `tensorlist` whose entries are themselves lists (preferably
`Vector{Int}`) where every integer entry provides a label for corresponding index/dimension
of the corresponding tensor in `tensorlist`. Positive integers are used to label indices
that need to be contracted, and such thus appear in two different entries within
`indexlist`, whereas negative integers are used to label indices of the output tensor, and
should appear only once.

Optional arguments in another list with the same length, `conjlist`, whose entries are of
type `Bool` and indicate whether the corresponding tensor object should be conjugated
(`true`) or not (`false`). Finally, a `Symbol` can be provided which provides a hook into
the global LRU cache of TensorOperations.jl to store temporaries. This symbol should be
unique for every `ncon` call withing a running application.

By default, contractions are performed in the order such that the indices being contracted
over are labelled by increasing integers, i.e. first the contraction corresponding to label
`1` is performed. The output tensor had an index order corresponding to decreasing
(negative, so increasing in absolute value) index labels. The keyword arguments `order` and
`output` allow to change these defaults.

See also the macro version [`@ncon`](@ref).
"""
function ncon(tensors, network,
                conjlist = fill(false, length(tensors)), sym = nothing ;
                order = nothing, output = nothing)
    length(tensors) >= 2 ||
        throw(ArgumentError("do not use `ncon` for less than two tensors"))
    length(tensors) == length(network) == length(conjlist) ||
        throw(ArgumentError("number of tensors and of index lists should be the same"))
    isnconstyle(network) || throw(ArgumentError("invalid NCON network: $network"))
    outputindices = Vector{Int}()
    for n in network
        for k in n
            if k < 0
                push!(outputindices, k)
            end
        end
    end
    if output === nothing
        output = sort(outputindices; rev = true)
    else
        for a in output
            a in outputindices ||
                throw(ArgumentError("invalid NCON network: $network -> $output"))
        end
        for a in outputindices
            a in output ||
                throw(ArgumentError("invalid NCON network: $network -> $output"))
        end
    end
    tree = order === nothing ? ncontree(network) : indexordertree(network, order)

    if sym !== nothing
        syma = Symbol(sym, "_a")
        symb = Symbol(sym, "_b")
    else
        syma = symb = nothing
    end
    A, IA, CA = contracttree(tensors, network, conjlist, tree[1], syma)
    B, IB, CB = contracttree(tensors, network, conjlist, tree[2], symb)
    IC = tuple(output...)

    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)
    T = promote_type(eltype(A), eltype(B))
    # end result: don't use cache
    C = similar_from_indices(T, oindA, oindB, indCinoAB, (), A, B, CA, CB)
    if sym !== nothing
        symcontract = (Symbol(sym, "_a′"), Symbol(sym, "_b′"), Symbol(sym, "_c′"))
    else
        symcontract = nothing
    end
    contract!(true, A, CA, B, CB, false, C,
                oindA, cindA, oindB, cindB, indCinoAB, (), symcontract)
    return C
end

function contracttree(tensors, network, conjlist, tree, sym)
    @nospecialize
    if tree isa Int
        return tensors[tree], tuple(network[tree]...), (conjlist[tree] ? :C : :N)
    end

    if sym !== nothing
        syma = Symbol(sym, "_a")
        symb = Symbol(sym, "_b")
    else
        syma = nothing
        symb = nothing
    end
    A, IA, CA = contracttree(tensors, network, conjlist, tree[1], syma)
    B, IB, CB = contracttree(tensors, network, conjlist, tree[2], symb)
    IC = tuple(symdiff(IA, IB)...)
    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)
    T = promote_type(eltype(A), eltype(B))
    if sym !== nothing
        symc = Symbol(sym, "_c")
        C = cached_similar_from_indices(symc, T, oindA, oindB, indCinoAB, (), A, B, CA, CB)
    else
        C = similar_from_indices(T, oindA, oindB, indCinoAB, (), A, B, CA, CB)
    end
    if sym !== nothing
        symcontract = (Symbol(sym, "_a′"), Symbol(sym, "_b′"), Symbol(sym, "_c′"))
    else
        symcontract = nothing
    end
    contract!(true, A, CA, B, CB, false, C,
                oindA, cindA, oindB, cindB, indCinoAB, (), symcontract)
    return C, IC, :N
end


text:
"""
    horner(c,x)

Evaluate a polynomial whose coefficients are given in descending order in `c`,
at the point `x`, using Horner's rule.
"""
function horner(c,x)

n = length(c)
y = c[1]
for k = 2:n
    y = x*y + c[k]
end

return y

end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
using Mimi

_default_years = 2000:2100
_default_regions = [:A, :B]

function run_getindex(; years = collect(_default_years), regions = _default_regions)
    # Test with one scalar parameter, one 1-D timestep array, and one 2-D timestep array
    types = [Mimi.ScalarModelParameter{Float64}, _get_timesteparray_type(years, 1), _get_timesteparray_type(years, 2)]

    names = [:d1, :d2, :d3]
    values = [
        types[1](4.),
        types[2](Array{Union{Missing, Float64}, 1}(rand(length(years)))),
        types[3](Array{Union{Missing, Float64}, 2}(rand(length(years), length(regions))))
        ]

    datum = NamedTuple{(names...,), Tuple{types...,}}(values)

    clock = _get_clock(years)

    while ! Mimi.finished(clock)
        ts = Mimi.timestep(clock)
        _run_timestep(datum, ts)
        Mimi.advance(clock)
    end
end

function _run_timestep(datum::NamedTuple, ts)
    datum.d1
    datum.d2[ts]
    datum.d3[ts, :]
    nothing
end

function _get_timesteparray_type(years, num_dims, dtype=Float64)
    if Mimi.isuniform(years)
        first, stepsize = Mimi.first_and_step(years)
        last = years[last]
        T = Mimi.TimestepArray{Mimi.FixedTimestep{first, stepsize, last}, Union{dtype, Missing}, num_dims}
    else
        T = Mimi.TimestepArray{Mimi.VariableTimestep{(years...,)}, Union{dtype, Missing}, num_dims}
    end
    return T
end

function _get_clock(years)
    if Mimi.isuniform(years)
        last = years[end]
        first, stepsize = Mimi.first_and_step(years)
        return Mimi.Clock{Mimi.FixedTimestep}(first, stepsize, last)
    else
        last_index = findfirst(isequal(last), years)
        times = (years[1:last_index]...,)
        return Mimi.Clock{Mimi.VariableTimestep}(times)
    end
end


text:
"""
Sources data
"""

inf_gen_case1 = StaticSource(1, #number
                :InfBus, #name
                nodes_case1[1], #bus
                1.05, #VR
                0.0, #VI
                0.000001) #Xth

inf_gen_case234 = StaticSource(1, #number
                :InfBus, #name
                nodes_case234[1], #bus
                1.02, #VR
                0.0, #VI
                0.000001) #Xth

inf_gen_case5 = StaticSource(1, #number
                :InfBus, #name
                nodes_case5[1], #bus
                1.00, #VR
                0.0, #VI
                0.000001) #Xth

"""
Generators data
"""

######## Machine Data #########

### Case 1: Classical machine against infinite bus ###
case1_machine = BaseMachine(0.0, #R
                            0.2995, #Xd_p
                            0.7087, #eq_p
                            100.0)  #MVABase

### Case 2: 4th Order Model with AVR (3-bus case) ###
case2_machine2 =  OneDOneQMachine(0.0, #R
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  100.0)   #MVABase

case2_machine3 =  OneDOneQMachine(0.0, #R
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  100.0)   #MVABase

### Case 3: 6th Order Model with AVR (3-bus case) ###
case3_machine2 = SimpleMarconatoMachine(0.0,
                                        1.3125, #Xd
                                        1.2578, #Xq
                                        0.1813, #Xd_p
                                        0.25, #Xq_p
                                        0.14, #Xd_pp
                                        0.18, #Xq_pp
                                        5.89, #Td0_p
                                        0.6, #Tq0_p
                                        0.5, #Td0_pp
                                        0.023, #Tq0_pp
                                        0.0, #T_AA
                                        100.0) #MVABase

case3_machine3 = SimpleMarconatoMachine(0.0,
                                        1.3125, #Xd
                                        1.2578, #Xq
                                        0.1813, #Xd_p
                                        0.25, #Xq_p
                                        0.14, #Xd_pp
                                        0.18, #Xq_pp
                                        5.89, #Td0_p
                                        0.6, #Tq0_p
                                        0.5, #Td0_pp
                                        0.023, #Tq0_pp
                                        0.0, #T_AA
                                        100.0) #MVABase

### Case 4: 8th Order Model with AVR (3-bus case) ###
case4_machine2 = MarconatoMachine(0.0,
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  0.14, #Xd_pp
                                  0.18, #Xq_pp
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  0.5, #Td0_pp
                                  0.023, #Tq0_pp
                                  0.0, #T_AA
                                  100.0) #MVABase

case4_machine3 = MarconatoMachine(0.0,
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  0.14, #Xd_pp
                                  0.18, #Xq_pp
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  0.5, #Td0_pp
                                  0.023, #Tq0_pp
                                  0.0, #T_AA
                                  100.0) #MVABase

### Case 5: 4th Order Model with AVR + TG + Multishaft ###
case5_machine =  OneDOneQMachine(0.0, #R
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  100.0)   #MVABase

######## Shaft Data #########

### Shaft for Case 1 ###
case1_shaft = SingleMass(3.148, #H
                         2.0) #D

### Shafts for Gens 2 and 3: Cases 2, 3 and 4 ###
case234_shaft2 = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                            0.0) #D

case234_shaft3 = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                            0.0) #D

case5_shaft = FiveMassShaft(3.01, #5.148, #H
                            0.3348, #H_hp
                            0.7306, #H_ip
                            0.8154, #H_lp
                            0.0452, #H_ex,
                            0.0, #2.0, #D
                            0.5180, #D_hp
                            0.2240, #D_ip
                            0.2240, #D_lp
                            0.1450, #D_ex
                            0.0518, #D_12
                            0.0224, #D_23
                            0.0224, #D_34
                            0.0145, #D_45
                            33.07, #K_hp
                            28.59, #K_ip
                            44.68, #K_lp
                            21.984) #K_ex

######## PSS Data #########
cases_no_pss = PSSFixed(0.0)

######## TG Data #########

### No TG for Cases 1, 2, 3, 4 ###
case1234_no_tg = TGFixed(1.0) #eff

### TG for Case 5 ###
case5_tg = TGTypeII(0.05, #R
                    1.0, #T1
                    2.0, #T2
                    1.5, #τ_max
                    0.1) #τ_min

######## AVR Data #########

### AVR for Case 1 ###
case1_avr = AVRFixed(0.0) #Vf not applicable in Classic Machines

### AVRs for Case 2, 3, 4 and 5 ###
case2345_avr2 = AVRTypeI(20.0, #Ka - Gain
                        0.01, #Ke
                        0.063, #Kf
                        0.2, #Ta
                        0.314, #Te
                        0.35, #Tf
                        0.001, #Tr
                        5.0, #Vrmax
                        -5.0, #Vrmin
                        0.0039, #Ae - 1st ceiling coefficient
                        1.555) #Be - 2nd ceiling coefficient

case2345_avr3 = AVRTypeI(20.0, #Ka - Gain
                        0.01, #Ke
                        0.063, #Kf
                        0.2, #Ta
                        0.314, #Te
                        0.35, #Tf
                        0.001, #Tr
                        5.0, #Vrmax
                        -5.0, #Vrmin
                        0.0039, #Ae - 1st ceiling coefficient
                        1.555) #Be - 2nd ceiling coefficient

######################### Generators ########################

### Case 1 Generator ###
case1_gen = DynGenerator(1, #Number
                         :Case1Gen,
                         nodes_case1[2], #bus
                         1.0, # ω_ref,
                         1.0, #V_ref
                         0.5, #P_ref
                         case1_machine, #machine
                         case1_shaft, #shaft
                         case1_avr, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

### Case 2 Generators ###
case2_gen2 = DynGenerator(1, #Number
                         :Case2Gen2,
                         nodes_case234[2], #bus
                         1.0, # ω_ref,
                         1.0142, #V_ref
                         1.0, #P_ref
                         case2_machine2, #machine
                         case234_shaft2, #shaft
                         case2345_avr2, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

case2_gen3 = DynGenerator(2, #Number
                         :Case2Gen3,
                         nodes_case234[3], #bus
                         1.0, # ω_ref,
                         1.0059, #V_ref
                         1.0, #P_ref
                         case2_machine3, #machine
                         case234_shaft3, #shaft
                         case2345_avr3, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

### Case 3 Generators ###
case3_gen2 = DynGenerator(1, #Number
                         :Case3Gen2,
                         nodes_case234[2], #bus
                         1.0, # ω_ref,
                         1.0142, #V_ref
                         1.0, #P_ref
                         case3_machine2, #machine
                         case234_shaft2, #shaft
                         case2345_avr2, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

case3_gen3 = DynGenerator(2, #Number
                         :Case3Gen3,
                         nodes_case234[3], #bus
                         1.0, # ω_ref,
                         1.0059, #V_ref
                         1.0, #P_ref
                         case3_machine3, #machine
                         case234_shaft3, #shaft
                         case2345_avr3, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

### Case 4 Generators ###
case4_gen2 = DynGenerator(1, #Number
                         :Case4Gen2,
                         nodes_case234[2], #bus
                         1.0, # ω_ref,
                         1.0142, #V_ref
                         1.0, #P_ref
                         case4_machine2, #machine
                         case234_shaft2, #shaft
                         case2345_avr2, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

case4_gen3 = DynGenerator(2, #Number
                         :Case4Gen3,
                         nodes_case234[3], #bus
                         1.0, # ω_ref,
                         1.0059, #V_ref
                         1.0, #P_ref
                         case4_machine3, #machine
                         case234_shaft3, #shaft
                         case2345_avr3, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

### Case 5 Generator ###
case5_gen = DynGenerator(1, #Number
                         :Case5Gen,
                         nodes_case5[2], #bus
                         1.0, # ω_ref,
                         1.0155, #V_ref
                         0.5, #P_ref
                         case5_machine, #machine
                         case5_shaft, #shaft
                         case2345_avr2, #avr
                         case5_tg, #tg
                         #case1234_no_tg,
                         cases_no_pss) #pss

######################### Dynamical System ########################

case1_DynSystem = DynamicSystem(nodes_case1,
                                branch_case1,
                                [case1_gen],
                                vcat(inf_gen_case1,loads_case1),
                                100.0,
                                60.0)

case2_DynSystem = DynamicSystem(nodes_case234,
                                branch_case234,
                                [case2_gen2, case2_gen3],
                                vcat(inf_gen_case234,loads_case234),
                                100.0,
                                60.0)

case3_DynSystem = DynamicSystem(nodes_case234,
                                branch_case234,
                                [case3_gen2, case3_gen3],
                                vcat(inf_gen_case234,loads_case234),
                                100.0,
                                60.0)

case4_DynSystem = DynamicSystem(nodes_case234,
                                branch_case234,
                                [case4_gen2, case4_gen3],
                                vcat(inf_gen_case234,loads_case234),
                                100.0,
                                60.0)

case5_DynSystem = DynamicSystem(nodes_case5,
                                branch_case5,
                                [case5_gen],
                                vcat(inf_gen_case5,loads_case5),
                                100.0,
                                60.0)


text:
using JuMP
using BenchmarkTools

function sum_iterate(con_refs)
    x = 0.0
    for con_ref in con_refs
        x += dual(con_ref)
    end
    return x
end

function sum_index(con_refs)
    x = 0.0
    for i in eachindex(con_refs)
        x += dual(con_refs[i])
    end
    return x
end

function dense_axis_constraints(n)
    model = Model()
    mock = MOIU.MockOptimizer(
        MOIU.Model{Float64}(),
        eval_variable_constraint_dual = false,
    )
    MOIU.set_mock_optimize!(
        mock,
        mock -> MOIU.mock_optimize!(
            mock,
            zeros(n),
            (MOI.SingleVariable, MOI.EqualTo{Float64}) => ones(n - 1),
        ),
    )
    MOIU.reset_optimizer(model, mock)

    @variable(model, x[1:n])
    set = MOI.EqualTo(0.0)
    con_refs = @time @constraint(model, [i = 2:n], x[i] in set)
    optimize!(model)
    @assert sum_iterate(con_refs) == n - 1
    @btime sum_iterate($con_refs)
    @assert sum_index(con_refs) == n - 1
    @btime sum_index($con_refs)
    return
end
function sparse_axis_constraints(n)
    model = Model()
    mock = MOIU.MockOptimizer(
        MOIU.Model{Float64}(),
        eval_variable_constraint_dual = false,
    )
    MOIU.set_mock_optimize!(
        mock,
        mock -> MOIU.mock_optimize!(
            mock,
            zeros(n),
            (MOI.SingleVariable, MOI.EqualTo{Float64}) => ones(div(n, 2)),
        ),
    )
    MOIU.reset_optimizer(model, mock)

    @variable(model, x[1:n])
    set = MOI.EqualTo(0.0)
    con_refs = @time @constraint(model, [i = 1:n; iseven(i)], x[i] in set)
    optimize!(model)
    @assert sum_iterate(con_refs) == div(n, 2)
    @btime sum_iterate($con_refs)
    @assert sum_index(con_refs) == div(n, 2)
    @btime sum_index($con_refs)
    return
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestCompositesCore

using Test
using MLJBase
using Tables
import MLJBase
using ..Models
using ..TestUtilities
using CategoricalArrays
using OrderedCollections
import Random.seed!
seed!(1234)

mutable struct Rubbish <: DeterministicComposite
    model_in_network
    model_not_in_network
    some_other_variable
end

knn = KNNRegressor()
model = Rubbish(knn, OneHotEncoder(), 42)
X, y = make_regression(10, 2)

@testset "logic for composite model update - fallback()" begin
    Xs = source(X)
    ys = source(y)
    mach0 = machine(Standardizer(), Xs)
    W = transform(mach0, Xs)
    mach1 = machine(model.model_in_network, W, ys)
    yhat = predict(mach1, W)
    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    fitresult, cache, _ = return!(mach, model, 0)
    network_model_names = getfield(fitresult, :network_model_names)
    @test network_model_names == [:model_in_network, nothing]
    old_model = cache.old_model
    glb_node = MLJBase.glb(mach)
    @test !MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # don't fallback if mutating field for a network model:
    model.model_in_network.K = 24
    @test !MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # do fallback if replacing field for a network model:
    model.model_in_network = KNNRegressor()
    @test MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # return to original state:
    model.model_in_network = knn
    @test !MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # do fallback if a non-network field changes:
    model.model_not_in_network.features = [:x1,]
    @test MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # return to original state:
    model.model_not_in_network = OneHotEncoder()
    @test !MLJBase.fallback(model, old_model, network_model_names, glb_node)

    # do fallback if any non-model changes:
    model.some_other_variable = 123412
    @test MLJBase.fallback(model, old_model, network_model_names, glb_node)

end

model = Rubbish(KNNRegressor(), Standardizer(), 42)

function MLJBase.fit(model::Rubbish, verbosity, X, y)
    Xs = source(X)
    ys = source(y)
    mach1 = machine(model.model_in_network, Xs, ys)
    yhat = predict(mach1, Xs)
    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    return!(mach, model, verbosity)
end

# `model` is instance of `Rubbish`
mach = fit!(machine(model, X, y), verbosity=0)

@testset "logic for composite model update - fit!" begin

    # immediately refit:
    @test_model_sequence(fit!(mach), [(:skip, model), ])

    # mutate a field for a network model:
    model.model_in_network.K = 24
    @test_model_sequence(fit!(mach),
                         [(:update, model), (:update, model.model_in_network)])

    # immediately refit:
    @test_model_sequence(fit!(mach), [(:skip, model), ])

    # replace a field for a network model:
    model.model_in_network = KNNRegressor()
    @test_model_sequence(fit!(mach),
                         [(:update, model), (:train, model.model_in_network)])

    # immediately refit:
    @test_model_sequence(fit!(mach), [(:skip, model), ])

    # mutate a field for a model not in network:
    model.model_not_in_network.features = [:x1,]
    @test_model_sequence(fit!(mach),
                         [(:update, model), (:train, model.model_in_network)])

    # immediately refit:
    @test_model_sequence(fit!(mach), [(:skip, model), ])

    # mutate some field that is not a model:
    model.some_other_variable = 123412
    @test_model_sequence(fit!(mach),
                         [(:update, model), (:train, model.model_in_network)])
end

N = 50
Xin = (a=rand(N), b=rand(N), c=rand(N));
yin = rand(N);

train, test = partition(eachindex(yin), 0.7);
Xtrain = MLJBase.selectrows(Xin, train);
ytrain = yin[train];

ridge_model = FooBarRegressor(lambda=0.1)
selector_model = FeatureSelector()

@testset "first test of hand-exported network" begin
    composite = SimpleDeterministicCompositeModel(model=ridge_model,
                                                  transformer=selector_model)

    fitresult, cache, rep = MLJBase.fit(composite, 0, Xtrain, ytrain);

    # test data anonymity:
    ss = sources(glb(values(MLJBase.signature(fitresult))...))
    @test all(isempty, ss)

    # to check internals:
    ridge = MLJBase.machines(fitresult.predict)[1]
    selector = MLJBase.machines(fitresult.predict)[2]
    ridge_old = deepcopy(ridge)
    selector_old = deepcopy(selector)

    # this should trigger no retraining:
    fitresult, cache, rep =
        @test_logs(
            (:info, r"^Not"),
            (:info, r"^Not"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain));
    @test ridge.fitresult == ridge_old.fitresult
    @test selector.fitresult == selector_old.fitresult

    # this should trigger update of selector and training of ridge:
    selector_model.features = [:a, :b]
    fitresult, cache, rep =
        @test_logs(
            (:info, r"^Updating"),
            (:info, r"^Training"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain));
    @test ridge.fitresult != ridge_old.fitresult
    @test selector.fitresult != selector_old.fitresult
    ridge_old = deepcopy(ridge)
    selector_old = deepcopy(selector)

    # this should trigger updating of ridge only:
    ridge_model.lambda = 1.0
    fitresult, cache, rep =
        @test_logs(
            (:info, r"^Not"),
            (:info, r"^Updating"),
            MLJBase.update(composite, 2, fitresult, cache, Xtrain, ytrain));
    @test ridge.fitresult != ridge_old.fitresult
    @test selector.fitresult == selector_old.fitresult

    predict(composite, fitresult, MLJBase.selectrows(Xin, test));

    Xs = source(Xtrain)
    ys = source(ytrain)

    mach = machine(composite, Xs, ys)
    yhat = predict(mach, Xs)
    fit!(yhat, verbosity=0)
    composite.transformer.features = [:b, :c]
    fit!(yhat, verbosity=0)
    fit!(yhat, rows=1:20, verbosity=0)
    yhat(MLJBase.selectrows(Xin, test));

end

mutable struct WrappedRidge <: DeterministicComposite
    ridge
end

# julia bug? If I return the following test to a @testset block, then
# the test marked with ******* fails (bizarre!)
#@testset "second test of hand-exported network" begin
function MLJBase.fit(model::WrappedRidge, verbosity::Integer, X, y)
    Xs = source(X)
    ys = source(y)

    stand = Standardizer()
    standM = machine(stand, Xs)
    W = transform(standM, Xs)

    boxcox = UnivariateBoxCoxTransformer()
    boxcoxM = machine(boxcox, ys)
    z = transform(boxcoxM, ys)

    ridgeM = machine(model.ridge, W, z)
    zhat = predict(ridgeM, W)
    yhat = inverse_transform(boxcoxM, zhat)

    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    return!(mach, model, verbosity)
end

MLJBase.input_scitype(::Type{<:WrappedRidge}) =
    Table(Continuous)
MLJBase.target_scitype(::Type{<:WrappedRidge}) =
    AbstractVector{<:Continuous}

ridge = FooBarRegressor(lambda=0.1)
model_ = WrappedRidge(ridge)
mach = machine(model_, Xin, yin)
id = objectid(mach)
fit!(mach, verbosity=0)
@test  objectid(mach) == id  # *********
yhat=predict(mach, Xin);
ridge.lambda = 1.0
fit!(mach, verbosity=0)
@test predict(mach, Xin) != yhat

#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))...]

# A wrap of above model:
mutable struct WrappedDummyClusterer <: UnsupervisedComposite
    model
end
WrappedDummyClusterer(; model=DummyClusterer()) =
    WrappedDummyClusterer(model)

@testset "third test of hand-exported network" begin
    function MLJBase.fit(model::WrappedDummyClusterer, verbosity::Int, X)
        Xs = source(X)
        W = transform(machine(OneHotEncoder(), Xs), Xs)
        m = machine(model.model, W)
        yhat = predict(m, W)
        Wout = transform(m, W)
        foo = node(η -> first(η), yhat)
        mach = machine(Unsupervised(),
                       Xs;
                       predict=yhat,
                       transform=Wout,
                       report=(foo=foo,))
        return!(mach, model, verbosity)
    end
    X, _ = make_regression(10, 5);
    model = WrappedDummyClusterer(model=DummyClusterer(n=2))
    mach = fit!(machine(model, X), verbosity=0)
    model.model.n = 3
    fit!(mach, verbosity=0)
    @test transform(mach, X) == selectcols(X, 1:3)
    r = report(mach)
    @test r.model.centres == MLJBase.matrix(X)[1:3,:]
    @test r.foo == predict(mach, rows=:)[1]
    fp = fitted_params(mach)
    @test :model in keys(fp)
    levs = fp.model.fitresult
    @test predict(mach, X) == fill(levs[1], 10)
end

## NETWORK WITH MULTIPLE NODES REPORTING STATE/ REFIT

mutable struct TwoStages <: DeterministicComposite
    model1
    model2
    model3
end

function MLJBase.fit(m::TwoStages, verbosity, X, y)
    Xs = source(X)
    ys = source(y)
    mach1 = machine(m.model1, Xs, ys)
    mach2 = machine(m.model2, Xs, ys)
    ypred1 = MLJBase.predict(mach1, Xs)
    ypred2 = MLJBase.predict(mach2, Xs)
    Y = MLJBase.table(hcat(ypred1, ypred2))
    mach3 = machine(m.model3, Y, ys)
    ypred3 = MLJBase.predict(mach3, Y)
    μpred = node(x->mean(x), ypred3)
    σpred = node((x, μ)->mean((x.-μ).^2), ypred3, μpred)
    mach = machine(Deterministic(),
                   Xs,
                   ys;
                   predict=ypred3,
                   report=(μpred=μpred,
                           σpred=σpred))
    return!(mach, m, verbosity)
end

@testset "Test exported-network with multiple saved nodes and refit" begin
    X, y = make_regression(100, 3)
    model3 = FooBarRegressor(lambda=1)
    twostages = TwoStages(FooBarRegressor(lambda=0.1),
                          FooBarRegressor(lambda=10), model3)
    mach = machine(twostages, X, y)
    fit!(mach, verbosity=0)
    rep = report(mach)
    # All machines have been fitted once
    @test rep.machines[1].state ==
        rep.machines[2].state ==
        rep.machines[3].state == 1
    # Retrieve current values of interest
    μpred = rep.μpred
    σpred = rep.σpred
    # Change model3 and refit
    model3.lambda = 10
    fit!(mach, verbosity=0)
    rep = report(mach)
    # Machines 1,2 have been fitted once and machine 3 twice
    @test rep.machines[1].state == rep.machines[2].state == 1
    @test rep.machines[3].state == 2
    # The new values have been updated
    @test rep.μpred != μpred
    @test rep.σpred != σpred
end

## COMPOSITE WITH COMPONENT MODELS STORED IN NTUPLE

# `modelnames` is a tuple of `Symbol`s, one for each `model` in `models`:
mutable struct Averager{modelnames} <: DeterministicComposite
    models::NTuple{<:Any,Deterministic}
    weights::Vector{Float64}
    Averager(modelnames, models, weights) =
        new{modelnames}(models, weights)
end

# special kw constructor, allowing one to specify the property names
# to be attributed to each component model (see below):
function Averager(; weights=Float64[], named_models...)
        nt = NamedTuple(named_models)
    modelnames = keys(nt)
    models = values(nt)
    return Averager(modelnames, models, weights)
end

# for example:
averager = Averager(weights=[1, 1],
                    model1=KNNRegressor(K=3),
                    model2=RidgeRegressor())

# so we can do `averager.model1` and `averager.model2`:
Base.propertynames(::Averager{modelnames}) where modelnames =
        tuple(:weights, modelnames...)
function Base.getproperty(averager::Averager{modelnames},
                          name::Symbol) where modelnames
    name === :weights && return getfield(averager, :weights)
    models = getfield(averager, :models)
    for j in eachindex(modelnames)
        name === modelnames[j] && return models[j]
    end
    error("type Averager has no field $name")
end

# overload multiplication of a node by a matrix:
import Base.*
*(preds::Node, weights) = node(p->p*weights, preds)

# learning network wrapped in a fit method:
function MLJBase.fit(averager::Averager{modelnames},
                     verbosity,
                     X,
                     y) where modelnames

    Xs = source(X)
    ys = source(y)

    weights = averager.weights

    machines = [machine(getproperty(averager, name), Xs, ys) for
                name in modelnames]
    predictions = hcat([predict(mach, Xs) for mach in machines]...)
    yhat = (1/sum(weights))*(predictions*weights)

    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    return!(mach, averager, verbosity)
end

@testset "composite with component models stored in ntuple" begin
    X, y = make_regression(10, 3);
    mach = machine(averager, X, y)
    fit!(mach, verbosity=0)
    fp = fitted_params(mach)
    @test keys(fp.model1) == (:tree, )
    @test keys(fp.model2) == (:coefficients, :intercept)
    r = report(mach)
    @test r.model1 == NamedTuple()
    @test r.model2 == NamedTuple()
    range(averager, :(model1.K), lower=2, upper=3)
end

## DATA FRONT-END IN AN EXPORTED LEARNING NETWORK

mutable struct Scale <: MLJBase.Static
    scaling::Float64
end

function MLJBase.transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X * s.scaling
    MLJBase.table(s.scaling * MLJBase.matrix(X), prototype=X)
end

function MLJBase.inverse_transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X / s.scaling
    MLJBase.table(MLJBase.matrix(X) / s.scaling, prototype=X)
end

mutable struct ElephantModel <: ProbabilisticComposite
    scaler
    clf
    cache::Bool
end

function MLJBase.fit(model::ElephantModel, verbosity, X, y)

    Xs = source(X)
    ys = source(y)

    scaler = model.scaler
    mach1 = machine(scaler, cache=model.cache)
    W = transform(mach1, Xs)

    # a classifier with reformat front-end:
    clf = model.clf
    mach2 = machine(clf, W, ys, cache=model.cache)
    yhat = predict(mach2, W)

    mach = machine(Probabilistic(), Xs, ys, predict=yhat)
    return!(mach, model, verbosity)
end

@testset "reformat/selectrows logic in composite model" begin

    X = (x1=ones(5), x2=ones(5))
    y = categorical(collect("abaaa"))
    model = ElephantModel(Scale(2.0),
                        ConstantClassifier(testing=true, bogus=1.0),
                        true)
    mach = machine(model, X, y, cache=false)

    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(mach, verbosity=0, rows=1:3)
               )
    @test mach.state == 1

    # new clf hyperparmater (same rows) means no reformatting or resampling:
    model.clf.bogus = 10
    @test_logs fit!(mach, verbosity=0, rows=1:3)
    @test mach.state == 2

    # however changing an upstream hyperparameter forces reformatting
    # and resampling:
    model.scaler.scaling = 3.1
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(mach, verbosity=0, rows=1:3))

end

@testset "operation nodes that are source nodes" begin

    mutable struct BananaComposite <: UnsupervisedComposite
        stand
    end
    BananaComposite(; stand=Standardizer()) = BananaComposite(stand)

    function MLJBase.fit(model::BananaComposite, verbosity, X)

        Xs = source(X)
        mach1 = machine(model.stand, Xs)
        X2 = transform(mach1, Xs)

        # node for the inverse_transform:

        network_mach = machine(Unsupervised(), Xs, transform=X2, inverse_transform=Xs)
        return!(network_mach, model, verbosity)

    end

    X = (x = Float64[1, 2, 3],)
    mach = machine(BananaComposite(), X)
    fit!(mach, verbosity=0, force=true)
    @test transform(mach, X).x ≈ Float64[-1, 0, 1]
    @test inverse_transform(mach, X) == X

end

end # module
true


text:
using Test

include("overt_to_mip_test.jl")
include("acc_feasibility_test.jl")
include("acc_reachability_test.jl")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDAFastRegistrationLib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDAFastRegistrationLib")
JLLWrappers.@generate_main_file("CUDAFastRegistrationLib", UUID("b2992ab9-cfb1-5a7d-8be5-46e2b4874ed9"))
end  # module CUDAFastRegistrationLib_jll


text:
abstract type AbstractPenalty{T} end
abstract type AbstractCPUPenalty{T} <: AbstractPenalty{T} end
abstract type AbstractProjection end

mutable struct PowerPenalty{T} <: AbstractCPUPenalty{T}
    p::T
end
@inline (P::PowerPenalty)(x) = x^(P.p)

mutable struct RationalPenalty{T} <: AbstractCPUPenalty{T}
    p::T
end
@inline (R::RationalPenalty)(x) = x / (1 + R.p * (1 - x))

mutable struct SinhPenalty{T} <: AbstractCPUPenalty{T}
    p::T
end
@inline (R::SinhPenalty)(x) = sinh(R.p * x) / sinh(R.p)

struct ProjectedPenalty{T,Tpen<:AbstractPenalty{T},Tproj} <: AbstractCPUPenalty{T}
    penalty::Tpen
    proj::Tproj
end
function ProjectedPenalty(penalty::AbstractPenalty{T}) where {T}
    return ProjectedPenalty(penalty, HeavisideProjection(10 * one(T)))
end
@inline (P::ProjectedPenalty)(x) = P.penalty(P.proj(x))
@forward_property ProjectedPenalty penalty

mutable struct HeavisideProjection{T} <: AbstractProjection
    β::T
end
@inline (P::HeavisideProjection)(x) = 1 - exp(-P.β * x) + x * exp(-P.β)
mutable struct SigmoidProjection{T} <: AbstractProjection
    β::T
end
@inline (P::SigmoidProjection)(x) = 1 / (1 + exp((P.β + 1) * (-x + 0.5)))

import Base: copy
copy(p::TP) where {TP<:AbstractPenalty} = TP(p.p)
copy(p::HeavisideProjection) = HeavisideProjection(p.β)
copy(p::SigmoidProjection) = SigmoidProjection(p.β)
copy(p::ProjectedPenalty) = ProjectedPenalty(copy(p.penalty), copy(p.proj))

function Utilities.setpenalty!(P::AbstractPenalty, p)
    P.p = p
    return P
end
function Utilities.setpenalty!(P::ProjectedPenalty, p)
    P.penalty.p = p
    return P
end


text:
using SafeTestsets

@safetestset "SBM tests" begin include("test_sbm.jl") end
@safetestset "Dataset tests" begin include("test_datasets.jl") end
@safetestset "Opt Methods tests" begin include("test_estimate.jl") end
@safetestset "Opt Results tests" begin include("test_results.jl") end
@safetestset "Heuristics tests" begin include("test_heuristic.jl") end
@safetestset "Exact tests" begin include("test_exact.jl") end


text:
#NOTE deprecating an substituting with Workspace
"""
    initialise(::ObsScheme, P, m, y::StartingPtPrior{T}, ::S, fpt)

Initialise the workspace for MCMC algorithm. Initialises containers for driving
Wiener processes `WWᵒ` & `WW`, for diffusion processes `XXᵒ` & `XX`, for
diffusion Law `Pᵒ` (parametetrised by proposal parameters) and defines the type
of Wiener process `Wnr`.
"""
function initialise(::ObsScheme, P, m, yPr::StartingPtPrior{T}, ::S,
                    fpt) where {ObsScheme <: AbstractObsScheme,T,S}
    y = startPt(yPr)
    Pᵒ = deepcopy(P)
    TW = typeof(sample([0], Wiener{S}()))
    TX = typeof(SamplePath([], zeros(T, 0)))
    XXᵒ = Vector{TX}(undef,m)
    WWᵒ = Vector{TW}(undef,m)
    Wnr = Wiener{S}()
    for i in 1:m
        WWᵒ[i] = Bridge.samplepath(P[i].tt, zero(S))
        sample!(WWᵒ[i], Wnr)
        WWᵒ[i], XXᵒ[i] = forcedSolve(Euler(), y, WWᵒ[i], P[i])    # this will enforce adherence to domain
        while !checkFpt(ObsScheme(), XXᵒ[i], fpt[i])
            sample!(WWᵒ[i], Wnr)
            forcedSolve!(Euler(), XXᵒ[i], y, WWᵒ[i], P[i])    # this will enforce adherence to domain
        end
        y = XXᵒ[i].yy[end]
    end
    y = startPt(yPr)
    ll = logpdf(yPr, y)
    ll += pathLogLikhd(ObsScheme(), XXᵒ, P, 1:m, fpt, skipFPT=true)
    ll += lobslikelihood(P[1], y)

    XX = deepcopy(XXᵒ)
    WW = deepcopy(WWᵒ)
    # needed for proper initialisation of the Crank-Nicolson scheme
    yPr = invStartPt(y, yPr, P[1])

    Wnr, WWᵒ, WW, XXᵒ, XX, Pᵒ, ll, yPr
end

#NOTE deprecated
"""
    swapXX!(𝔅::ChequeredBlocking, XX)

Swap containers between `XX` and `𝔅.XX`
"""
function swapXX!(𝔅::BlockingSchedule, XX)
    for block in 𝔅.blocks[𝔅.idx]
        swap!(XX, 𝔅.XX, block)
    end
end

#NOTE deprecated
"""
    swapXX!(𝔅::NoBlocking, XX)

nothing to do
"""
swapXX!(𝔅::NoBlocking, XX) = nothing

#NOTE deprecated
"""
    impute!(::ObsScheme, 𝔅::ChequeredBlocking, Wnr, y, WWᵒ, WW, XXᵒ, XX, P, ll,
            fpt; ρ=0.0, verbose=false, it=NaN, headStart=false) where
            ObsScheme <: AbstractObsScheme -> acceptedLogLikhd, acceptDecision

Imputation step of the MCMC scheme (without blocking).
...
# Arguments
- `::ObsScheme`: observation scheme---first-passage time or partial observations
- `𝔅`: object with relevant information about blocking
- `Wnr`: type of the Wiener process
- `yPr`: prior over the starting point of the diffusion path
- `WWᵒ`: containers for proposal Wiener paths
- `WW`: containers with old Wiener paths
- `XXᵒ`: containers for proposal diffusion paths
- `XX`: containers with old diffusion paths
- `P`: laws of the diffusion path (proposal and target)
- `11`: log-likelihood of the old (previously accepted) diffusion path
- `fpt`: info about first-passage time conditioning
- `ρ`: memory parameter for the Crank-Nicolson scheme
- `verbose`: whether to print updates info while sampling
- `it`: iteration index of the MCMC algorithm
- `headStart`: flag for whether to 'ease into' fpt conditions
...
"""
function impute!_deprecated(::ObsScheme, 𝔅::ChequeredBlocking, Wnr, yPr, WWᵒ, WW, XXᵒ, XX,
                 P, ll, fpt; ρ=0.0, verbose=false, it=NaN, headStart=false,
                 solver::ST=Ralston3()) where
                 {ObsScheme <: AbstractObsScheme, ST}
    θ = params(P[1].Target)             # current parameter
    𝔅 = next(𝔅, XX, θ)
    solveBackRec!(𝔅, 𝔅.P, ST())         # compute (H, Hν, c) for given blocks

    swapXX!(𝔅, XX)                      # move current path to object 𝔅
    noiseFromPath!(𝔅, 𝔅.XX, 𝔅.WW, 𝔅.P) # find noise WW that generates XX under 𝔅

    # compute white noise generating starting point under 𝔅
    yPr𝔅 = invStartPt(𝔅.XX[1].yy[1], yPr, 𝔅.P[1])

    for (blockIdx, block) in enumerate(𝔅.blocks[𝔅.idx])
        blockFlag = Val{block[1]}()
        y = 𝔅.XX[block[1]].yy[1]       # current starting point

        # set the starting point for the block
        yᵒ, yPrᵒ = proposalStartPt(𝔅, blockFlag, y, yPr𝔅, 𝔅.P[block[1]], ρ)

        # sample path in block
        sampleSegments!(block, Wnr, 𝔅.WW, 𝔅.WWᵒ, 𝔅.P, yᵒ, 𝔅.XXᵒ, ρ)
        setEndPtManually!(𝔅, blockIdx, block)

        # loglikelihoods
        llᵒ = startPtLogPdf(blockFlag, yPrᵒ, yᵒ)
        llᵒ += pathLogLikhd(ObsScheme(), 𝔅.XXᵒ, 𝔅.P, block, fpt)
        llᵒ += lobslikelihood(𝔅.P[block[1]], yᵒ)

        llPrev = startPtLogPdf(blockFlag, yPr𝔅, y)
        llPrev += pathLogLikhd(ObsScheme(), 𝔅.XX, 𝔅.P, block, fpt; skipFPT=true)
        llPrev += lobslikelihood(𝔅.P[block[1]], y)

        printInfo(verbose, it, value(llPrev), value(llᵒ), "impute")
        if acceptSample(llᵒ-llPrev, verbose)
            swap!(𝔅.XX, 𝔅.XXᵒ, block)
            registerAccpt!(𝔅, blockIdx, true)
            yPr𝔅 = yPrᵒ # can do something non-trivial only for the first block
        else
            registerAccpt!(𝔅, blockIdx, false)
        end
    end
    swapXX!(𝔅, XX) # move accepted path from object 𝔅 to general container XX
    noiseFromPath!(𝔅, XX, WW, P) # compute noise WW that generated XX under law P
    # compute white noise generating starting point under P
    y = XX[1].yy[1]
    yPr = invStartPt(y, yPr𝔅, P[1])

    ll = logpdf(yPr, y) # starting point contribution
    ll += pathLogLikhd(ObsScheme(), XX, P, 1:length(P), fpt; skipFPT=true)
    ll += lobslikelihood(P[1], y)

    # acceptance indicator does not matter for sampling with blocking
    return ll, true, 𝔅, yPr
end

struct Workspace{ObsScheme,S,TX,TW,R,Q}
    Wnr::Wiener{S}
    XXᵒ::Vector{TX}
    XX::Vector{TX}
    WWᵒ::Vector{TW}
    WW::Vector{TW}
    Pᵒ::Vector{R}
    P::Vector{R}
    fpt::Vector
    ρ::Float64 #TODO use vector instead for blocking
    result::Vector{Q}
    resultᵒ::Vector{Q}

    function Workspace(::ObsScheme, P::Vector{R}, m, yPr::StartingPtPrior{T},
                       ::S, fpt, ρ, updtCoord) where {ObsScheme <: AbstractObsScheme,R,T,S}
        y = startPt(yPr)
        Pᵒ = deepcopy(P)
        TW = typeof(sample([0], Wiener{S}()))
        TX = typeof(SamplePath([], zeros(T, 0)))
        XXᵒ = Vector{TX}(undef,m)
        WWᵒ = Vector{TW}(undef,m)
        Wnr = Wiener{S}()
        for i in 1:m
            WWᵒ[i] = Bridge.samplepath(P[i].tt, zero(S))
            sample!(WWᵒ[i], Wnr)
            WWᵒ[i], XXᵒ[i] = forcedSolve(Euler(), y, WWᵒ[i], P[i])    # this will enforce adherence to domain
            while !checkFpt(ObsScheme(), XXᵒ[i], fpt[i])
                sample!(WWᵒ[i], Wnr)
                forcedSolve!(Euler(), XXᵒ[i], y, WWᵒ[i], P[i])    # this will enforce adherence to domain
            end
            y = XXᵒ[i].yy[end]
        end
        y = startPt(yPr)
        ll = logpdf(yPr, y)
        ll += pathLogLikhd(ObsScheme(), XXᵒ, P, 1:m, fpt, skipFPT=true)
        ll += lobslikelihood(P[1], y)

        XX = deepcopy(XXᵒ)
        WW = deepcopy(WWᵒ)
        # needed for proper initialisation of the Crank-Nicolson scheme
        yPr = invStartPt(y, yPr, P[1])

        θ = params(P[1].Target)
        ϑs = [[θ[j] for j in idx(uc)] for uc in updtCoord]

        result = [DiffResults.GradientResult(ϑ) for ϑ in ϑs]
        resultᵒ = [DiffResults.GradientResult(ϑ) for ϑ in ϑs]
        Q = eltype(result)

        new{ObsScheme,S,TX,TW,R,Q}(Wnr, XXᵒ, XX, WWᵒ, WW, Pᵒ, P, fpt, ρ, result, resultᵒ), ll, yPr
    end

    function Workspace(𝓦𝓢::Workspace{ObsScheme,S,TX,TW,R,Q}, new_ρ::Float64
                       ) where {ObsScheme,S,TX,TW,R,Q}
        new{ObsScheme,S,TX,TW,R,Q}(𝓦𝓢.Wnr, 𝓦𝓢.XXᵒ, 𝓦𝓢.XX, 𝓦𝓢.WWᵒ, 𝓦𝓢.WW,
                                   𝓦𝓢.Pᵒ, 𝓦𝓢.P, 𝓦𝓢.fpt, new_ρ, 𝓦𝓢.result,
                                   𝓦𝓢.resultᵒ)
    end
end

"""
    mcmc(::ObsScheme, obs, obsTimes, yPr::StartingPtPrior, w, P˟, P̃, Ls, Σs,
         numSteps, tKernel, priors; fpt=fill(NaN, length(obsTimes)-1), ρ=0.0,
         dt=1/5000, timeChange=true, saveIter=NaN, verbIter=NaN,
         updtCoord=(Val((true,)),), paramUpdt=true, skipForSave=1,
         updtType=(MetropolisHastingsUpdt(),), solver::ST=Ralston3(), warmUp=0)

Gibbs sampler alternately imputing unobserved parts of the path and updating
unknown coordinates of the parameter vector (the latter only if paramUpdt==true)
...
# Arguments
- `::ObsScheme`: observation scheme---first-passage time or partial observations
- `obs`: vector with observations
- `obsTimes`: times of the observations
- `yPr`: prior over the starting point of the diffusion path
- `w`: dummy variable whose type must agree with the type of the Wiener process
- `P˟`: law of the target diffusion (with initial θ₀ set)
- `P̃`: law of the auxiliary process (with initial θ₀ set)
- `Ls`: vector of observation operators (one per each observation)
- `Σs`: vector of covariance matrices of the noise (one per each observaiton)
- `numSteps`: number of mcmc iterations
- `tKernel`: transition kernel (also with initial θ₀ set)
- `priors`: a list of lists of priors
- `τ`: time-change transformation
- `fpt`: info about first-passage time conditioning
- `ρ`: memory parameter for the Crank-Nicolson scheme
- `dt`: time-distance for the path imputation
- `saveIter`: save path `XX` once every `saveIter` many iterations
- `verbIter`: print out progress info once every `verbIter` many iterations
- `updtCoord`: list of objects declaring indices of to-be-updated parameters
- `paramUpdt`: flag for whether to update parameters at all
- `skipForSave`: when saving paths, save only one in every `skipForSave` points
- `updtType`: list of types of updates to cycle through
- `solver`: numerical solver used for computing backward ODEs
- `warmUp`: number of steps for which no parameter update is to be made
...
"""
function mcmc(::Type{K}, ::ObsScheme, obs, obsTimes, yPr::StartingPtPrior, w,
              P˟, P̃, Ls, Σs, numSteps, tKernel, priors, τ;
              fpt=fill(NaN, length(obsTimes)-1), ρ=0.0, dt=1/5000, saveIter=NaN,
              verbIter=NaN, updtCoord=(Val((true,)),), paramUpdt=true,
              skipForSave=1, updtType=(MetropolisHastingsUpdt(),),
              blocking::Blocking=NoBlocking(),
              blockingParams=([], 0.1, NoChangePt()),
              solver::ST=Ralston3(), changePt::CP=NoChangePt(), warmUp=0,
              adaptiveProp=NoAdaptation()
              ) where {K, ObsScheme <: AbstractObsScheme, ST, Blocking, CP}
    P = findProposalLaw( updtType[1], K, obs, obsTimes, P˟, P̃, Ls, Σs, τ; dt=dt, solver=ST(),
                         changePt=CP(getChangePt(blockingParams[3])) )
    m = length(obs)-1
    updtLen = length(updtCoord)
    𝓦𝓢, ll, yPr = Workspace(ObsScheme(), P, m, yPr, w, fpt, ρ, updtCoord)
    init_adaptation!(adaptiveProp, 𝓦𝓢)

    Paths = []
    accImpCounter = 0
    accUpdtCounter = [0 for i in 1:updtLen]
    θ = params(P˟)
    θchain = Vector{typeof(θ)}(undef, (numSteps-warmUp)*updtLen+1)
    θchain[1] = copy(θ)
    recomputeODEs = [any([e in dependsOnParams(P[1].Pt) for e
                         in idx(uc)]) for uc in updtCoord]

    updtStepCounter = 1
    𝔅 = setBlocking(blocking, blockingParams, 𝓦𝓢)
    display(𝔅)
    for i in 1:numSteps
        #print(i, ", ")
        verbose = (i % verbIter == 0)
        i > warmUp && savePath!(Paths, blocking == NoBlocking() ? 𝓦𝓢.XX : 𝔅.XX,
                                (i % saveIter == 0), skipForSave)
        ll, acc, 𝔅, yPr = impute!(𝔅, yPr, 𝓦𝓢, ll; verbose=verbose, it=i,
                                  solver=ST())

        accImpCounter += 1*acc

        if paramUpdt && i > warmUp
            for j in 1:updtLen
                (ll, acc, θ,
                 yPr) = updateParam!(updtType[j], 𝔅, tKernel, θ, updtCoord[j],
                                     yPr, 𝓦𝓢, ll, priors[j], recomputeODEs[j];
                                     solver=ST(), verbose=verbose, it=i, uidx=j)
                accUpdtCounter[j] += 1*acc
                updtStepCounter += 1
                θchain[updtStepCounter] = copy(θ)
                verbose && print("\n")
            end
            verbose && print("------------------------------------------------",
                             "------\n")
        end
        addPath!(adaptiveProp, 𝓦𝓢.XX, i)
        print_adaptation_info(adaptiveProp, accImpCounter, accUpdtCounter, i)
        adaptiveProp, 𝓦𝓢, yPr, ll = adaptationUpdt!(adaptiveProp, 𝓦𝓢, yPr, i,
                                                     ll, ObsScheme(), ST())
        adaptiveProp = still_adapting(adaptiveProp)
    end
    displayAcceptanceRate(𝔅)
    Time = collect(Iterators.flatten(p.tt[1:skipForSave:end-1] for p in P))
    θchain, accImpCounter/numSteps, accUpdtCounter./numSteps, Paths, Time
end

#NOTE deprecated, will be removed once blocking uses containers in ws
function sample_segments!(iRange, Wnr, WW, WWᵒ, P, y, XXᵒ, ρ)
    for i in iRange
        y = sample_segment!(i, Wnr, WW, WWᵒ, P, y, XXᵒ, ρ)
    end
end

#NOTE deprecated, will be removed once blocking uses containers in ws
function sample_segment!(i, Wnr, WW, WWᵒ, P, y, XXᵒ, ρ)
    sample!(WWᵒ[i], Wnr)
    crank_nicolson!(WWᵒ[i].yy, WW[i].yy, ρ)
    solve!(Euler(), XXᵒ[i], y, WWᵒ[i], P[i])
    XXᵒ[i].yy[end]
end

#NOTE deprecated
"""
    update_target_laws!(𝔅::NoBlocking, θᵒ)

Nothing to do
"""
update_target_laws!(𝔅::NoBlocking, θᵒ) = nothing

#NOTE deprecated
"""
    update_target_laws!(𝔅::BlockingSchedule, θᵒ)

Set new parameter `θᵒ` for the target laws in blocking object `𝔅`
"""
function update_target_laws!(𝔅::BlockingSchedule, θᵒ)
    for block in 𝔅.blocks[𝔅.idx]
        for i in block
            𝔅.P[i] = GuidPropBridge(𝔅.P[i], θᵒ)
        end
    end
end

#NOTE deprecated
"""
    update_proposal_laws!(𝔅::BlockingSchedule, θᵒ)

Set new parameter `θᵒ` for the proposal laws inside blocking object `𝔅`
"""
function update_proposal_laws!(𝔅::BlockingSchedule, θᵒ)
    for block in 𝔅.blocks[𝔅.idx]
        for i in block
            𝔅.Pᵒ[i] = GuidPropBridge(𝔅.Pᵒ[i], θᵒ)
        end
    end
end

#NOTE deprecated
#fetchTargetLaw(𝔅::NoBlocking, P) = P[1].Target

#NOTE deprecated
#fetchTargetLaw(𝔅::BlockingSchedule, P) = 𝔅.P[1].Target

function save_path!(ws, wsXX, bXX) #TODO deprecate bXX
    XX = ws.no_blocking_used ? wsXX : bXX
    skip = ws.skip_for_save
    push!(ws.paths, collect(Iterators.flatten(XX[i].yy[1:skip:end-1]
                                               for i in 1:length(XX))))
end

# remember to remove ws.no_blocking_used

#NOTE deprecated
function update_laws!(Ps, θᵒ, ws)
    for block in ws.blocking.blocks[ws.blidx]
        for i in block
            Ps[i] = GuidPropBridge(Ps[i], θᵒ)
        end
    end
end


text:
@testset "matern" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @testset "MaternKernel" begin
        ν = 2.1
        k = MaternKernel(; ν=ν)
        matern(x, ν) = 2^(1 - ν) / gamma(ν) * (sqrt(2ν) * x)^ν * besselk(ν, sqrt(2ν) * x)
        @test MaternKernel(; nu=ν).ν == [ν]
        @test kappa(k, x) ≈ matern(x, ν)
        @test kappa(k, 0.0) == 1.0
        @test metric(MaternKernel()) == Euclidean()
        @test metric(MaternKernel(; ν=2.0)) == Euclidean()
        @test repr(k) == "Matern Kernel (ν = $(ν), metric = Euclidean(0.0))"

        k2 = MaternKernel(; ν=ν, metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(() -> MaternKernel(; nu=ν))

        test_params(k, ([ν],))
    end
    @testset "Matern32Kernel" begin
        k = Matern32Kernel()
        @test kappa(k, x) ≈ (1 + sqrt(3) * x)exp(-sqrt(3) * x)
        @test k(v1, v2) ≈ (1 + sqrt(3) * norm(v1 - v2))exp(-sqrt(3) * norm(v1 - v2))
        @test kappa(Matern32Kernel(), x) == kappa(k, x)
        @test metric(Matern32Kernel()) == Euclidean()
        @test repr(k) == "Matern 3/2 Kernel (metric = Euclidean(0.0))"

        k2 = Matern32Kernel(; metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(Matern32Kernel)
    end
    @testset "Matern52Kernel" begin
        k = Matern52Kernel()
        @test kappa(k, x) ≈ (1 + sqrt(5) * x + 5 / 3 * x^2)exp(-sqrt(5) * x)
        @test k(v1, v2) ≈
            (
            1 + sqrt(5) * norm(v1 - v2) + 5 / 3 * norm(v1 - v2)^2
        )exp(-sqrt(5) * norm(v1 - v2))
        @test kappa(Matern52Kernel(), x) == kappa(k, x)
        @test metric(Matern52Kernel()) == Euclidean()
        @test repr(k) == "Matern 5/2 Kernel (metric = Euclidean(0.0))"

        k2 = Matern52Kernel(; metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(Matern52Kernel)
    end
    @testset "Coherence Materns" begin
        @test kappa(MaternKernel(; ν=0.5), x) ≈ kappa(ExponentialKernel(), x)
        @test kappa(MaternKernel(; ν=1.5), x) ≈ kappa(Matern32Kernel(), x)
        @test kappa(MaternKernel(; ν=2.5), x) ≈ kappa(Matern52Kernel(), x)
    end
end