text:
# Functions supporting trendline identification (support/resistance, zigzag, elliot waves, etc.)

function maxima{Float64}(x::Vector{Float64}; threshold::Float64=0.0, order::Int=1)
    @assert threshold >= 0.0 "Threshold must be positive."
    @assert order > 0 "Order must be a positive integer."
    n = size(x,1)
    crit = falses(n)
    @inbounds for i=2:n-1
        if (x[i]-x[i-1] >= threshold) && (x[i]-x[i+1] >= threshold)
            crit[i] = true
        end
    end
    while order > 1
        idx = find(crit)
        crit[idx[!maxima(x[crit], threshold=threshold)]] = false
        order -= 1
    end
    return crit
end

function minima{Float64}(x::Vector{Float64}; threshold::Float64=0.0, order::Int=1)
    @assert threshold <= 0.0 "Threshold must be negative."
    @assert order > 0 "Order must be a positive integer."
    n = size(x,1)
    crit = falses(n)
    @inbounds for i=2:n-1
        if (x[i]-x[i-1] <= threshold) && (x[i]-x[i+1] <= threshold)
            crit[i] = true
        end
    end
    while order > 1
        idx = find(crit)
        crit[idx[!minima(x[crit], threshold=threshold)]] = false
        order -= 1
    end
    return crit
end

function interpolate(x1::Int, x2::Int, y1::Float64, y2::Float64)
	m = (y2-y1)/(x2-x1)
	b = y1 - m*x1
	x = collect(x1:1.0:x2)
	y = m*x + b
	return y
end

function resistance{Float64}(x::Vector{Float64}; order::Int=1, threshold::Float64=0.0)
    out = zeros(x)
    crit = maxima(x, threshold=threshold, order=order)
    out[!crit] = NaN
    idx = find(crit)
    @inbounds for i=2:length(idx)
        out[idx[i-1]:idx[i]] = interpolate(idx[i-1], idx[i], x[i-1], x[i])
    end
    return out
end

function support{Float64}(x::Vector{Float64}; order::Int=1, threshold::Float64=0.0)
    out = zeros(x)
    crit = minima(x, threshold=threshold, order=order)
    out[!crit] = NaN
    idx = find(crit)
    @inbounds for i=2:length(idx)
        out[idx[i-1]:idx[i]] = interpolate(idx[i-1], idx[i], x[i-1], x[i])
    end
    return out
end


text:
using Revise
using VCFTools
using MendelImpute
using GeneticVariation
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLD2, FileIO, JLSO
using BenchmarkTools
using GroupSlices
using LinearAlgebra
using Profile

chr = 20
maf = 0.01
width = 64

BLAS.set_num_threads(1)
Threads.nthreads() != 10 && error("not 10 threads!")
Profile.init(n = 10^7, delay = 0.01)

Random.seed!(2020)
tgtfile = "target.chr$chr.typedOnly.maf$maf.masked.vcf.gz"
reffile = "./single_precision/ref.chr$chr.w$width.maf$maf.excludeTarget.jlso"
outfile = "mendel.imputed.dp$width.maf$maf.vcf.gz"

phase(tgtfile, reffile, outfile=outfile, impute=true, width=width,
    dynamic_programming = false)
Profile.clear()
@profile phase(tgtfile, reffile, outfile=outfile, impute=true, width=width,
    dynamic_programming = false)

open("profile_result.txt", "w") do s
    Profile.print(IOContext(s, :displaysize => (24, 500)), format=:flat)
end

rm(outfile, force=true)

ext: jl

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using HTTP
using BufferedStreams

function load_first_datapoint(filename; set_phase=:real)
    data = load_Data(filename)
    M = data[1,:]

    if set_phase == :real
        phase = angle(M[end])
        M = M .* exp(-1im * phase)
        M = real.(M)
    elseif set_phase == :abs
        M = abs.(M)
    end
    return M
end

function load_Data(filename)
    if filename[1:4] == "http"
        io = Base.BufferStream()
        mytask = @async readdata(io)
        HTTP.request("GET", filename, response_stream=io)
        data = fetch(mytask)
    else
        filename = expanduser(filename)
        io = open(filename)
        data = readdata(io)
    end

    close(io)
    return data
end

function readdata(io)
    for _=1:12
        Char.(read(io, Char))
    end
    if read(io, Int32) != 501; error(); end # 500 real; 501 complex; 502 double real; 503 xy_real; 504 xy_complex

    size1 = read(io, Int32)
    size2 = read(io, Int32)
    if read(io, Int32) != 1; error(); end
    if read(io, Int32) != 1; error(); end

    data = zeros(ComplexF32, size1, size2)
    for i = 1:length(data)
        data[i] = read(io, Float32) + 1im * read(io, Float32)
    end

    return data[7:end - 7,:]
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    FractionPartition(fraction, shuffle=true)

A method for partitioning spatial objects according to a given `fraction`.
Optionally `shuffle` elements before partitioning.
"""
struct FractionPartition <: PartitionMethod
  fraction::Float64
  shuffle::Bool

  function FractionPartition(fraction, shuffle)
    @assert 0 < fraction < 1 "fraction must be in interval (0,1)"
    new(fraction, shuffle)
  end
end

FractionPartition(fraction) = FractionPartition(fraction, true)

function partition(object, method::FractionPartition)
  n = nelms(object)
  f = round(Int, method.fraction * n)

  locs = method.shuffle ? randperm(n) : 1:n
  subsets = [locs[1:f], locs[f+1:n]]

  SpatialPartition(object, subsets)
end


text:
@symbol_func function K_kappa(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n

  cur_nu_T = cur_reactor.nu_T

  cur_K = 0.4744

  cur_K *= ( 1 + cur_reactor.f_D )

  cur_K *= ( 1 + cur_nu_n )

  cur_K *= ( 1 + cur_nu_T )

  cur_K /= 1 + cur_nu_n + cur_nu_T

  cur_K *= cur_reactor.epsilon ^ 2

  cur_K *= cur_reactor.kappa_95

  cur_K *= g(cur_reactor)

  cur_K
end


text:
#==============================================================================#
# ACM.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/acm-2015-12-08.normal.json
#==============================================================================#

__precompile__()

module ACM

using AWSCore

"""
    using AWSSDK.ACM.add_tags_to_certificate
    add_tags_to_certificate([::AWSConfig], arguments::Dict)
    add_tags_to_certificate([::AWSConfig]; CertificateArn=, Tags=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "AddTagsToCertificate", arguments::Dict)
    acm([::AWSConfig], "AddTagsToCertificate", CertificateArn=, Tags=)

# AddTagsToCertificate Operation

Adds one or more tags to an ACM certificate. Tags are labels that you can use to identify and organize your AWS resources. Each tag consists of a `key` and an optional `value`. You specify the certificate on input by its Amazon Resource Name (ARN). You specify the tag by using a key-value pair.

You can apply a tag to just one certificate if you want to identify a specific characteristic of that certificate, or you can apply the same tag to multiple certificates if you want to filter for a common relationship among those certificates. Similarly, you can apply the same tag to multiple resources if you want to specify a relationship among those resources. For example, you can add the same tag to an ACM certificate and an Elastic Load Balancing load balancer to indicate that they are both used by the same website. For more information, see [Tagging ACM certificates](http://docs.aws.amazon.com/acm/latest/userguide/tags.html).

To remove one or more tags, use the [RemoveTagsFromCertificate](@ref) action. To view all of the tags that have been applied to the certificate, use the [ListTagsForCertificate](@ref) action.

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains the ARN of the ACM certificate to which the tag is to be applied. This must be of the form:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

## `Tags = [[ ... ], ...]` -- *Required*
The key-value pair that defines the tag. The tag value is optional.
```
 Tags = [[
        "Key" => <required> ::String,
        "Value" =>  ::String
    ], ...]
```

# Exceptions

`ResourceNotFoundException`, `InvalidArnException`, `InvalidTagException` or `TooManyTagsException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/AddTagsToCertificate)
"""
@inline add_tags_to_certificate(aws::AWSConfig=default_aws_config(); args...) = add_tags_to_certificate(aws, args)

@inline add_tags_to_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "AddTagsToCertificate", args)

@inline add_tags_to_certificate(args) = add_tags_to_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.delete_certificate
    delete_certificate([::AWSConfig], arguments::Dict)
    delete_certificate([::AWSConfig]; CertificateArn=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "DeleteCertificate", arguments::Dict)
    acm([::AWSConfig], "DeleteCertificate", CertificateArn=)

# DeleteCertificate Operation

Deletes a certificate and its associated private key. If this action succeeds, the certificate no longer appears in the list that can be displayed by calling the [ListCertificates](@ref) action or be retrieved by calling the [GetCertificate](@ref) action. The certificate will not be available for use by AWS services integrated with ACM.

**Note**
> You cannot delete an ACM certificate that is being used by another AWS service. To delete a certificate that is in use, the certificate association must first be removed.

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains the ARN of the ACM certificate to be deleted. This must be of the form:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

# Exceptions

`ResourceNotFoundException`, `ResourceInUseException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/DeleteCertificate)
"""
@inline delete_certificate(aws::AWSConfig=default_aws_config(); args...) = delete_certificate(aws, args)

@inline delete_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "DeleteCertificate", args)

@inline delete_certificate(args) = delete_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.describe_certificate
    describe_certificate([::AWSConfig], arguments::Dict)
    describe_certificate([::AWSConfig]; CertificateArn=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "DescribeCertificate", arguments::Dict)
    acm([::AWSConfig], "DescribeCertificate", CertificateArn=)

# DescribeCertificate Operation

Returns detailed metadata about the specified ACM certificate.

# Arguments

## `CertificateArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the ACM certificate. The ARN must have the following form:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

# Returns

`DescribeCertificateResponse`

# Exceptions

`ResourceNotFoundException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/DescribeCertificate)
"""
@inline describe_certificate(aws::AWSConfig=default_aws_config(); args...) = describe_certificate(aws, args)

@inline describe_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "DescribeCertificate", args)

@inline describe_certificate(args) = describe_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.export_certificate
    export_certificate([::AWSConfig], arguments::Dict)
    export_certificate([::AWSConfig]; CertificateArn=, Passphrase=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "ExportCertificate", arguments::Dict)
    acm([::AWSConfig], "ExportCertificate", CertificateArn=, Passphrase=)

# ExportCertificate Operation

Exports a private certificate issued by a private certificate authority (CA) for use anywhere. You can export the certificate, the certificate chain, and the encrypted private key associated with the public key embedded in the certificate. You must store the private key securely. The private key is a 2048 bit RSA key. You must provide a passphrase for the private key when exporting it. You can use the following OpenSSL command to decrypt it later. Provide the passphrase when prompted.

`openssl rsa -in encrypted_key.pem -out decrypted_key.pem`

# Arguments

## `CertificateArn = ::String` -- *Required*
An Amazon Resource Name (ARN) of the issued certificate. This must be of the form:

`arn:aws:acm:region:account:certificate/12345678-1234-1234-1234-123456789012`

## `Passphrase = blob` -- *Required*
Passphrase to associate with the encrypted exported private key. If you want to later decrypt the private key, you must have the passphrase. You can use the following OpenSSL command to decrypt a private key:

`openssl rsa -in encrypted_key.pem -out decrypted_key.pem`

# Returns

`ExportCertificateResponse`

# Exceptions

