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

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

Create a CLI token to use Airflow CLI.

# Arguments
- `name`: Create a CLI token request for a MWAA environment.

"""
function create_cli_token(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "POST", "/clitoken/$(Name)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function create_cli_token(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "POST",
        "/clitoken/$(Name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_environment(dag_s3_path, execution_role_arn, name, network_configuration, source_bucket_arn)
    create_environment(dag_s3_path, execution_role_arn, name, network_configuration, source_bucket_arn, params::Dict{String,<:Any})

Creates an Amazon Managed Workflows for Apache Airflow (MWAA) environment.

# Arguments
- `dag_s3_path`: The relative path to the DAGs folder on your Amazon S3 bucket. For
  example, dags. To learn more, see Adding or updating DAGs.
- `execution_role_arn`: The Amazon Resource Name (ARN) of the execution role for your
  environment. An execution role is an AWS Identity and Access Management (IAM) role that
  grants MWAA permission to access AWS services and resources used by your environment. For
  example, arn:aws:iam::123456789:role/my-execution-role. To learn more, see Amazon MWAA
  Execution role.
- `name`: The name of the Amazon MWAA environment. For example, MyMWAAEnvironment.
- `network_configuration`: The VPC networking components used to secure and enable network
  traffic between the AWS resources for your environment. To learn more, see About networking
  on Amazon MWAA.
- `source_bucket_arn`: The Amazon Resource Name (ARN) of the Amazon S3 bucket where your
  DAG code and supporting files are stored. For example,
  arn:aws:s3:::my-airflow-bucket-unique-name. To learn more, see Create an Amazon S3 bucket
  for Amazon MWAA.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AirflowConfigurationOptions"`: A list of key-value pairs containing the Apache Airflow
  configuration options you want to attach to your environment. To learn more, see Apache
  Airflow configuration options.
- `"AirflowVersion"`: The Apache Airflow version for your environment. For example,
  v1.10.12. If no value is specified, defaults to the latest version. Valid values: v1.10.12.
- `"EnvironmentClass"`: The environment class type. Valid values: mw1.small, mw1.medium,
  mw1.large. To learn more, see Amazon MWAA environment class.
- `"KmsKey"`: The AWS Key Management Service (KMS) key to encrypt the data in your
  environment. You can use an AWS owned CMK, or a Customer managed CMK (advanced). To learn
  more, see Get started with Amazon Managed Workflows for Apache Airflow.
- `"LoggingConfiguration"`: Defines the Apache Airflow logs to send to CloudWatch Logs:
  DagProcessingLogs, SchedulerLogs, TaskLogs, WebserverLogs, WorkerLogs.
- `"MaxWorkers"`: The maximum number of workers that you want to run in your environment.
  MWAA scales the number of Apache Airflow workers up to the number you specify in the
  MaxWorkers field. For example, 20. When there are no more tasks running, and no more in the
  queue, MWAA disposes of the extra workers leaving the one worker that is included with your
  environment, or the number you specify in MinWorkers.
- `"MinWorkers"`: The minimum number of workers that you want to run in your environment.
  MWAA scales the number of Apache Airflow workers up to the number you specify in the
  MaxWorkers field. When there are no more tasks running, and no more in the queue, MWAA
  disposes of the extra workers leaving the worker count you specify in the MinWorkers field.
  For example, 2.
- `"PluginsS3ObjectVersion"`: The version of the plugins.zip file on your Amazon S3 bucket.
  A version must be specified each time a plugins.zip file is updated. To learn more, see How
  S3 Versioning works.
- `"PluginsS3Path"`: The relative path to the plugins.zip file on your Amazon S3 bucket.
  For example, plugins.zip. If specified, then the plugins.zip version is required. To learn
  more, see Installing custom plugins.
- `"RequirementsS3ObjectVersion"`: The version of the requirements.txt file on your Amazon
  S3 bucket. A version must be specified each time a requirements.txt file is updated. To
  learn more, see How S3 Versioning works.
- `"RequirementsS3Path"`: The relative path to the requirements.txt file on your Amazon S3
  bucket. For example, requirements.txt. If specified, then a file version is required. To
  learn more, see Installing Python dependencies.
- `"Schedulers"`: The number of Apache Airflow schedulers to run in your environment.
- `"Tags"`: The key-value tag pairs you want to associate to your environment. For example,
  \"Environment\": \"Staging\". To learn more, see Tagging AWS resources.
- `"WebserverAccessMode"`: The Apache Airflow Web server access mode. To learn more, see
  Apache Airflow access modes.
- `"WeeklyMaintenanceWindowStart"`: The day and time of the week to start weekly
  maintenance updates of your environment in the following format: DAY:HH:MM. For example:
  TUE:03:30. You can specify a start time in 30 minute increments only. Supported input
  includes the following:   MON|TUE|WED|THU|FRI|SAT|SUN:([01]d|2[0-3]):(00|30)
"""
function create_environment(
    DagS3Path,
    ExecutionRoleArn,
    Name,
    NetworkConfiguration,
    SourceBucketArn;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "PUT",
        "/environments/$(Name)",
        Dict{String,Any}(
            "DagS3Path" => DagS3Path,
            "ExecutionRoleArn" => ExecutionRoleArn,
            "NetworkConfiguration" => NetworkConfiguration,
            "SourceBucketArn" => SourceBucketArn,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_environment(
    DagS3Path,
    ExecutionRoleArn,
    Name,
    NetworkConfiguration,
    SourceBucketArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "PUT",
        "/environments/$(Name)",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "DagS3Path" => DagS3Path,
                    "ExecutionRoleArn" => ExecutionRoleArn,
                    "NetworkConfiguration" => NetworkConfiguration,
                    "SourceBucketArn" => SourceBucketArn,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Create a JWT token to be used to login to Airflow Web UI with claims based Authentication.

# Arguments
- `name`: Create an Airflow Web UI login token request for a MWAA environment.

"""
function create_web_login_token(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "POST", "/webtoken/$(Name)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function create_web_login_token(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "POST",
        "/webtoken/$(Name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Deletes an Amazon Managed Workflows for Apache Airflow (MWAA) environment.

# Arguments
- `name`: The name of the Amazon MWAA environment. For example, MyMWAAEnvironment.

"""
function delete_environment(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "DELETE",
        "/environments/$(Name)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_environment(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "DELETE",
        "/environments/$(Name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the details of an Amazon Managed Workflows for Apache Airflow (MWAA) environment.

# Arguments
- `name`: The name of the Amazon MWAA environment. For example, MyMWAAEnvironment.

"""
function get_environment(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "GET",
        "/environments/$(Name)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_environment(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "GET",
        "/environments/$(Name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the Amazon Managed Workflows for Apache Airflow (MWAA) environments.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to retrieve per page. For example, 5
  environments per page.
- `"NextToken"`: Retrieves the next page of the results.
"""
function list_environments(; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "GET", "/environments"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_environments(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "GET",
        "/environments",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Lists the key-value tag pairs associated to the Amazon Managed Workflows for Apache Airflow
(MWAA) environment. For example, \"Environment\": \"Staging\".

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Amazon MWAA environment. For
  example, arn:aws:airflow:us-east-1:123456789012:environment/MyMWAAEnvironment.

"""
function list_tags_for_resource(
    ResourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "GET",
        "/tags/$(ResourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    ResourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "GET",
        "/tags/$(ResourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    publish_metrics(environment_name, metric_data)
    publish_metrics(environment_name, metric_data, params::Dict{String,<:Any})

An operation for publishing metrics from the customers to the Ops plane.

# Arguments
- `environment_name`: Publishes environment metric data to Amazon CloudWatch.
- `metric_data`: Publishes metric data points to Amazon CloudWatch. CloudWatch associates
  the data points with the specified metrica.

"""
function publish_metrics(
    EnvironmentName, MetricData; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "POST",
        "/metrics/environments/$(EnvironmentName)",
        Dict{String,Any}("MetricData" => MetricData);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function publish_metrics(
    EnvironmentName,
    MetricData,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "POST",
        "/metrics/environments/$(EnvironmentName)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("MetricData" => MetricData), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Associates key-value tag pairs to your Amazon Managed Workflows for Apache Airflow (MWAA)
environment.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Amazon MWAA environment. For
  example, arn:aws:airflow:us-east-1:123456789012:environment/MyMWAAEnvironment.
- `tags`: The key-value tag pairs you want to associate to your environment. For example,
  \"Environment\": \"Staging\". To learn more, see Tagging AWS resources.

"""
function tag_resource(ResourceArn, Tags; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "POST",
        "/tags/$(ResourceArn)",
        Dict{String,Any}("Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    ResourceArn,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "POST",
        "/tags/$(ResourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Tags" => Tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes key-value tag pairs associated to your Amazon Managed Workflows for Apache Airflow
(MWAA) environment. For example, \"Environment\": \"Staging\".

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Amazon MWAA environment. For
  example, arn:aws:airflow:us-east-1:123456789012:environment/MyMWAAEnvironment.
- `tag_keys`: The key-value tag pair you want to remove. For example, \"Environment\":
  \"Staging\".

"""
function untag_resource(
    ResourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "DELETE",
        "/tags/$(ResourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    ResourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mwaa(
        "DELETE",
        "/tags/$(ResourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Updates an Amazon Managed Workflows for Apache Airflow (MWAA) environment.

# Arguments
- `name`: The name of your Amazon MWAA environment. For example, MyMWAAEnvironment.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AirflowConfigurationOptions"`: A list of key-value pairs containing the Apache Airflow
  configuration options you want to attach to your environment. To learn more, see Apache
  Airflow configuration options.
- `"AirflowVersion"`: The Apache Airflow version for your environment. For example,
  v1.10.12. If no value is specified, defaults to the latest version. Valid values: v1.10.12.
- `"DagS3Path"`: The relative path to the DAGs folder on your Amazon S3 bucket. For
  example, dags. To learn more, see Adding or updating DAGs.
- `"EnvironmentClass"`: The environment class type. Valid values: mw1.small, mw1.medium,
  mw1.large. To learn more, see Amazon MWAA environment class.
- `"ExecutionRoleArn"`: The Amazon Resource Name (ARN) of the execution role in IAM that
  allows MWAA to access AWS resources in your environment. For example,
  arn:aws:iam::123456789:role/my-execution-role. To learn more, see Amazon MWAA Execution
  role.
- `"LoggingConfiguration"`: Defines the Apache Airflow logs to send to CloudWatch Logs:
  DagProcessingLogs, SchedulerLogs, TaskLogs, WebserverLogs, WorkerLogs.
- `"MaxWorkers"`: The maximum number of workers that you want to run in your environment.
  MWAA scales the number of Apache Airflow workers up to the number you specify in the
  MaxWorkers field. For example, 20. When there are no more tasks running, and no more in the
  queue, MWAA disposes of the extra workers leaving the one worker that is included with your
  environment, or the number you specify in MinWorkers.
- `"MinWorkers"`: The minimum number of workers that you want to run in your environment.
  MWAA scales the number of Apache Airflow workers up to the number you specify in the
  MaxWorkers field. When there are no more tasks running, and no more in the queue, MWAA
  disposes of the extra workers leaving the worker count you specify in the MinWorkers field.
  For example, 2.
- `"NetworkConfiguration"`: The VPC networking components used to secure and enable network
  traffic between the AWS resources for your environment. To learn more, see About networking
  on Amazon MWAA.
- `"PluginsS3ObjectVersion"`: The version of the plugins.zip file on your Amazon S3 bucket.
  A version must be specified each time a plugins.zip file is updated. To learn more, see How
  S3 Versioning works.
- `"PluginsS3Path"`: The relative path to the plugins.zip file on your Amazon S3 bucket.
  For example, plugins.zip. If specified, then the plugins.zip version is required. To learn
  more, see Installing custom plugins.
- `"RequirementsS3ObjectVersion"`: The version of the requirements.txt file on your Amazon
  S3 bucket. A version must be specified each time a requirements.txt file is updated. To
  learn more, see How S3 Versioning works.
- `"RequirementsS3Path"`: The relative path to the requirements.txt file on your Amazon S3
  bucket. For example, requirements.txt. If specified, then a file version is required. To
  learn more, see Installing Python dependencies.
- `"Schedulers"`: The number of Apache Airflow schedulers to run in your Amazon MWAA
  environment.
- `"SourceBucketArn"`: The Amazon Resource Name (ARN) of the Amazon S3 bucket where your
  DAG code and supporting files are stored. For example,
  arn:aws:s3:::my-airflow-bucket-unique-name. To learn more, see Create an Amazon S3 bucket
  for Amazon MWAA.
- `"WebserverAccessMode"`: The Apache Airflow Web server access mode. To learn more, see
  Apache Airflow access modes.
- `"WeeklyMaintenanceWindowStart"`: The day and time of the week to start weekly
  maintenance updates of your environment in the following format: DAY:HH:MM. For example:
  TUE:03:30. You can specify a start time in 30 minute increments only. Supported input
  includes the following:   MON|TUE|WED|THU|FRI|SAT|SUN:([01]d|2[0-3]):(00|30)
"""
function update_environment(Name; aws_config::AbstractAWSConfig=global_aws_config())
    return mwaa(
        "PATCH",
        "/environments/$(Name)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_environment(
    Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mwaa(
        "PATCH",
        "/environments/$(Name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using Test
using TimeseriesSurrogates

include("all_method_tests.jl")
include("reproducibility.jl")

ext: jl


text:
module SciPyDiffEq

using Reexport
@reexport using DiffEqBase
using PyCall

abstract type SciPyAlgoritm <: DiffEqBase.AbstractODEAlgorithm end
struct RK45 <: SciPyAlgoritm end
struct RK23 <: SciPyAlgoritm end
struct Radau <: SciPyAlgoritm end
struct BDF <: SciPyAlgoritm end
struct LSODA <: SciPyAlgoritm end
struct odeint <: SciPyAlgoritm end

const integrate = PyNULL()

function __init__()
  copy!(integrate, pyimport_conda("scipy.integrate", "scipy", "conda-forge"))
end

function DiffEqBase.__solve(
    prob::DiffEqBase.AbstractODEProblem,
    alg::SciPyAlgoritm,timeseries=[],ts=[],ks=[];
    dense = true, dt = nothing, dtmax = abs(prob.tspan[2] - prob.tspan[1]),
    dtmin = eps(eltype(prob.tspan)),save_everystep = false,
    saveat=eltype(prob.tspan)[],timeseries_errors=true,
    reltol = 1e-3, abstol = 1e-6, maxiters = 10_000,
    kwargs...)

    p = prob.p
    tspan = prob.tspan
    u0 = prob.u0

    if DiffEqBase.isinplace(prob)
        f = function (t,u)
            du = similar(u)
            prob.f(du,u,p,t)
            du
        end
    else
        f = (t,u) -> prob.f(u,p,t)
    end

    _saveat = isempty(saveat) ? nothing : saveat
    if _saveat isa Array
        __saveat = _saveat
    elseif _saveat isa Number
        __saveat = Array(tspan[1]:_saveat:tspan[2])
    elseif _saveat isa Nothing
        if save_everystep
            __saveat = nothing
        else
            __saveat = [tspan[1],tspan[2]]
        end
    else
        __saveat = Array(_saveat)
    end

    if alg isa odeint
        __saveat === nothing && error("saveat is required for odeint!")
        sol,fullout = integrate.odeint(f,u0,__saveat,
                                       hmax = dtmax,
                                       rtol = reltol, atol = abstol,
                                       full_output=1, tfirst = true,
                                       mxstep = maxiters)
        tcur = fullout["tcur"]
        retcode = fullout["tcur"] == __saveat[end] ? :Success : :Failure
        ts = __saveat
        y = sol

        if typeof(u0) <: AbstractArray
            timeseries = Vector{typeof(u0)}(undef,length(ts))
            for i=1:length(ts)
                timeseries[i] = @view y[i,:]
            end
        else
            timeseries = y
        end

    else
        sol = integrate.solve_ivp(f,tspan,u0,
                                  first_step = dt,
                                  max_step = dtmax,
                                  rtol = reltol, atol = abstol,
                                  t_eval = __saveat,
                                  dense_output=dense,
                                  method = string(alg)[13:end-2])
        ts = sol["t"]
        y = sol["y"]
        retcode = sol["success"] == false ? :Failure : :Success

        if typeof(u0) <: AbstractArray
            timeseries = Vector{typeof(u0)}(undef,length(ts))
            for i=1:length(ts)
                timeseries[i] = @view y[:,i]
            end
        else
            timeseries = y
        end

    end

    if !(alg isa odeint) && dense
        _interp = PyInterpolation(sol["sol"])
    else
        _interp = DiffEqBase.LinearInterpolation(ts,timeseries)
    end

    DiffEqBase.build_solution(prob,alg,ts,timeseries,
                              interp = _interp,
                              dense = dense,
                              retcode = retcode,
                              timeseries_errors = timeseries_errors)
end

struct PyInterpolation{T} <: DiffEqBase.AbstractDiffEqInterpolation
    pydense::T
end
function (PI::PyInterpolation)(t,idxs,deriv,p,continuity)
    if idxs !== nothing
        return PI.pydense(t)[idxs]
    else
        return PI.pydense(t)
    end
end
DiffEqBase.interp_summary(::PyInterpolation) = "Interpolation from SciPy"

end # module


text:
@testset "ParticleSwarm Tuning Parameters" begin
    warning = "ParticleSwarm requires at least 3 particles. Resetting n_particles=3. " *
              "ParticleSwarm requires w ≥ 0. Resetting w=1. " *
              "ParticleSwarm requires c1 ≥ 0. Resetting c1=2. " *
              "ParticleSwarm requires c2 ≥ 0. Resetting c2=2. " *
              "ParticleSwarm requires 0 ≤ prob_shift < 1. Resetting prob_shift=0.25. "
    ps = @test_logs (:warn, warning) ParticleSwarm(
        n_particles=2, w=-1, c1=-1, c2=-1, prob_shift=2
    )
    @test ps.n_particles == 3
    @test ps.w == 1.0
    @test ps.c1 == 2.0
    @test ps.c2 == 2.0
    @test ps.prob_shift == 0.25
end

for acceleration in (CPU1(), CPUProcesses(), CPUThreads())
    @testset "EvoTree Tuning with ParticleSwarm and $(typeof(acceleration))" begin
        rng = StableRNG(123)
        features = rand(rng, 10_000) .* 5 .- 2
        X = MLJBase.table(reshape(features, (size(features)[1], 1)))
        y = sin.(features) .* 0.5 .+ 0.5
        y = EvoTrees.logit(y) + randn(rng, size(y))
        y = EvoTrees.sigmoid(y)

        tree = EvoTreeRegressor(rng=rng)
        r1 = range(tree, :max_depth; values=[3:7;])
        r2 = range(tree, :η; lower=-2, upper=0, scale=exp10)

        baseline_self_tuning_tree = TunedModel(
            model=tree,
            tuning=RandomSearch(rng=StableRNG(1234)),
            resampling=CV(nfolds=5, rng=StableRNG(8888)),
            range=[r1, r2],
            measure=mae,
            n=15,
            acceleration=acceleration
        )
        baseline_mach = machine(baseline_self_tuning_tree, X, y)
        fit!(baseline_mach, verbosity=0)
        baseline_rep = report(baseline_mach)
        baseline_best_loss = baseline_rep.best_history_entry.measurement[1]

        self_tuning_tree = TunedModel(
            model=tree,
            tuning=ParticleSwarm(rng=StableRNG(1234)),
            resampling=CV(nfolds=5, rng=StableRNG(8888)),
            range=[r1, r2],
            measure=mae,
            n=15,
            acceleration=acceleration
        )
        mach = machine(self_tuning_tree, X, y)
        fit!(mach, verbosity=0)
        rep = report(mach)
        best_loss = rep.best_history_entry.measurement[1]

        # Compare with random search result with the same settings
        @test best_loss < baseline_best_loss ||
              isapprox(best_loss, baseline_best_loss; atol=1e-3)
    end
end


text:
# - LOAD SYSTEM PACKAGES ------------------------------------------------------- #
using HTTP
# ------------------------------------------------------------------------------ #

# constants -
const path_to_package = dirname(pathof(@__MODULE__))

# - LOCAL FILES ---------------------------------------------------------------- #
include("Types.jl")
include("Network.jl")
include("Reactions.jl")
include("Sequence.jl")
include("Lists.jl")
# ------------------------------------------------------------------------------ #


text:
include("AbstractTypes/AbstractType.jl")
mutable struct PowerFlow <: PhysicalEdge
    from::Node #generator
    to::Node #MO
    quantity::Float64 # transmitted power (MWh)
    marketPrice::Float64 # EUR/MWh
end


text:
import Gmsh: gmsh

function generate_mesh(x0,radius; Lx=80e3, Ly=30e3, H_bdt = 10e3, lowres=0.02, highres=0.001, file)

    @assert length(x0) == length(radius)

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("model")

    xmin, xmax = 0.0, Lx
    ymin, ymax = 0.0, Ly
    L_bdt = Ly-H_bdt

    p1 = gmsh.model.geo.addPoint(xmin, ymin, 0, lowres)
    p2 = gmsh.model.geo.addPoint(xmax, ymin,  0, lowres)
    p3 = gmsh.model.geo.addPoint(xmax, ymax, 0, lowres)
    p4 = gmsh.model.geo.addPoint(xmin, ymax, 0, lowres)
    pbdt1 = gmsh.model.geo.addPoint(xmin, L_bdt, 0, lowres)
    pbdt2 = gmsh.model.geo.addPoint(xmax, L_bdt, 0, lowres)
    l_bot = gmsh.model.geo.addLine(p1, p2)
    lr = gmsh.model.geo.addLine(p2, p3)
    lr_bot = gmsh.model.geo.addLine(p2, pbdt2)
    lr_top = gmsh.model.geo.addLine(pbdt2,p3)
    l_top = gmsh.model.geo.addLine(p3, p4)
    ll = gmsh.model.geo.addLine(p4, p1)
    ll_top = gmsh.model.geo.addLine(p4, pbdt1)
    ll_bot = gmsh.model.geo.addLine(pbdt1,p1)
    l_bdt1 = gmsh.model.geo.addLine(pbdt1, pbdt2)
    l_bdt2 = gmsh.model.geo.addLine(pbdt2, pbdt1)
    cl_lcrust = gmsh.model.geo.addCurveLoop([l_bot, lr_bot, l_bdt2, ll_bot])
    cl_ucrust = gmsh.model.geo.addCurveLoop([l_bdt1, lr_top, l_top, ll_top])
    #cl_domain = gmsh.model.geo.addCurveLoop([l_bot,lr,l_top,ll])
    # Add two points that will be used as fixed in ux :
    # ptop = gmsh.model.geo.addPoint(0.5, 2, 0, lowres)
    # pbot = gmsh.model.geo.addPoint(0.5, 0, 0, lowres)
    # embed these points in the surface to force the mesh to define nodes at their location.
    # pbot = gmsh.model.geo.addPoint(xmax/2, ymin, 0, lowres)
    # ptop = gmsh.model.geo.addPoint(xmax/2, ymax, 0, lowres)

    #pmid = gmsh.model.geo.addPoint(xmax/2, ymax/2, 0, highres)
    #gmsh.finalize()

    #circ = gmsh.model.occ.addCircle(xmax/2, ymax/2, 0.0, xmax/20)
    T = typeof(p1)
    curved_loops = T[]
    curved_loops_phys = T[]
    surfaces = T[]
    surfaces_phys = T[]
    for i in eachindex(radius)
        x,y = x0[i]
        r = radius[i]
        pmid = gmsh.model.geo.addPoint(x, y, 0, 2highres)
        pc1 = gmsh.model.geo.addPoint(x, y-r, 0, 2highres)
        pc2 = gmsh.model.geo.addPoint(x+r, y, 0, 2highres)
        pc3 = gmsh.model.geo.addPoint(x, y+r, 0, 2highres)
        pc4 = gmsh.model.geo.addPoint(x-r, y, 0, 2highres)
        arc1 = gmsh.model.geo.addCircleArc(pc1,pmid,pc2)
        arc2 = gmsh.model.geo.addCircleArc(pc2,pmid,pc3)
        arc3 = gmsh.model.geo.addCircleArc(pc3,pmid,pc4)
        arc4 = gmsh.model.geo.addCircleArc(pc4,pmid,pc1)
        cl_seed = gmsh.model.geo.addCurveLoop([arc1,arc2,arc3,arc4])
        seed = gmsh.model.geo.addPlaneSurface([cl_seed])

        cl_seed_phys = gmsh.model.addPhysicalGroup(1, [cl_seed])
        seed_phys = gmsh.model.addPhysicalGroup(2, [seed])

        push!(curved_loops, cl_seed)
        push!(surfaces, seed)
        push!(curved_loops_phys, cl_seed_phys)
        push!(surfaces_phys, seed_phys)
    end
    # domain excluding seeds
    upper_crust = gmsh.model.geo.addPlaneSurface([cl_ucrust,curved_loops...])
    upper_crust_phys = gmsh.model.addPhysicalGroup(2, [upper_crust])

    lower_crust = gmsh.model.geo.addPlaneSurface([cl_lcrust])
    lower_crust_phys = gmsh.model.addPhysicalGroup(2, [lower_crust])

    # domain = gmsh.model.geo.addPlaneSurface([cl_domain,curved_loops...])
    # domain_phys = gmsh.model.addPhysicalGroup(2, [domain])

    # gmsh.model.geo.synchronize()
    # gmsh.model.mesh.embed(0, [l_bdt], 2, domain)
    # gmsh.model.mesh.embed(0, [ptop], 1, l3)
    #gmsh.model.mesh.embed(2, [s2], 2, s1)

    bottom = gmsh.model.addPhysicalGroup(1, [l_bot])
    right = gmsh.model.addPhysicalGroup(1, [lr_bot,lr_top])
    top = gmsh.model.addPhysicalGroup(1, [l_top])
    left = gmsh.model.addPhysicalGroup(1, [ll_bot,ll_top])
    #l_bdt = gmsh.model.addPhysicalGroup(1, [l_bdt1])

    #gmsh.model.setPhysicalName(0, clamped, "clamped")
    gmsh.model.setPhysicalName(1, bottom, "bottom")
    gmsh.model.setPhysicalName(1, right, "right")
    gmsh.model.setPhysicalName(1, top, "top")
    gmsh.model.setPhysicalName(1, left, "left")
    #gmsh.model.setPhysicalName(1, l_bdt, "bdt")
    gmsh.model.setPhysicalName(2, lower_crust_phys, "lower_crust")
    gmsh.model.setPhysicalName(2, upper_crust_phys, "upper_crust")
    #gmsh.model.setPhysicalName(2, domain_phys, "domain")
    for i in eachindex(surfaces_phys)
        gmsh.model.setPhysicalName(1, curved_loops_phys[i], "circle_seed_$i")
        gmsh.model.setPhysicalName(2, surfaces_phys[i], "seed_$i")
    end
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    gmsh.write(file)
    gmsh.finalize()

    return nothing
end

function generate_mesh2(x0,radius; Lx=80e3, Ly=30e3, H_bdt = 10e3, lowres=0.02, highres=0.001, file)

    @assert length(x0) == length(radius)

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("model")

    xmin, xmax = 0.0, Lx
    ymin, ymax = 0.0, Ly
    L_bdt = Ly-H_bdt

    p1 = gmsh.model.geo.addPoint(xmin, ymin, 0, lowres)
    p2 = gmsh.model.geo.addPoint(xmax, ymin,  0, lowres)
    p3 = gmsh.model.geo.addPoint(xmax, ymax, 0, lowres)
    p4 = gmsh.model.geo.addPoint(xmin, ymax, 0, lowres)
    pbdt1 = gmsh.model.geo.addPoint(xmin, L_bdt, 0, lowres)
    pbdt2 = gmsh.model.geo.addPoint(xmax, L_bdt, 0, lowres)
    l_bot = gmsh.model.geo.addLine(p2,p1)
    lr = gmsh.model.geo.addLine(p2, p3)
    lr_bot = gmsh.model.geo.addLine(pbdt2,p2)
    lr_top = gmsh.model.geo.addLine(pbdt2,p3)
    l_top = gmsh.model.geo.addLine(p3, p4)
    ll = gmsh.model.geo.addLine(p4, p1)
    ll_top = gmsh.model.geo.addLine(p4, pbdt1)
    ll_bot = gmsh.model.geo.addLine(p1, pbdt1)
    l_bdt1 = gmsh.model.geo.addLine(pbdt1, pbdt2)
    l_bdt2 = gmsh.model.geo.addLine(pbdt2, pbdt1)
    # cl_lcrust = gmsh.model.geo.addCurveLoop([l_bdt1, lr_bot, l_bot, ll_bot])
    # cl_ucrust = gmsh.model.geo.addCurveLoop([l_bdt1, lr_top, l_top, ll_top])
    cl_domain = gmsh.model.geo.addCurveLoop([l_bot,lr,l_top,ll])
    # Add two points that will be used as fixed in ux :
    # ptop = gmsh.model.geo.addPoint(0.5, 2, 0, lowres)
    # pbot = gmsh.model.geo.addPoint(0.5, 0, 0, lowres)
    # embed these points in the surface to force the mesh to define nodes at their location.
    # pbot = gmsh.model.geo.addPoint(xmax/2, ymin, 0, lowres)
    # ptop = gmsh.model.geo.addPoint(xmax/2, ymax, 0, lowres)

    #pmid = gmsh.model.geo.addPoint(xmax/2, ymax/2, 0, highres)
    #gmsh.finalize()

    #circ = gmsh.model.occ.addCircle(xmax/2, ymax/2, 0.0, xmax/20)
    T = typeof(p1)
    curved_loops = T[]
    curved_loops_phys = T[]
    surfaces = T[]
    surfaces_phys = T[]
    for i in eachindex(radius)
        x,y = x0[i]
        r = radius[i]
        pmid = gmsh.model.geo.addPoint(x, y, 0, 2highres)
        pc1 = gmsh.model.geo.addPoint(x, y-r, 0, 2highres)
        pc2 = gmsh.model.geo.addPoint(x+r, y, 0, 2highres)
        pc3 = gmsh.model.geo.addPoint(x, y+r, 0, 2highres)
        pc4 = gmsh.model.geo.addPoint(x-r, y, 0, 2highres)
        arc1 = gmsh.model.geo.addCircleArc(pc1,pmid,pc2)
        arc2 = gmsh.model.geo.addCircleArc(pc2,pmid,pc3)
        arc3 = gmsh.model.geo.addCircleArc(pc3,pmid,pc4)
        arc4 = gmsh.model.geo.addCircleArc(pc4,pmid,pc1)
        cl_seed = gmsh.model.geo.addCurveLoop([arc1,arc2,arc3,arc4])
        seed = gmsh.model.geo.addPlaneSurface([cl_seed])

        cl_seed_phys = gmsh.model.addPhysicalGroup(1, [cl_seed])
        seed_phys = gmsh.model.addPhysicalGroup(2, [seed])

        push!(curved_loops, cl_seed)
        push!(surfaces, seed)
        push!(curved_loops_phys, cl_seed_phys)
        push!(surfaces_phys, seed_phys)
    end
    # domain excluding seeds
    upper_crust = gmsh.model.geo.addPlaneSurface([cl_ucrust,curved_loops...])
    upper_crust_phys = gmsh.model.addPhysicalGroup(2, [upper_crust])

    lower_crust = gmsh.model.geo.addPlaneSurface([cl_lcrust])
    lower_crust_phys = gmsh.model.addPhysicalGroup(2, [lower_crust])

    # domain = gmsh.model.geo.addPlaneSurface([cl_domain,curved_loops...])
    # domain_phys = gmsh.model.addPhysicalGroup(2, [domain])

    # gmsh.model.geo.synchronize()
    # gmsh.model.mesh.embed(0, [l_bdt], 2, domain)
    # gmsh.model.mesh.embed(0, [ptop], 1, l3)
    #gmsh.model.mesh.embed(2, [s2], 2, s1)

    bottom = gmsh.model.addPhysicalGroup(1, [l_bot])
    right = gmsh.model.addPhysicalGroup(1, [lr_bot,lr_top])
    top = gmsh.model.addPhysicalGroup(1, [l_top])
    left = gmsh.model.addPhysicalGroup(1, [ll_bot,ll_top])
    #l_bdt = gmsh.model.addPhysicalGroup(1, [l_bdt1])

    #gmsh.model.setPhysicalName(0, clamped, "clamped")
    gmsh.model.setPhysicalName(1, bottom, "bottom")
    gmsh.model.setPhysicalName(1, right, "right")
    gmsh.model.setPhysicalName(1, top, "top")
    gmsh.model.setPhysicalName(1, left, "left")
    #gmsh.model.setPhysicalName(1, l_bdt, "bdt")
    gmsh.model.setPhysicalName(2, lower_crust_phys, "lower_crust")
    gmsh.model.setPhysicalName(2, upper_crust_phys, "upper_crust")
    #gmsh.model.setPhysicalName(2, domain_phys, "domain")
    for i in eachindex(surfaces_phys)
        gmsh.model.setPhysicalName(1, curved_loops_phys[i], "circle_seed_$i")
        gmsh.model.setPhysicalName(2, surfaces_phys[i], "seed_$i")
    end
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    gmsh.write(file)
    gmsh.finalize()

    return nothing
end

function generate_mesh3(x0,radius; Lx=80e3, Ly=30e3, H_bdt = 10e3, lowres=0.02, highres=0.001, el_type = Triangle, file)

    @assert length(x0) == length(radius)

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("model")

    xmin, xmax = 0.0, Lx
    ymin, ymax = 0.0, Ly
    L_bdt = Ly-H_bdt

    p1 = gmsh.model.geo.addPoint(xmin, ymin, 0, lowres)
    p2 = gmsh.model.geo.addPoint(xmax, ymin,  0, lowres)
    p3 = gmsh.model.geo.addPoint(xmax, ymax, 0, lowres)
    p4 = gmsh.model.geo.addPoint(xmin, ymax, 0, lowres)
    pbdt1 = gmsh.model.geo.addPoint(xmin-1, L_bdt, 0, lowres)
    pbdt2 = gmsh.model.geo.addPoint(xmax+1, L_bdt, 0, lowres)
    l_bot = gmsh.model.geo.addLine(p1, p2)
    lr = gmsh.model.geo.addLine(p2, p3)
    l_top = gmsh.model.geo.addLine(p3, p4)
    ll = gmsh.model.geo.addLine(p4, p1)
    l_bdt = gmsh.model.geo.addLine(pbdt1, pbdt2)
    # cl_lcrust = gmsh.model.geo.addCurveLoop([l_bot, lr_bot, l_bdt2, ll_bot])
    # cl_ucrust = gmsh.model.geo.addCurveLoop([l_bdt1, lr_top, l_top, ll_top])
    cl_domain = gmsh.model.geo.addCurveLoop([l_bot,lr,l_top,ll])
    domain = gmsh.model.geo.addPlaneSurface([cl_domain])
    # Add two points that will be used as fixed in ux :
    # ptop = gmsh.model.geo.addPoint(0.5, 2, 0, lowres)
    # pbot = gmsh.model.geo.addPoint(0.5, 0, 0, lowres)
    # embed these points in the surface to force the mesh to define nodes at their location.
    # pbot = gmsh.model.geo.addPoint(xmax/2, ymin, 0, lowres)
    # ptop = gmsh.model.geo.addPoint(xmax/2, ymax, 0, lowres)

    #pmid = gmsh.model.geo.addPoint(xmax/2, ymax/2, 0, highres)
    #gmsh.finalize()

    #circ = gmsh.model.occ.addCircle(xmax/2, ymax/2, 0.0, xmax/20)
    T = typeof(p1)
    curved_loops = T[]
    curved_loops_phys = T[]
    surfaces = T[]
    surfaces_phys = T[]
    for i in eachindex(radius)
        x,y = x0[i]
        r = radius[i]
        pmid = gmsh.model.geo.addPoint(x, y, 0, 2highres)
        pc1 = gmsh.model.geo.addPoint(x, y-r, 0, 2highres)
        pc2 = gmsh.model.geo.addPoint(x+r, y, 0, 2highres)
        pc3 = gmsh.model.geo.addPoint(x, y+r, 0, 2highres)
        pc4 = gmsh.model.geo.addPoint(x-r, y, 0, 2highres)
        arc1 = gmsh.model.geo.addCircleArc(pc1,pmid,pc2)
        arc2 = gmsh.model.geo.addCircleArc(pc2,pmid,pc3)
        arc3 = gmsh.model.geo.addCircleArc(pc3,pmid,pc4)
        arc4 = gmsh.model.geo.addCircleArc(pc4,pmid,pc1)
        cl_seed = gmsh.model.geo.addCurveLoop([arc1,arc2,arc3,arc4])
        seed = gmsh.model.geo.addPlaneSurface([cl_seed])

        gmsh.model.geo.synchronize()
        gmsh.model.mesh.embed(0, [pc1,pc2,pc3,pc4], 2, domain)
        gmsh.model.mesh.embed(1, [arc1,arc2,arc3,arc4], 2, domain)
        #cl_seed_phys = gmsh.model.addPhysicalGroup(1, [cl_seed])
        #seed_phys = gmsh.model.addPhysicalGroup(2, [seed])

        push!(curved_loops, cl_seed)
        push!(surfaces, seed)
        #push!(curved_loops_phys, cl_seed_phys)
        #push!(surfaces_phys, seed_phys)
    end
    # domain excluding seeds
    # upper_crust = gmsh.model.geo.addPlaneSurface([cl_ucrust,curved_loops...])
    # upper_crust_phys = gmsh.model.addPhysicalGroup(2, [upper_crust])
    #
    # lower_crust = gmsh.model.geo.addPlaneSurface([cl_lcrust])
    # lower_crust_phys = gmsh.model.addPhysicalGroup(2, [lower_crust])

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(1, [l_bdt], 2, domain)

    # gmsh.model.geo.synchronize()
    #gmsh.model.mesh.embed(2, [s2], 2, s1)

    bottom = gmsh.model.addPhysicalGroup(1, [l_bot])
    right = gmsh.model.addPhysicalGroup(1, [lr])
    top = gmsh.model.addPhysicalGroup(1, [l_top])
    left = gmsh.model.addPhysicalGroup(1, [ll])
    domain_phys = gmsh.model.addPhysicalGroup(2, [domain])
    #l_bdt = gmsh.model.addPhysicalGroup(1, [l_bdt1])

    #gmsh.model.setPhysicalName(0, clamped, "clamped")
    gmsh.model.setPhysicalName(1, bottom, "bottom")
    gmsh.model.setPhysicalName(1, right, "right")
    gmsh.model.setPhysicalName(1, top, "top")
    gmsh.model.setPhysicalName(1, left, "left")
    #gmsh.model.setPhysicalName(1, l_bdt, "bdt")
    # gmsh.model.setPhysicalName(2, lower_crust_phys, "lower_crust")
    # gmsh.model.setPhysicalName(2, upper_crust_phys, "upper_crust")
    gmsh.model.setPhysicalName(2, domain_phys, "domain")
    for i in eachindex(surfaces_phys)
        gmsh.model.setPhysicalName(1, curved_loops_phys[i], "circle_seed_$i")
        gmsh.model.setPhysicalName(2, surfaces_phys[i], "seed_$i")
    end
    gmsh.model.geo.synchronize()
    #gmsh.option.setNumber("Mesh.Algorithm", 2)
    #gmsh.option.setNumber("Mesh.RecombineAll", 1)
    if el_type == Quadrilateral
        gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 3)
        gmsh.model.mesh.setRecombine(2, domain)
    end
    gmsh.model.mesh.generate(2)
    gmsh.write(file)
    gmsh.finalize()

    return nothing
end

function t15()
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)

    # Copied from t1.py...
    lc = 1e-2
    gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
    gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
    gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
    gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(3, 2, 2)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.addLine(4, 1, 4)
    gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
    gmsh.model.geo.addPlaneSurface([1], 1)

    # We change the mesh size to generate a coarser mesh
    lc = lc * 4
    gmsh.model.geo.mesh.setSize([(0, 1), (0, 2), (0, 3), (0, 4)], lc)

    # We define a new point
    gmsh.model.geo.addPoint(0.02, 0.02, 0., lc, 5)

    # We have to synchronize before embedding entites:
    gmsh.model.geo.synchronize()

    # One can force this point to be included ("embedded") in the 2D mesh, using the
    # `embed()' function:
    gmsh.model.mesh.embed(0, [5], 2, 1)

    # In the same way, one can use `embed()' to force a curve to be embedded in the
    # 2D mesh:
    gmsh.model.geo.addPoint(0.02, 0.12, 0., lc, 6)
    gmsh.model.geo.addPoint(0.04, 0.18, 0., lc, 7)
    gmsh.model.geo.addLine(6, 7, 5)

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(1, [5], 2, 1)

    # Points and curves can also be embedded in volumes
    gmsh.model.geo.extrude([(2, 1)], 0, 0, 0.1)

    p = gmsh.model.geo.addPoint(0.07, 0.15, 0.025, lc)

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(0, [p], 3, 1)

    gmsh.model.geo.addPoint(0.025, 0.15, 0.025, lc, p + 1)
    l = gmsh.model.geo.addLine(7, p + 1)

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(1, [l], 3, 1)

    # Finally, we can also embed a surface in a volume:
    gmsh.model.geo.addPoint(0.02, 0.12, 0.05, lc, p + 2)
    gmsh.model.geo.addPoint(0.04, 0.12, 0.05, lc, p + 3)
    gmsh.model.geo.addPoint(0.04, 0.18, 0.05, lc, p + 4)
    gmsh.model.geo.addPoint(0.02, 0.18, 0.05, lc, p + 5)

    gmsh.model.geo.addLine(p + 2, p + 3, l + 1)
    gmsh.model.geo.addLine(p + 3, p + 4, l + 2)
    gmsh.model.geo.addLine(p + 4, p + 5, l + 3)
    gmsh.model.geo.addLine(p + 5, p + 2, l + 4)

    ll = gmsh.model.geo.addCurveLoop([l + 1, l + 2, l + 3, l + 4])
    s = gmsh.model.geo.addPlaneSurface([ll])

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(2, [s], 3, 1)

    # Note that with the OpenCASCADE kernel (see `t16.py'), when the `fragment()'
    # function is applied to entities of different dimensions, the lower dimensional
    # entities will be autmatically embedded in the higher dimensional entities if
    # necessary.

    gmsh.model.mesh.generate(3)

    gmsh.write("t15.msh")

    gmsh.finalize()

end

function test_embed()
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)

    # Copied from t1.py...
    lc = 1e-2
    xmax, ymax = 1.0, 2.0
    gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
    gmsh.model.geo.addPoint(xmax, 0, 0, lc, 2)
    gmsh.model.geo.addPoint(xmax, ymax, 0, lc, 3)
    gmsh.model.geo.addPoint(0, ymax, 0, lc, 4)
    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(3, 2, 2)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.addLine(4, 1, 4)
    gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
    gmsh.model.geo.addPlaneSurface([1], 1)

    # In the same way, one can use `embed()' to force a curve to be embedded in the
    # 2D mesh:
    gmsh.model.geo.addPoint(-0.0001, ymax/2, 0., lc, 5)
    gmsh.model.geo.addPoint(xmax+0.0001, ymax/2, 0., lc, 6)
    gmsh.model.geo.addLine(5, 6, 5)

    gmsh.model.geo.synchronize()
    gmsh.model.mesh.embed(1, [5], 2, 1)

    gmsh.model.mesh.generate(2)

    gmsh.write("test_embed.msh")

    gmsh.finalize()

end

function test_embed2()
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)

    # Copied from t1.py...
    lc = 1e-2
    xmax, ymax = 1.0, 2.0
    gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
    gmsh.model.geo.addPoint(xmax, 0, 0, lc, 2)
    gmsh.model.geo.addPoint(xmax, ymax, 0, lc, 3)
    gmsh.model.geo.addPoint(0, ymax, 0, lc, 4)
    gmsh.model.geo.addPoint(0.0001, ymax/2, 0., lc, 5)
    gmsh.model.geo.addPoint(xmax-0.0001, ymax/2, 0., lc, 6)
    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(2, 6, 2)
    gmsh.model.geo.addLine(6, 5, 3)
    gmsh.model.geo.addLine(5, 1, 4)

    gmsh.model.geo.addLine(6, 3, 5)
    gmsh.model.geo.addLine(3, 4, 6)
    gmsh.model.geo.addLine(4, 5, 7)

    gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
    gmsh.model.geo.addCurveLoop([-3, 5, 6, 7], 2)
    gmsh.model.geo.addPlaneSurface([1,2], 1)

    # In the same way, one can use `embed()' to force a curve to be embedded in the
    # 2D mesh:
    # gmsh.model.geo.addPoint(0.0001, ymax/2, 0., lc, 5)
    # gmsh.model.geo.addPoint(xmax-0.0001, ymax/2, 0., lc, 6)
    # gmsh.model.geo.addLine(5, 6, 5)
    #
    # gmsh.model.geo.synchronize()
    # gmsh.model.mesh.embed(1, [5], 2, 1)

    gmsh.model.mesh.generate(2)

    gmsh.write("test_embed2.msh")

    gmsh.finalize()

end


text:
module panel
import Term:
    split_lines,
    reshape_text,
    do_by_line,
    join_lines,
    truncate,
    textlen,
    fillin,
    do_by_line

import ..console: console_width, console_height
import ..measure: Measure
import ..renderables: AbstractRenderable, RenderablesUnion, Renderable, RenderableText
import ..segment: Segment
using ..box
import ..layout: pad, vstack, Padding
import ..style: apply_style

export Panel, TextBox

abstract type AbstractPanel <: AbstractRenderable end

# ---------------------------------------------------------------------------- #
#                                     PANEL                                    #
# ---------------------------------------------------------------------------- #
"""
    Panel

`Renderable` with a panel surrounding some content:
        ╭──────────╮
        │ my panel │
        ╰──────────╯
"""
mutable struct Panel <: AbstractPanel
    segments::Vector
    measure::Measure
    function Panel(x1, x2; kwargs...)
        # this is necessary to handle the special case in which 2 objs are passed
        if x1 isa Vector
            return new(x1, x2)
        else
            return Panel(vstack(x1, x2); kwargs...)
        end
    end
end

"""
    Panel(
            style::String,
            title_style::String,
            title_justify::Symbol,
            subtitle::String,
            subtitle_style::String,
            subtitle_justify::String,
            panel_measure::Measure,
            content_measure::Measure,
            Δw::Int,
            Δh::Int,
            padding::Padding
        )

Construct a `Panel` given all required info.
"""
function render(
                content;
                box::Symbol=:ROUNDED,
                style::String="default",
                title::Union{String,Nothing} = nothing,
                title_style::Union{Nothing, String} = nothing,
                title_justify::Symbol = :left,
                subtitle::Union{String,Nothing} = nothing,
                subtitle_style::Union{Nothing, String} = nothing,
                subtitle_justify::Symbol = :left,
                justify::Symbol = :left,
                panel_measure::Measure,
                content_measure::Measure,
                Δw::Int,
                Δh::Int,
                padding::Padding
    )::Panel

    # create top/bottom rows with titles
    box = eval(box)  # get box object from symbol
    top = get_title_row(
        :top,
        box,
        title;
        width = panel_measure.w,
        style = style,
        title_style = title_style,
        justify = title_justify,
    )

    bottom = get_title_row(
        :bottom,
        box,
        subtitle;
        width = panel_measure.w,
        style = style,
        title_style = subtitle_style,
        justify = subtitle_justify,
    )

    # get left/right vertical lines
    σ(s) = apply_style("[" * style * "]" * s * "[/" * style * "]")
    left, right = σ(box.mid.left), "\e[0m"*σ(box.mid.right)*"\e[0m"

    # get an empty padding line
    empty = [Segment(left * " "^(panel_measure.w-2) * right)]

    # add lines with content fn
    function makecontent_line(cline)::Segment
        line = pad(apply_style(cline), panel_measure.w-Δw, justify)
        line = pad(line, padding.left, padding.right)

        # make line
        return Segment(left * line * right)
    end

    # check if we need extra lines at the bottom to reach target height
    if content_measure.h < panel_measure.h - 2 - Δh
        n_extra = panel_measure.h-2 - content_measure.h
    else
        n_extra = 0
    end

    # create segments
    initial_segments::Vector{Segment} = [
        top,                                        # top border
        repeat(empty, padding.top)...,              # top padding
    ]

    # content
    content_sgs::Vector{Segment} = content.measure.w > 0 ?  map(s -> makecontent_line(s.text), content.segments) : []

    final_segments::Vector{Segment} = [
        repeat(empty, n_extra)...,                  # lines to reach target height
        repeat(empty, padding.bottom)...,           # bottom padding
        bottom,                                     # bottom border
    ]

    segments = vcat(initial_segments, content_sgs, final_segments)

    return Panel(
        segments, panel_measure
    )
end

"""
    Panel(;
        fit::Symbol=:nofit,
        width::Int = 88,
        height::Union{Nothing,Int} = nothing,
        kwargs...
    )

Construct a `Panel` with no content
"""
function Panel(;
            fit::Bool=false,
            width::Int = 88,
            height::Int = 2,
            padding::Union{Vector, Padding, NTuple} = Padding(0, 0, 0, 0),
            kwargs...
    )
    # get panel measure
    if fit
        # hardcoded size of empty 'fitted' panel
        panel_measure = Measure(3, 2)
    else
        panel_measure = Measure(width, height)
    end

    # get empty content measure
    content = ""
    content_measure = Measure(0, 0)

    # get padding
    padding = padding isa Padding ? padding : Padding(padding...)

    # make panel
    return Panel(
        content;
        panel_measure=panel_measure,
        content_measure=content_measure,
        Δw=padding.left + padding.right + 2,
        Δh=padding.top + padding.bottom,
        padding=padding,
        kwargs...
    )
end

"""
    Panel(
        content::AbstractRenderable;
        width::Int = 88,
        height::Union{Nothing,Int} = nothing,
        fit::Bool=false,
        padding::Union{Padding, NTuple} = Padding(2, 2, 0, 0),
        kwargs...
    )

Construct a `Panel` around of a `AbstractRenderable`
"""
function Panel(
        content::Union{AbstractString, AbstractRenderable};
        width::Int = 88,
        height::Union{Nothing,Int} = nothing,
        fit::Bool=false,
        padding::Union{Padding, NTuple} = Padding(2, 2, 0, 0),
        kwargs...
    )

    content = content isa RenderableText ? string(content) : content

    # get measure
    WIDTH = console_width(stdout)
    padding = padding isa Padding ? padding : Padding(padding...)
    Δw = padding.left + padding.right + 2
    Δh = padding.top + padding.bottom

    # define convenience function
    function resize_content(content, _width)
        if content isa AbstractString
            content = reshape_text(content, _width)
        end
        return content, Measure(content)
    end

    # get measure of panel's box and optionally resize text.
    content_measure = Measure(content)
    if fit
        # if content width too large, resize content if its text renderable
        if content_measure.w > WIDTH - Δw
            content, content_measure = resize_content(content, WIDTH - Δw )
            panel_measure = Measure(WIDTH, content_measure.h+Δh+2)
        else
            panel_measure = Measure(content_measure.w+Δw, content_measure.h+Δh+2)
        end
    end

    if !fit
        # check that the content fits within the given width
        if content isa AbstractString
            width = min(width, WIDTH)
            if content_measure.w > width-Δw
                content, content_measure = resize_content(content, width-Δw)
            end
        else
            # if width too small for content, try to enlarge
            width = width < content_measure.w+Δw ?  min(content_measure.w+Δw, WIDTH-Δw) : width
        end

        # get target height
        _h = content_measure.h + Δh + 2
        height = isnothing(height) ? _h  : max(height, _h)
        panel_measure = Measure(width, height)
    end

    return render(
        Renderable(content);
        panel_measure=panel_measure,
        content_measure=content_measure,
        Δw=Δw,
        Δh=Δh,
        padding=padding,
        kwargs...
    )
end

"""
    Panel(renderables; kwargs...)

`Panel` constructor for creating a panel out of multiple renderables at once.
"""
Panel(renderables::Vector{RenderablesUnion}; kwargs...) = Panel(vstack(renderables...); kwargs...)
Panel(texts::Vector{AbstractString}; kwargs...) = Panel(join_lines(texts); kwargs...)
Panel(renderables...; kwargs...) = Panel(vstack(renderables...); kwargs...)
# ---------------------------------------------------------------------------- #
#                                    TextBox                                   #
# ---------------------------------------------------------------------------- #

"""
    TextBox

Creates a `Panel` and fits input text to it.
The pannel is hidden so that the result is just a text box.
"""
mutable struct TextBox <: AbstractPanel
    segments::Vector
    measure::Measure
end

"""
    TextBox(
        text::Union{Vector, AbstractString};
        width::Union{Nothing, Int}=nothing,
        title::Union{Nothing, String}=nothing,
        title_style::Union{String, Nothing}="default",
        title_justify::Symbol=:left,
        subtitle::Union{String, Nothing}=nothing,
        subtitle_style::Union{String, Nothing}="default",
        subtitle_justify::Symbol=:left,
        justify::Symbol=:left,
        fit::Symbol=:fit,
        )

Creates an hidden `Panel` with `text` in it.

If a `width` is passed, the input `text` is reshaped to have
that size, unless `fit=:truncate` in which case it's cut to size.
If no `width` is passed and `fit=true` the `TextBox`'s size
matches the size of the input `text`.
Other arguments behave like `Panel`.

See also [`Panel`](@ref).
"""
function TextBox(
    text::Union{AbstractString};
    width::Int = 88,
    title::Union{Nothing,String} = nothing,
    title_style::Union{String,Nothing} = "default",
    title_justify::Symbol = :left,
    subtitle::Union{String,Nothing} = nothing,
    subtitle_style::Union{String,Nothing} = "default",
    subtitle_justify::Symbol = :left,
    justify::Symbol = :left,
    fit::Symbol = :nofit,
    padding::Union{Padding, NTuple} = Padding(2, 2, 0, 0),
)
    padding = padding isa Padding ? padding : Padding(padding...)
    Δw =  padding.left + padding.right

    text = apply_style(text)

    # fit text or get width
    if fit==:fit
        # the box's size depends on the text's size
        width = Measure(text).w + Δw

        # too large, fit to console
        if width > console_width(stdout)
            width=console_width(stdout)
            fit = true
        end
    else
        width = width > console_width(stdout) ? console_width(stdout) - Δw : width
    end

    # truncate or reshape text
    if fit == :truncate
        # truncate the text to fit the given width
        text = do_by_line(ln -> truncate(ln, width - 7), apply_style(text))
    else
        text = reshape_text(text, width-Δw-2)
    end

    # create panel with text inside
    panel = Panel(
        text;
        style = "hidden",
        title = title,
        title_style = title_style,
        title_justify = title_justify,
        subtitle = subtitle,
        subtitle_style = subtitle_style,
        subtitle_justify = subtitle_justify,
        justify = justify,
        width = width,
        fit=false,
        padding=padding
    )

    return TextBox(panel.segments, Measure(panel.measure.w, panel.measure.h))
end

TextBox(texts...; kwargs...) = TextBox(join_lines(texts...); kwargs...)

end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.cluster;
#
# import typedb.client.api.database.Database;
# import typedb.client.api.database.DatabaseManager;
# import typedb.client.common.exception.TypeDBClientException;
# import typedb.client.core.CoreDatabaseManager;
# import typedb.common.collection.Pair;
# import typedb.protocol.ClusterDatabaseProto;
#
# import java.util.List;
# import java.util.Map;
#
# import static typedb.client.common.exception.ErrorMessage.Client.CLUSTER_ALL_NODES_FAILED;
# import static typedb.client.common.rpc.RequestBuilder.Cluster.DatabaseManager.allReq;
# import static typedb.client.common.rpc.RequestBuilder.Cluster.DatabaseManager.getReq;
# import static typedb.common.collection.Collections.pair;
# import static java.util.stream.Collectors.toList;
# import static java.util.stream.Collectors.toMap;
#
# public class ClusterDatabaseManager implements DatabaseManager.Cluster {
#     private final Map<String, CoreDatabaseManager> databaseMgrs;
#     private final ClusterClient client;
#
#     public ClusterDatabaseManager(ClusterClient client) {
#         this.client = client;
#         this.databaseMgrs = client.coreClients().entrySet().stream()
#                 .map(c -> pair(c.getKey(), c.getValue().databases()))
#                 .collect(toMap(Pair::first, Pair::second));
#         ;
#     }
#
#     @Override
#     public boolean contains(String name) {
#         StringBuilder errors = new StringBuilder();
#         for (String address : databaseMgrs.keySet()) {
#             try {
#                 return databaseMgrs.get(address).contains(name);
#             } catch (TypeDBClientException e) {
#                 errors.append("- ").append(address).append(": ").append(e).append("\n");
#             }
#         }
#         throw new TypeDBClientException(CLUSTER_ALL_NODES_FAILED, errors.toString());
#     }
#
#     @Override
#     public void create(String name) {
#         for (CoreDatabaseManager databaseMgr : databaseMgrs.values()) {
#             if (!databaseMgr.contains(name)) {
#                 databaseMgr.create(name);
#             }
#         }
#     }
#
#     @Override
#     public Database.Cluster get(String name) {
#         StringBuilder errors = new StringBuilder();
#         for (String address : databaseMgrs.keySet()) {
#             try {
#                 ClusterDatabaseProto.ClusterDatabaseManager.Get.Res res = client.stub(address).databasesGet(getReq(name));
#                 return ClusterDatabase.of(res.getDatabase(), this);
#             } catch (TypeDBClientException e) {
#                 errors.append("- ").append(address).append(": ").append(e).append("\n");
#             }
#         }
#         throw new TypeDBClientException(CLUSTER_ALL_NODES_FAILED, errors.toString());
#     }
#
#     @Override
#     public List<Database.Cluster> all() {
#         StringBuilder errors = new StringBuilder();
#         for (String address : databaseMgrs.keySet()) {
#             try {
#                 ClusterDatabaseProto.ClusterDatabaseManager.All.Res res = client.stub(address).databasesAll(allReq());
#                 return res.getDatabasesList().stream().map(db -> ClusterDatabase.of(db, this)).collect(toList());
#             } catch (TypeDBClientException e) {
#                 errors.append("- ").append(address).append(": ").append(e).append("\n");
#             }
#         }
#         throw new TypeDBClientException(CLUSTER_ALL_NODES_FAILED, errors.toString());
#     }
#
#     Map<String, CoreDatabaseManager> databaseMgrs() {
#         return databaseMgrs;
#     }
# }


text:
#=
  These routines return +/-Inf when given +/-Inf.
  That differs from nextfloat(-Inf) == -realmax(), prevfloat(Inf) == realmax()
  (prevfloat(Inf)==Inf makes more sense to me, and likely is more helpful).
  And steps twice when given values of very small magnitude (see paper).
  The alternative implementation, converting to [U]Int and adding/subtracting 1,
  returns NaN when given +/-Inf -- and checking for Inf adds branching.
  ref:
  'Predecessor and Successor In Rounding To Nearest' by Siegried M. Rump
  "The routines deliver the exact answer except for a small range near underflow,
   in which case the true result is overestimated by eta [the value added/subtracted below]."
=#
# exact for |x| > 8.900295434028806e-308
nextNearerToZero(x::Float64)   = (x-1.1102230246251568e-16*x)-5.0e-324 # do not simplify
nextAwayFromZero(x::Float64)   = (x+1.1102230246251568e-16*x)+5.0e-324 # do not simplify
# exact for |x| > 4.814825f-35
nextNearerToZero(x::Float32)   = (x-5.960465f-8*x)-1.435f-42 # do not simplify
nextAwayFromZero(x::Float32)   = (x+5.960465f-8*x)+1.435f-42 # do not simplify
# the multiplicative formulation for Float16 is exact for |x| > Float16(0.25)
# which is quite coarse, we do not use that here
nextNearerToZero(x::Float16) = (x < 0) ? nextfloat(x) : prevfloat(x)
nextAwayFromZero(x::Float16) = (x < 0) ? prevfloat(x) : nextfloat(x)
nextFloat(x::Float16) = nextfloat(x)
prevFloat(x::Float16) = prevfloat(x)

nextFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextNearerToZero(x) : nextAwayFromZero(x)
prevFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextAwayFromZero(x) : nextNearerToZero(x)


text:
using Poirot, Test, Random
using Poirot: Singleton

Random.seed!(0)

@test infer(() -> 1) == Singleton(1)
@test infer(() -> "foo") == Singleton("foo")

@test infer(() -> rand()) == Uniform(0, 1)
@test infer(() -> randn()) == Normal(0, 1)
@test infer(() -> rand(Bool)) == Bernoulli(1//2)

coin() = rand(Bool)

@test infer(() -> coin() & coin()) == Bernoulli(0.25)
@test infer(() -> coin() | coin()) == Bernoulli(0.75)

@test infer() do
  a = coin()
  b = coin()
  observe(a | b)
  a & b
end == Bernoulli(1/3)

infer() do
  a = coin()
  b = coin()
  observe(a | b)
  observe(!b)
  a
end == Singleton(true)

@test infer() do
  x = randn()
  observe(x > 0)
  x
end isa Poirot.Empirical

d = infer() do
  x = 10*rand()
  x^2 > 50
end

@test d isa Bernoulli
@test mean(d) ≈ (10-sqrt(50))/10

d = infer() do
  x = 10*rand()
  observe(x > 5)
  x^2 > 50
end

@test d isa Bernoulli
@test mean(d) ≈ (10-sqrt(50))/5

d = infer() do
  identical = rand(Bernoulli(1/3))
  boy1 = rand(Bernoulli(1/2))
  # TODO first version gives a different answer
  # boy2 = identical ? boy1 : rand(Bernoulli(1/2))
  boy2 = rand(Bernoulli(identical ? boy1*1.0 : 1/2))
  observe(boy1 & boy2)
  identical
end

@test mean(d) == 0.5


text:
using LinearAlgebra: lmul!, norm
using Parameters: Parameters
using NQCDynamics.NonadiabaticDistributions: ElectronicDistribution
using NQCDynamics: nbeads
using NQCModels: nstates

function NQCDynamics.RingPolymerSimulation{eCMM}(atoms::Atoms{S,T}, model::Model, n_beads::Integer; γ=0, kwargs...) where {S,T}
    NQCDynamics.RingPolymerSimulation(atoms, model, eCMM{T}(nstates(model), γ), n_beads; kwargs...)
end

function DynamicsMethods.DynamicsVariables(sim::RingPolymerSimulation{<:eCMM}, v, r, ::ElectronicDistribution)

    F = nstates(sim)
    n_beads = length(sim.beads)

    radius = sqrt(2 + 2F*sim.method.γ*n_beads)
    points = generate_random_points_on_nsphere(2F*n_beads, radius)
    half = length(points) ÷ 2
    qmap = points[begin:half]
    pmap = points[half+1:end]
    qmap = reshape(qmap, (F, n_beads))
    pmap = reshape(pmap, (F, n_beads))

    return ComponentVector(v=v, r=r, pmap=pmap, qmap=qmap)
end

function acceleration!(dv, u, sim::RingPolymerSimulation{<:eCMM})

    Parameters.@unpack γ, temp_q, temp_p = sim.method

    Calculators.evaluate_derivative!(sim.calculator, DynamicsUtils.get_positions(u))
    for I in CartesianIndices(dv)
        qmap = get_mapping_positions(u, I[3])
        pmap = get_mapping_momenta(u, I[3])
        D = sim.calculator.derivative[I]
        D̄ = tr(D) / nstates(sim)
        Dtraceless = D - Diagonal(fill(D̄, nstates(sim)))
        mul!(temp_q, Dtraceless, qmap)
        mul!(temp_p, Dtraceless, pmap)
        dv[I] = dot(qmap, temp_q)
        dv[I] += dot(pmap, temp_p)
        dv[I] += 2D̄
    end
    lmul!(-1/2, dv)
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masses)
    DynamicsUtils.apply_interbead_coupling!(dv, DynamicsUtils.get_positions(u), sim)
end

function set_mapping_force!(du, u, sim::RingPolymerSimulation{<:eCMM})
    Calculators.evaluate_potential!(sim.calculator, DynamicsUtils.get_positions(u))
    for i in range(sim.beads)
        V = sim.calculator.potential[i]
        V̄ = tr(V) / nstates(sim)
        Vtraceless = V - Diagonal(fill(V̄, nstates(sim)))
        mul!(get_mapping_positions(du, i), Vtraceless, get_mapping_momenta(u, i))
        mul!(get_mapping_momenta(du, i), Vtraceless, get_mapping_positions(u, i))
        lmul!(-1, get_mapping_momenta(du, i))
    end
end

function DynamicsUtils.classical_hamiltonian(sim::RingPolymerSimulation{<:eCMM}, u)
    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)

    spring = RingPolymers.get_spring_energy(sim.beads, sim.atoms.masses, r)
    kinetic = DynamicsUtils.classical_kinetic_energy(sim, v)

    Calculators.evaluate_potential!(sim.calculator, r)
    potential = zero(eltype(u))
    for i=1:nbeads(sim)
        V = sim.calculator.potential[i]
        V̄ = tr(V) / nstates(sim)
        Vtraceless = V - Diagonal(fill(V̄, nstates(sim)))
        qmap = get_mapping_positions(u, i)
        pmap = get_mapping_momenta(u, i)
        potential += 0.5 * (pmap'Vtraceless*pmap + qmap'Vtraceless*qmap) + V̄
    end

    return kinetic + spring + potential
end

function inverse_mapping_kernel(qmap, pmap, γ, N)
    F = length(qmap)
    prefactor = (1 + N*F) / (1 + N*F*γ)^2
    subtractor = (1 - γ) / (1 + N*F*γ)
    return prefactor .* (qmap.^2 .+ pmap.^2)./2 .- subtractor
end

function Estimators.diabatic_population(sim::RingPolymerSimulation{<:eCMM}, u)
    population = zeros(nstates(sim))
    for i=1:nbeads(sim)
        qmap = get_mapping_positions(u, i)
        pmap = get_mapping_momenta(u, i)
        population .+= inverse_mapping_kernel(qmap, pmap, sim.method.γ, nbeads(sim))
    end
    return population
end

function Estimators.initial_diabatic_population(sim::RingPolymerSimulation{<:eCMM}, u)
    population = zeros(nstates(sim))
    for i=1:nbeads(sim)
        qmap = get_mapping_positions(u, i)
        pmap = get_mapping_momenta(u, i)
        population .+= mapping_kernel(qmap, pmap, sim.method.γ)
    end
    return population
end


text:
import TAMode
using Test
using Turing

f = open("log.txt", "w")
redirect_stdout(f)
redirect_stderr(f)

@testset "Test sampling." begin
    samp = sample(TAMode.A549model, NUTS(0.65), 500)
    write("chain-file-10_1.jls", samp)
end

close(f)


text:
import InteractiveUtils
import Pkg
import Test

@debug(string("Julia depot paths: "), Base.DEPOT_PATH)
@debug(string("Julia load paths: "), Base.LOAD_PATH)

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
@debug(string("Julia version info: ",))
if !isnothing(logger)
    InteractiveUtils.versioninfo(logger.stream; verbose=true)
end

@debug(string("Attempting to import PredictMD...",))
import PredictMD
@debug(string("Successfully imported PredictMD.",))
@debug(string("PredictMD version: "),PredictMD.version(),)
@debug(string("PredictMD package directory: "),PredictMD.package_directory(),)

@debug(string("Julia depot paths: "), Base.DEPOT_PATH)
@debug(string("Julia load paths: "), Base.LOAD_PATH)

PredictMD.import_all()


text:
struct Vern6Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern6Cache) = (c.utilde,c.atmp)
du_cache(c::Vern6Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9)

function alg_cache(alg::Vern6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Vern6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype);
  k5 = zero(rate_prototype); k6 = zero(rate_prototype); k7 = zero(rate_prototype);
  k8 = zero(rate_prototype); k9 = zero(rate_prototype);
  utilde = similar(u,axes(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u));
  Vern6Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,utilde,tmp,atmp,tab)
end

alg_cache(alg::Vern6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Vern6ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct Vern7Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern7Cache) = (c.utilde,c.atmp)
du_cache(c::Vern7Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10)

function alg_cache(alg::Vern7,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Vern7ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype);
  k5 = zero(rate_prototype); k6 = zero(rate_prototype); k7 = zero(rate_prototype); k8 = zero(rate_prototype);
  k9 = zero(rate_prototype); k10 = zero(rate_prototype); utilde = similar(u,axes(u))
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u))
  Vern7Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,tmp,atmp,tab)
end

alg_cache(alg::Vern7,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Vern7ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct Vern8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern8Cache) = (c.utilde,c.atmp)
du_cache(c::Vern8Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13)

function alg_cache(alg::Vern8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Vern8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype);
  k4 = zero(rate_prototype);
  k5 = zero(rate_prototype); k6 = zero(rate_prototype); k7 = zero(rate_prototype);
  k8 = zero(rate_prototype); tmp = similar(u)
  k9 = zero(rate_prototype); k10 = zero(rate_prototype); k11 = zero(rate_prototype);
  k12 = zero(rate_prototype); k13 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits,axes(u))
  Vern8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,tab)
end

alg_cache(alg::Vern8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Vern8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct Vern9Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern9Cache) = (c.utilde,c.atmp)
du_cache(c::Vern9Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16)

function alg_cache(alg::Vern9,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Vern9ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype);k3 = zero(rate_prototype);
  k4 = zero(rate_prototype);
  k5 = zero(rate_prototype); k6 = zero(rate_prototype);k7 = zero(rate_prototype);
  k8 = zero(rate_prototype);
  k9 = zero(rate_prototype); k10 = zero(rate_prototype); k11 = zero(rate_prototype);
  k12 = zero(rate_prototype);
  k13 = zero(rate_prototype); k14 = zero(rate_prototype); k15 = zero(rate_prototype);
  k16 =zero(rate_prototype);
  utilde = similar(u,axes(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u));
  Vern9Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,utilde,tmp,atmp,tab)
end

alg_cache(alg::Vern9,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Vern9ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))


text:
using Documenter # , ArbNumerics

makedocs(
    # modules = [ArbNumerics],
    sitename = "ArbNumerics",
    pages  = Any[
        "Overview"                 => "index.md",
        "The Types"                => "thetypes.md",
        "Setting Values"           => "settingvalues.md",
        "Seeing Values"            => "seeingvalues.md",
        "Precisions"               => "precisions.md",
        "Rounding"                 => "rounding.md",
        "Elementary Functions"     => "mathfunctions.md",
        "Elliptic Integrals"       => "ellipticintegrals.md",
        "Elliptic Functions"       => "ellipticfunctions.md",
        "Matrix Functions"         => "matrixfunctions.md",
        "Interval Functions"       => "intervalfunctions.md",
        "LibArb Differences"       => "libarb_differences.md",
        "References"               => "references.md",
        "Index"                    => "documentindex.md"
        ]
    )

deploydocs(
    repo = "github.com/JeffreySarnoff/ArbNumerics.jl.git",
    target = "build"
)


text:
# input: nz, nx, dz, dx, nSteps, nPoints_pml, nPad, dt, f0, survey_fname, data_dir_name, scratch_dir_name, isAc
using JSON
using DataStructures
using Dierckx

function paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml, nPad, para_fname, survey_fname, data_dir_name; if_win=false, filter_para=nothing, if_src_update=false, scratch_dir_name::String="")

  para = OrderedDict()
  para["nz"] = nz
	para["nx"] = nx
	para["dz"] = dz
	para["dx"] = dx
	para["nSteps"] = nSteps
	para["dt"] = dt
	para["f0"] = f0
	para["nPoints_pml"] = nPml
	para["nPad"] = nPad
  if if_win != false
    para["if_win"] = true
  end
  if filter_para != nothing
    para["filter"] = filter_para
  end
  if if_src_update != false
    para["if_src_update"] = true
  end
	para["survey_fname"] = survey_fname
  para["data_dir_name"] = data_dir_name
  if !isdir(data_dir_name)
    mkdir(data_dir_name)
  end
  # if nStepsWrap != nothing
  #   para["nStepsWrap"] = nStepsWrap
  # end

	if(scratch_dir_name != "")
      para["scratch_dir_name"] = scratch_dir_name
      if !isdir(scratch_dir_name)
        mkdir(scratch_dir_name)
      end
  end
  para_string = JSON.json(para)

  open(para_fname,"w") do f
    write(f, para_string)
  end
end

# all shots share the same number of receivers
function surveyGen(z_src, x_src, z_rec, x_rec, survey_fname; Windows=nothing, Weights=nothing)
  nsrc = length(x_src)
  nrec = length(x_rec)
  survey = OrderedDict()
  survey["nShots"] = nsrc
  for i = 1:nsrc
    shot = OrderedDict()
    shot["z_src"] = z_src[i]
    shot["x_src"] = x_src[i]
    shot["nrec"] = nrec
    shot["z_rec"] = z_rec
    shot["x_rec"] = x_rec
    if Windows != nothing
      shot["win_start"] = Windows["shot$(i-1)"][:start]
      shot["win_end"] = Windows["shot$(i-1)"][:end]
    end
    if Weights != nothing
      # shot["weights"] = Int64.(Weights["shot$(i-1)"][:weights])
      shot["weights"] = Weights["shot$(i-1)"][:weights]
    end
    survey["shot$(i-1)"] = shot
  end

  survey_string = JSON.json(survey)
  open(survey_fname,"w") do f
    write(f, survey_string)
  end

end

function sourceGene(f, nStep, delta_t)
#  Ricker wavelet generation and integration for source
#  Dongzhuo Li @ Stanford
#  May, 2015

  e = pi*pi*f*f;
  t_delay = 1.2/f;
  source = Matrix{Float64}(undef, 1, nStep)
  for it = 1:nStep
      source[it] = (1-2*e*(delta_t*(it-1)-t_delay)^2)*exp(-e*(delta_t*(it-1)-t_delay)^2);
  end

  for it = 2:nStep
      source[it] = source[it] + source[it-1];
  end
  source = source * delta_t;
end

"""
  cs_bounds_cloud(cpImg, Bounds)

Get `vs` high and low bounds from log point cloud
1st row of Bounds: vp ref line
2nd row of Bounds: vs high ref line
3rd row of Bounds: vs low ref line
"""
function cs_bounds_cloud(cpImg, Bounds)
  cs_high_itp = Spline1D(Bounds[1,:], Bounds[2,:]; k=1)
  cs_low_itp = Spline1D(Bounds[1,:], Bounds[3,:]; k=1)
  csHigh = zeros(size(cpImg))
  csLow = zeros(size(cpImg))
  for i = 1:size(cpImg, 1)
    for j = 1:size(cpImg, 2)
      csHigh[i,j] = cs_high_itp(cpImg[i,j])
      csLow[i,j] = cs_low_itp(cpImg[i,j])
    end
  end
  return csHigh, csLow
end

"""
  klauderWave(fmin, fmax, t_sweep, nStepTotal, nStepDelay, delta_t)

Generates Klauder wavelet.
"""
function klauderWave(fmin, fmax, t_sweep, nStepTotal, nStepDelay, delta_t)
  nStep = nStepTotal - nStepDelay
  source = Matrix{Float64}(undef, 1, nStep+nStep-1)
  source_half = Matrix{Float64}(undef, 1, nStep-1)
  K = (fmax - fmin) / t_sweep
  f0 = (fmin + fmax) / 2.0
  t_axis = delta_t:delta_t:(nStep-1)*delta_t
  source_half = sin.(pi * K .* t_axis .* (t_sweep .- t_axis)) .* cos.(2.0 * pi * f0 .* t_axis) ./ (pi*K.*t_axis*t_sweep)
  for i = 1:nStep-1
    source[i] = source_half[end-i+1]
  end
  for i = nStep+1:2*nStep-1
    source[i] = source_half[i-nStep]
  end
  source[nStep] = 1.0
  source_crop = source[:,nStep-nStepDelay:end]
  return source_crop
end

# function klauderWave(fmin, fmax, t_sweep, nStep, delta_t)
# #  Klauder wavelet
# #  Dongzhuo Li @ Stanford
# #  August, 2019
#   source = Matrix{Float64}(undef, 1, nStep)
#   K = (fmax - fmin) / t_sweep
#   f0 = (fmin + fmax) / 2.0
#   t_axis = delta_t:delta_t:(nStep-1)*delta_t
#   source_part = sin.(pi * K .* t_axis .* (t_sweep .- t_axis)) .* cos.(2.0 * pi * f0 .* t_axis) ./ (pi*K.*t_axis*t_sweep)
#   for i = 2:nStep
#     source[i] = source_part[i-1]
#   end
#   source[1] = 1.0
#   return source
# end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm
# using PyCall
# Io3D = @pyimport("io3d")

function randomImage(shape, structure, noiseFraction=0.0)
   """ Generation of random image of given shape and structure.
      Return scipy.ndarray(shape)
   """
   println("noiseFraction = $noiseFraction")
   ranges = [shape[k]/structure[k] for k=1:length(shape)]
   println("ranges = $ranges")
   random_array = rand(Int8, structure)
end

shape = [4,4]
structure = [8,8]
random_array = randomImage(shape, structure, 0.0)

function boundaryOps(skeletons)

	characteristicMatrices = map(Lar.characteristicMatrix,skeletons)
	couples = zip(characteristicMatrices[1:end-1], characteristicMatrices[2:end])
	operators = [ div(,2) for (A,B) in couples ]

end

function imageChainComplex (shape)
	_,skeletons = Lar.cuboidGrid(shape,true)

	imageLAR = (shape, skeletons, operators)
	dump(imageLAR,filename)
	print "filename =",filename
	return shape, skeletons, operators
end

def imageChainComplex (shape):
   tokens = str(shape)[1:-1].split(',')
   tokens = [token.strip() for token in tokens]
   filename = "tmp/larimage-" + "-".join(tokens) + ".pickle"

   if os.path.isfile(filename):
      shape, skeletons, operators = loadImageLAR(filename)
   else:
      skeletons = gridSkeletons(list(shape))
      operators = boundaryOps(skeletons)
      imageLAR = (shape, skeletons, operators)
      dump(imageLAR,filename)
      print "filename =",filename
   return shape, skeletons, operators


text:
module BluseBeamformerRecipes

using BeamformerRecipes
using FringeExplorer
using RadioInterferometry
using EarthOrientation
using LinearAlgebra
using Dates
using Redis
using Blio
using YAML

# Re-export all names from BeamformerRecipes
for name in names(BeamformerRecipes)
    @eval export $name
end

include("telinfo.jl")
include("cal_solution.jl")
include("guppiraw.jl")
include("redisguppi.jl")
include("targets2bfr.jl")

export run_targets_listener

"""
    BeamformerRecipe(grh, beam_names, beam_positions; redis, telinfo_file)

Create a BeamformerRecipe for the scan described by `grh` (a `GuppiRaw.Header`
object) with beams whose names are given by `beam_names` and whose positions are
given as `beam_positions`.  Ancillary info from `redis` and `telinfo_file` will
be used.

The following fields will be
populated:

- `DimInfo`
- `TelInfo`
- `ObsInfo`
- `CalInfo`
"""
function BeamformerRecipes.BeamformerRecipe(
    grh::GuppiRaw.Header, beam_names, beam_positions;
    redis::RedisConnection=RedisConnection(host="redishost"),
    telinfo_file::AbstractString=joinpath(ENV["HOME"], "telinfo.yml")
)
    # Get number of beams
    @assert length(beam_names) == size(beam_positions, 2) "beam names and positions size mismatch"
    nbeams = length(beam_names)

    subname = subarray_name(grh)
    tstart = starttime(grh)

    # Read cal_solution from redis
    nants, nchan, ants, cals = cal_solution(redis, subname, tstart)

    # Load telinfo and then filter/reorder according to nants
    telinfo=TelInfo(telinfo_file)
    # Get indices of `ants` in telinfo.antenna_names
    idxs = indexin(ants, telinfo.antenna_names)
    # Filter/reorder antenna fields in telinfo
    telinfo.antenna_positions = telinfo.antenna_positions[:,idxs]
    telinfo.antenna_names = telinfo.antenna_names[idxs]
    telinfo.antenna_numbers = telinfo.antenna_numbers[idxs]
    telinfo.antenna_diameters = telinfo.antenna_diameters[idxs]

    # Create frequency array
    freqs=collect(range(fch1(grh), step=foff(grh), length=nchan))

    # Time step is 1 second (for now)
    Δt = 1.0
    Δjd = Δt/86400

    # Get delays and delay rates each Δt for "DWELL" (or 300) seconds
    ntimes = ceil(Int, get(grh, "DWELL", 300) / Δt)
    delays = Array{Float64}(undef, nants, nbeams, ntimes)
    rates = Array{Float64}(undef, nants, nbeams, ntimes)
    times = collect(range(tstart, step=Δt, length=ntimes))
    jdstart = datetime2julian(unix2datetime(tstart))
    jds = collect(range(jdstart, step=Δjd, length=ntimes))

    # Use a single dut1 value for the entire scan.  Use value from RAW header,
    # if present, otherwise use value from EarthOrientation for nearest midnight
    # to jdstart.
    jdmid = floor(jdstart - 0.5) + 0.5
    dut1 = get(grh, "UT1_UTC", EarthOrientation.getΔUT1(jdmid))

    # Get boresight RA/Dec
    α = ra(grh)
    δ = dec(grh)

    # Add bore sight position to "front" of beam_positions
    beam_positions = hcat([α, δ], beam_positions)

    # Get hour angles and declinations for each beam at jdstart
    hdobs = radec2hadec(beam_positions, jdstart,
        latitude=telinfo.latitude,
        longitude=telinfo.longitude,
        altitude=telinfo.altitude,
        dut1=dut1, xp=0.0, yp=0.0
    )

    # Convert observed hours angles to "observed transit times" (transit at tobs=0)
    tobs = ha2t.(hdobs[1,:])
    dobs =       hdobs[2,:]

    # Create wdw_m and wdw matrices
    wdw_m = zeros(Float64, 2, 3)
    wdw = zeros(Float64, 2, nants)

    # Convert antpos into topocentric XYZ
    antpos = antpos_topo_xyz(telinfo)

    # Convert antpos into nanoseconds
    antpos ./= RadioInterferometry.CMPNS

    # For each time step
    for ti in 1:ntimes
        # For each beam
        for bi in 1:nbeams
            # Get dwd matrix for current time/beam
            td2wdw!(wdw_m, tobs[bi], dobs[bi])
            # Get delay and rate for each antenna
            mul!(wdw, wdw_m, antpos)
            # Store in delays and rates
            delays[:, bi, ti] .= wdw[1, :]
            rates[ :, bi, ti] .= wdw[2, :]
        end

        # Step tobs
        tobs .+= Δt
    end

    # Subtract boresight delay and rate from each antenna,beam
    delays = delays[:, 2:end, :] .- delays[:,1:1,:]
    rates  = rates[ :, 2:end, :] .- rates[ :,1:1,:]

    BeamformerRecipe(
        diminfo=DimInfo(
            nants=nants,
            npol=2,
            nchan=nchan,
            nbeams=nbeams,
            ntimes=ntimes),
        telinfo=telinfo,
        obsinfo=ObsInfo(
            obsid=obsid(grh, telinfo.telescope_name),
            freq_array=freqs,
            phase_center_ra=α,
            phase_center_dec=δ,
            instrument_name=instrument(grh)
        ),
        calinfo=CalInfo(; cals...),
        beaminfo=BeamInfo(
            ras = beam_positions[1, 2:end], # Skip boresight position
            decs= beam_positions[2, 2:end], # Skip boresight position
            src_names=beam_names
        ),
        delayinfo=DelayInfo(
            delays=delays,
            rates=rates,
            time_array=times,
            jds=jds,
            dut1=dut1
        )
    )
end

"""
    BeamformerRecipes.BeamformerRecipe(rawname; redis, telinfo_file, ring_arcsec)

Create a BeamformerRecipe for the scan corresponding to GUPPI RAW file `rawname`
using ancillary info from `redis` and `telinfo_file`.  One beam will be at bore
sight and 60 other beams will be arranged in four concentric rings with radii
`1:4 .* ring_arcsec`.  The following fields will be populated:

- `DimInfo` (all but `nbeams` and `ntimes`)
- `TelInfo`
- `ObsInfo`
- `CalInfo`
"""
function BeamformerRecipes.BeamformerRecipe(
    rawname::AbstractString;
    redis::RedisConnection=RedisConnection(host="redishost"),
    telinfo_file::AbstractString=joinpath(ENV["HOME"], "telinfo.yml"),
    ring_arcsec=10
)
    # Read GuppiRaw.Header from `rawname`
    @info "reading header from $rawname"
    grh = open(io->read(io, GuppiRaw.Header), rawname)

    # Get boresight RA/Dec
    α = ra(grh)
    δ = dec(grh)

    # Get beam positions
    nrings = 4
    beam_positions = beamrings(α, δ, nrings=nrings, dϕ=deg2rad(ring_arcsec/3600))

    # Get beam names
    src_name = get(grh, "SRC_NAME", "UNKNOWN")
    beam_names = ["$(src_name)_R$(r)B$(b)" for r=1:nrings for b=1:6r]
    pushfirst!(beam_names, src_name)

    BeamformerRecipe(grh, beam_names, beam_positions; redis, telinfo_file)
end

end


text:
using NearestNeighbors
include("./preprocess.jl")

dimensions = 3
numOfPoints = 1<<19
numOfQueries = 1<<19
Points = rand(Float32, numOfPoints ,dimensions)
Queries = rand(Float32, numOfQueries, dimensions)

numOfGrids = 1<<4
BlockOfPoint = AssignPointsToBlock(Points, numOfGrids, dimensions)
BlockOfQuery = AssignPointsToBlock(Queries, numOfGrids, dimensions)

PointsPerBlock, IntegralPointsPerBlock = CountPointsPerBlock(Points, numOfGrids, dimensions)
QueriesPerBlock, IntegralQueriesPerBlock = CountPointsPerBlock(Queries, numOfGrids, dimensions)

OrderedPoints = ReorderPointsByBlock(Points, BlockOfPoint)
OrderedQueries = ReorderPointsByBlock(Queries, BlockOfQuery)

############################################################
# Uncomment only the implementation you want to validate.  #
############################################################
include("kernel_simple.jl")

println("Running for simple version: ")
gpu_idxs, gpu_dists = cuda_knn_simple(OrderedPoints, OrderedQueries,PointsPerBlock,
    QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
    numOfQueries, numOfGrids, dimensions)
##############################################
# include("kernel_simple_view_function.jl")
#
# println("Running for simple version with @view and distance calcs in function: ")
# gpu_idxs, gpu_dists = cuda_knn_simple_view_function(OrderedPoints, OrderedQueries,PointsPerBlock,
#     QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
#     numOfQueries, numOfGrids, dimensions)
###############################################
# include("kernel_with_skip.jl")
#
# println("Running for version with skip: ")
# gpu_idxs, gpu_dists = cuda_knn_with_skip(OrderedPoints, OrderedQueries,PointsPerBlock,
#     QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
#     numOfQueries, numOfGrids, dimensions)
###############################################
# include("kernel_with_skip_view_function.jl")
# println("Running for version with skip and @view and distance calcs in function: ")
# gpu_idxs, gpu_dists = cuda_knn_with_skip_view_function(OrderedPoints, OrderedQueries,PointsPerBlock,
#     QueriesPerBlock, IntegralPointsPerBlock, IntegralQueriesPerBlock,numOfPoints,
#     numOfQueries, numOfGrids, dimensions)

tr = NearestNeighbors.KDTree(OrderedPoints')
cpu_idxs, cpu_dists = NearestNeighbors.knn(tr, OrderedQueries', 1)

cpu_idxs = vcat(cpu_idxs...)

for i in 1:length(cpu_idxs)
    res = gpu_idxs[i]
    if cpu_idxs[i] != res
        println(i , " : " ,cpu_idxs[i], "     ", res)
        println("\tblock pt", BlockOfPoint[cpu_idxs[i]], " wrong pt ", BlockOfPoint[ res] ," block query: ", BlockOfQuery[i])
        end
end


text:
#lmmdata.jl

struct LMMData{T <: AbstractFloat}
    # Fixed effect matrix
    xv::Matrix{T}
    # Responce vector
    yv::Vector{T}
    function LMMData(xa::Matrix{T}, ya::Vector{T}) where T <: AbstractFloat
        new{T}(xa, ya)
    end
end

struct LMMDataViews{T1, T2} <: AbstractLMMDataBlocks
    # Fixed effect matrix views
    xv::T1
    # Responce vector views
    yv::T2
    function LMMDataViews(xv, yv, vcovblock)
        x1 = view(xv, vcovblock[1],:)
        y1 = view(yv, vcovblock[1])
        x = Vector{typeof(x1)}(undef, length(vcovblock))
        y = Vector{typeof(y1)}(undef, length(vcovblock))
        x[1] = x1
        y[1] = y1
        for i = 2:length(vcovblock)
            x[i] = view(xv, vcovblock[i],:)
            y[i] = view(yv, vcovblock[i])
        end
        new{typeof(x), typeof(y)}(x, y)
    end
    function LMMDataViews(lmm)
        return LMMDataViews(lmm.data.xv, lmm.data.yv, lmm.covstr.vcovblock)
    end
end


text:
module Included
using Hijack, Test

@testset "included testset" for _=1:1
    @test true
    push!(Hijack.RUN, 2)
    # test that `testset=true` is forwarded
    include("included_testset2.jl")
end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# caller_name.jl

using MIRTjim: caller_name

using Test: @test, @testset

@testset "caller_name" begin

    function f2()
        caller_name()
    end

    line = 2 + @__LINE__
    function f1()
        f2() # this is two lines below @__LINE__ above
    end

    @test isa(f1(), String)
    @test f1()[end-12:end] == ".jl $line f1(): "
end