text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: translate
using AWS.Compat
using AWS.UUIDs

"""
    create_parallel_data(client_token, name, parallel_data_config)
    create_parallel_data(client_token, name, parallel_data_config, params::Dict{String,<:Any})

Creates a parallel data resource in Amazon Translate by importing an input file from Amazon
S3. Parallel data files contain examples that show how you want segments of text to be
translated. By adding parallel data, you can influence the style, tone, and word choice in
your translation output.

# Arguments
- `client_token`: A unique identifier for the request. This token is automatically
  generated when you use Amazon Translate through an AWS SDK.
- `name`: A custom name for the parallel data resource in Amazon Translate. You must assign
  a name that is unique in the account and region.
- `parallel_data_config`: Specifies the format and S3 location of the parallel data input
  file.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: A custom description for the parallel data resource in Amazon Translate.
- `"EncryptionKey"`:
"""
function create_parallel_data(
    ClientToken, Name, ParallelDataConfig; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "CreateParallelData",
        Dict{String,Any}(
            "ClientToken" => ClientToken,
            "Name" => Name,
            "ParallelDataConfig" => ParallelDataConfig,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_parallel_data(
    ClientToken,
    Name,
    ParallelDataConfig,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "CreateParallelData",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClientToken" => ClientToken,
                    "Name" => Name,
                    "ParallelDataConfig" => ParallelDataConfig,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_parallel_data(name)
    delete_parallel_data(name, params::Dict{String,<:Any})

Deletes a parallel data resource in Amazon Translate.

# Arguments
- `name`: The name of the parallel data resource that is being deleted.

"""
function delete_parallel_data(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "DeleteParallelData",
        Dict{String,Any}("Name" => Name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_parallel_data(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "DeleteParallelData",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Name" => Name), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_terminology(name)
    delete_terminology(name, params::Dict{String,<:Any})

A synchronous action that deletes a custom terminology.

# Arguments
- `name`: The name of the custom terminology being deleted.

"""
function delete_terminology(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "DeleteTerminology",
        Dict{String,Any}("Name" => Name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_terminology(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "DeleteTerminology",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Name" => Name), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_text_translation_job(job_id)
    describe_text_translation_job(job_id, params::Dict{String,<:Any})

Gets the properties associated with an asynchronous batch translation job including name,
ID, status, source and target languages, input/output S3 buckets, and so on.

# Arguments
- `job_id`: The identifier that Amazon Translate generated for the job. The
  StartTextTranslationJob operation returns this identifier in its response.

"""
function describe_text_translation_job(
    JobId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "DescribeTextTranslationJob",
        Dict{String,Any}("JobId" => JobId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_text_translation_job(
    JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "DescribeTextTranslationJob",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("JobId" => JobId), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_parallel_data(name)
    get_parallel_data(name, params::Dict{String,<:Any})

Provides information about a parallel data resource.

# Arguments
- `name`: The name of the parallel data resource that is being retrieved.

"""
function get_parallel_data(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "GetParallelData",
        Dict{String,Any}("Name" => Name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_parallel_data(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "GetParallelData",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Name" => Name), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_terminology(name)
    get_terminology(name, params::Dict{String,<:Any})

Retrieves a custom terminology.

# Arguments
- `name`: The name of the custom terminology being retrieved.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"TerminologyDataFormat"`: The data format of the custom terminology being retrieved. If
  you don't specify this parameter, Amazon Translate returns a file that has the same format
  as the file that was imported to create the terminology.  If you specify this parameter
  when you retrieve a multi-directional terminology resource, you must specify the same
  format as that of the input file that was imported to create it. Otherwise, Amazon
  Translate throws an error.
"""
function get_terminology(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "GetTerminology",
        Dict{String,Any}("Name" => Name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_terminology(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "GetTerminology",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Name" => Name), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    import_terminology(merge_strategy, name, terminology_data)
    import_terminology(merge_strategy, name, terminology_data, params::Dict{String,<:Any})

Creates or updates a custom terminology, depending on whether or not one already exists for
the given terminology name. Importing a terminology with the same name as an existing one
will merge the terminologies based on the chosen merge strategy. Currently, the only
supported merge strategy is OVERWRITE, and so the imported terminology will overwrite an
existing terminology of the same name. If you import a terminology that overwrites an
existing one, the new terminology take up to 10 minutes to fully propagate and be available
for use in a translation due to cache policies with the DataPlane service that performs the
translations.

# Arguments
- `merge_strategy`: The merge strategy of the custom terminology being imported. Currently,
  only the OVERWRITE merge strategy is supported. In this case, the imported terminology will
  overwrite an existing terminology of the same name.
- `name`: The name of the custom terminology being imported.
- `terminology_data`: The terminology data for the custom terminology being imported.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: The description of the custom terminology being imported.
- `"EncryptionKey"`: The encryption key for the custom terminology being imported.
"""
function import_terminology(
    MergeStrategy, Name, TerminologyData; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "ImportTerminology",
        Dict{String,Any}(
            "MergeStrategy" => MergeStrategy,
            "Name" => Name,
            "TerminologyData" => TerminologyData,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function import_terminology(
    MergeStrategy,
    Name,
    TerminologyData,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "ImportTerminology",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "MergeStrategy" => MergeStrategy,
                    "Name" => Name,
                    "TerminologyData" => TerminologyData,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_parallel_data()
    list_parallel_data(params::Dict{String,<:Any})

Provides a list of your parallel data resources in Amazon Translate.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of parallel data resources returned for each request.
- `"NextToken"`: A string that specifies the next page of results to return in a paginated
  response.
"""
function list_parallel_data(; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "ListParallelData"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_parallel_data(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "ListParallelData", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_terminologies()
    list_terminologies(params::Dict{String,<:Any})

Provides a list of custom terminologies associated with your account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of custom terminologies returned per list request.
- `"NextToken"`: If the result of the request to ListTerminologies was truncated, include
  the NextToken to fetch the next group of custom terminologies.
"""
function list_terminologies(; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "ListTerminologies"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_terminologies(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "ListTerminologies", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_text_translation_jobs()
    list_text_translation_jobs(params::Dict{String,<:Any})

Gets a list of the batch translation jobs that you have submitted.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filter"`: The parameters that specify which batch translation jobs to retrieve. Filters
  include job name, job status, and submission time. You can only set one filter at a time.
- `"MaxResults"`: The maximum number of results to return in each page. The default value
  is 100.
- `"NextToken"`: The token to request the next page of results.
"""
function list_text_translation_jobs(; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "ListTextTranslationJobs"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_text_translation_jobs(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "ListTextTranslationJobs",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_text_translation_job(client_token, data_access_role_arn, input_data_config, output_data_config, source_language_code, target_language_codes)
    start_text_translation_job(client_token, data_access_role_arn, input_data_config, output_data_config, source_language_code, target_language_codes, params::Dict{String,<:Any})

Starts an asynchronous batch translation job. Batch translation jobs can be used to
translate large volumes of text across multiple documents at once. For more information,
see async. Batch translation jobs can be described with the DescribeTextTranslationJob
operation, listed with the ListTextTranslationJobs operation, and stopped with the
StopTextTranslationJob operation.  Amazon Translate does not support batch translation of
multiple source languages at once.

# Arguments
- `client_token`: A unique identifier for the request. This token is auto-generated when
  using the Amazon Translate SDK.
- `data_access_role_arn`: The Amazon Resource Name (ARN) of an AWS Identity Access and
  Management (IAM) role that grants Amazon Translate read access to your input data. For more
  information, see identity-and-access-management.
- `input_data_config`: Specifies the format and S3 location of the input documents for the
  translation job.
- `output_data_config`: Specifies the S3 folder to which your job output will be saved.
- `source_language_code`: The language code of the input language. For a list of language
  codes, see what-is-languages. Amazon Translate does not automatically detect a source
  language during batch translation jobs.
- `target_language_codes`: The language code of the output language.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"JobName"`: The name of the batch translation job to be performed.
- `"ParallelDataNames"`: The name of a parallel data resource to add to the translation
  job. This resource consists of examples that show how you want segments of text to be
  translated. When you add parallel data to a translation job, you create an Active Custom
  Translation job.  This parameter accepts only one parallel data resource.  Active Custom
  Translation jobs are priced at a higher rate than other jobs that don't use parallel data.
  For more information, see Amazon Translate pricing.  For a list of available parallel data
  resources, use the ListParallelData operation. For more information, see
  customizing-translations-parallel-data.
- `"Settings"`: Settings to configure your translation output, including the option to mask
  profane words and phrases.
- `"TerminologyNames"`: The name of a custom terminology resource to add to the translation
  job. This resource lists examples source terms and the desired translation for each term.
  This parameter accepts only one custom terminology resource. For a list of available custom
  terminology resources, use the ListTerminologies operation. For more information, see
  how-custom-terminology.
"""
function start_text_translation_job(
    ClientToken,
    DataAccessRoleArn,
    InputDataConfig,
    OutputDataConfig,
    SourceLanguageCode,
    TargetLanguageCodes;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "StartTextTranslationJob",
        Dict{String,Any}(
            "ClientToken" => ClientToken,
            "DataAccessRoleArn" => DataAccessRoleArn,
            "InputDataConfig" => InputDataConfig,
            "OutputDataConfig" => OutputDataConfig,
            "SourceLanguageCode" => SourceLanguageCode,
            "TargetLanguageCodes" => TargetLanguageCodes,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_text_translation_job(
    ClientToken,
    DataAccessRoleArn,
    InputDataConfig,
    OutputDataConfig,
    SourceLanguageCode,
    TargetLanguageCodes,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "StartTextTranslationJob",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClientToken" => ClientToken,
                    "DataAccessRoleArn" => DataAccessRoleArn,
                    "InputDataConfig" => InputDataConfig,
                    "OutputDataConfig" => OutputDataConfig,
                    "SourceLanguageCode" => SourceLanguageCode,
                    "TargetLanguageCodes" => TargetLanguageCodes,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    stop_text_translation_job(job_id)
    stop_text_translation_job(job_id, params::Dict{String,<:Any})

Stops an asynchronous batch translation job that is in progress. If the job's state is
IN_PROGRESS, the job will be marked for termination and put into the STOP_REQUESTED state.
If the job completes before it can be stopped, it is put into the COMPLETED state.
Otherwise, the job is put into the STOPPED state. Asynchronous batch translation jobs are
started with the StartTextTranslationJob operation. You can use the
DescribeTextTranslationJob or ListTextTranslationJobs operations to get a batch translation
job's JobId.

# Arguments
- `job_id`: The job ID of the job to be stopped.

"""
function stop_text_translation_job(JobId; aws_config::AbstractAWSConfig=global_aws_config())
    return translate(
        "StopTextTranslationJob",
        Dict{String,Any}("JobId" => JobId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_text_translation_job(
    JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "StopTextTranslationJob",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("JobId" => JobId), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    translate_text(source_language_code, target_language_code, text)
    translate_text(source_language_code, target_language_code, text, params::Dict{String,<:Any})

Translates input text from the source language to the target language. For a list of
available languages and language codes, see what-is-languages.

# Arguments
- `source_language_code`: The language code for the language of the source text. The
  language must be a language supported by Amazon Translate. For a list of language codes,
  see what-is-languages. To have Amazon Translate determine the source language of your text,
  you can specify auto in the SourceLanguageCode field. If you specify auto, Amazon Translate
  will call Amazon Comprehend to determine the source language.
- `target_language_code`: The language code requested for the language of the target text.
  The language must be a language supported by Amazon Translate.
- `text`: The text to translate. The text string can be a maximum of 5,000 bytes long.
  Depending on your character set, this may be fewer than 5,000 characters.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Settings"`: Settings to configure your translation output, including the option to mask
  profane words and phrases.
- `"TerminologyNames"`: The name of the terminology list file to be used in the
  TranslateText request. You can use 1 terminology list at most in a TranslateText request.
  Terminology lists can contain a maximum of 256 terms.
"""
function translate_text(
    SourceLanguageCode,
    TargetLanguageCode,
    Text;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "TranslateText",
        Dict{String,Any}(
            "SourceLanguageCode" => SourceLanguageCode,
            "TargetLanguageCode" => TargetLanguageCode,
            "Text" => Text,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function translate_text(
    SourceLanguageCode,
    TargetLanguageCode,
    Text,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "TranslateText",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "SourceLanguageCode" => SourceLanguageCode,
                    "TargetLanguageCode" => TargetLanguageCode,
                    "Text" => Text,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_parallel_data(client_token, name, parallel_data_config)
    update_parallel_data(client_token, name, parallel_data_config, params::Dict{String,<:Any})

Updates a previously created parallel data resource by importing a new input file from
Amazon S3.

# Arguments
- `client_token`: A unique identifier for the request. This token is automatically
  generated when you use Amazon Translate through an AWS SDK.
- `name`: The name of the parallel data resource being updated.
- `parallel_data_config`: Specifies the format and S3 location of the parallel data input
  file.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: A custom description for the parallel data resource in Amazon Translate.
"""
function update_parallel_data(
    ClientToken, Name, ParallelDataConfig; aws_config::AbstractAWSConfig=global_aws_config()
)
    return translate(
        "UpdateParallelData",
        Dict{String,Any}(
            "ClientToken" => ClientToken,
            "Name" => Name,
            "ParallelDataConfig" => ParallelDataConfig,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_parallel_data(
    ClientToken,
    Name,
    ParallelDataConfig,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return translate(
        "UpdateParallelData",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClientToken" => ClientToken,
                    "Name" => Name,
                    "ParallelDataConfig" => ParallelDataConfig,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
# Julia wrapper for header: /opt/ffmpeg/include/libavutil/stereo3d.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_stereo3d_alloc,
    av_stereo3d_create_side_data

function av_stereo3d_alloc()
    ccall((:av_stereo3d_alloc,libavutil),Ptr{AVStereo3D},())
end

function av_stereo3d_create_side_data(frame)
    ccall((:av_stereo3d_create_side_data,libavutil),Ptr{AVStereo3D},(Ptr{AVFrame},),frame)
end

ext: jl


text:
"""
    A = JopReghost(spc, zt, zm, dt, dx)

Build a 2D receiver-side reghosting and redatuming operator in the FK domain using
simple phase-shift operations (e.g. Posthumus, 1993).  `spc::JetSpaceSymmetric` is
the domain and range of `A`. `zt` is target depth where the receivers will be
re-dataumed to, and `zm` is the measurement depth where the recording is made.
`dt` is the time sampling interval and `dx` is the receiver sampling interval.

Since this operator is built in the **FK** domain, `spc` has symmetry in the
frequency dimension.  Hence `spc::JotSpaceSymmetric`.  We provide a convenience
method:

    sp = symspace(JopReghost, T, nw, nk)

where `T` is either `Float32` or `Float64`, `nw` are the number of frequency samples,
and `nk` are the number of wavenumber samples.

# Example:
```julia
P = JopPad(JetSpace(Float64,nt,nx), 1:2nt, 1:nx)
F = JopFft(range(P))
G = JopReghost(range(F), 20.0, 20.0, .004, 10.0)
A = P'*F'*G*F*P
data_reghost = A * data_noghost
```
Note that in the above example,
```julia
range(F) == symspace(JopReghost_df!, Float64, size(range(P))...)
```
"""
function JopReghost(spc::JetSSpace, zt, zm, dt, dx)
    nw = size(spc,1)
    nk = size(spc,2)
    nkpos = div(nk,2)+1
    nkneg = nk - nkpos
    wn = pi/dt
    dw = wn/nw
    kn = pi/dx
    dk = 2*kn/nk
    w = collect(range(0.0, stop=dw*(nw-1), length=nw))
    kpos = range(0.0, stop=dk*(nkpos-1), length=nkpos)
    kneg = range(-dk*nkneg, stop=-dk, length=nkneg)
    k = [kpos;kneg]

    @assert Base.length(w) == nw
    @assert Base.length(k) == nk
    @assert Base.length(kpos) == nkpos
    @assert Base.length(kneg) == nkneg

    JopLn(dom = spc, rng = spc, df! = JopReghost_df!, df′! = JopReghost_df′!, s = (zt=zt, zm=zm, dt=dt, dx=dx, k=k, w=w))
end

export JopReghost

JopReghost_df!(d::AbstractArray{Complex{T},2}, m::AbstractArray{Complex{T},2}; zt, zm, dt, dx, k, w, kwargs...) where {T} = reghost(d, m, -one(T), zt, zm, dt, dx, k, w)
JopReghost_df′!(m::AbstractArray{Complex{T},2}, d::AbstractArray{Complex{T},2}; zt, zm, dt, dx, k, w, kwargs...) where {T} = reghost(m, d, one(T), zt, zm, dt, dx, k, w)

function reghost(d::AbstractArray{Complex{T}}, m, sgn, zt, zm, dt, dx, k, w) where {T}
    nw, nk = size(d)
    @assert size(m) == (nw,nk)

    dzs = zt - zm
    dzg = zt + zm
    for iw = 1:nw
        wc = w[iw]^2 / T(1500)^2
        for ik = 1:nk
            kz2 = wc - k[ik]^2
            if kz2 >= 0.0
                kz = sgn*sqrt(kz2)
                d[iw,ik] = (exp(im*kz*dzs) - exp(im*kz*dzg)) * m[iw,ik]  # todo -- obliquity factor
            else
                d[iw,ik] = zero(T)
            end
        end
    end
    d
end

function JopReghost_symspace_map(I, n, _n)
    if I[1] > _n[1]
        return CartesianIndex((2 + n[1] - I[1], I[2]))
    else
        return CartesianIndex(I)
    end
end

function Jets.symspace(JopReghost, ::Type{T}, nw, nk) where {T}
    n = (nw,nk)
    _n = (div(nw,2)+1, nk)
    JetSSpace(Complex{T}, n, _n, I->JopReghost_symspace_map(I, n, _n))
end


text:
const colours = Dict(0=>:black,1=>:red,2=>:green,3=>:blue,4=>:cyan,5=>:magenta,6=>:yellow,7=>:orange)

## HierarchicalDomain

using Plots

Plots.plot(H::HierarchicalDomain;grid=true,kwds...) = plot!(plot(grid=grid),H;kwds...)
Plots.plot!(H::HierarchicalDomain;kwds...) = plot!(current(),H;kwds...)

function Plots.plot!(plt::Plots.Plot,H::HierarchicalDomain;kwds...)
    H1,H2 = partition(H)
    plot!(plt,H1;kwds...)
    plot!(plt,H2;kwds...)
end


text:
using Sockets
import AssetRegistry, JSON
using WebIO
using .WebSockets: is_upgrade, upgrade, writeguarded
using .WebSockets: HTTP

struct WSConnection{T} <: WebIO.AbstractConnection
    sock::T
end

Sockets.send(p::WSConnection, data) = writeguarded(p.sock, JSON.json(data))
Base.isopen(p::WSConnection) = isopen(p.sock)

const bundle_key = AssetRegistry.register(normpath(joinpath(
    WebIO.packagepath, "generic-http-provider", "dist", "generic-http.js"
)))

"""
Serve an asset from the asset registry.
"""
function serve_assets(req)
    if haskey(AssetRegistry.registry, req.target)
        filepath = AssetRegistry.registry[req.target]
        if isfile(filepath)
            return HTTP.Response(
                200,
                [],
                body = read(filepath)
            )
        end
    end
    return HTTP.Response(404)
end

function websocket_handler(ws)
    conn = WSConnection(ws)
    while isopen(ws)
        data, success = WebSockets.readguarded(ws)
        !success && break
        msg = JSON.parse(String(data))
        WebIO.dispatch(conn, msg)
    end
end

struct WebIOServer{S}
    server::S
    serve_task::Task
end

kill!(server::WebIOServer) = put!(server.server.in, HTTP.Servers.KILL)

const singleton_instance = Ref{WebIOServer}()

const routing_callback = Ref{Any}((req)-> missing)

"""
Run the WebIO server.

# Examples
```
using WebSockets, WebIO
app = Ref{Any}(node(:div, "hi"))
function serve_app(req)
    req.target != "/" && return missing
    return sprint() do io
        print(io, \"\"\"
            <!doctype html><html><head>
            <meta charset="UTF-8"></head><body>
        \"\"\")
        show(io, MIME"application/webio"(), app[])
        print(io, "</body></html>")
    end
end
server = WebIO.WebIOServer(serve_app, logger = stdout, verbose = true)
```
"""
function WebIOServer(
        default_response::Function = (req)-> missing;
        baseurl::String = "127.0.0.1", http_port::Int = 8081,
        verbose = false, singleton = true,
        websocket_route = "/webio_websocket/",
        server_kw_args...
    )
    # TODO test if actually still running, otherwise restart even if singleton
    if !singleton || !isassigned(singleton_instance)
        function handler(req)
            response = default_response(req)
            response !== missing && return response
            return serve_assets(req)
        end
        function wshandler(req, sock)
            req.target == websocket_route && websocket_handler(sock)
        end
        server = WebSockets.ServerWS(handler, wshandler; server_kw_args...)
        server_task = @async WebSockets.serve(server, baseurl, http_port, verbose)
        singleton_instance[] = WebIOServer(server, server_task)
    end
    return singleton_instance[]
end

const webio_server_config = Ref{typeof((url = "", bundle_url = "", http_port = 0, ws_url = ""))}()

"""
Fetches the global configuration for our http + websocket server from environment
variables. It will memoise the result, so after a first call, any update to
the environment will get ignored.
"""
function global_server_config()
    if !isassigned(webio_server_config)

        setbaseurl!(get(ENV, "JULIA_WEBIO_BASEURL", ""))

        url = get(ENV, "WEBIO_SERVER_HOST_URL", "127.0.0.1")
        http_port = parse(Int, get(ENV, "WEBIO_HTTP_PORT", "8081"))
        ws_default = string("ws://", url, ":", http_port, "/webio_websocket/")
        ws_url = get(ENV, "WEBIO_WEBSOCKT_URL", ws_default)
        # make it possible, to e.g. host the bundle online
        bundle_url = get(ENV, "WEBIO_BUNDLE_URL", string(WebIO.baseurl[], bundle_key))
        webio_server_config[] = (
            url = url, bundle_url = bundle_url,
            http_port = http_port, ws_url = ws_url
        )
    end
    webio_server_config[]
end

"""
Show a WebIO application.

This method ensures that the requisite asset server and WebSocket endpoints
are running and prints the required HTML scripts and WebIO mounting code.

# Examples
```
function Base.display(d::MyWebDisplay, m::MIME"application/webio", app)
    println(d.io, "outer html")
    show(io, m, app)
    println(d.io, "close outer html")
end
```
"""
function Base.show(io::IO, m::WEBIO_APPLICATION_MIME, app::Application)
    c = global_server_config()
    WebIOServer(routing_callback[], baseurl = c.url, http_port = c.http_port)
    println(io, "<script>window._webIOWebSocketURL = $(repr(c.ws_url));</script>")
    println(io, "<script src=$(repr(c.bundle_url))></script>")
    show(io, "text/html", app)
    return
end


text:
struct RootNode <: Node
    children::Tuple{Any}
    function RootNode(node)
        new((node,))
    end
end

function _traversal(node, visitfn::Function, visitfnargs::Union{Any, Nothing}, visited)
    if haskey(visited, node)
        return visited[node]
    else
        new_children = [_traversal(child, visitfn, visitfnargs, visited) for child in node.children]
        result = visitfnargs !== nothing ?
            visitfn(visited, node, visitfnargs, new_children...) :
                visitfn(visited, node, new_children...)
        visited[node] = result
        return result
    end
end

function _posttraversal(node, visitfn::Function, visitfnargs::Union{Any, Nothing}, visited)
    if haskey(visited, node)
        return visited[node]
    else
        result = visitfnargs !== nothing ?
            visitfn(node, visitfnargs, (x -> _posttraversal(x, visitfn, visitfnargs, visited))) :
                visitfn(node, (x -> _posttraversal(x, visitfn, visitfnargs, visited)))
        visited[node] = result
        return result
    end
end

"""
    traversal(node, pretraversalfn::Function, visitfn::Function,
        pretraversalfnargs::Union{Any, Nothing}, visitfnargs::Union{Any, Nothing},
        visited::Union{Dict, Nothing}, posttraversal=false)

Traverses an AST, calling visitfn on the nodes of the tree.

# Arguments
- `node`: the root of the tree to traverse. Assumes that each node has children
attribute.
- `pretraversalfn::Function` A function to be executed before the traversal
- `visitfn::Function` The function that is executed on each visit. If
posttraversal is false, and visitfnargs nothing, the function must have signature
`visitfn(visited, node, new_children...)`. If visitfnargs is not nothing, the
signature is `visitfn(visited, node, visitfnargs, new_children...). If
posttraversal is true, the signature must be `visitfn(node, fn)`,
or `visitfn(node, visitfn, fn)`, where `fn` is a function that calls traversal
on a node.
- `pretraversalfnargs::Union{Any, Nothing}` Arguments for the pretraversalfn.
- `visitfnargs::Union{Any, Nothing}` Arguments for traversalfn.
- `visited::Union{Dict, Nothing}` If not nothing, a dictionary with the results
of calling traversalfn on different nothing. If nothing, this dictionary is
created in traversal.
- `posttraversal=false` Determines the way traversal traverses each node. If
`posttraversal=false`, the traversal first processes the children of a node,
and then the parent node (example: [`togem(node::Node)`](@ref)). If
`posttraversal=true`, the visitfn first processes the parent and then calls
visitfn on its children if needed (example:
[`differentiate(A::AbstractTensor{0}, x::VariableTensor{0})`](@ref)).
"""
function traversal(node, pretraversalfn::Function, visitfn::Function,
        pretraversalfnargs::Union{Any, Nothing}, visitfnargs::Union{Any, Nothing},
        visited::Union{Dict, Nothing}, posttraversal::Bool)
    if visited === nothing
        visited = Dict{Any, Any}()
    end
    root = RootNode(node)
    root = pretraversalfnargs !== nothing ? pretraversalfn(root, pretraversalfnargs...) : pretraversalfn(root)
    root = posttraversal ? _posttraversal(
        root, visitfn, visitfnargs, visited) : _traversal(
                root, visitfn, visitfnargs, visited)
    return root.children[1]
end

function traversal(node, pretraversalfn::Function, visitfn::Function,
        pretraversalfnargs::Union{Any, Nothing}, visitfnargs::Union{Any, Nothing},
        visited::Union{Dict, Nothing})
    return traversal(node, pretraversalfn, visitfn, pretraversalfnargs,
            visitfnargs, visited, false)
end

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
"""
Author: Michael Palumbo
Created: December 2019
Contact: mlp95@psu.edu
Based on code by Alex Wise (aw@psu.edu)
Refactored and optimized by Eric Ford
"""

"""
Module for estimating the radial velocity based on the CCF
"""
module RVFromCCF

import Statistics: mean
import PDMats: PDiagMat
#import Polynomials
import ..RvSpectMLBase: searchsortednearest, speed_of_light_mps

include("common.jl")
export AbstractMeasureRvFromCCF
export measure_rv_from_ccf, measure_rvs_from_ccf

include("find_min_ccf.jl")
export MeasureRvFromMinCCF

include("calc_ccf_centroid.jl")
export MeasureRvFromCCFCentroid

include("fit_quadratic_to_ccf.jl")
export MeasureRvFromCCFQuadratic

include("fit_gaussian_to_ccf.jl")
#import .FitGaussianToCCF
#import .FitGaussianToCCF: MeasureRvFromCCFGaussian
export MeasureRvFromCCFGaussian

import ..EchelleCCFs: calc_ccf_template, calc_normalized_ccfs
include("ccf_sample_covar.jl")
export calc_ccf_sample_covar

include("ccf_covar_model.jl")
export calc_ccf_sample_covar_reduced_rank, calc_ccf_sample_covar_and_near_diag_model, est_covar_for_obs

include("fit_template_plus_deriv_to_ccf.jl")
#import .FitTemplateToCCF
#import .FitTemplateToCCF: MeasureRvFromCCFTemplate
export MeasureRvFromCCFTemplate

include("fit_template_plus_deriv_to_ccf_non_diag_covar.jl")
#import .FitTemplateToCCF
#import .FitTemplateToCCF: MeasureRvFromCCFTemplateNonDiagCovar
export MeasureRvFromCCFTemplateNonDiagCovar

end # module

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
module EvoTrees

export init_evotree, grow_evotree!, grow_tree, fit_evotree, predict,
    fit_evotree_gpu, predict_gpu,
    EvoTreeRegressor, EvoTreeCount, EvoTreeClassifier, EvoTreeGaussian,
    EvoTreeRModels, importance, Random

using Base.Threads: @threads
using Statistics
using StatsBase: sample, quantile
using Random
using Distributions
using StaticArrays
using CategoricalArrays
using CUDA
using BSON: @save, @load
using NetworkLayout
using RecipesBase
import MLJModelInterface
import MLJModelInterface: fit, update, predict

include("models.jl")
include("structs.jl")
include("loss.jl")
include("eval.jl")
include("predict.jl")
include("find_split.jl")
include("fit.jl")
include("importance.jl")
include("plot.jl")
include("MLJ.jl")

include("gpu/structs_gpu.jl")
include("gpu/loss_gpu.jl")
include("gpu/eval_gpu.jl")
include("gpu/predict_gpu.jl")
include("gpu/find_split_gpu.jl")
include("gpu/fit_gpu.jl")

function save(model::GBTree, path)
    @save path model
end

function load(path)
    @load path model
    return model
end

end # module


text:
"""
$(TYPEDEF)

PlutoVistaPlot is a wrapper struct around different possible backends.
As default, for 1D plots, [`PlutoPlotlyPlot`](@ref) and for 2D and 3D plots,
[`PlutoVTKPlot`](@ref) are chosen. The backend is automatically chosen
when the first plot is invoked.

"""
mutable struct PlutoVistaPlot
    backend::Union{PlutoPlotlyPlot,PlutoCanvasPlot,PlutoVTKPlot,Nothing}
    resolution
    args
    PlutoVistaPlot(::Nothing)=new(nothing, (0,0), nothing)
end

"""
$(SIGNATURES)

Constructor for PlutoVistaPlot.

Relevant keyword arguments:
- `resolution`: Screen resolution in pixels
- `dim`: Space dimension of subsequent plots

"""
function PlutoVistaPlot(;resolution=(300,300),kwargs...)
    p=PlutoVistaPlot(nothing)
    p.resolution=resolution
    p.args=kwargs
    p
end

"""
    backend!(p::PlutoVistaPlot; datadim=1, backend=:default, clear=false)

Overwrite backend if clear is true.
"""
function backend!(p::PlutoVistaPlot; datadim=1, backend=:default, clear=false, kwargs...)
    if backend==:default
        if datadim == 1
            backend=:plotly
        else
            backend=:vtk
        end
    end

    if backend == :plotly
        if isnothing(p.backend)  || typeof(p.backend)!=PlutoPlotlyPlot || clear
            p.backend=PlutoPlotlyPlot(;resolution=p.resolution, merge(p.args,kwargs)...)
        end
    elseif backend == :vtk
        if isnothing(p.backend)  || typeof(p.backend)!=PlutoVTKPlot || clear
            p.backend=PlutoVTKPlot(;resolution=p.resolution,merge(p.args,kwargs)...)
        end
    else
        error("No valid backend: $(backend). Valid backends: :default, :vtk, :plotly")
    end
    p.backend
end

Base.show(io::IO, mime::MIME"text/html", p::Nothing)=nothing
Base.show(io::IO, mime::MIME"text/html", p::PlutoVistaPlot)=Base.show(io,mime,p.backend)

####################################################################################
"""
$(SIGNATURES)

Plot `y` over `y`.

Arguments:
- `x` vector of x coordinates
- `y` vector of y coordinates (function values)

Keyword  arguments:
- `limits=(1,-1)`: function limits
- `xlimits=(1,-1)` : x axis limits
- `xlabel=""`: x axis label
- `ylabel=""`: y axis label
- `title="" : plot title
- `xscale=:linear`: linear or log scale for x axis
- `yscale=:linear`: linear or log scale for y axis
- `legend=:none`: legend placement (:lt,:rt ...)
- `clear=false`: clear plot contents
- `label=""`: label of particular plot
- `color=:auto`
- `titlefontsize=12`
- `axisfontsize=10`
- `tickfontsize=10`
- `legendfontsize=10`
- `linewidth=2`
- `linestyle=:solid`:
- `markersize=6`
- `markercount=10`
- `markertype=:none`, possible values: `:dtriangle`, `:utriangle`, `:rtriangle`, `:ltriangle`, `:circle`, `:square`, `:cross`, `:+ `, `:xcross`, `:x `,`:diamond`, `:star5`, `:pentagon`, `:hexagon

"""
plot(x,y; kwargs...)=plot!(PlutoVistaPlot(;kwargs...),x,y;kwargs...)

"""
$(SIGNATURES)

Add additional x-y plot to `p`
"""
plot!(p::PlutoVistaPlot,x,y; backend=:plotly, clear=false, kwargs...) = plot!(backend!(p; datadim=1,backend=backend,clear=false, kwargs...),
                                                                              x,y;
                                                                              clear=clear,kwargs...)

"""
$(SIGNATURES)

Create empty 1D plot.
"""
function plot(;datadim=1,backend=:default,kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=datadim,backend=backend, kwargs...)
end

##########################################################################
"""
$(SIGNATURES)

Filled colored tricontour with isolines. By default, a vtk.js based backend is used.
Plots piecewise linear function on  triangular grid.

Arguments:
- `pts`:  `2 x n_points` array of point coordinates
- `tris`: `3 x n_tris` array of point indices describing triangles
- `f`: `n`-vector of function values

Keyword arguments:
- `title=""`
- `clear=false`
- `xlabel="x"`: x axis label
- `ylabel="y"`: y axis label
- `colormap=:viridis`
- `levels=0`: either number of interior isolevels, or vector of isolevel values
- `colorbarticks=:default` : colorbar ticks. Default: levels
- `limits=:auto`: function limits
- `aspect=1`: xy aspect ratio (plotly backend)

"""
tricontour(pts,tris,f; kwargs...)=tricontour!(PlutoVistaPlot(;kwargs...),pts,tris,f; kwargs...)

"""
$(SIGNATURES)

Add tricontour to existing plot
"""
tricontour!(p::PlutoVistaPlot,pts,tris,f;backend=:vtk, kwargs...)=tricontour!(backend!(p; datadim=2,backend=backend, kwargs...),
                                                                              pts,tris,f; kwargs...)

"""
$(SIGNATURES)

Create empty tricontour plot
"""
function tricontour(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=2,kwargs...)
end

##########################################################################
"""
$(SIGNATURES)

Create triangular mesh from x and y coordinates, and call tricontour.
"""
contour(X,Y,f; kwargs...)=contour!(PlutoVistaPlot(;kwargs...),X,Y,f; kwargs...)

"""
$(SIGNATURES)

Add contour to existing plot
"""
contour!(p::PlutoVistaPlot,X,Y,f; backend=:vtk, kwargs...)=contour!(backend!(p;datadim=2,backend=backend, kwargs...),
                                                                    X,Y,f; kwargs...)
"""
$(SIGNATURES)

Create empty contour plot
"""
function contour(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=2,kwargs...)
end

##########################################################################
"""
$(SIGNATURES)
Quiver plot, using vtk.js backend.

Arguments:

- `pts`: `2 x n_pts` array of points
- `qvec`: `2 x n_pts` array of vector values
"""
quiver2d(pts,qvec; kwargs...)=quiver2d!(PlutoVistaPlot(;kwargs...),pts,qvec; kwargs...)

"""
$(SIGNATURES)

Add quiver to existing plot
"""
quiver2d!(p::PlutoVistaPlot,pts,qvec;backend=:vtk, kwargs...)=quiver2d!(backend!(p; datadim=2,backend=backend, kwargs...),
                                                                              pts,qvec; kwargs...)
"""
$(SIGNATURES)

Create empty quiver plot
"""
function quiver2d(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=3,kwargs...)
end

##########################################################################
"""
$(SIGNATURES)

Plot triangular mesh, showing triangle boundaries and outer boundaries, using vtk.js backend.

Arguments:
- `pts`:  `2 x n_points` array of point coordinates
- `tris`: `3 x n_tris` array of point indices describing triangles

Keyword arguments:

- `markers=nothing`: Optional `n_tris` integer vector of triangle markers
- `colormap`: optional colormap for triangle markers
- `edges`: `2 x n_edges` optional  array of point indices describing edges
- `edgemarkers=nothing`: optional `n_edges` vector of integer edge markers
- `edgecolormap=nothing`: optional colormap for edge markers

"""
trimesh(pts,tris; kwargs...)=trimesh!(PlutoVistaPlot(;kwargs...),pts,tris; kwargs...)

"""
$(SIGNATURES)

Add trimesh to plot.
"""
trimesh!(p::PlutoVistaPlot,pts,tris; backend=:vtk, kwargs...)=trimesh!(backend!(p;datadim=2,backend=backend, kwargs...),
                                                                        pts,tris; kwargs...)

"""
$(SIGNATURES)

Create empty trimesh plot.
"""
function trimesh(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=2,kwargs...)
end

##########################################################################
"""
$(SIGNATURES)

Plot piecewise linear function on tetrahedral grid using vtk backend.
The plot consists of three parts:

- Transparent isosurfaces of function values according to `levels`
- Transparent plane cuts with function "hetmapes" according to `xplanes`,`yplanes`,`zplanes`
- Transparent domain outline

Arguments:
- `pts`:  `3 x n_points` array of point coordinates
- `tets`: `4 x n_tets` array of point indices describing tetrahedra
- `f`: `n_points`-vector of function values

Keyword arguments:
- `title=""`
- `clear=false`
- `colormap=:viridis`
- `faces=nothing`: optional `3 x n_faces` array of boundary faces
- `facemarkers=nothing`: optional `n_faces` integer vector of face markers
- `facecolormap=nothing`: optional colormap of facemarkers
- `levels=0`: either number of interior isolevels, or vector of isolevel values
- `limits=:auto`: function limits
- `xplanes`: either number or array of x coordinate values of x-orthogonal plane sections
- `yplanes`: either number or array of y coordinate values of y-orthogonal plane sections
- `zplanes`: either number or array of z coordinate values of z-orthogonal plane sections
- `levelalpha=0.25`: alpha value for isosurfaces and plane cuts
- `outlinealpha=0.1`: alpha value for outline. Outliene is for value 0.0

"""
tetcontour(pts,tets,f; kwargs...)=tetcontour!(PlutoVistaPlot(;kwargs...),pts,tets,f; kwargs...)

"""
$(SIGNATURES)

Add tetcontour to plot
"""
tetcontour!(p::PlutoVistaPlot,pts,tets,f;backend=:vtk, kwargs...)=tetcontour!(backend!(p; datadim=2,backend=backend, kwargs...),
                                                                              pts,tets,f; kwargs...)
"""
$(SIGNATURES)

Create empty tetcontour plot.
"""
function tetcontour(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=3,kwargs...)
end

##########################################################################
"""
$(SIGNATURES)

Plot tetrahedral mesh, showing tetrahedron boundaries and outer boundaries, using vtk.js backend.
The plot consists of two parts:

- A subset of tetrahedra visible after cutting of all tets with are on the positive side of the respective x, y, z planes
- A transpatemt outline of the boundary.

Arguments:
- `pts`:  `3 x n_points` array of point coordinates
- `tets`: `4 x n_tris` array of point indices describing tetrahedra

Keyword arguments:

- `markers=nothing`: Optional `n_tets` integer vector of tetrahedron markers
- `colormap`: optional colormap for tetrahedron markers
- `faces=nothing`: optional `3 x n_faces` array of boundary faces
- `facemarkers=nothing`: optional `n_faces` integer vector of face markers
- `facecolormap=nothing`: optional colormap of facemarkers
- `xplanes`: array of x coordinate values for cut-off in x-direction
- `yplanes`: array of y coordinate values for cut-off in y-direction
- `zplanes`: array of z coordinate values for cut-off in z-direction
- `outlinealpha=0.1`: alpha value for outline. Outliene is for value 0.0
"""
tetmesh(pts,tets; kwargs...)=tetmesh!(PlutoVistaPlot(;kwargs...),pts,tets; kwargs...)

"""
$(SIGNATURES)

Add tetmesh to plot
"""
tetmesh!(p::PlutoVistaPlot,pts,tets; backend=:vtk, kwargs...)=tetmesh!(backend!(p;datadim=3,backend=backend, kwargs...),
                                                                        pts,tets; kwargs...)

"""
$(SIGNATURES)

Create empty tetmesh plot.
"""
function tetmesh(;kwargs...)
    p=PlutoVistaPlot(;kwargs...)
    backend!(p;datadim=2,kwargs...)
end

#########################################################
"""
$(SIGNATURES)

Enter new command named `cmd`.

The idea is to pass one single Dict with all plot data to javascript
using `publish_to_js`.

For this purpose we need a some "language" hidden in
behind the dict.

Plot   elements  are   described  by   commands  executed   one  after
another. For  this purpose,  we use integers  converted to  strings as
dict keys and pass the command name as the corresponding entry.

The entry `cmdcount` keeps track of  the number of commands. So we can
parse all  commands from 1  to jsdict[:cmdcount] in javascript  in the
same sequence as they have been entered.

Parameters are named and entered into the dictionary with the
command number as prefix.

E.g. for a polyline as command number 5, we create the entries

```
"5" => "polyline"
"5x" => Vector of x coordinates in canvas coordinate system
"5y" => Vector of y coordinates in canvas coordinate system
```
"""
function command!(p::T,cmd) where {T <: AbstractPlutoVistaBackend}
    p.jsdict["cmdcount"]=p.jsdict["cmdcount"]+1
    pfx=string(p.jsdict["cmdcount"])
    p.jsdict[pfx]=cmd
    p
end

"""
$(SIGNATURES)

After [`command!`](@ref), create a parameter entry
"""
function parameter!(p::T,name,value) where {T <: AbstractPlutoVistaBackend}
    pfx=string(p.jsdict["cmdcount"])
    p.jsdict[pfx*name]=value
    p
end

"""
$(SIGNATURES)

Create triangulation data from X and Y coordinate vectors
"""
function triang(X,Y)
    nx=length(X)
    ny=length(Y)
    num_pts=nx*ny
    num_tris=2*(nx-1)*(ny-1)
    pts=zeros(Float32,2,num_pts)
    tris=zeros(Int32,3,num_tris)

    ipoint=0
    for iy=1:ny
        for ix=1:nx
            ipoint=ipoint+1
            pts[1,ipoint]=X[ix]
            pts[2,ipoint]=Y[iy]
        end
    end

    @assert(ipoint==num_pts)

    itri=0
    for iy=1:ny-1
        for ix=1:nx-1
	    ip=ix+(iy-1)*nx
	    p00 = ip
	    p10 = ip+1
	    p01 = ip  +nx
	    p11 = ip+1+nx

            itri=itri+1
            tris[1,itri]=p00
            tris[2,itri]=p10
            tris[3,itri]=p11

            itri=itri+1
            tris[1,itri]=p11
            tris[2,itri]=p01
            tris[3,itri]=p00
        end
    end
    @assert(itri==num_tris)
    pts, tris
end


text:
export Database, execute, commit
mutable struct Database
    conn::PyObject
    c::PyObject
    sqlite3::PyObject
    commit_after_execute::Bool
    filename::String
end

"""
    Database(filename::Union{Missing, String} = missing;
        commit_after_execute::Bool = true)

Creates a database from `filename`. If `filename` is not provided, an in-memory database is created.
If `commit_after_execute` is false, no commit operation is performed after each [`execute`](@ref).

- do block syntax:
```julia
Database() do db
    execute(db, "create table mytable (a real, b real)")
end
```
The database is automatically closed after execution. Therefore, if execute is a query operation,
users need to store the results in a global variable.

- Query meta information
```julia
keys(db) # all tables
keys(db, "mytable") # all column names in `db.mytable`
```
"""
function Database(filename::Union{Missing, String} = missing;
    commit_after_execute::Bool = true)
    filename = coalesce(filename, ":memory:")
    sqlite3 = pyimport("sqlite3")
    db = Database(PyNULL(), PyNULL(), sqlite3, commit_after_execute, filename)
    db
end

function connect(db::Database)
    sqlite3 = pyimport("sqlite3")
    db.conn = sqlite3.connect(db.filename)
    db.c = db.conn.cursor()
end

function Database(f::Function, args...; kwargs...)
    db = Database(args...;kwargs...)
    out = f(db)
    close(db)
    out
end

"""
    execute(db::Database, sql::String, args...)

Executes the SQL statement `sql` in `db`. Users can also use the do block syntax.
```julia
execute(db) do
    "create table mytable (a real, b real)"
end
```

`execute` can also be used to insert a batch of records
```julia
t1 = rand(10)
t2 = rand(10)
param = collect(zip(t1, t2))
execute(db, "INSERT TO mytable VALUES (?,?)", param)
```
or
```julia
execute(db, "INSERT TO mytable VALUES (?,?)", t1, t2)
```
"""
function execute(db::Database, sql::String, args...)
    connect(db)
    if length(args)>=1
        if length(args)>1
            param = collect(zip(args...))
        else
            param = args[1]
        end
        db.c.executemany(sql, param)
    else
        db.c.execute(sql)
    end
    if db.commit_after_execute
        commit(db)
    end
    db.c
end

function execute(f::Function, db::Database; kwargs...)
    sql = f()
    execute(db, sql; kwargs...)
end

"""
    commit(db::Database)

Commits changes to `db`.
"""
function commit(db::Database)
    db.conn.commit()
end

function Base.:close(db::Database)
    if !(db.conn==PyNULL())
        commit(db)
        db.c.close()
        db.conn.close()
        db.c = PyNULL()
        db.conn = PyNULL()
    end
    nothing
end

function Base.:keys(db::Database)
    ret = execute(db, "select name from sqlite_master")|>collect
    close(db)
    tables = String[]
    for k = 1:length(ret)
        push!(tables, ret[k][1])
    end
    tables
end

function Base.:keys(db::Database, table::String)
    execute(db, "select * from $table limit 1")
    ret = db.c.description
    close(db)
    columns = String[]
    for r in ret
        push!(columns, r[1])
    end
    columns
end

function Base.:push!(db::Database, table::String, nts::NamedTuple...)
    v1 = []
    v2 = []
    for nt in nts
        cols = propertynames(nt)
        push!(v1, join(["\"$c\"" for c in cols], ","))
        push!(v2, join([isnothing(nt[i]) ? "NULL" : "\"$(nt[i])\"" for i = 1:length(nt)], ","))
    end
    v1 = join(v1, ",")
    v2 = join(v2, ",")
   execute(db,
"""
    INSERT OR REPLACE INTO $table ($v1) VALUES ($v2)
"""
   )
   close(db)
end

function Base.:delete!(db::Database, table::String)
    execute(db, "drop table $table")
    close(db)
end

function Base.:getindex(db::Database, table::String)
    c = execute(db, "select * from $table")
    out = collect(c)
    close(db)
    out
end


text:
module BitwiseKoans

function not(x)
end

function and(x, y)
end

function or(x, y)
end

function xor_koan(x, y)
end

function right_shift(x, y)
end

function left_shift(x, y)
end

end


text:
module FromFileTestPack

export MyInt, add_same_rand1, add_same_rand2, int_square, int_cube, int_unwrap, int_square_unwrap

using FromFile
@from "types.jl" import MyInt
@from "add_same_rand/entrypoint.jl" import SameRand: add_same_rand1
@from "add_same_rand/another_entrypoint.jl" import add_same_rand2
@from "myint_ops/ops.jl" import Square.int_square, Cube.int_cube
@from "myint_ops/unwrap.jl" import int_unwrap, int_square_unwrap
end


text:
# Basic control flow
#
# # Fibonnacci numbers
#
# The Fibonnacci sequence is: 1, 1, 2, 3, 5, 8, 13..., where each subsequent
# number is the sum of the prior two numbers. Create a function that computes
# the N-th Fibonnacci number.
#   Options:
#       * Iterative: use a for loop and keep track of prior two values to add them
#           What happens if you ask for the 95th number? The 96th? 97th?
#       * Recursive: work backwards and add the prior two numbers together
#           How long does it take to get the 40th number? The 41st? 42nd?
#       * Matrix magic

# # Monte-carlo pi
#
# Create a function that estimates pi by repeatedly throwing a dart at a square
# board. It's easy to test if the particular dart landed inside the circle —
# just ask if x^2 + y^2 < 1! Then comparing the number inside vs. the total
# gives you the ratio of the circle's area to the square's area.

using Images
load(download("http://corysimon.github.io/images/julia/myplot.png"))

# Plot the error as the number of iterations increases

# # Packages
#
# Add the Primes package. What does it do? How do you compute the 100th prime?
# A very large prime number?
#
#

max_stars_repo_path: Courses/Introduction to Julia/Day One Exercises.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
function compute!(I::IntegralHelper, entry::String)

    if entry == "S"
        compute_S!(I)
    elseif entry == "T"
        compute_T!(I)
    elseif entry == "V"
        compute_V!(I)
    elseif entry == "ERI"
        compute_ERI!(I)
    else
        throw(FermiException("Invalid key for IntegralHelper: $(entry)."))
    end
end

function compute_S!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["S"] = FermiMDArray(ao_1e(bs, "overlap", T))
end

function compute_T!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["T"] = FermiMDArray(ao_1e(bs, "kinetic", T))
end

function compute_V!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["V"] = FermiMDArray(ao_1e(bs, "nuclear", T))
end

function compute_ERI!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractDFERI}

    bs = I.orbitals.basisset
    auxbs = I.eri_type.basisset
    J = FermiMDArray(ao_2e2c(auxbs, T))
    Pqp = FermiMDArray(ao_2e3c(bs, auxbs, T))
    Jh = similar(J)
    Jh .= Array(real(J^(-1/2)))
    @tensor b[Q,p,q] := Pqp[p,q,P]*Jh[P,Q]
    I.cache["ERI"] = b
end

function compute_ERI!(I::IntegralHelper{T, Chonky, AtomicOrbitals}) where T<:AbstractFloat
    bs = I.orbitals.basisset
    I.cache["ERI"] = FermiMDArray(ao_2e4c(bs, T))
end

function compute_ERI!(I::IntegralHelper{T, SparseERI, AtomicOrbitals}) where T<:AbstractFloat
    bs = I.orbitals.basisset
    I.cache["ERI"] = Fermi.FermiSparse(sparse_ao_2e4c(bs, T)...)
end

function index2(i::Signed, j::Signed)::Signed
    if i < j
        return (j * (j + 1)) >> 1 + i
    else
        return (i * (i + 1)) >> 1 + j
    end
end

# Produces all unique indices ijkl for the two-electron integral
# Note, these indexes start from zero, because we use it in a ccall
function find_indices(nbf::Signed)

    out = NTuple{4,Int16}[]
    N = Int16(nbf - 1)
    ZERO = zero(Int16)

    for i = ZERO:N
        for j = i:N
            for k = ZERO:N
                for l = k:N
                    if index2(i,j) < index2(k,l)
                        continue
                    end
                    push!(out, (i,j,k,l))
                end
            end
        end
    end

    return out
end

include("AtomicIntegrals/OneElectron.jl")
include("AtomicIntegrals/TwoElectronTwoCenter.jl")
include("AtomicIntegrals/TwoElectronThreeCenter.jl")
include("AtomicIntegrals/TwoElectronFourCenter.jl")


text:
# @testset "Similarity depth cache" begin
#     db = DB()

#     function similarity_depth(s1::Synset, s2::Synset, max_depth=5)
#         if s1 == s2
#             return 0
#         end
#         explored = Set{Synset}([s1])
#         active_set = Vector{Synset}([s1])
#         frontier = Vector{Synset}()
#         for depth in 1:max_depth
#             for synset in active_set
#                 for related in SemanticSimilarity.similar_to(db, synset)
#                     if related == s2
#                         return depth
#                     end
#                     if !(related in explored)
#                         push!(frontier, related)
#                         push!(explored, related)
#                     end
#                 end
#             end
#             active_set, frontier = frontier, active_set
#             if isempty(active_set)
#                 break
#             end
#             empty!(frontier)
#         end
#         return nothing
#     end

#     @showprogress "Similarity depths" for groups_for_word in values(SemanticSimilarity.CACHE.similarity_groups)
#         for groups in groups_for_word
#             synset = first(first(groups))
#             for i in 1:length(groups)
#                 for element in groups[i]
#                     @test similarity_depth(synset, element) == i - 1
#                 end
#             end
#         end
#     end
# end


text:
function converge(cur_T_bar::AbstractFloat, tmp_dict::Dict)
  cur_reactor = Reactor(cur_T_bar, tmp_dict)
  cur_reactor.T_bar <= 0 && return []

  cur_reactor.sigma_v = calc_sigma_v(cur_reactor)
  isnan(cur_reactor.sigma_v) && return []

  converge(cur_reactor)
end

function converge(cur_reactor::AbstractReactor; no_pts=no_pts_eta_CD)
  cur_func = function (cur_eta_CD::Number)
    cur_eta_CD = real(cur_eta_CD)

    tmp_reactor = deepcopy(cur_reactor)

    tmp_reactor.eta_CD = cur_eta_CD

    solve!(tmp_reactor)

    tmp_reactor.is_good || return Inf

    tmp_eta_CD = calc_eta_CD(tmp_reactor)

    isnan(tmp_eta_CD) && return Inf

    cur_equation = cur_eta_CD - tmp_eta_CD

    cur_equation
  end

  cur_root_list = find_roots(
    cur_func, min_eta_CD, max_eta_CD,
    no_pts = no_pts, abstol = 1e-4
  )

  cur_root_list
end


text:
module Output

export PlanktonOutputWriter
export write_output!
export humanize_filesize

using LinearAlgebra, Statistics, JLD2, Serialization, Printf

using PlanktonIndividuals.Grids
using PlanktonIndividuals.Fields
using PlanktonIndividuals.Diagnostics
using PlanktonIndividuals.Model

using PlanktonIndividuals: AbstractMode, CarbonMode, QuotaMode, MacroMolecularMode

import Base: show

include("output_writers.jl")
include("write_outputs.jl")

end


text:
# homogeneous velocity and density model
nz = 101; nx = 201;
vel = 3000 * ones(nz, nx);  # m/s
vel[51:end,:] .= 3500;  # m/s
rho = 2000 * ones(nz, nx);  # kg/m^3

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 1.0;  # use second as unit

# source location
isx1 = floor(Int64, nx/2); isz1 = 2;
isx2 = collect(5:40:nx-5); ns=length(isx2); isz2 = 2*ones(ns);

test_type = [(isz1, isx1, true , Float32, 1.0e-5 ),
             (isz1, isx1, true , Float64, 1.0e-10),
             (isz1, isx1, false, Float32, 1.0e-5 ),
             (isz1, isx1, false, Float64, 1.0e-10),
             (isz2, isx2, true , Float32, 1.0e-5 ),
             (isz2, isx2, true , Float64, 1.0e-10),
             (isz2, isx2, false, Float32, 1.0e-5 ),
             (isz2, isx2, false, Float64, 1.0e-10),
            ]
@testset "forward modelling and pressure field reconstruction (forward and backward)" begin

    for (isz, isx, free_surface, data_format, tol) in test_type

        # generate finite difference stencil
        params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                 data_format=data_format, fd_flag="taylor", order=3, npml=20, apml=900.);

        # create source term
        ns = length(isz)
        if ns == 1
           source_term = Source(isz, isx, params; ot=0.0, fdom=15.0, amp=100000);

        elseif ns > 1
           source_term = get_multi_sources(isz, isx, params; ot=0.3*rand(ns), fdom=15.0, amp=100000);
        end

        # save the pressure field and boundary value
        path_pre  = joinpath(homedir(), "Desktop/pressure.rsf");
        path_bnd  = joinpath(homedir(), "Desktop/boundary_value.rsf");
        path_lwfd = joinpath(homedir(), "Desktop/last_wfd.rsf");

        # forward modeling of simultaneous sources
        multi_step_forward(source_term, params; path_pre=path_pre, path_bnd=path_bnd, path_lwfd=path_lwfd);
        (hdr, pre1) = read_RSdata(path_pre);

        # reconstruct pressure field using boundary infomation
        pre2 = pressure_reconstruct_forward(path_bnd, source_term, params);

        # reconstruct pressure field using boundary infomation
        pre3 = pressure_reconstruct_backward(path_bnd, path_lwfd, source_term, params);

        # compute the difference
        @test sqrt(sum((pre1 .- pre2) .^ 2)) / sqrt(sum(pre1 .^ 2)) < tol
        @test sqrt(sum((pre1 .- pre3) .^ 2)) / sqrt(sum(pre1 .^ 2)) < tol
    end
end


text:
abstract flowable
abstract flowblock
#abstract flowtable{T}
abstract buy_controller
abstract sell_controller

#all structures have built in storage for twice their maximum daily output,
#x24 cycle output, selling for commerce

type population <: flowable #functions as both resident and consumer
  amount::Real
end

type popblock
  pop::population
end

type workblock
  work::population
end

max_stars_repo_name: JKSano/netgameproj

max_issues_repo_name: JKSano/netgameproj

max_forks_repo_name: JKSano/netgameproj


text:
"""
# Stipple
Stipple is a reactive UI library for Julia. It provides a rich API for building rich web UIs with 2-way bindings between HTML UI elements and Julia.
It requires minimum configuration, automatically setting up the WebSockets communication channels and automatically keeping the data in sync.

Stipple allows creating powerful reactive web data dashboards using only Julia coding. It employs a declarative programming model, the framework
taking care of the full data sync workflow.
"""
module Stipple

using Logging, Reexport

@reexport using Observables
@reexport using Genie
@reexport using Genie.Renderer.Html
import Genie.Renderer.Json.JSONParser.JSONText

const Reactive = Observables.Observable
const R = Reactive

WEB_TRANSPORT = Genie.WebChannels

export R, Reactive, ReactiveModel, @R_str
export newapp

#===#

function __init__()
  Genie.config.websockets_server = true
end

#===#

abstract type ReactiveModel end

#===#

const JS_SCRIPT_NAME = "stipple.js"
const JS_DEBOUNCE_TIME = 300 #ms

#===#

function render end
function update! end
function watch end

function js_methods(m::Any)
  ""
end

#===#

const COMPONENTS = Dict()

function register_components(model::Type{M}, keysvals::Vector{Pair{K,V}}) where {M<:ReactiveModel, K, V}
  haskey(COMPONENTS, model) || (COMPONENTS[model] = Pair{K,V}[])
  push!(COMPONENTS[model], keysvals...)
end

function components(m::Type{M}) where {M<:ReactiveModel}
  haskey(COMPONENTS, m) || return ""

  response = Dict(COMPONENTS[m]...) |> Genie.Renderer.Json.JSONParser.json
  replace(response, "\""=>"")
end

#===#

function Observables.setindex!(observable::Observable, val, keys...; notify=(x)->true)
  count = 1
  observable.val = val

  for f in Observables.listeners(observable)
    if in(count, keys)
      count += 1
      continue
    end

    if notify(f)
      if f isa Observables.InternalFunction
        f(val)
      else
        Base.invokelatest(f, val)
      end
    end

    count += 1
  end

end

#===#

include("Typography.jl")
include("Elements.jl")
include("Layout.jl")
include("Generator.jl")

@reexport using .Typography
@reexport using .Elements
@reexport using .Layout
using .Generator

const newapp = Generator.newapp

#===#

function update!(model::M, field::Symbol, newval::T, oldval::T)::M where {T,M<:ReactiveModel}
  update!(model, getfield(model, field), newval, oldval)
end

function update!(model::M, field::Reactive, newval::T, oldval::T)::M where {T,M<:ReactiveModel}
  field[1] = newval

  model
end

function update!(model::M, field::Any, newval::T, oldval::T)::M where {T,M<:ReactiveModel}
  setfield!(model, field, newval)

  model
end

#===#

function watch(vue_app_name::String, fieldtype::Any, fieldname::Symbol, channel::String, debounce::Int, model::M)::String where {M<:ReactiveModel}
  js_channel = channel == "" ? "window.Genie.Settings.webchannels_default_route" : "'$channel'"
  output = """
  $vue_app_name.\$watch(function () {return this.$fieldname}, _.debounce(function(newVal, oldVal){
    Genie.WebChannels.sendMessageTo($js_channel, 'watchers', {'payload': {'field':'$fieldname', 'newval': newVal, 'oldval': oldVal}});
  }, $debounce));
  """
  # in production mode vue does not fill `this.expression` in the watcher, so we do it manually
  if Genie.Configuration.isprod()
    output *= "$vue_app_name._watchers[$vue_app_name._watchers.length - 1].expression = 'function () {return this.$fieldname}'"
  end
  output *= "\n\n"
  return output
end

#===#

function Base.parse(::Type{T}, v::T) where {T}
  v::T
end

function init(model::M, ui::Union{String,Vector} = ""; vue_app_name::String = Stipple.Elements.root(model),
              endpoint::String = JS_SCRIPT_NAME, channel::String = Genie.config.webchannels_default_route,
              debounce::Int = JS_DEBOUNCE_TIME, transport::Module = Genie.WebChannels)::M where {M<:ReactiveModel}

  global WEB_TRANSPORT = transport
  transport == Genie.WebChannels || (Genie.config.websockets_server = false)

  deps_routes(channel)

  Genie.Router.channel("/$(channel)/watchers") do
    payload = Genie.Router.@params(:payload)["payload"]
    client = Genie.Router.@params(:WS_CLIENT)

    payload["newval"] == payload["oldval"] && return "OK"

    field = Symbol(payload["field"])
    val = getfield(model, field)

    valtype = isa(val, Reactive) ? typeof(val[]) : typeof(val)

    newval = try
      if AbstractFloat >: valtype && Integer >: typeof(payload["newval"])
        convert(valtype, payload["newval"])
      else
        Base.parse(valtype, payload["newval"])
      end
    catch ex
      @error ex
      payload["newval"]
    end

    oldval = try
      if AbstractFloat >: valtype && Integer >: typeof(payload["oldval"])
        convert(valtype, payload["oldval"])
      else
        Base.parse(valtype, payload["oldval"])
      end
    catch ex
      @error ex
      payload["oldval"]
    end

    push!(model, field => newval, channel = channel, except = client)
    update!(model, field, newval, oldval)

    "OK"
  end

  ep = channel == Genie.config.webchannels_default_route ? endpoint : "js/$channel/$endpoint"
  Genie.Router.route("/$(ep)") do
    Stipple.Elements.vue_integration(model, vue_app_name = vue_app_name, endpoint = ep, channel = "", debounce = debounce) |> Genie.Renderer.Js.js
  end

  setup(model, channel)
end

function setup(model::M, channel = Genie.config.webchannels_default_route)::M where {M<:ReactiveModel}
  for f in fieldnames(typeof(model))
    isa(getproperty(model, f), Reactive) || continue

    on(getproperty(model, f)) do v
      push!(model, f => v, channel = channel)
    end
  end

  model
end

#===#

function Base.push!(app::M, vals::Pair{Symbol,T};
                    channel::String = Genie.config.webchannels_default_route,
                    except::Union{Genie.WebChannels.HTTP.WebSockets.WebSocket,Nothing,UInt} = nothing) where {T,M<:ReactiveModel}
  WEB_TRANSPORT.broadcast(channel,
                          Genie.Renderer.Json.JSONParser.json(Dict( "key" => julia_to_vue(vals[1]),
                                                                    "value" => Stipple.render(vals[2], vals[1]))),
                          except = except)
end

function Base.push!(app::M, vals::Pair{Symbol,Reactive{T}};
                    channel::String = Genie.config.webchannels_default_route,
                    except::Union{Genie.WebChannels.HTTP.WebSockets.WebSocket,Nothing,UInt} = nothing) where {T,M<:ReactiveModel}
  push!(app, Symbol(julia_to_vue(vals[1])) => vals[2][], channel = channel, except = except)
end

#===#

RENDERING_MAPPINGS = Dict{String,String}()
mapping_keys() = collect(keys(RENDERING_MAPPINGS))

function rendering_mappings(mappings = Dict{String,String})
  merge!(RENDERING_MAPPINGS, mappings)
end

function julia_to_vue(field, mapping_keys = mapping_keys())
  if in(string(field), mapping_keys)
    parts = split(RENDERING_MAPPINGS[string(field)], "-")

    if length(parts) > 1
      extraparts = map((x) -> uppercasefirst(string(x)), parts[2:end])
      string(parts[1], join(extraparts))
    else
      parts |> string
    end
  else
    field |> string
  end
end

function Stipple.render(app::M, fieldname::Union{Symbol,Nothing} = nothing)::Dict{Symbol,Any} where {M<:ReactiveModel}
  result = Dict{String,Any}()

  for field in fieldnames(typeof(app))
    result[julia_to_vue(field)] = Stipple.render(getfield(app, field), field)
  end

  Dict(:el => Elements.elem(app), :data => result, :components => components(typeof(app)), :methods => JSONText("{ $(js_methods(app)) }"), :mixins =>JSONText("[watcherMixin]"))
end

function Stipple.render(val::T, fieldname::Union{Symbol,Nothing} = nothing) where {T}
  val
end

function Stipple.render(o::Reactive{T}, fieldname::Union{Symbol,Nothing} = nothing) where {T}
  Stipple.render(o[], fieldname)
end

#===#

const DEPS = Function[]

vuejs() = Genie.Configuration.isprod() ? "vue.min.js" : "vue.js"

function deps_routes(channel::String = Genie.config.webchannels_default_route) :: Nothing
  Genie.Router.route("/js/stipple/$(vuejs())") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", vuejs()), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Router.route("/js/stipple/vue_filters.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "vue_filters.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Router.route("/js/stipple/underscore-min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "underscore-min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Router.route("/js/stipple/stipplecore.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "stipplecore.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  (WEB_TRANSPORT == Genie.WebChannels ? Genie.Assets.channels_support(channel) : Genie.Assets.webthreads_support(channel))

  nothing
end

function deps(channel::String = Genie.config.webchannels_default_route) :: String

  endpoint = (channel == Genie.config.webchannels_default_route) ?
              Stipple.JS_SCRIPT_NAME :
              "js/$(channel)/$(Stipple.JS_SCRIPT_NAME)"

  string(
    (WEB_TRANSPORT == Genie.WebChannels ? Genie.Assets.channels_support(channel) : Genie.Assets.webthreads_support(channel)),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipple/underscore-min.js"),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipple/$(vuejs())"),
    join([f() for f in DEPS], "\n"),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipple/stipplecore.js"),
    Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipple/vue_filters.js"),

    # if the model is not configured and we don't generate the stipple.js file, no point in requesting it
    in(Symbol("get_$(replace(endpoint, '/' => '_'))"), Genie.Router.named_routes() |> keys |> collect) ?
      string(
        Genie.Renderer.Html.script("Stipple.init({theme: 'stipple-blue'});"),
        Genie.Renderer.Html.script(src="$(Genie.config.base_path)$(endpoint)?v=$(Genie.Configuration.isdev() ? rand() : 1)")
      ) :
      @warn "The Reactive Model is not initialized - make sure you call Stipple.init(YourModel()) to initialize it"
  )
end

#===#

function camelcase(s::String) :: String
  replacements = [replace(s, r.match=>uppercase(r.match[2:end])) for r in eachmatch(r"_.", s) |> collect] |> unique
  isempty(replacements) ? s : first(replacements)
end

function Core.NamedTuple(kwargs::Dict) :: NamedTuple
  NamedTuple{Tuple(keys(kwargs))}(collect(values(kwargs)))
end

function Core.NamedTuple(kwargs::Dict, property::Symbol, value::String) :: NamedTuple
  value = "$value $(get!(kwargs, property, ""))" |> strip
  kwargs = delete!(kwargs, property)
  kwargs[property] = value

  NamedTuple(kwargs)
end

macro R_str(s)
  :(Symbol($s))
end

function set_multi_user_mode(value)
  global MULTI_USER_MODE = value
end

function jsonify(val; escape_untitled::Bool = true) :: String
  escape_untitled ?
    replace(Genie.Renderer.Json.JSONParser.json(val), "\"undefined\""=>"undefined") :
    Genie.Renderer.Json.JSONParser.json(val)
end

end


text:
using Base.Test
using Mimi

# Try directly using unitcheck
@test unitcheck("kg", "kg")
@test !unitcheck("kg", "MT")
@test !unitcheck("kg", "")

# Create a model with some matching and some mis-matching units
@defcomp FooUnits begin
    output = Variable(unit="kg")
end

@defcomp BarUnits begin
    input = Parameter(unit="MT")
end

@defcomp BazUnits begin
    input = Parameter(unit="kg")
end

m = Model()
setindex(m, :time, 1)
foo = addcomponent(m, FooUnits)
bar = addcomponent(m, BarUnits)
baz = addcomponent(m, BazUnits)

# Check that we cannot connect foo and bar...
@test_throws ErrorException bar[:input] = foo[:output]
# ...unless we pass ignoreunits=true
connectparameter(m, :BarUnits, :input,  :FooUnits, :output, ignoreunits=true)
# But we can connect foo and baz
baz[:input] = foo[:output]


text:
export LinearLayer, product, Primitive, NeuralKernelNetwork

using .Flux
using .Flux: softplus, @functor

# Linear layer, perform linear transformation to input array
# x₁ = W * x₀
struct LinearLayer{T, MT<:AbstractArray{T}}
    W::MT
end
@functor LinearLayer
LinearLayer(in_dim, out_dim) = LinearLayer(randn(out_dim, in_dim))
(lin::LinearLayer)(x) = softplus.(lin.W) * x

function Base.show(io::IO, layer::LinearLayer)
    print(io, "LinearLayer(", size(layer.W, 2), ", ", size(layer.W, 1), ")")
end

# Product function, given an 2d array whose size is M×N, product layer will
# multiply every m neighboring rows of the array elementwisely to obtain
# an new array of size (M÷m)×N
function product(x, step=2)
    m, n = size(x)
    m%step == 0 || error("the first dimension of inputs must be multiple of step")
    new_x = reshape(x, step, m÷step, n)
    .*([new_x[i, :, :] for i in 1:step]...)
end

# Primitive layer, mainly act as a container to hold basic kernels for the neural kernel network
struct Primitive{T}
    kernels::T
    Primitive(ks...) = new{typeof(ks)}(ks)
end
@functor Primitive

# flatten k kernel matrices of size Mk×Nk, and concatenate these 1d array into a k×(Mk*Nk) 2d array
_cat_kernel_array(x) = vcat([reshape(x[i], 1, :) for i in 1:length(x)]...)

# NOTE, though we implement `ew` & `pw` function for Primitive, it isn't a subtype of Kernel type,
# I do this because it will facilitate writing NeuralKernelNetwork
ew(p::Primitive, x) = _cat_kernel_array(map(k->ew(k, x), p.kernels))
pw(p::Primitive, x) = _cat_kernel_array(map(k->pw(k, x), p.kernels))

ew(p::Primitive, x, x′) = _cat_kernel_array(map(k->ew(k, x, x′), p.kernels))
pw(p::Primitive, x, x′) = _cat_kernel_array(map(k->pw(k, x, x′), p.kernels))

function Base.show(io::IO, layer::Primitive)
    print(io, "Primitive(")
    join(io, layer.kernels, ", ")
    print(io, ")")
end

# Neural Kernel Network, since kernel space ( stationary kernel ) is closed under linear combination
# ( with positive coefficient ) and element-wise multiplication, we can use a neural network like structure
# to build composite kernels. This type contains a `Primitive` layer which holds basic kerenls and a specialised
# nerual network architecture to perform kernel composition. It should function like a normal `Stheno` kernel.
struct NeuralKernelNetwork{PT, NNT} <: Kernel
    primitives::PT
    nn::NNT
end
@functor NeuralKernelNetwork

# use this function to reshape the 1d array back to kernel matrix
_rebuild_kernel(x, n, m) = reshape(x, n, m)
_rebuild_diag(x) = reshape(x, :)

ew(nkn::NeuralKernelNetwork, x) = _rebuild_diag(nkn.nn(ew(nkn.primitives, x)))
function pw(nkn::NeuralKernelNetwork, x)
    return _rebuild_kernel(nkn.nn(pw(nkn.primitives, x)), length(x), length(x))
end

ew(nkn::NeuralKernelNetwork, x, x′) = _rebuild_diag(nkn.nn(ew(nkn.primitives, x, x′)))
function pw(nkn::NeuralKernelNetwork, x, x′)
    return _rebuild_kernel(nkn.nn(pw(nkn.primitives, x, x′)), length(x), length(x′))
end

function Base.show(io::IO, kernel::NeuralKernelNetwork)
    print(io, "NeuralKernelNetwork(")
    join(io, [kernel.primitives, kernel.nn], ", ")
    print(io, ")")
end


text:
include("TR1DModels.jl")
include("TR1DModels3.jl")
include("driversTR.jl")

export TR1D

""" `TR1D(ϕ, a, b, stp, α, β, η₁, η₂, red. aug, M)`
   - ϕ is a linemodel
   - minimization in the interval [a,b], default [0,∞].
      - ϕ'(a) is assumed negative
   - stp : stopping object; atol and rtol are bypassed and the conditions below apply
   - α<0<β is the (possibly asymmetric) stopping tolerance, default β=-α=1e-8
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<h(a)
   - η₁ and η₂ : trust region treshold for unsuccessesful and very successful iterations, default η₁=0.25, η₂ = 0.75
   - red, aug : trust region reduction and augmentation ratio, default red=0.5, aug=5
   - M a model in the sense of trust region, usually quadratic Taylor model.

`TR1D` is a scalar trust region implementation
"""
function TR1D(ϕ       :: OneDModel;
              a       :: T = 0.0, b :: T = Inf,
              stp     :: AbstractStopping =
              NLPStopping(ϕ, OneDAtX(a, zeros(2)),
                          tol_check=(a,b,c) -> ones(2)),
              #stp     :: AbstractStopping = LS_Stopping(ϕ, LSAtT(a)),
              α       :: T = -1e-8,  β  :: T = 1e-8,
              η₁      :: Float64 = 0.25,
              η₂      :: Float64 = 0.75,
              red     :: Float64 = 0.5,
              aug     :: Float64 = 5.0,
              M       :: TR1DModel = Taylor2(0.0,0.0,0.0) )  where T
    # Specialized TR for handling non-negativity constraint on t, i.e. a <= t
    # Trust region parameters
    Δp = 1.0  # >=0
    Δn = 0.0  # <=0

    t=a
    q(d) = m(M, d)

    ϕt = obj(ϕ, t)
    dϕt = grad(ϕ,t)
    dϕtd = NaN
    ####################
    # à revoir, mais ça semble OK, on tente de trouver un meilleur minimiseur dans l'algo,
    # mais le critère d'arrêt ne teste que la descente par rapport à t=a.
    # on stoppe donc lorsque ϕ(t)<ϕ(a)   et   α ≤ dϕ(t) ≤  β
    #
    # Ce serait plus propre que ϕa=0 toujours...autrement, avec les erreurs d'arrondi,
    # le truc de tricher la valeur de ϕ courante ne marche possiblement pas... voir(*) plus bas
    ####################
    stp.meta.optimality_check = (p::OneDModel, s::AbstractState) -> [s.fx; s.gx] #
    stp.meta.atol = T(0)
    stp.meta.rtol = T(0)
    stp.meta.tol_check_neg = (a,b,c) -> [-Inf; α]
    stp.meta.tol_check     = (a,b,c) -> [  ϕt; β]

    M = initModel(M, ϕ, t, ϕt, dϕt)
    @info log_header([:iter, :f, :dual, :step, :slope, :success], [Int, T, T, T, T, String],
                     hdr_override=Dict(:f=>"h(t)", :dual=>"h'(t)", :step=>"Δn", :slope=>"Δp"))

    #OK = update_and_start!(stp, x = a, fx=ϕt, gx = dϕt)
    ϕ1  = obj(ϕ, T(1))
    dϕ1 = grad(ϕ, T(1))
    OK = update_and_start!(stp, x = T(1), fx=ϕ1, gx = dϕ1) # try one for line search

    while !OK
        dN = solveModel(M, t)
        #if (q(Δp)<q(Δn)) #| (Δn==0.0)
        if (Δm(M, Δp, Δn) < 0.0) #| (Δn==0.0)
            d=Δp
            @debug "dp"
        else
            d=Δn
            @debug "dn"
        end

        #if  (dN < Δp) & (dN > Δn) & (q(d)>q(dN))
        if  (dN < Δp) & (dN > Δn) & (Δm(M, d, dN) > 0.0)
#        if  (dN < Δp) & (dN > Δn) & ( (dϕt+0.5*H*(d+dN))*(d-dN) >0.0)
#        if  (dN < Δp) & (dN > Δn) & ( Δm(M, d-dN) >0.0)
            d=dN
            #@show "dN = ", d
        end

        ϕtestTR = obj(ϕ,t+d)
        pred = Δm(M, d)
        if pred > 0   @show pred  end
        @debug "d = ", d
        @debug "H = ", M.H
        @debug "dϕt = ", M.g
        @debug "pred = ", pred

        if (pred > 0) #|| (Δp < 1.0)
            @show pred
            @show M
            @show d, dN, Δp, Δn
            @show Δm(M, dN)
            @show q(d), q(dN), q(d)-q(dN)

            disc = M.H^2 - 2.0*M.C*M.g
            dC = disc < 0 ? -sign(M.C)*Inf : (-M.H + sqrt(disc)) / M.C
            @show dC

            tg = t+2.0*Δn:0.01*Δp:t+1.0*Δp
            n, = size(tg)
            fg=zeros(n)
            qg = zeros(n)
            dqg = zeros(n)
            for i = 1:n
                fg[i] = obj(ϕ,tg[i])
                #qg[i] = q(tg[i]-t)
                qg[i] = m(M, tg[i]-t)
                dqg[i] = dm(M, tg[i]-t)
            end
            @show "t = ", t, "d = ", d
            @show "Δp = ", Δp, "Δn = ", Δn
            plt = plot(tg,[fg, qg, dqg], label = ["f"  "p"  "dp"])
            return plt
            @assert pred<=0
        end

        goodgrad = false
        #ared = ϕtestTR-ϕt # inclure ici le truc numérique de H & Z
        if abs(d) < 1e-8  dϕtd = grad(ϕ,t+d); ared = (dϕtd + dϕt)*d/2.0; goodgrad = true
        else  ared = ϕtestTR-ϕt
        end

        ratio = ared / pred

        if ratio < η₁ # Unsuccessful
            @debug ared,  pred,  t,  d, dϕt
            Δp = red*Δp
            Δn = red*Δn
            @info log_row(Any[stp.meta.nb_of_stop, ϕt, dϕt, Δn, Δp, "U"])
        else             # Successful
            t = t + d
            if goodgrad
                dϕt = dϕtd
            else
                dϕt = grad(ϕ,t)
            end

            ϕt = ϕtestTR

            M = updateModel!(M, ϕ, t, ϕt, dϕt)

            if ratio > η₂
                Δp = aug * Δp
                Δn = max(a, aug * Δn) - t
            else
                Δn = max(a, Δn) - t
            end
            OK = update_and_stop!(stp, x=t, fx=ϕt, gx=dϕt)
            @info log_row(Any[stp.meta.nb_of_stop, ϕt, dϕt, Δn, Δp,  "S"])
        end;
    end;

    return stp

end


text:
# auxiliary function to compute necessary contributions for refined Simpson rule
function residual!(
    f!   :: Function,
    buff :: Matrix{Float64},
    a    :: Float64,
    b    :: Float64,
    n    :: Int64
    )    :: Nothing

    # compute refined step size
    h = 0.5 * (b - a) / n

    # compute integral contributions for new nodes
    for i in 1 : n
        f!(buff, a + (2.0 * i - 1.0) * h, +4.0 * h / 3.0)
    end

    # compute integral contributions for reweighted nodes
    for i in 1 : n ÷ 2
        f!(buff, a + (4.0 * i - 2.0) * h, -2.0 * h / 3.0)
    end

    return nothing
end

# integrate inplace function over an interval [a, b], with b >= a, using an adaptive Simpson rule with Richardson extrapolation on the converged result
function simps!(
    f!    :: Function,
    buff1 :: Matrix{Float64},
    buff2 :: Matrix{Float64},
    a     :: Float64,
    b     :: Float64,
    atol  :: Float64,
    rtol  :: Float64
    )     :: Nothing

    # reset result buffer
    @turbo buff1 .= 0.0

    # compute initial approximation
    m = 0.5 * (a + b)
    f!(buff1, a, 1.0 * (b - a) / 6.0)
    f!(buff1, m, 2.0 * (b - a) / 3.0)
    f!(buff1, b, 1.0 * (b - a) / 6.0)

    # compute improved approximation
    @turbo buff2 .= 0.5 .* buff1
    residual!((b, x, dx) -> f!(b, x, dx), buff2, a, b, 2)
    n = 4

    # continue improving the integral approximation until convergence is reached
    while true
        # perform Richardson extrapolation
        @turbo buff1 .*= -1.0 / 15.0
        @turbo buff1 .+= 16.0 / 15.0 .* buff2

        # compute errors
        norm1 = norm(buff1)
        norm2 = norm(buff2)
        @turbo buff1 .-= buff2
        adiff = norm(buff1)
        rdiff = adiff / max(norm1, norm2)

        # if result has converged or maximum depth (= 20) is reached, terminate
        if adiff < atol || rdiff < rtol || n > 2^20
            @turbo buff1 .+= buff2
            break
        end

        # update current solution
        @turbo buff1 .= buff2

        # compute improved approximation
        @turbo buff2 .= 0.5 .* buff1
        residual!((b, x, dx) -> f!(b, x, dx), buff2, a, b, n)
        n *= 2
    end

    return nothing
end

# integrate inplace function over an interval [a, b], with b >= a, by pre-discretizing the integration domain linearly and applying an adaptive Simpson rule to each subdomain
# note: tbuff[1] is not reset for convenience in flow integration
function integrate_lin!(
    f!    :: Function,
    tbuff :: NTuple{3, Matrix{Float64}},
    a     :: Float64,
    b     :: Float64,
    eval  :: Int64,
    atol  :: Float64,
    rtol  :: Float64
    )     :: Nothing

    @assert b >= a "Upper integration bound must be larger than or equal to lower bound."

    if b > a
        # split integration domain in subdomains of equal length
        h = (b - a) / eval

        # iterate over subdomains and apply adaptive Simpson rule
        for i in 1 : eval
            simps!((b, x, dx) -> f!(b, x, dx), tbuff[2], tbuff[3], a + (i - 1) * h, a + i * h, atol, rtol)
            tbuff[1] .+= tbuff[2]
        end
    end

    return nothing
end

# integrate inplace function over an interval [a, b], with b >= a and a > 0, by pre-discretizing the integration domain logarithmically and applying an adaptive Simpson rule to each subdomain
# note: tbuff[1] is not reset for convenience in flow integration
# note: the sgn keyword can be usend to map [a, b] -> [-b, -a]
function integrate_log!(
    f!    :: Function,
    tbuff :: NTuple{3, Matrix{Float64}},
    a     :: Float64,
    b     :: Float64,
    eval  :: Int64,
    atol  :: Float64,
    rtol  :: Float64
    ;
    sgn   :: Float64 = 1.0
    )     :: Nothing

    @assert b >= a  "Upper integration bound must be larger than or equal to lower bound."
    @assert a > 0.0 "Lower bound must be larger zero."

    if b > a
        # determine logarithmic factor
        ξ = (b / a)^(1.0 / eval)

        # iterate over subdomains and apply adaptive Simpson rule
        for i in 1 : eval
            simps!((b, x, dx) -> f!(b, sgn * x, dx), tbuff[2], tbuff[3], ξ^(i - 1) * a, ξ^i * a, atol, rtol)
            tbuff[1] .+= tbuff[2]
        end
    end

    return nothing
end


text:
module SolveDiffEq

const VecI = AbstractVector # Input  Vector
const MatI = AbstractMatrix # Input  Matrix
const MatB = AbstractMatrix # Buffer Matrix
const ArrI = AbstractArray  # Input  Array
const ArrO = AbstractArray  # Output Array

const ATOL = 1.0e-12
const RTOL = 1.0e-12

import LinearAlgebra: BLAS

if haskey(ENV, "BLAS_THREAD_NUM")
    BLAS.set_num_threads(parse(Int, ENV["BLAS_THREAD_NUM"]))
else
    BLAS.set_num_threads(4)
end

export BulirschStoer, apply_step!

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

macro def(genre::Symbol, ex::Union{Expr,Symbol}, vars::Symbol...)
    o = genre == :prop ? :(::) : genre == :vars ? :(=) : error("@def(genre = $genre, ...) is invalid.")
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        e[i] = Expr(o, vars[i], ex)
    end
    return Expr(:escape, Expr(:block, e...))
end

macro get(obj::Symbol, vars::Symbol...)
    n = length(vars)
    e = Vector{Expr}(undef, n)
    @inbounds for i in 1:n
        vari = vars[i]
        e[i] = :($vari = $obj.$vari)
    end
    return Expr(:escape, Expr(:block, e...))
end

function unsafe_cpy!(des::AbstractArray, src::AbstractArray)
    @simd for i in eachindex(des)
        @inbounds des[i] = src[i]
    end
    return nothing
end

derivative!(ηn::ArrO, ηm::ArrI, tn::Real, g!::Function) = g!(ηn, ηm, tn)

# @code_warntype ✓
# Don't use `:tuple` here, NTuple will cause
# unexpected allocations during the extrapolation.
function get_stepinv2(ns::NTuple{LMAX,Int}) where LMAX
    if @generated
        a = Vector{Expr}(undef, LMAX-1)
        for k in 1:LMAX-1
            b = Vector{Expr}(undef, LMAX-k)
            @inbounds for i in 1:LMAX-k
                b[i] = :(abs2(ns[$i] / ns[$(i+k)]) - 1.0)
            end
            @inbounds a[k] = Expr(:vect, b...)
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:vect, a...))
        end
    else
        xs = Vector{Vector{Float64}}(undef, LMAX-1)
        for k in 1:LMAX-1
            a = Vector{Float64}(undef, LMAX-k)
            @inbounds for i in 1:LMAX-k
                a[i] = abs2(ns[i] / ns[i+k]) - 1.0
            end
            xs[k] = a
        end
        return xs
    end
end

# @code_warntype ✓
function get_stepsize(ns::NTuple{LMAX,Int}, h::Real) where LMAX
    if @generated
        a = Vector{Expr}(undef, LMAX)
        @inbounds for k in 1:LMAX
            a[i] = :(h / ns[$i])
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:vect, a...))
        end
    else
        hs = Vector{Float64}(undef, LMAX)
        @simd for i in eachindex(hs)
            @inbounds hs[i] = h / ns[i]
        end
        return hs
    end
end

#=
Extrapolate a polynomial τ0 + τ1 * h² + τ2 * h⁴ + ...
by h[i] = h / n[i] where n[1] < n[2] < n[3] < ...
as
    f(h[1]²) = T[11] → T[12] → T[13]
    f(h[2]²) = T[22] → T[23] /
    f(h[3]²) = T[33] /
=#
function extrap2zero!(ys::ArrO, yr::ArrO, ym::MatI, bm::MatB, xs::VecI, ℓ::Int)
    one2n = eachindex(ys)

    @inbounds for i in one2n
        yr[i] = ym[i,ℓ]
    end

    @inbounds for j in 1:ℓ, i in one2n
        bm[i,j] = ym[i,j] - yr[i]
    end

    for k in 1:ℓ-1
        @inbounds xk = xs[k]
        for j in 1:ℓ-k
            @inbounds xkj = xk[j]
            @inbounds for i in one2n
                bm[i,j] += (bm[i,j] - bm[i,j+1]) / xkj
            end
        end
    end

    @inbounds for i in one2n
        ys[i]  = bm[i,1] + yr[i]
        yr[i] += bm[i,2]
    end

    return nothing
end

# @code_warntype ✓
function estimate_err(ys::ArrI, yr::ArrI, y0::ArrI)
    err = 0.0 # error
    @inbounds for i in eachindex(y0)
        scal = ATOL + RTOL * max(abs(y0[i]), abs(yr[i])) # scale
        err += abs2((ys[i] - yr[i]) / scal)
    end
    return sqrt(err / length(y0))
end

struct BulirschStoer{LMAX}
    @def prop Vector{Float64} hs ηt ηm ηn yr # vector of reference
    @def prop Matrix{Float64} ym bm # matrix of sequences, buffer of matrix
    xs::Vector{Vector{Float64}}
    ns::NTuple{LMAX,Int}

    BulirschStoer{13}(H::Real, ydim::Int) = BulirschStoer(H, ydim, (2,  6, 10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50))
    BulirschStoer{ 9}(H::Real, ydim::Int) = BulirschStoer(H, ydim, (2,  8, 14, 20, 26, 32, 38, 44, 50))
    BulirschStoer{ 7}(H::Real, ydim::Int) = BulirschStoer(H, ydim, (2, 10, 18, 26, 34, 42, 50))

    function BulirschStoer(H::Real, ydim::Int, ns::NTuple{LMAX,Int}) where LMAX
        @def vars Vector{Float64}(undef, ydim)       ηt ηm ηn yr
        @def vars Matrix{Float64}(undef, ydim, LMAX) ym bm
        xs = get_stepinv2(ns)
        hs = get_stepsize(ns, H)
        return new{LMAX}(hs, ηt, ηm, ηn, yr, ym, bm, xs, ns)
    end
end

# @code_warntype ✓
function apply_step!(ys::ArrO, y0::ArrI, t0::Real, g!::Function, o::BulirschStoer{LMAX}) where LMAX
    @get o ns xs hs ηt ηm ηn ym yr bm

    prevE = Inf # previous error
    one2n = eachindex(ηt)

    for ℓ in 1:LMAX
        @inbounds hℓ = hs[ℓ]
        @inbounds nℓ = ns[ℓ]
        h2 = 2.0 * hℓ

        unsafe_cpy!(ηm, y0)
        derivative!(ηn, ηm, t0, g!)
        BLAS.axpby!(1., ηm, hℓ, ηn)

        tj = t0 + hℓ; j = 1
        while j < nℓ
            derivative!(ηt, ηn, tj, g!)
            BLAS.axpby!(1., ηm, h2, ηt)
            @inbounds for i in one2n
                ηm[i], ηn[i] = ηn[i], ηt[i]
            end
            tj += hℓ; j += 1
        end

        derivative!(ηt, ηn, tj, g!)
        @inbounds for i in one2n
            ηt[i] = 0.5 * (ηn[i] + ηm[i] + hℓ * ηt[i])
        end

        if isone(ℓ)
            @inbounds for i in one2n
                ys[i] = ym[i,1] = ηt[i]
                yr[i] = ym[i,2] = 0.0
            end
        else
            @inbounds for i in one2n
                ym[i,ℓ] = ηt[i]
            end
            extrap2zero!(ys, yr, ym, bm, xs, ℓ)
        end

        thisE = estimate_err(ys, yr, y0)
        if thisE > prevE || iszero(thisE)
            break
        end
        prevE = thisE
    end

    return nothing
end

end # module


text:
@testset "polysegment" begin
    @info "Test: Polysegment"

    @testset "API" begin
        img_gray = imresize(testimage("lena_gray_256"); ratio=0.25)
        img = copy(img_gray)

        # binarize
        f = Polysegment()
        binarized_img_1 = binarize(img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_1) == Gray{N0f8}

        binarized_img_2 = binarize(Gray{Bool}, img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_2) == Gray{Bool}

        binarized_img_3 = similar(img, Bool)
        binarize!(binarized_img_3, img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_3) == Bool

        binarized_img_4 = copy(img_gray)
        binarize!(binarized_img_4, f)
        @test eltype(binarized_img_4) == Gray{N0f8}

        @test binarized_img_1 == binarized_img_2
        @test binarized_img_1 == binarized_img_3
        @test binarized_img_1 == binarized_img_4

        for T in generate_test_types([Float32, N0f8, Bool], [Gray])
            @test eltype(binarize(T, img, f)) == T
        end
    end

    @testset "Types" begin
        # Gray
        img_gray = imresize(testimage("lena_gray_256"); ratio=0.25)
        f = Polysegment()

        type_list = generate_test_types([Float32, N0f8], [Gray])
        for T in type_list
            img = T.(img_gray)
            @test_reference "References/Polysegment_Gray.png" Gray.(binarize(img, f)) by=binarization_equality()
        end

        # Color3
        img_color = imresize(testimage("lena_color_256"); ratio=0.25)
        f = Polysegment()

        type_list = generate_test_types([Float32, N0f8], [RGB, Lab])
        for T in type_list
            img = T.(img_gray)
            @test_reference "References/Polysegment_Color3.png" Gray.(binarize(img, f)) by=binarization_equality()
        end

        @testset "explicit Bool type" begin
            img_ori = imresize(testimage("cameraman"); ratio=0.25)
            for T in [RGB, Gray]
                img = T.(img_ori)
                f = Polysegment()
                for out_T in [Gray{Bool}, Bool]
                    @test binarize(out_T, img, f) == out_T.(binarize(img, f))
                end
            end
        end
    end

    @testset "Numerical" begin
        # Check that the image only has ones or zeros.
        img = imresize(testimage("lena_gray_256"); ratio=0.25)
        f = Polysegment()
        img₀₁ = binarize(img, f)
        non_zeros = findall(x -> x != 0.0 && x != 1.0, img₀₁)
        @test length(non_zeros) == 0

        # Check that ones and zeros have been assigned to the correct side of the threshold.
        maxval, maxpos = findmax(Gray.(img))
        @test img₀₁[maxpos] == 1
        minval, minpos = findmin(Gray.(img))
        @test img₀₁[minpos] == 0
    end

end