`ResourceNotFoundException`, `RequestInProgressException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/ExportCertificate)
"""
@inline export_certificate(aws::AWSConfig=default_aws_config(); args...) = export_certificate(aws, args)

@inline export_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "ExportCertificate", args)

@inline export_certificate(args) = export_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.get_certificate
    get_certificate([::AWSConfig], arguments::Dict)
    get_certificate([::AWSConfig]; CertificateArn=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "GetCertificate", arguments::Dict)
    acm([::AWSConfig], "GetCertificate", CertificateArn=)

# GetCertificate Operation

Retrieves a certificate specified by an ARN and its certificate chain . The chain is an ordered list of certificates that contains the end entity certificate, intermediate certificates of subordinate CAs, and the root certificate in that order. The certificate and certificate chain are base64 encoded. If you want to decode the certificate to see the individual fields, you can use OpenSSL.

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains a certificate ARN in the following format:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

# Returns

`GetCertificateResponse`

# Exceptions

`ResourceNotFoundException`, `RequestInProgressException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/GetCertificate)
"""
@inline get_certificate(aws::AWSConfig=default_aws_config(); args...) = get_certificate(aws, args)

@inline get_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "GetCertificate", args)

@inline get_certificate(args) = get_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.import_certificate
    import_certificate([::AWSConfig], arguments::Dict)
    import_certificate([::AWSConfig]; Certificate=, PrivateKey=, <keyword arguments>)

    using AWSCore.Services.acm
    acm([::AWSConfig], "ImportCertificate", arguments::Dict)
    acm([::AWSConfig], "ImportCertificate", Certificate=, PrivateKey=, <keyword arguments>)

# ImportCertificate Operation

Imports a certificate into AWS Certificate Manager (ACM) to use with services that are integrated with ACM. Note that [integrated services](http://docs.aws.amazon.com/acm/latest/userguide/acm-services.html) allow only certificate types and keys they support to be associated with their resources. Further, their support differs depending on whether the certificate is imported into IAM or into ACM. For more information, see the documentation for each service. For more information about importing certificates into ACM, see [Importing Certificates](http://docs.aws.amazon.com/acm/latest/userguide/import-certificate.html) in the *AWS Certificate Manager User Guide*.

**Note**
> ACM does not provide [managed renewal](http://docs.aws.amazon.com/acm/latest/userguide/acm-renewal.html) for certificates that you import.

Note the following guidelines when importing third party certificates:

*   You must enter the private key that matches the certificate you are importing.

*   The private key must be unencrypted. You cannot import a private key that is protected by a password or a passphrase.

*   If the certificate you are importing is not self-signed, you must enter its certificate chain.

*   If a certificate chain is included, the issuer must be the subject of one of the certificates in the chain.

*   The certificate, private key, and certificate chain must be PEM-encoded.

*   The current time must be between the `Not Before` and `Not After` certificate fields.

*   The `Issuer` field must not be empty.

*   The OCSP authority URL, if present, must not exceed 1000 characters.

*   To import a new certificate, omit the `CertificateArn` argument. Include this argument only when you want to replace a previously imported certificate.

*   When you import a certificate by using the CLI, you must specify the certificate, the certificate chain, and the private key by their file names preceded by `file://`. For example, you can specify a certificate saved in the `C:\\temp` folder as `file://C:\\temp\\certificate_to_import.pem`. If you are making an HTTP or HTTPS Query request, include these arguments as BLOBs.

*   When you import a certificate by using an SDK, you must specify the certificate, the certificate chain, and the private key files in the manner required by the programming language you're using.

This operation returns the [Amazon Resource Name (ARN)](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html) of the imported certificate.

# Arguments

## `CertificateArn = ::String`
The [Amazon Resource Name (ARN)](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html) of an imported certificate to replace. To import a new certificate, omit this field.

## `Certificate = blob` -- *Required*
The certificate to import.

## `PrivateKey = blob` -- *Required*
The private key that matches the public key in the certificate.

## `CertificateChain = blob`
The PEM encoded certificate chain.

# Returns

`ImportCertificateResponse`

# Exceptions

`ResourceNotFoundException` or `LimitExceededException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/ImportCertificate)
"""
@inline import_certificate(aws::AWSConfig=default_aws_config(); args...) = import_certificate(aws, args)

@inline import_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "ImportCertificate", args)

@inline import_certificate(args) = import_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.list_certificates
    list_certificates([::AWSConfig], arguments::Dict)
    list_certificates([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.acm
    acm([::AWSConfig], "ListCertificates", arguments::Dict)
    acm([::AWSConfig], "ListCertificates", <keyword arguments>)

# ListCertificates Operation

Retrieves a list of certificate ARNs and domain names. You can request that only certificates that match a specific status be listed. You can also filter by specific attributes of the certificate.

# Arguments

## `CertificateStatuses = ["PENDING_VALIDATION", "ISSUED", "INACTIVE", "EXPIRED", "VALIDATION_TIMED_OUT", "REVOKED" or "FAILED", ...]`
Filter the certificate list by status value.

## `Includes = [ ... ]`
Filter the certificate list. For more information, see the [Filters](@ref) structure.
```
 Includes = [
        "extendedKeyUsage" =>  ["TLS_WEB_SERVER_AUTHENTICATION", "TLS_WEB_CLIENT_AUTHENTICATION", "CODE_SIGNING", "EMAIL_PROTECTION", "TIME_STAMPING", "OCSP_SIGNING", "IPSEC_END_SYSTEM", "IPSEC_TUNNEL", "IPSEC_USER", "ANY", "NONE" or "CUSTOM", ...],
        "keyUsage" =>  ["DIGITAL_SIGNATURE", "NON_REPUDIATION", "KEY_ENCIPHERMENT", "DATA_ENCIPHERMENT", "KEY_AGREEMENT", "CERTIFICATE_SIGNING", "CRL_SIGNING", "ENCIPHER_ONLY", "DECIPHER_ONLY", "ANY" or "CUSTOM", ...],
        "keyTypes" =>  ["RSA_2048", "RSA_1024", "RSA_4096", "EC_prime256v1", "EC_secp384r1" or "EC_secp521r1", ...]
    ]
```

## `NextToken = ::String`
Use this parameter only when paginating results and only in a subsequent request after you receive a response with truncated results. Set it to the value of `NextToken` from the response you just received.

## `MaxItems = ::Int`
Use this parameter when paginating results to specify the maximum number of items to return in the response. If additional items exist beyond the number you specify, the `NextToken` element is sent in the response. Use this `NextToken` value in a subsequent request to retrieve additional items.

# Returns

`ListCertificatesResponse`

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/ListCertificates)
"""
@inline list_certificates(aws::AWSConfig=default_aws_config(); args...) = list_certificates(aws, args)

@inline list_certificates(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "ListCertificates", args)

@inline list_certificates(args) = list_certificates(default_aws_config(), args)

"""
    using AWSSDK.ACM.list_tags_for_certificate
    list_tags_for_certificate([::AWSConfig], arguments::Dict)
    list_tags_for_certificate([::AWSConfig]; CertificateArn=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "ListTagsForCertificate", arguments::Dict)
    acm([::AWSConfig], "ListTagsForCertificate", CertificateArn=)

# ListTagsForCertificate Operation

Lists the tags that have been applied to the ACM certificate. Use the certificate's Amazon Resource Name (ARN) to specify the certificate. To add a tag to an ACM certificate, use the [AddTagsToCertificate](@ref) action. To delete a tag, use the [RemoveTagsFromCertificate](@ref) action.

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains the ARN of the ACM certificate for which you want to list the tags. This must have the following form:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

# Returns

`ListTagsForCertificateResponse`

# Exceptions

`ResourceNotFoundException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/ListTagsForCertificate)
"""
@inline list_tags_for_certificate(aws::AWSConfig=default_aws_config(); args...) = list_tags_for_certificate(aws, args)

@inline list_tags_for_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "ListTagsForCertificate", args)

@inline list_tags_for_certificate(args) = list_tags_for_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.remove_tags_from_certificate
    remove_tags_from_certificate([::AWSConfig], arguments::Dict)
    remove_tags_from_certificate([::AWSConfig]; CertificateArn=, Tags=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "RemoveTagsFromCertificate", arguments::Dict)
    acm([::AWSConfig], "RemoveTagsFromCertificate", CertificateArn=, Tags=)

# RemoveTagsFromCertificate Operation

Remove one or more tags from an ACM certificate. A tag consists of a key-value pair. If you do not specify the value portion of the tag when calling this function, the tag will be removed regardless of value. If you specify a value, the tag is removed only if it is associated with the specified value.

To add tags to a certificate, use the [AddTagsToCertificate](@ref) action. To view all of the tags that have been applied to a specific ACM certificate, use the [ListTagsForCertificate](@ref) action.

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains the ARN of the ACM Certificate with one or more tags that you want to remove. This must be of the form:

`arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012`

For more information about ARNs, see [Amazon Resource Names (ARNs) and AWS Service Namespaces](http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html).

## `Tags = [[ ... ], ...]` -- *Required*
The key-value pair that defines the tag to remove.
```
 Tags = [[
        "Key" => <required> ::String,
        "Value" =>  ::String
    ], ...]
```

# Exceptions

`ResourceNotFoundException`, `InvalidArnException` or `InvalidTagException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/RemoveTagsFromCertificate)
"""
@inline remove_tags_from_certificate(aws::AWSConfig=default_aws_config(); args...) = remove_tags_from_certificate(aws, args)

@inline remove_tags_from_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "RemoveTagsFromCertificate", args)

@inline remove_tags_from_certificate(args) = remove_tags_from_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.request_certificate
    request_certificate([::AWSConfig], arguments::Dict)
    request_certificate([::AWSConfig]; DomainName=, <keyword arguments>)

    using AWSCore.Services.acm
    acm([::AWSConfig], "RequestCertificate", arguments::Dict)
    acm([::AWSConfig], "RequestCertificate", DomainName=, <keyword arguments>)

# RequestCertificate Operation

Requests an ACM certificate for use with other AWS services. To request an ACM certificate, you must specify a fully qualified domain name (FQDN) in the `DomainName` parameter. You can also specify additional FQDNs in the `SubjectAlternativeNames` parameter.

If you are requesting a private certificate, domain validation is not required. If you are requesting a public certificate, each domain name that you specify must be validated to verify that you own or control the domain. You can use [DNS validation](http://docs.aws.amazon.com/acm/latest/userguide/gs-acm-validate-dns.html) or [email validation](http://docs.aws.amazon.com/acm/latest/userguide/gs-acm-validate-email.html). We recommend that you use DNS validation. ACM issues public certificates after receiving approval from the domain owner.

# Arguments

## `DomainName = ::String` -- *Required*
Fully qualified domain name (FQDN), such as www.example.com, that you want to secure with an ACM certificate. Use an asterisk (*) to create a wildcard certificate that protects several sites in the same domain. For example, *.example.com protects www.example.com, site.example.com, and images.example.com.

The first domain name you enter cannot exceed 63 octets, including periods. Each subsequent Subject Alternative Name (SAN), however, can be up to 253 octets in length.

## `ValidationMethod = "EMAIL" or "DNS"`
The method you want to use if you are requesting a public certificate to validate that you own or control domain. You can [validate with DNS](http://docs.aws.amazon.com/acm/latest/userguide/gs-acm-validate-dns.html) or [validate with email](http://docs.aws.amazon.com/acm/latest/userguide/gs-acm-validate-email.html). We recommend that you use DNS validation.

## `SubjectAlternativeNames = [::String, ...]`
Additional FQDNs to be included in the Subject Alternative Name extension of the ACM certificate. For example, add the name www.example.net to a certificate for which the `DomainName` field is www.example.com if users can reach your site by using either name. The maximum number of domain names that you can add to an ACM certificate is 100. However, the initial limit is 10 domain names. If you need more than 10 names, you must request a limit increase. For more information, see [Limits](http://docs.aws.amazon.com/acm/latest/userguide/acm-limits.html).

The maximum length of a SAN DNS name is 253 octets. The name is made up of multiple labels separated by periods. No label can be longer than 63 octets. Consider the following examples:

*   `(63 octets).(63 octets).(63 octets).(61 octets)` is legal because the total length is 253 octets (63+1+63+1+63+1+61) and no label exceeds 63 octets.

*   `(64 octets).(63 octets).(63 octets).(61 octets)` is not legal because the total length exceeds 253 octets (64+1+63+1+63+1+61) and the first label exceeds 63 octets.

*   `(63 octets).(63 octets).(63 octets).(62 octets)` is not legal because the total length of the DNS name (63+1+63+1+63+1+62) exceeds 253 octets.

## `IdempotencyToken = ::String`
Customer chosen string that can be used to distinguish between calls to `RequestCertificate`. Idempotency tokens time out after one hour. Therefore, if you call `RequestCertificate` multiple times with the same idempotency token within one hour, ACM recognizes that you are requesting only one certificate and will issue only one. If you change the idempotency token for each call, ACM recognizes that you are requesting multiple certificates.

## `DomainValidationOptions = [[ ... ], ...]`
The domain name that you want ACM to use to send you emails so that you can validate domain ownership.
```
 DomainValidationOptions = [[
        "DomainName" => <required> ::String,
        "ValidationDomain" => <required> ::String
    ], ...]
```

## `Options = ["CertificateTransparencyLoggingPreference" =>  "ENABLED" or "DISABLED"]`
Currently, you can use this parameter to specify whether to add the certificate to a certificate transparency log. Certificate transparency makes it possible to detect SSL/TLS certificates that have been mistakenly or maliciously issued. Certificates that have not been logged typically produce an error message in a browser. For more information, see [Opting Out of Certificate Transparency Logging](http://docs.aws.amazon.com/acm/latest/userguide/acm-bestpractices.html#best-practices-transparency).

## `CertificateAuthorityArn = ::String`
The Amazon Resource Name (ARN) of the private certificate authority (CA) that will be used to issue the certificate. If you do not provide an ARN and you are trying to request a private certificate, ACM will attempt to issue a public certificate. For more information about private CAs, see the [AWS Certificate Manager Private Certificate Authority (PCA)](http://docs.aws.amazon.com/acm-pca/latest/userguide/PcaWelcome.html) user guide. The ARN must have the following form:

`arn:aws:acm-pca:region:account:certificate-authority/12345678-1234-1234-1234-123456789012`

# Returns

`RequestCertificateResponse`

# Exceptions

`LimitExceededException`, `InvalidDomainValidationOptionsException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/RequestCertificate)
"""
@inline request_certificate(aws::AWSConfig=default_aws_config(); args...) = request_certificate(aws, args)

@inline request_certificate(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "RequestCertificate", args)

@inline request_certificate(args) = request_certificate(default_aws_config(), args)

"""
    using AWSSDK.ACM.resend_validation_email
    resend_validation_email([::AWSConfig], arguments::Dict)
    resend_validation_email([::AWSConfig]; CertificateArn=, Domain=, ValidationDomain=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "ResendValidationEmail", arguments::Dict)
    acm([::AWSConfig], "ResendValidationEmail", CertificateArn=, Domain=, ValidationDomain=)

# ResendValidationEmail Operation

Resends the email that requests domain ownership validation. The domain owner or an authorized representative must approve the ACM certificate before it can be issued. The certificate can be approved by clicking a link in the mail to navigate to the Amazon certificate approval website and then clicking **I Approve**. However, the validation email can be blocked by spam filters. Therefore, if you do not receive the original mail, you can request that the mail be resent within 72 hours of requesting the ACM certificate. If more than 72 hours have elapsed since your original request or since your last attempt to resend validation mail, you must request a new certificate. For more information about setting up your contact email addresses, see [Configure Email for your Domain](http://docs.aws.amazon.com/acm/latest/userguide/setup-email.html).

# Arguments

## `CertificateArn = ::String` -- *Required*
String that contains the ARN of the requested certificate. The certificate ARN is generated and returned by the [RequestCertificate](@ref) action as soon as the request is made. By default, using this parameter causes email to be sent to all top-level domains you specified in the certificate request. The ARN must be of the form:

`arn:aws:acm:us-east-1:123456789012:certificate/12345678-1234-1234-1234-123456789012`

## `Domain = ::String` -- *Required*
The fully qualified domain name (FQDN) of the certificate that needs to be validated.

## `ValidationDomain = ::String` -- *Required*
The base validation domain that will act as the suffix of the email addresses that are used to send the emails. This must be the same as the `Domain` value or a superdomain of the `Domain` value. For example, if you requested a certificate for `site.subdomain.example.com` and specify a **ValidationDomain** of `subdomain.example.com`, ACM sends email to the domain registrant, technical contact, and administrative contact in WHOIS and the following five addresses:

*   admin@subdomain.example.com

*   administrator@subdomain.example.com

*   hostmaster@subdomain.example.com

*   postmaster@subdomain.example.com

*   webmaster@subdomain.example.com

# Exceptions

`ResourceNotFoundException`, `InvalidStateException`, `InvalidArnException` or `InvalidDomainValidationOptionsException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/ResendValidationEmail)
"""
@inline resend_validation_email(aws::AWSConfig=default_aws_config(); args...) = resend_validation_email(aws, args)

@inline resend_validation_email(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "ResendValidationEmail", args)

@inline resend_validation_email(args) = resend_validation_email(default_aws_config(), args)

"""
    using AWSSDK.ACM.update_certificate_options
    update_certificate_options([::AWSConfig], arguments::Dict)
    update_certificate_options([::AWSConfig]; CertificateArn=, Options=)

    using AWSCore.Services.acm
    acm([::AWSConfig], "UpdateCertificateOptions", arguments::Dict)
    acm([::AWSConfig], "UpdateCertificateOptions", CertificateArn=, Options=)

# UpdateCertificateOptions Operation

Updates a certificate. Currently, you can use this function to specify whether to opt in to or out of recording your certificate in a certificate transparency log. For more information, see [Opting Out of Certificate Transparency Logging](http://docs.aws.amazon.com/acm/latest/userguide/acm-bestpractices.html#best-practices-transparency).

# Arguments

## `CertificateArn = ::String` -- *Required*
ARN of the requested certificate to update. This must be of the form:

`arn:aws:acm:us-east-1:*account*:certificate/*12345678-1234-1234-1234-123456789012*`

## `Options = ["CertificateTransparencyLoggingPreference" =>  "ENABLED" or "DISABLED"]` -- *Required*
Use to update the options for your certificate. Currently, you can specify whether to add your certificate to a transparency log. Certificate transparency makes it possible to detect SSL/TLS certificates that have been mistakenly or maliciously issued. Certificates that have not been logged typically produce an error message in a browser.

# Exceptions

`ResourceNotFoundException`, `LimitExceededException`, `InvalidStateException` or `InvalidArnException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/acm-2015-12-08/UpdateCertificateOptions)
"""
@inline update_certificate_options(aws::AWSConfig=default_aws_config(); args...) = update_certificate_options(aws, args)

@inline update_certificate_options(aws::AWSConfig, args) = AWSCore.Services.acm(aws, "UpdateCertificateOptions", args)

@inline update_certificate_options(args) = update_certificate_options(default_aws_config(), args)

end # module ACM

#==============================================================================#
# End of file
#==============================================================================#


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = 1.0;
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    f = [mu[1]*x[2]*u[1] -mu[2]*x[1]*u[1]];
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = 0.0;
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = mu[2];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[1]*n[1] + f[2]*n[2] + tau[1]*(u[1]-0.0);
    return fb;
end
function initu(x, mu, eta)
    u0 = exp(-(x[1]-0.25)*(x[1]-0.25)/0.01-(x[2]-0.0)*(x[2]-0.0)/0.01);
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# utilities

using Random
using DataFrames
using CSV

# random number generator for system:
seed = 10
rng = MersenneTwister(seed);

global dt = 0.1
#global tf = 400000.0
#global tf = 20000.0
tf = 1000000.0
global t_burn = 1000.0

global linelist_W = DataFrame(exposure_days = [0.0],
                       t_latent = [0.0],
                       t_incubation = [0.0],
                       t_post_incubation = [0.0],
                       tested_positive = [0],
                       expressed_symptoms = [0],
                       time_discharged = [0.0],
                       symptomatic = [0],
                       vaccinated = [0],
                       FoI_max = [0.0],
                       FoI_community = [0.0]
                      )

global linelist_W_entry = Dict(:exposure_days => 0.0,
                       :t_latent => 0.0,
                       :t_incubation => 0.0,
                       :t_post_incubation => 0.0,
                       :tested_positive => 0,
                       :expressed_symptoms => 0,
                       :time_discharged => 0.0,
                       :symptomatic => 0,
                       :vaccinated => 0,
                       :FoI_max => 0.0,
                       :FoI_community => 0.0
                      )

global linelist_T = DataFrame(exposure_days = [0.0],
                       days_in_quar = [0.0],
                       days_in_extended_quar = [0.0],
                       days_in_iso = [0.0],
                       t_latent = [0.0],
                       t_incubation = [0.0],
                       t_post_incubation = [0.0],
                       time_discharged = [0.0],
                       index_case = [0],
                       symptomatic = [0],
                       vaccinated = [0],
                       compliance = [0.0],
                       FoI_max = [0.0],
                       FoI_community = [0.0]
                      )

global linelist_T_entry = Dict(:exposure_days => 0.0,
                        :days_in_quar => 0.0,
                        :days_in_extended_quar => 0.0,
                        :days_in_iso => 0.0,
                        :t_latent => 0.0,
                        :t_incubation => 0.0,
                        :t_post_incubation => 0.0,
                        :time_discharged => 0.0,
                        :index_case => 0,
                        :symptomatic => 0,
                        :vaccinated => 0,
                        :compliance => 0.0,
                        :FoI_max => 0.0,
                        :FoI_community => 0.0
                      )

function write_config_file(
  output_filename::String,
  R0::Float64,
  t_final::Float64,
  delta_t::Float64,
  t_relax::Float64,
  imporation_rate::Float64,
  vac_importation_rate::Float64,
  quar_duration::Float64,
  ext_duration::Float64,
  iso_duration::Float64,
  quar_environment::String,
  medi_hotel_flag::Bool,
  test_days::Vector{Float64},
  ext_test_days::Vector{Float64},
  test_report_delay::Float64,
  clinical_detection::Bool,
  groups_size::Int64,
  n_travellers::Int64,
  n_workers::Int64,
  work_roster::AbstractVector,
  compliance_prob::Float64,
  pVac_Workers::Float64,
  pVac_Travellers::Float64,
  Vaccinated_index_cases::Bool,
  Vac_efficacy_infection::Float64,
  Vac_efficacy_transmission::Float64,
  transmission_Fac_TTsame::Float64,
  transmissionFac_TTdiff::Float64,
  transmission_Fac_TW::Float64,
  transmission_Fac_WW::Float64)

  open(output_filename, "w") do io
    write(io, join(["configuration input parameters: \n
             ******* \n
             Reproductive Ratio (homogeneous calib.): $R0  \n
             simulation time: $t_final  \n
             time step: $delta_t  \n
             burn-in time (arb): $t_relax \n
             baseline importation rate p(I | unvac., arrival): $imporation_rate \n
             vaccinated importation rate p(I | vac., arrival): $vac_importation_rate \n
             quarantine duration: $quar_duration \n
             contact extension duration: $ext_duration \n
             case isolation duration (extended from symptom onset): $iso_duration \n
             quarantine environment: $quar_environment \n
             contacts transferred to MediHotel: $medi_hotel_flag \n
             test days: $test_days \n
             extension test days: $ext_test_days \n
             test report delay: $test_report_delay \n
             clinical detection (via symptoms): $clinical_detection \n
             group size (close contacts): $groups_size \n
             n travellers: $n_travellers \n
             n workers: $n_workers \n
             work roster: $work_roster \n
             compliance probability: $compliance_prob \n
             proportion of workers vaccinated: $pVac_Workers \n
             proportion of travellers vaccinated: $pVac_Travellers \n
             all index cases vaccinated: $Vaccinated_index_cases \n
             vaccine efficacy against infection: $Vac_efficacy_infection \n
             vaccine efficacy against transmission: $Vac_efficacy_transmission \n
             Hotel transmission factor (Traveller, Traveller, same group): $transmission_Fac_TTsame \n
             Hotel transmission factor (Traveller, Traveller, different group): $transmissionFac_TTdiff \n
             Hotel transmission factor (Traveller, Worker): $transmission_Fac_TW \n
             Hotel transmission factor (Worker, Worker): $transmission_Fac_WW \n
             *******"]))
           end
end

function write_output_summary(
    output_filename::String,
    edays_pw_T::Float64,
    edays_pw_T_uv::Float64,
    edays_pw_T_vac::Float64,
    edays_pw_W::Float64,
    edays_pw_W_uv::Float64,
    edays_pw_W_vac::Float64,
    edays_per_inf_arrival_T::Float64,
    edays_per_inf_arrival_W::Float64,
    infected_arrivals_pw_uvac::Float64,
    infected_arrivals_pw_vac::Float64,
    detections_pw_tests::Float64,
    detections_pw_clin::Float64,
    travellers_discharged_pw::Float64,
    I_travellers_discharged_pw::Float64,
    secondary_cases_pw::Float64)

  open(output_filename, "w") do io
    write(io, join(["simulation output summary: \n
             ******* \n
             exposure days per week (Travellers): $edays_pw_T
             exposure days per week (Travellers, unvaccinated): $edays_pw_T_uv
             exposure days per week (Travellers, vaccinated): $edays_pw_T_vac \n
             exposure days per week (Workers): $edays_pw_W \n
             exposure days per week (Workers, unvaccinated): $edays_pw_W_uv \n
             exposure days per week (Workers, vaccinated): $edays_pw_W_vac \n
             exposure days per infected arrival (Travellers): $edays_per_inf_arrival_T \n
             exposure days per infected arrival (Workers): $edays_per_inf_arrival_W \n
             infected arrivals per week (unvaccinated): $infected_arrivals_pw_uvac \n
             infected arrivals per week (vaccinated): $infected_arrivals_pw_vac \n
             detected cases per week (PCR test): $detections_pw_tests \n
             detected cases per week (symptom onset): $detections_pw_clin \n
             travellers discharged per week: $travellers_discharged_pw \n
             infected travellers discharged per week: $I_travellers_discharged_pw \n
             secondary cases per week: $secondary_cases_pw \n
             *******"]))
           end

end


text:
import Base: showerror

struct IGCParseException <: Exception
    line_id::Int
    line::String
    e::Exception
end
IGCParseException(line, e; line_id=-1) = IGCParseException(line_id, line, e)

function showerror(io::IO, e::IGCParseException)
    msg = "Can't parse"
    if e.line_id > 0
        msg *= " line $(e.line_id)"
    end
    msg *= " \"$(e.line)\""
    print(io, msg)
    # print(io, " ")
    # print(io, e.e)
    throw(e.e)
end

struct IGCWriteException <: Exception
    msg::String
end

struct IGCNonUniqueRecordException <: Exception
    msg::String
end


text:
## add in simple derivative operator using the ForwardDiff package
## D(f), D(f,k) returns the derivative to power k (k=1 is default)

"""

Take derivative of order `k` of a function.

Arguments:
* `f::Function`: a mathematical function from R to R.
* `k::Int=1`: A non-negative integer specifying the order of the
  derivative. Values larger than 8 can be slow to compute.

Wrapper around `derivative` function in `ForwardDiff`

"""
function D(f::Function, k::Int=1)
    k < 0 && error("The order of the derivative must be non-negative")
    k == 0 && return(x -> f(float(x)))
    D(x -> ForwardDiff.derivative(f, x), k-1)
end

D2(f::Function) = D(f, 2)

## This conflicts with a definition in Calculus, but is better.
## Base.ctranspose(f::Function) = x -> ForwardDiff.derivative(f, x)


text:
using Test
import XSBenchThreads
import ArgParse

include("../../src/Main.jl")

@testset "test_main" begin
    @test main(["-g", "32000", "-m", "event"]) == 0
end


text:
using Test, StaticArrays, EllZeroTrendFiltering
using EllZeroTrendFiltering: minimize_wrt_x2

# Addition of quaratic forms
Q = QuadraticForm(@SMatrix([2.0 0; 0 1]), @SVector([0.0,0]), 1.0) +
    QuadraticForm(@SMatrix([2.0 1; 1 2]), @SVector([1.0,0]), 0.0)

@test Q == QuadraticForm(@SMatrix([4.0 1; 1 3]), @SVector([1.0,0]), 1.0)

# Maybe somewhat surprising that the tests pass considering equalities between
# Floats, or is the aritmhetic exact

# x^2 + y^2 + 1
Q1 = QuadraticForm(@SMatrix([1.0 0; 0 0]), @SVector([0.0, 0]), 0.5)
p1 = QuadraticPolynomial(1.0, 0.0, 0.5)
@test minimize_wrt_x2(Q1, p1) == QuadraticPolynomial(1,0,1)
@test Q1(1, 2) == 1.5

# (2x - y + 1)^2 + x^2 - 2x + 4
# = 5x^2 + y^2 + 2 + 2x - 2y - 4xy
Q2 = QuadraticForm(@SMatrix([5.0 -2; -2 0.5]), @SVector([2.0, -2]), 1.0)
p2 = QuadraticPolynomial(0.5, 0.0, 5.0)
@test minimize_wrt_x2(Q2, p2) == QuadraticPolynomial(1,-2,5)

# x^2 + 2x + 2
Q3 = QuadraticForm(@SMatrix([1.0 0; 0 0]), @SVector([2.0, 0]), 2.0)
p3 = QuadraticPolynomial(0.0, 0.0, 0.0)
@test minimize_wrt_x2(Q3, p3) == QuadraticPolynomial(1.0,2.0,2.0)

# y^2 + 2y + 2
Q4 = QuadraticForm(@SMatrix([0.0 0; 0 1]), @SVector([0.0, 2]),2.0)
p4 = QuadraticPolynomial(0.0, 0.0, 0.0)
@test minimize_wrt_x2(Q4, p4) == QuadraticPolynomial(0,0,1)


text:
using Test
using GraphIO.Graph6

@testset "Graph6" begin
    n1 = (30, UInt8.([93]))
    n2 = (12345, UInt8.([126; 66; 63; 120]))
    n3 = (460175067, UInt8.([126; 126; 63; 90; 90; 90; 90; 90]))
    ns = [n1; n2; n3]
    for n in ns
        @test Graph6._g6_N(n[1]) == n[2]
        @test Graph6._g6_Np(n[2])[1] == n[1]
    end

    for g in values(graphs)
        readback_test(Graph6Format(), g, "graph1")
    end

    f = gettempname()
    write_test(Graph6Format(), graphs, f; remove=false, silent=true)
    read_test_mult(Graph6Format(), graphs, f)
    rm(f)
end


text:
# import stuff
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using Statistics
@everywhere using GRASS
@everywhere using SharedArrays
@everywhere using EchelleCCFs
using CSV
using DataFrames
using LaTeXStrings
using LsqFit

# define some functions
include(GRASS.moddir * "figures/fig_functions.jl")

# some global stuff
const N = round.(Int, 2 .^ range(6, 10, step=0.5))
const Nt = 100
const Nloop = 24

# get command line args and output directories
run, plot = parse_args(ARGS)
grassdir, plotdir, datadir = check_plot_dirs()

function main()
    # set up parameters for lines
    lines = [5434.5]
    depths = [0.8]
    res = 700000.0
    top = NaN
    contiguous_only=true
    spec = SpecParams(lines=lines, depths=depths, resolution=res,
                      extrapolate=true, contiguous_only=contiguous_only)

    # allocate shared arrays
    avg_avg_res = SharedArray{Float64}(length(N))
    std_avg_res = SharedArray{Float64}(length(N))
    avg_rms_res = SharedArray{Float64}(length(N))
    std_rms_res = SharedArray{Float64}(length(N))

    # calculate
    @sync @distributed for i in 1:length(N)
    	println("running resolution N = " * string(N[i]))
        disk = DiskParams(N=N[i], Nt=Nt)
    	avg_avg1, std_avg1, avg_rms1, std_rms1 = spec_loop(spec, disk, Nloop, top=top)
        avg_avg_res[i] = avg_avg1
        std_avg_res[i] = std_avg1
    	avg_rms_res[i] = avg_rms1
    	std_rms_res[i] = std_rms1
    end

    # make data frame
    df = DataFrame()
    df[!,:res] = N
    df[!,:avg_avg_res] = avg_avg_res
    df[!,:std_avg_res] = std_avg_res
    df[!,:avg_rms_res] = avg_rms_res
    df[!,:std_rms_res] = std_rms_res

    # write to CSV
    fname = datadir * "rms_vs_res.csv"
    CSV.write(fname, df)
    return nothing
end

# run the simulation
if run
    main()
end

# plotting code block
if plot
    # plotting imports
    import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff()
    using PyCall; animation = pyimport("matplotlib.animation")
    mpl.style.use(GRASS.moddir * "figures/fig.mplstyle")

    # read in the data
    fname = datadir * "rms_vs_res.csv"
    df = CSV.read(fname, DataFrame)
    res = df.res
    rms_res = df.avg_rms_res
    std_res = df.std_rms_res

    # get the error
    err_res = std_res ./ sqrt(Nloop)

    # fit the data
    @. power_law(x, p) = p[1] * x^(-p[2])
    fit = curve_fit(power_law, res, rms_res, [1.0, 1.0])
    res_fit = range(6, 2400, length=1000)
    println(">>> Best fit power law index = " * string(fit.param[2]))
    x = string(round(fit.param[2], sigdigits=2))

    # initialize figure and set log scales
    fig = plt.figure()
    ax1 = fig.add_subplot()
    ax1.set_xscale("log", base=2)
    # ax1.set_yscale("log", base=10)

    # plot the data
    ax1.errorbar(res, rms_res, yerr=err_res, capsize=3.0, color="black", fmt=".")
    ax1.plot(res_fit, power_law(res_fit, fit.param), "k--", alpha=0.4,
             label = L"{\rm Power\ law\ index\ } \approx\ %$x ")

    # add a line corresponding to average footprint area
    ax1.axvline(132, ymin=0, ymax=1, color="black", lw=1.5, ls="--")

    # plot the literature values
    xrng = ax1.get_xlim()
    xs = range(xrng[1], xrng[2], length=2)
    Elsworth = [repeat([0.319-0.09], 2), repeat([0.319+0.09], 2)]
    Palle = [repeat([0.461-0.1], 2), repeat([0.461+0.1], 2)]
    ax1.fill_between(xs, Elsworth..., alpha=0.5, fc="tab:orange", ec="white",
                     hatch="/", label=L"{\rm Elsworth\ et\ al.\ (1994)}")
    ax1.fill_between(xs, Palle..., alpha=0.5, fc="tab:green", ec="white",
                     hatch="\\", label=L"\textnormal{Pall\'{e}}\ {\rm et\ al.\ (1999)}")

    # set the axes limits + labels
    ax1.set_xlim(2^5.9, 2^10.1)
    ax1.set_ylim(0.01, 1.15)
    ax1.yaxis.set_major_formatter(mpl.ticker.ScalarFormatter())
    ax1.set_xlabel(L"N")
    ax1.set_ylabel(L"{\rm RMS}_{\rm RV}\ {\rm (m s}^{-1})")
    ax1.legend()

    # save the figure
    fig.savefig(plotdir * "fig4.pdf")
    plt.clf(); plt.close()
    println(">>> Figure written to: " * plotdir * "fig4.pdf")
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using GraphPlot
#using SimpleGraphs
#using Graphs
using OffsetArrays
#using Compose, Cairo, Fontconfig

function readPRP(filename)

	allLines = open(filename) do f
		readlines(f)
	end

	type = parse(Int, split(allLines[1], " ")[2])
	n = parse(Int, split(allLines[2], " ")[2])
	l = parse(Int, split(allLines[3], " ")[2])
	u = parse(Int, split(allLines[4], " ")[2])
	f = parse(Int, split(allLines[5], " ")[2])

	#To handle scientific notation (1e+10 for exemple)
	C = Int(parse(Float64, split(allLines[6], " ")[2]))

	Q = parse(Int, split(allLines[7], " ")[2])
	k = parse(Int, split(allLines[8], " ")[2])

	#=
		n = Nombre de revendeur (Ne comprend pas le fournisseur)
		l = Horizon de planification ( t dans {1,...,l} )
			(Les revendeurs seront fournis l fois (une fois par t))

		f = Cout fixe (Cout de setup), a chaque production à un temps t

		u = Cout par unité produite

		C = Capacité de production max du fournisseur ??? --> oui c'est M_t dans le PL du LSP, dans la contrainte p_t≤M_ty_t> quand y_t=1 on peut produire au max M_t=C, en fait, on considère que la capacité de production du fournisseur est infini

		??? k = Nombre de véhicule de transports disponible ??? ou nombre de tournées
		Q = Capacité maximale d'un véhicule de transport
	=#
	params = Dict("type" => type, "n" => n, "l" => l, "u" => u, "f" => f, "C" => C, "Q" => Q, "k" => k)

	nextIndex = 9

	if type == 2
		mc = parse(Int, split(allLines[9], " ")[2])

		#Qu'es ce que le type 2 (B) et c'est quoi ce "mc" en plus ?
		#Les coûts sont calculés selon une formule, sur les instances de type 2, le calcul est pondérée par ce mc
		params["mc"] = mc

		nextIndex += 1
	end

	#=
	Informations du fournisseur et des revendeurs: [Fourni{}, reven_1{}, ..., reven_n{}]
		h = Cout de stockage
		L = Capacité de stockage max (capacité de stockage du fournisseur est CONSIDÉRÉ COMME INFINI)
		L0 = Stock initial

		x et y (2 premières valeurs) = Coordonnées sur la "carte" des revendeurs et du fournisseur
	=#
	nodes = Vector{Dict}() #élement à l'indice i = noeud n°i

	for line in allLines[nextIndex:nextIndex+n]

		#Read line
		allElems = split(line, " ")

		newNode = Dict{String, Int}()

		newNode["initial_index"]=Int(parse(Float64, allElems[1]))
		newNode["x"] = Int(parse(Float64, allElems[2]))
		newNode["y"] = Int(parse(Float64, allElems[3]))

		#Tous les éléments d'indice impair de la ligne line en commençant par l'élément d'indice 5
		for i in filter(isodd, eachindex(allElems[5:end])) .+ 4 #[5, 7, 9]
			newNode[allElems[i]] = Int(parse(Float64, allElems[i+1]))
		end
		#=
		Equivalent à :
			newNode["h"]=Int(parse(Float64, allElems[6])) #h
			newNode["L"]]=Int(parse(Float64, allElems[7])) #L
			newNode["L0"]=Int(parse(Float64, allElems[8])) #L0
		=#

		push!(nodes, newNode)

	end

	#0 based index nodes, more coherent with the project
	nodes = OffsetVector(nodes, 0:(size(nodes)[1] - 1))

	nextIndex = nextIndex+n+2

	#=
	d_it:
		Demandes des revendeurs i au temps t.
	=#
	demands = Array{Int, 2}(undef, n, l)#élement à l'indice i,t = demande du revendeur n°i au pas de temps t

	for (i, line) in enumerate(allLines[nextIndex:nextIndex+n-1])

		cost = Vector{Int}()

		for elem in split(line, " ")[2:end]

			#Pour eviter une erreur avec les fins de lignes
			if cmp("", elem) == 0
				continue
			end

			push!(cost, Int(parse(Float64, elem)) )
		end

		demands[i, 1:end] = cost

	end

	#=
	c_ij:
		Cout de transport du noeud i vers le noeud j
		(exemple: transportFee[(1, 2)] -> Cout du noeud 1 vers le 2)
	=#
	transportFee = Dict{Tuple{Int, Int}, Float64}()

	#Le type determine la fonction cout
	if type == 1
		costFunc = (x1, y1, x2, y2) -> floor(1/2 + sqrt( (x1 - x2)^2 + (y1 - y2)^2 ))
	elseif type == 2
		costFunc = (x1, y1, x2, y2) -> mc * sqrt( (x1 - x2)^2 + (y1 - y2)^2 )
	end

	for i in 0:(n - 1)
		for j in (i + 1):n

			xi=nodes[i]["x"]
			yi=nodes[i]["y"]
			xj=nodes[j]["x"]
			yj=nodes[j]["y"]
			transportFee[(i, j)] = costFunc(xi, yi, xj, yj)

			#Graphe Orienté Complet
			transportFee[(j, i)] = transportFee[(i,j)]

		end
	end

	return params, nodes, demands, transportFee

end

# function generateGraphComplet(nodes)

#     g=SimpleDiGraph(length(nodes))	 # Creation of a directed graph with 1 node

# 	for i in 0:length(nodes)-1
# 		for j in 0:length(nodes)-1

# 			if(i!=j)
# 				add_edge!(g,i,j)
# 			end

# 		end
# 	end

#     return g

# end

# function generateGraph(nodes,edges)

#     g=SimpleDiGraph(length(nodes))	 # Creation of a directed graph with 1 node

# 	for (node1,node2) in edges
# 		add_edge!(g,node1,node2)
# 	end

#     return g

# end

# function generateGraphPDF(destination_file_name, file_name, graph)
# 	nodesize = [Graphs.outdegree(graph, v) for v in Graphs.vertices(graph)] #Gérer la taille des noeuds
# 	layout = (args...)->spring_layout(args...; C=20) #Gérer l'espacement des noeuds
# 	draw(PDF(destination_file_name*"/"*file_name, 16, 16), gplot(graph, nodesize=nodesize, layout=layout, EDGELINEWIDTH=0.01))
# end


text:
function load_exp(filename)
    exp_data = readdlm(filename)  # [t, T, m]
    index = indexin(unique(exp_data[:, 1]), exp_data[:, 1])
    exp_data = exp_data[index, :]
    exp_data[:, 3] = exp_data[:, 3] / maximum(exp_data[:, 3])
    return exp_data
end

l_exp_data = [];
l_exp_info = zeros(Float64, length(l_exp), 3);
for (i_exp, value) in enumerate(l_exp)
    filename = string("exp_data/expdata_no", string(value), ".txt")
    exp_data = Float64.(load_exp(filename))

    if value == 4
        exp_data = exp_data[1:60, :]
    elseif value == 5
        exp_data = exp_data[1:58, :]
    elseif value == 6
        exp_data = exp_data[1:60, :]
    elseif value == 7
        exp_data = exp_data[1:71, :]
    end

    push!(l_exp_data, exp_data)
    l_exp_info[i_exp, 1] = exp_data[1, 2] # initial temperature, K
end
l_exp_info[:, 2] = readdlm("exp_data/beta.txt")[l_exp];
l_exp_info[:, 3] = log.(readdlm("exp_data/ocen.txt")[l_exp] .+ llb);


text:
using LinearAlgebra: norm_sqr

function polar_to_cartesian(ρ, θ)
    x = ρ * cos(θ)
    y = ρ * sin(θ)
    Vector2(x, y)
end
polar_to_cartesian(polar::Vector2) = polar_to_cartesian(polar[1], polar[2])

function cartesian_to_polar(x, y)
    ρ = √(x ^ 2 + y ^ 2)
    θ = acos(x / ρ)
    Vector2(ρ, θ)
end
cartesian_to_polar(xy::Vector3) = cartesian_to_polar(xy[1], xy[2])

@_with_default_rng function random_point_on_circle(radius::Real)
    θ = 2π * rand(rng)
    polar_to_cartesian(radius, θ)
end

@_with_default_rng random_2d_direction() = random_point_on_circle(rng, 1.0)

@_with_default_rng function random_point_on_disk(radius::Real)
    ρ = √(rand(rng) * radius ^ 2)
    random_point_on_circle(rng, ρ)
end


text:
#WireObject.jl
# much of the magic happens here as this object can be passed around,
# representing Wire object state.

type WireObject{R};
  lexical_representation::String
end

Base.range{R}(::WireObject{R}) = R
Base.length{R}(::WireObject{R}) = length(R)

lr(w) = w.lexical_representation

#lexical transformations for logical operators.
#unary logical operators
Base.:~{R}(w::WireObject{R}) = WireObject{R}(string("~(", w.lexical_representation, ")"))
Base.:&{R}(w::WireObject{R}) = WireObject{0:0v}(string("&(", w.lexical_representation, ")"))
Base.:|{R}(w::WireObject{R}) = WireObject{0:0v}(string("|(", w.lexical_representation, ")"))
Base.:^{R}(w::WireObject{R}) = WireObject{0:0v}(string("^(", w.lexical_representation, ")"))
#make single-wire objects pass through without doing any collection.
Base.:&(w::WireObject{0:0v}) = w
Base.:|(w::WireObject{0:0v}) = w
Base.:^(w::WireObject{0:0v}) = w

Base.:&(w1::WireObject, w2::WireObject, w3::WireObject, ws::WireObject...) = WireObject{0:0v}(string("&({",join([lr(w) for w in [w1, w2, w3, ws...]], ", "),"})"))
Base.:|(w1::WireObject, w2::WireObject, w3::WireObject, ws::WireObject...) = WireObject{0:0v}(string("|({",join([lr(w) for w in [w1, w2, w3, ws...]], ", "),"})"))
Base.:^(w1::WireObject, w2::WireObject, w3::WireObject, ws::WireObject...) = WireObject{0:0v}(string("^({",join([lr(w) for w in [w1, w2, w3, ws...]], ", "),"})"))

#presume that R & S have been checked to be the same size (for now).
#binary logical operators
Base.:&{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) & $(lr(rhs)))")
Base.:|{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) | $(lr(rhs)))")
Base.:^{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) ^ $(lr(rhs)))")

xorxnor{R}(w::WireObject{R}) = WireObject{1:0v}("{^($(lr(w))), ~^($(lr(w)))}")

#lexical transformations for arithmetic operatiors
Base.:*{R}(lhs::Int, rhs::WireObject{R}) = WireObject{range(length(R) * lhs)}("{$lhs{$(lr(rhs))}}")

#other arithmetic operators
Base.:+{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) + $(lr(rhs)))")
Base.:-{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) - $(lr(rhs)))")
Base.:*{R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R) + length(S))}("($(lr(lhs)) * $(lr(rhs)))")

Base.:-{R}(lhs::WireObject{R}) = WireObject{range(length(R))}("-($(lr(lhs)))")

#shifters
Base.:(<<){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) << $(lr(rhs)))")
Base.:(>>){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) >> $(lr(rhs)))")
Base.:(>>>){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("(\$signed($(lr(lhs))) >>> $(lr(rhs)))")

#comparative operators
Base.:>{R,S}(lhs::WireObject{R}, rhs::WireObject{S})  = WireObject{range(length(R))}("($(lr(lhs)) > $(lr(rhs)))")
Base.:<{R,S}(lhs::WireObject{R}, rhs::WireObject{S})  = WireObject{range(length(R))}("($(lr(lhs)) < $(lr(rhs)))")
Base.:(>=){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) >= $(lr(rhs)))")
Base.:(<=){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) <= $(lr(rhs)))")
Base.:(==){R,S}(lhs::WireObject{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(lr(lhs)) == $(lr(rhs)))")

function Base.getindex{R}(w::WireObject{R}, r::VerilogRange)
  if (r.stop == r.start)
    WireObject{range(length(r))}(string("$(lr(w))[$(r.stop)]"))
  elseif (r.stop > r.start)
    WireObject{range(length(r))}(string("$(lr(w))[$(r.stop):$(r.start)]"))
  else
    WireObject{range(length(r))}(string("{",join(["$(lr(w))[$idx]" for idx in reverse(r)],", "),"}"))
  end
end

Base.getindex{R}(w::WireObject{R}, r::RelativeRange) = getindex(w, parse_msb(r, R))
Base.getindex{R}(w::WireObject{R}, ::Type{msb}) = getindex(w, R.stop)
Base.getindex{R}(w::WireObject{R}, m::msb) = getindex(w, R.stop - m.value)

Base.getindex{R}(w::WireObject{R}, i::Int) = WireObject{0:0v}(string("$(lr(w))[$i]"))
#Concatenation with wires using the Wire() operator.  Make the assumption that
#any "wire" object that doesn't derive from an existing WireObject must be a
#constant.  We can then transparently overload the Wire() operator with no
#ambiguities.
typealias WOO Union{WireObject, Wire}

wo_concat(w::WireObject) = w.lexical_representation
wo_concat(w::Wire) = string("$(length(w))'b", join(reverse(w.values).*1))

function (::Type{Wire})(ws::WOO...)
  l = sum(length, ws)
  WireObject{range(l)}(string("{", join([wo_concat(w) for w in ws],",") ,"}"))
end

#binary logical operators with constants on the right:
Base.:&{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}(string("($(lr(lhs)) & {$(wo_concat(rhs))})"))
Base.:|{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}(string("($(lr(lhs)) | {$(wo_concat(rhs))})"))
Base.:^{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}(string("($(lr(lhs)) ^ {$(wo_concat(rhs))})"))
#binary logical operators with constants on the left:
Base.:&{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}(string("({$(wo_concat(lhs))} & $(lr(rhs)))"))
Base.:|{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}(string("({$(wo_concat(lhs))} | $(lr(rhs)))"))
Base.:^{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}(string("({$(wo_concat(lhs))} ^ $(lr(rhs)))"))

#arithmetic operators with constants on the left:
Base.:+{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) + $(lr(rhs)))")
Base.:-{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) - $(lr(rhs)))")
Base.:*{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R) + length(S))}("($(wo_concat(lhs)) * $(lr(rhs)))")
#arithmetic operators with constants on the right:
Base.:+{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) + $(wo_concat(rhs)))")
Base.:-{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) - $(wo_concat(rhs)))")
Base.:*{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R) + length(S))}("($(lr(lhs)) * $(wo_concat(rhs)))")

#comparative operators with constants on the left:
Base.:<{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) < $(lr(rhs)))")
Base.:>{R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) > $(lr(rhs)))")
Base.:(<=){R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) <= $(lr(rhs)))")
Base.:(>=){R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) >= $(lr(rhs)))")
Base.:(==){R,S}(lhs::Wire{R}, rhs::WireObject{S}) = WireObject{range(length(R))}("($(wo_concat(lhs)) == $(lr(rhs)))")
#comparative operators with constants on the right:
Base.:<{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) < $(wo_concat(rhs)))")
Base.:>{R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) > $(wo_concat(rhs)))")
Base.:(<=){R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) <= $(wo_concat(rhs)))")
Base.:(>=){R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) >= $(wo_concat(rhs)))")
Base.:(==){R,S}(lhs::WireObject{R}, rhs::Wire{S}) = WireObject{range(length(R))}("($(lr(lhs)) == $(wo_concat(rhs)))")


text:
# Verify solutions against matpower results
using Test
using ExaPF
using FiniteDiff
using ForwardDiff
using LinearAlgebra
using KernelAbstractions

@testset "RGM Optimal Power flow 9 bus case" begin
    datafile = joinpath(INSTANCES_DIR, "case9.m")

    nlp = ExaPF.ReducedSpaceEvaluator(datafile)
    uk = ExaPF.initial(nlp)

    # solve power flow
    ExaPF.update!(nlp, uk)

    # reduced gradient method
    iterations = 0
    iter_max = 100
    step = 0.0001
    norm_grad = 10000
    norm_tol = 1e-5

    iter = 1

    # initial gradient
    grad = similar(uk)
    wk = similar(uk)
    up = copy(uk)
    fill!(grad, 0)

    while norm_grad > norm_tol && iter < iter_max
        ExaPF.update!(nlp, uk)
        c = ExaPF.objective(nlp, uk)
        ExaPF.gradient!(nlp, grad, uk)
        # compute control step
        wk = uk - step*grad
        ExaPF.project!(uk, wk, nlp.u_min, nlp.u_max)
        norm_grad = norm(uk .- up, Inf)
        iter += 1
        up .= uk
    end
    @test iter == 39
    @test isapprox(uk, [1.1, 1.1, 1.1, 1.343109921105559, 0.9421135274454701], atol=1e-4)
end


text:
using DitherPunk
using Documenter
using Literate

EXAMPLE_DIR = joinpath(@__DIR__, "literate")
OUT_DIR = joinpath(@__DIR__, "src/generated")

# Use Literate.jl to generate docs and notebooks of examples
for example in readdir(EXAMPLE_DIR)
    EXAMPLE = joinpath(EXAMPLE_DIR, example)

    Literate.markdown(EXAMPLE, OUT_DIR; documenter=true) # markdown for Documenter.jl
    Literate.notebook(EXAMPLE, OUT_DIR) # .ipynb notebook
    Literate.script(EXAMPLE, OUT_DIR) # .jl script
end

DocMeta.setdocmeta!(DitherPunk, :DocTestSetup, :(using DitherPunk); recursive=true)

makedocs(;
    modules=[DitherPunk],
    authors="Adrian Hill",
    sitename="DitherPunk.jl",
    format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", assets=String[]),
    pages=[
        "Getting started" => "index.md",
        "Gallery" => Any[
            "Gradient" => "generated/gallery_gradient.md",
            "Images" => "generated/gallery_images.md",
        ],
        "Fun applications" => Any[
            "SDF halftoning" => "generated/sdf_halftoning.md",
            "Color cycling" => "generated/color_cycling.md",
            "Palette swapping" => "generated/palette_swapping.md",
        ],
        "API Reference" => "api.md",
    ],
)

deploydocs(; repo="github.com/JuliaImages/DitherPunk.jl")


text:
# Top level test script for Stan.jl
using Base.Test

println("Running tests for Stan-j0.5-v1.0.2:")

code_tests = ["test_env.jl",
              "test_utilities.jl",
              "test_cmdtype.jl"]

# Run execution_tests only if CmdStan is installed and CMDSTAN_HOME is set correctly.
execution_tests = [
  "test_bernoulli.jl",
  "test_bernoulli_optimize.jl",
  "test_bernoulli_diagnose.jl",
  "test_bernoulli_variational.jl",
  "test_binomial.jl",
  "test_binormal.jl",
  "test_schools8.jl",
  "test_dyes.jl",
  "test_kidscore.jl"
]

for my_test in code_tests
    println("\n  * $(my_test) *")
    include(my_test)
end

if CMDSTAN_HOME != ""
  println("CMDSTAN_HOME set. Try to run bernoulli.")
  try
    for my_test in execution_tests
        println("\n  * $(my_test) *")
        include(my_test)
    end
  catch e
     println("CMDSTAN_HOME initialized, but CmdStan not installed properly.")
     println(e)
     println("No simulation runs have been performed.")
  end
else
  println("\n\nCMDSTAN_HOME not initialized. Skipping all tests that depend on CmdStan!\n")
end

println("\n")


text:
using Test
using MLDataUtils
using StatsBase
using UnicodePlots
using DataFrames
using Statistics

tests = [
    "tst_datapattern.jl"
    "tst_noisy_function.jl"
    "tst_feature_scaling.jl"
    "tst_datasets.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
isLoadingPlottingEndpoints = false
try
    getfield(Main, :RoMEPlotting)
    # If that works, then we can enable plotting
    isLoadingPlottingEndpoints = true
catch ex
end

if isLoadingPlottingEndpoints
    @info "[ZMQ Server] Plotting endpoints are enabled!"

    function plotKDE(configDict, fg, requestDict)::Dict{String, Any}
        error("Not implemented yet!")
    end

    function plotPose(configDict, fg, requestDict)::Dict{String, Any}
        error("Not implemented yet!")
    end

    function drawPoses(configDict, fg, requestDict)::Dict{String, Any}
        plotRequest = Unmarshal.unmarshal(PlotRequest, requestDict["payload"])
        pl = RoMEPlotting.drawPoses(fg);
        io = IOBuffer();
        Gadfly.draw(PNG(io, (plotRequest.widthPx)px, (plotRequest.heightPx)px), pl);
        return Dict{String, Any}("status" => "ok", "mimetype" => "image/png", "encoding" => "none", "data" => io.data)
    end

    function drawPosesLandms(configDict, fg, requestDict)::Dict{String, Any}
        plotRequest = Unmarshal.unmarshal(PlotRequest, requestDict["payload"])
        pl = RoMEPlotting.drawPosesLandms(fg);
        io = IOBuffer();
        Gadfly.draw(PNG(io, (plotRequest.widthPx)px, (plotRequest.heightPx)px), pl);
        return Dict{String, Any}("status" => "ok", "mimetype" => "image/png", "encoding" => "none", "data" => io.data)
    end
else
    @info "[ZMQ Server] Plotting endpoints are disabled in ZMQ server. Please call 'using RoMEPlotting' before 'using Caesar' if you would like to enable them."
end


text:
"""
    from_toml(::Type{T}, filename::String; kw...) where T

Convert a given TOML file `filename` to an option type `T`. Valid fields
can be override by keyword arguments. See also [`from_dict`](@ref).
"""
function from_toml(::Type{T}, filename::String; kw...) where {T}
    is_option(T) || error("not an option type")
    d = TOML.parsefile(filename)

    filepath = normpath(filename)
    d["#metadata#"] = Dict{String, Any}(
        "file" => filepath,
        "dir" => dirname(filepath),
        "format" => "TOML",
    )
    return from_dict(T, d; kw...)
end

"""
    from_toml_if_exists(::Type{T}, filename::String; kw...) where T

Similar to [`from_toml`](@ref) but will create the option instance
via `from_kwargs(T;kw...)` instead of error if the file does not exist.
"""
function from_toml_if_exists(::Type{T}, filename::String; kw...) where {T}
    if isfile(filename)
        return from_toml(T, filename; kw...)
    else
        return from_kwargs(T; kw...)
    end
end


text:
module TestChecks
using Test
using ArgCheck
using ArgCheck: pretty_string
using Random: randstring

macro catch_exception_object(code)
    quote
        err = try
            $(esc(code))
            nothing
        catch e
            e
        end
        if err == nothing
            error("Expected exception, got $err.")
        end
        err
    end
end

@testset "Chained comparisons" begin
    #6
    x=y=z = 1
    @test x == y == z
    @argcheck x == y == z
    z = 2
    @test_throws ArgumentError @argcheck x == y == z

    @test_throws ArgumentError @argcheck 1 ≈ 2 == 2
    @argcheck 1 == 1 ≈ 1 < 2 > 1.2
    @test_throws DimensionMismatch @argcheck 1 < 2 ==3 DimensionMismatch
end

@testset "@argcheck" begin
    @test_throws ArgumentError @argcheck false
    @argcheck true

    x = 1
    @test_throws ArgumentError (@argcheck x > 1)
    @argcheck x>0 # does not throw

    n =2; m=3
    @test_throws DimensionMismatch (@argcheck n==m DimensionMismatch)
    @argcheck n==n DimensionMismatch

    denominator = 0
    @test_throws DivideError (@argcheck denominator != 0 DivideError())
    @argcheck 1 !=0 DivideError()
end

# exotic cases
struct MyError <: Exception
    msg::String
end
struct MyExoticError <: Exception
    a::Int
    b::Int
end

falsy(args...;kw...) = false
truthy(args...;kw...) = true

@testset "exotic cases" begin
    @argcheck truthy()
    @test_throws ArgumentError @argcheck falsy()

    @argcheck begin
        multi_line_true_is_no_problem = true
        multi_line_true_is_no_problem
    end
    @test_throws DimensionMismatch @argcheck let
        falsy(1,2)
    end DimensionMismatch

    op() = (x,y) -> x < y
    x = 1; y = 2
    @argcheck op()(x,y)
    @test_throws ArgumentError @argcheck op()(y,x)
    @test_throws ArgumentError @argcheck begin false end
    @test_throws DivideError @argcheck falsy() DivideError()
    err = @catch_exception_object @argcheck false MyExoticError(1,2)
    @test err === MyExoticError(1,2)
end

@testset "error messages" begin
    x = 1.23455475675
    y = 2.345345345
    # comparison
    err = @catch_exception_object @argcheck x == y MyError
    @test err isa MyError
    msg = err.msg
    @test occursin(string(x), msg)
    @test occursin(string(y), msg)
    @test occursin("x", msg)
    @test occursin("y", msg)
    @test occursin("==", msg)

    x = 1.2
    y = 1.34
    z = -345.234
    err = @catch_exception_object @argcheck x < y < z
    msg = err.msg
    @test occursin(string(z), msg)
    @test occursin(string(y), msg)
    @test occursin("y", msg)
    @test occursin("z", msg)
    @test occursin("<", msg)
    @test !occursin(string(x), msg)

    ≦(a,b) = false
    err = @catch_exception_object @argcheck x ≦ y ≦ z
    msg = err.msg
    @test occursin("x", msg)
    @test occursin("y", msg)
    @test occursin(string(x), msg)
    @test occursin(string(y), msg)
    @test occursin("≦", msg)

    s = randstring()
    arr = rand(1000:9999, 1000)
    x = randn()
    err = @catch_exception_object @check falsy(x, arr, s)
    @test typeof(err) == CheckError
    msg = err.msg
    @test length(msg) < 2000
    @test occursin(pretty_string(x), msg)
    @test occursin(pretty_string(arr), msg)
    @test occursin(pretty_string(s), msg)

    x = 1.2
    y = 1.34
    z = -345.234
    err = @catch_exception_object @argcheck falsy([x y; z z])
    msg = err.msg
    @test occursin(string(x), msg)
    @test occursin(string(z), msg)
    @test occursin(string(y), msg)
    @test occursin("y", msg)
    @test occursin("z", msg)
    @test occursin("x", msg)
    @test occursin("f", msg)

    fail_function(args...) = false
    err = @catch_exception_object @argcheck fail_function(x,y,z) DimensionMismatch
    msg = err.msg

    @test err isa DimensionMismatch
    @test occursin(string(x), msg)
    @test occursin(string(z), msg)
    @test occursin(string(y), msg)
    @test occursin("y", msg)
    @test occursin("z", msg)
    @test occursin("x", msg)
    @test occursin("Got", msg)
    @test occursin("fail_function", msg)

    err = @catch_exception_object @argcheck issorted([2,1])
    @test !occursin("Got", err.msg)

    z = "some_keyword_arg"
    err = @catch_exception_object @argcheck falsy(x=z)
    @test occursin("z", err.msg)
    @test occursin(z, err.msg)

    args = ["these", "are", "splatargs"]
    err = @catch_exception_object @argcheck falsy(args...)
    @test occursin("args", err.msg)
    @test occursin(args[1], err.msg)

    myatol = 0.1
    err = @catch_exception_object @argcheck isapprox(1,2,atol=myatol)
    @test occursin("atol", err.msg)
    @test occursin("0.1", err.msg)

    kw = (atol=1.34, rtol=0.02)
    err = @catch_exception_object @argcheck isapprox(10, 20; kw...)
    @test occursin("atol", err.msg)
    @test occursin("1.34", err.msg)
    @test occursin("rtol", err.msg)
    @test occursin("0.02", err.msg)

    args = (10, 20)
    kw = (atol=1.34, rtol=0.02)
    err = @catch_exception_object @argcheck isapprox(args...; kw...)
    @test occursin("atol", err.msg)
    @test occursin("1.34", err.msg)
    @test occursin("rtol", err.msg)
    @test occursin("0.02", err.msg)
    @test occursin("10", err.msg)
    @test occursin("20", err.msg)

    # check argument orders
    err = @catch_exception_object @argcheck let
        a = 1.0; b = 1.2; atol = 0.1; nvalue = false; rtol = 0.05;
        @check isapprox(a, b, atol=atol, nans=nvalue; rtol)
    end
    locations = map(["a", "b", "atol", "nvalue", "rtol"]) do name
        findfirst(name, err.msg)
    end
    @test all(x -> x isa UnitRange, locations)
    @test issorted(getindex.(locations, 1))

    x = 1.234
    err = @catch_exception_object @argcheck (!isfinite)(x)
    @test occursin(string(x), err.msg)
    err = @catch_exception_object @argcheck !isfinite(x)
    @test_broken occursin(string(x), err.msg)
end

@testset "complicated calls" begin
    @argcheck issorted([2,1], rev=true)
    @argcheck issorted([2,1]; rev=true)
    xs = [[1,2]]
    @argcheck issorted(xs...)

    @argcheck truthy(xs...,xs...)
    @test_throws MyError @argcheck falsy(xs...,xs...) MyError

    kw1 = Dict(:x =>1)
    kw2 = Dict(:y =>2)
    @argcheck truthy(;kw1...)
    @argcheck truthy(;kw1..., kw2...)
    @test_throws MyError @argcheck falsy(cos, xs...,xs...;kw1...,kw2...,foo=3) MyError

    if VERSION > v"1.4.0"
        atol = 2
        @argcheck isapprox(1, 2; atol)
        @argcheck isapprox(1, 2; atol, rtol=0.2)
        @argcheck isapprox(1, 2; :atol=>2)
        @argcheck isapprox(1, 2; :atol=>2, rtol=0.1)
        @argcheck isapprox(1, 2; (() -> :atol)()=>2, rtol=0.1)

        foo = 3
        foosym = :foo
        @argcheck truthy(cos, xs...,xs...;kw1...,kw2...,foo) MyError
        @test_throws MyError @argcheck falsy(cos, xs...,xs...;kw1...,kw2...,foo) MyError
        @argcheck truthy(cos, xs...,xs...;kw1...,kw2...,foosym=>1) MyError
        @test_throws MyError @argcheck falsy(cos, xs...,xs...;kw1...,kw2...,foosym=>1) MyError
    end
end

@testset "custom message" begin
    x = 0
    @test_throws ArgumentError @argcheck x > 1 "this should not happen"
    @argcheck true "this should not happen"
end

@testset "@check" begin
    @check true
    E = CheckError
    @test_throws E @check false
    @test_throws E @check false "oh no"
    @test_throws DimensionMismatch @check false DimensionMismatch
    s = randstring()
    msg = sprint(showerror, CheckError(s))
    @test msg == "CheckError: $s"
end

@testset "pretty_string" begin
    @test pretty_string("asd") == "\"asd\""

    data = rand(10000:99999, 1000)
    str = pretty_string(data)
    @test length(str) < 1000
    @test occursin(string(last(data)), str)
    @test occursin(string(first(data)),str)
    @test !occursin("\n",str)

    data = randn()
    @test parse(Float64,pretty_string(data)) === data
end

@testset "marker" begin
    for ex in [
            :(@argcheck some_expr),
            :(@check some_expr),
            :(@check A < b),
            :(@check A < b MyError),
        ]
        ex = macroexpand(TestChecks, ex)
        @test Meta.isexpr(ex, :block)
        @test first(ex.args) == ArgCheck.LABEL_BEGIN_CHECK
        @test last(ex.args) == ArgCheck.LABEL_END_CHECK
    end
    @test ArgCheck.LABEL_BEGIN_CHECK != ArgCheck.LABEL_END_CHECK
    @test Meta.isexpr(ArgCheck.LABEL_BEGIN_CHECK, :meta)
    @test Meta.isexpr(ArgCheck.LABEL_END_CHECK, :meta)
end

end#module


text:
function number_field_over_subfield(C::ClassField_pp; using_norm_relation::Bool = true, using_stark_units::Bool = false)
  K = base_field(C)
  subs = principal_subfields(K)
  subs = [x for x in subs if degree(x[1]) != degree(K)]
  if isempty(subs)
    if using_norm_relation
      return ray_class_field_cyclic_pp_Brauer(C)
    else
      return ray_class_field_cyclic_pp(C)
    end
  end
  fl = false
  fl, C1 = translate_extension(subs[1][2], C)
  i = 1
  while !fl && i <= length(subs)-1
    i += 1
    fl, C1 = translate_extension(subs[i][2], C)
  end
  if !fl
    if using_norm_relation
      return ray_class_field_cyclic_pp_Brauer(C)
    else
      return ray_class_field_cyclic_pp(C)
    end
  end
  number_field_over_subfield(C1)
  translate_up(subs[i][2], C, C1)
  return C
end

function translate_extension(mL::NfToNfMor, C::ClassField_pp)
  L = domain(mL)
  OL = maximal_order(L)
  K = codomain(mL)
  OK = maximal_order(K)
  n = Int(exponent(C))
  ordext = Int(degree(C))
  d = divexact(discriminant(OK), discriminant(OL)^(divexact(degree(K), degree(L))))
  f = factor(ideal(OL, d))
  F = factor(ideal(OK, d))
  mR = C.rayclassgroupmap
  fM0 = copy(factored_modulus(C))
  fm0 = Dict{NfOrdIdl, Int}()
  for (p, v) in fM0
    p1 = Hecke.intersect_prime(mL, p)
    if !haskey(fm0, p1)
      if iscoprime(minimum(p1, copy = false), n)
        fm0[p1] = 1
      else
        fm0[p1] = v
      end
    end
    ep = divexact(ramification_index(p), ramification_index(p1))
    fM0[p] = ep*v
  end
  #Now, I have problems, so I need to add the ramification of the other extension.
  for (p, v) in f
    if !haskey(fm0, p)
      if isone(gcd(minimum(p), n))
        fm0[p] = 1
      else
        fm0[p] = v
      end
    else
      if !isone(gcd(minimum(p), n))
        fm0[p] = max(v, fm0[p])
      end
    end
    lPP = prime_decomposition(mL, p)
    for (P, vP) in lPP
      if haskey(fM0, P)
        fM0[P] = max(fM0[P], 2*vP*fm0[p])
      else
        fM0[P] = vP*fm0[p]*2
      end
    end
  end
  infplc = InfPlc[]
  if iszero(mod(n, 2))
    infplc = real_places(L)
  end
  @vprint :ClassField 3 "Checking if I can compute the field over a subfield\n"
  r, mr = Hecke.ray_class_group(OL, fm0, infplc, n_quo = n)
  if exponent(r) < n || order(r) < degree(C)
    return false, C
  end
  #Now, the norm group of K over L
  @vtime :ClassField 3 ngL, mngL = Hecke.norm_group(mL, mr)
  @hassert :ClassField 1 divisible(divexact(fmpz(degree(codomain(mL))), degree(domain(mL))), divexact(order(r), order(ngL)))
  if !divisible(order(ngL), degree(C)) || !divisible(exponent(C), n)
    return false, C
  end
  #Finally, the norm group of C over L
  #I use the usual strategy, as in check_abelian_extension
  for (p, v) in F
    if iscoprime(minimum(p, copy = false), degree(C))
      fM0[p] = 1
    else
      if haskey(fM0, p)
        fM0[p] = max(v, fM0[p])
      else
        fM0[p] = v
      end
    end
  end
  inf_plc2 = InfPlc[]
  if !isempty(infplc)
    inf_plc2 = real_places(K)
  end
  @vtime :ClassField 3 RM, mRM = ray_class_group(OK, fM0, inf_plc2, n_quo = n)
  @vtime :ClassField 3 lP, gS = Hecke.find_gens(mRM, coprime_to = minimum(defining_modulus(mR)[1]))
  listn = NfOrdIdl[norm(mL, x) for x in lP]
  # Create the map between R and r by taking norms
  preimgs = Vector{GrpAbFinGenElem}(undef, length(listn))
  for i = 1:length(preimgs)
    preimgs[i] = mr\listn[i]
  end
  proj = hom(gS, preimgs)
  #compute the norm group of C in RM
  prms = Vector{GrpAbFinGenElem}(undef, length(lP))
  for i = 1:length(lP)
    prms[i] = C.quotientmap(mR\lP[i])
  end
  RMtoR = hom(gS, prms)
  k, mk = kernel(RMtoR, false)
  @hassert :ClassField 1 isisomorphic(cokernel(mk, false)[1], codomain(C.quotientmap))
  mp = mk*proj
  ck, mck = cokernel(mp, false)
  #If everything could work, then ck should be the direct product of the abelian extension I am searching for and
  #the maximal abelian subextension of K/L
  G1 = snf(cokernel(mngL, false)[1])[1]
  G2 = snf(codomain(C.quotientmap))[1]
  if !has_quotient(ck, map(Int, vcat(G1.snf, G2.snf)))
    return false, C
  end
  s, ms = image(mngL*mck, false)
  fl, ms1 = has_complement(ms)
  @assert fl
  mq1 = cokernel(ms1, false)[2]
  mqq = mck * mq1
  @hassert :ClassField 1 domain(mqq) == r
  C1 = ClassField_pp{MapRayClassGrp, GrpAbFinGenMap}()
  C1.quotientmap = mqq
  C1.rayclassgroupmap = mr
  return true, C1
end

function translate_up(mL::NfToNfMor, C::ClassField_pp, C1::ClassField_pp)
  K = base_field(C)
  Ky = PolynomialRing(K, "y", cached = false)[1]
  L = domain(mL)
  d = degree(C1)
  CEK = cyclotomic_extension(K, d)
  CEL = cyclotomic_extension(L, d)
  img = gen(CEK.Kr)
  if degree(CEK.Kr) != euler_phi(d)
    pp = map_coefficients(mL, CEL.Kr.pol, cached = false)
    while !iszero(pp(img))
      mul!(img, img, gen(CEK.Kr))
    end
  end
  mrel = hom(CEL.Kr, CEK.Kr, mL, img)
  #@hassert :Fields 1 isconsistent(mrel)
  g = mrel(CEL.mp[1](gen(CEL.Ka)))
  mp = hom(CEL.Ka, CEK.Ka, CEK.mp[1]\(g), check = false)
  #Then, the fac elem corresponding to the generator of the Kummer Extension
  C.a = FacElem(Dict{nf_elem, fmpz}(mp(x) => v for (x, v) in C1.a))
  #Now, the Kummer extension
  Lzeta = codomain(mp)
  Lt = PolynomialRing(Lzeta, "t", cached = false)[1]
  d1 = degree(C1.K)
  coeffs = Vector{nf_elem}(undef, d1 + 1)
  coeffs[1] = mp(coeff(C1.K.pol, 0))
  for s = 2:length(coeffs)-1
    coeffs[s] = zero(Lzeta)
  end
  coeffs[end] = one(Lzeta)
  C.K = number_field(Lt(coeffs), cached = false, check = false)[1]
  #The target extension
  fdef = map_coefficients(mL, C1.A.pol, parent = Ky, cached = false)
  C.A = number_field(fdef, cached = false, check = false)[1]
  #Now, the primitive element of the target extension seen in Cpp.K
  mrel2 = hom(C1.K, C.K, mp, gen(C.K))
  C.pe = mrel2(C1.pe)
  CEKK = cyclotomic_extension(K, d)
  @hassert :ClassField 1 iszero(map_coefficients(CEKK.mp[2], fdef, cached = false)(C.pe))
  C.o = d1
  return nothing
end


text:
# type StepReg{W<:Weight} <: OnlineStat
#     β0::Float64
#     β::VecF
#     intercept::Bool
#     c::CovMatrix{W}
#     s::MatF         # placeholder for swept version of c
#     set::IntSet
#     weight::W
#     n::Int
#     nups::Int
# end
# function StepReg(p::Integer, wgt::Weight = EqualWeight(); intercept::Bool = true)
#     StepReg(
#         0.0, zeros(p), intercept, CovMatrix(p), zeros(p, p), IntSet([]),
#         wgt, 0, 0
#     )
# end
# function StepReg(x::AMat, y::AVec, wgt::Weight = EqualWeight())
#     o = StepReg(size(x, 2), wgt)
#     fit!(o, x, y)
#     o
# end
#
# function StatsBase.fit!(o::StepReg, x::AVec, y::Real)
#     fit!(o.c, vcat(x, y))
# end
# function value(o::StepReg)
#     o.intercept ? vcat(o.β0, o.β) : o.β
# end
#
# ## TESTING
# n, p = 1000, 10
# x = randn(n, p)
# β = collect(1.:p)
# y = x * β + randn(n)
# o = StepReg(x, y)
# display(o)
#
# # # Online Stepwise Regression
# # #
# # # With each updated batch, there is the possibility of one variable
# # # entering or leaving the model based on Mallow's Cp
# #
# # #-------------------------------------------------------# Type and Constructors
# # """
# # Online stepwise regression.
# #
# # At each update, there is the possible of one variable entering or leaving the model.
# # Performs best using `fitbatch!()` with large batches.
# # """
# # type StepwiseReg{W <: Weight} <: OnlineStat
# #     C::CovarianceMatrix{W}  # Cov([X y])
# #     s::MatF                 # "Swept" version of [X y]' [X y] / n
# #     set::IntSet             # set of coefficients included in the model
# #     n::Int
# # end
# #
# # function StepwiseReg(x::AMatF, y::AVecF, wgt::Weight = EqualWeight())
# #     n, p = size(x)
# #     o = StepwiseReg(p, wgt)
# #     fitbatch!(o, x, y)
# #     o
# # end
# #
# # function StepwiseReg(p::Integer, wgt::Weighting = default(Weighting))
# #     c = CovarianceMatrix(p + 1, wgt)
# #     StepwiseReg(c, zeros(p + 1, p + 1), IntSet([]), 0)
# # end
# #
# #
# # #-----------------------------------------------------------------------# state
# # statenames(o::StepwiseReg) = [:β, :nobs]
# # state(o::StepwiseReg) = Any[coef(o), nobs(o)]
# #
# # function StatsBase.coef(o::StepwiseReg)
# #     β = vec(o.s[end, 1:end - 1])
# #     for i in setdiff(1:length(β), o.set)
# #         β[i] = 0.
# #     end
# #     β
# # end
# #
# # #----------------------------------------------------------------------# fit!
# # function fit!(o::StepwiseReg, x::AVecF, y::Float64)
# #     fitbatch!(o, x', collect(y))
# # end
# #
# # function fitbatch!(o::StepwiseReg, x::AMatF, y::AVecF)
# #     n, p = size(x)
# #     o.n += n
# #     fitbatch!(o.C, hcat(x, y))
# #     copy!(o.s, o.C.A)
# #
# #     # get average squared error using all predictors
# #     sweep!(o.s, 1:p)
# #     ase = o.s[end, end]
# #
# #     copy!(o.s, o.C.A)
# #     sweep!(o.s, collect(o.set))
# #
# #     # Find best index to add/remove
# #     s = 1
# #     val = o.s[end, 1] ^ 2 / (o.s[1, 1] * ase) * (o.n - in(1, o.set)) + 2.0 - in(1, o.set)
# #     for i in 2:p
# #         newval = o.s[end, i] ^ 2 / (o.s[i, i] * ase) * (o.n - in(i, o.set)) - 2.0 * in(i, o.set)
# #         if newval > val
# #             val = newval
# #             s = i
# #         end
# #     end
# #
# #     if s in o.set
# #         delete!(o.set, s)
# #     else
# #         push!(o.set, s)
# #     end
# #
# #     DEBUG("Active set: ", o.set)
# # end
# #
# #
# # ########## TEST code
# # # log_severity!(DebugSeverity)
# # #
# # # n,p = 10_000, 10
# # # x = randn(n,p)
# # # β = collect(1.:p)
# # # β[5] = 0.
# # # y = x*β + randn(n)
# # #
# # # o = StepwiseReg(p)
# # # @time fit!(o, x, y, 500)
# # # print(coef(o))


text:
"""
	$(SIGNATURES)

Container for a collection of `ModelStats`. May also contain nested `StatsCollection`s.
"""
struct StatsCollection{T}
    meta :: Dict{Symbol, Any}
    d :: Dict{Any, ModelStats}
end

## -----------  Constructor

StatsCollection{T}(meta :: Dict{Symbol, Any}) where T =
    StatsCollection{T}(meta,  Dict{Any, ModelStats}());

StatsCollection{T}() where T = StatsCollection{T}(Dict{Symbol,Any}());

get_meta(s :: StatsCollection{T}, mKey) where T = s.meta[mKey];

## -----------  Adding stats

"""
	$(SIGNATURES)

Add `ModelStats` or a nested `StatsCollection`.
"""
function add_mstats!(m :: StatsCollection{T}, s :: ModelStats{T2}) where {T, T2}
    @assert !has_mstats(m, T2) "$s already exists in $m";
    m.d[T2] = s;
    return nothing;
end

"""
	$(SIGNATURES)

Delete `ModelStats`.
"""
function delete_mstats!(m :: StatsCollection{T}, msName) where {T}
    @assert has_mstats(m, msName) "$msName not found in $m";
    delete!(m.d, msName);
    return nothing;
end

"""
	$(SIGNATURES)

Replace `ModelStats`. Add it if it does not exist.
"""
function replace_mstats!(m :: StatsCollection{T}, s :: ModelStats{T2}) where {T, T2}
    m.d[T2] = s;
    return nothing;
end

## -----------  Retrieve

"""
	$(SIGNATURES)

List of all the `groups`. That is, `m` contains `ModelStats{grp}` for each `grp in get_groups(m)`.
"""
get_groups(m :: StatsCollection{T}) where T = keys(m.d);

"""
	$(SIGNATURES)

Does the `ModelStats{groups}` object exist?
"""
function has_mstats(m :: StatsCollection{T}, groups) where T
    return haskey(m.d, groups);
end

"""
	$(SIGNATURES)

Retrieve the `ModelStats{groups}` object. Errors if `groups` not found.
"""
function get_mstats(m :: StatsCollection{T}, groups) where T
    @assert has_mstats(m, groups)  "$groups not found in $m"
    ms = m.d[groups];
    return ms
end

## ----------  Directly retrieve individual statistics

"""
	$(SIGNATURES)

Does `StatsCollection` contain variable `statName`?
"""
has_variable(m :: StatsCollection{T}, groups, statName :: Symbol) where T =
    has_variable(get_mstats(m, groups), statName);

"""
	$(SIGNATURES)

Retrieve a statistic. Groups are specified.
"""
get_stats(m :: StatsCollection{T}, groups, statName :: Symbol;
    defaultValue = :error) where T =
	get_values(get_mstats(m, groups), statName; defaultValue = defaultValue);

"""
	$(SIGNATURES)

Retrieve a statistic without providing groups.

# Example
```julia
get_stats(m, :workTime_gp) == get_stats(m, (:gpa, :parental), :workTime)
# Equivalently
m.workTime_gp == get_stats(m, :workTime_gp)
m.workTime_gp[2,3] == get_stats(m, :workTime_gp)[2,3]
```
"""
function get_stats(m :: StatsCollection{T}, sg :: Symbol;
    defaultValue = :error) where T
    return get_stats(
        m, groups_from_name(m, sg), statname_from_name(m, sg);
        defaultValue = defaultValue)
end

function Base.getproperty(m :: StatsCollection{T}, sg :: Symbol) where T
    if sg == :d
        return getfield(m, :d)
    elseif sg == :meta
        return getfield(m, :meta)
    else
        return get_stats(m, sg);
    end
end

Base.getindex(m :: StatsCollection{T}, sg :: Symbol, idx...) where T =
    getindex(get_stats(m, sg), idx...);

"""
	$(SIGNATURES)

Set values for a statistic.

# Example
```julia
set_stats!(m, (:gpa, :parental), [1,2,3])
set_stats!(m, :workTime_gp, [1,2,3])
# Also with indices
m.workTime_gp[2,3] = x
```
"""
function set_stats!(m :: StatsCollection{T}, groups, statName :: Symbol, x) where T
	setproperty!(get_mstats(m, groups), statName, x);
end

set_stats!(m :: StatsCollection{T}, sg :: Symbol, x) where T =
    set_stats!(m, groups_from_name(m, sg), statname_from_name(m, sg), x);

# `sc.x1_gp = newValues`
Base.setproperty!(m :: StatsCollection{T}, sg :: Symbol, x) where T =
    set_stats!(m, sg, x);
Base.setproperty!(m :: StatsCollection{T}, x :: ModelStats{T2}) where {T, T2} =
    replace_mstats!(m, x);

Base.setindex!(m :: StatsCollection{T}, sg :: Symbol, x, idx...) where T =
    setindex!(get_stats(m, sg), x, idx...);

"""
	$(SIGNATURES)

Mapping of name of a statistic into groups. User defines this for each parametric `StatsCollection{T}`.

# Example:
```julia
groups_from_name(m, :workTime_gp) == (:gpa, :parental)
```
"""
groups_from_name(m :: StatsCollection, statName :: Symbol) =
    error("User defines this for each parametric StatsCollection");

"""
	$(SIGNATURES)

Mapping of name of a statistic into its name without the group info. User defines this for each parametric `StatsCollection{T}`.

# Example:
```julia
groups_from_name(m, :workTime_gp) == :workTime
```
"""
statname_from_name(m :: StatsCollection, statName :: Symbol) =
    error("User defines this for each parametric StatsCollection");

## -------  Testing

function make_test_scollection(T, rng)
    sc = StatsCollection{T}();
    tmsV = (:TMS, :TMS2, :TMS3);
    nVarV = (3, 4, 5);
    for j = 1 : length(tmsV)
        add_mstats!(sc, make_test_model_stats(tmsV[j], nVarV[j], rng));
    end
    return sc
end

# ------------


text:
include("abstract_cfr_policy.jl")
include("tabular_cfr.jl")
include("outcome_sampling_mccfr.jl")
include("external_sampling_mccfr.jl")
include("best_response_policy.jl")
include("nash_conv.jl")
include("deep_cfr.jl")