text:
@testset "vectorization of 2-arg functions" begin
    binary_math_functions = [
        besselh, hankelh1, hankelh2, hankelh1x, hankelh2x,
        besseli, besselix, besselj, besseljx, besselk, besselkx, bessely, besselyx,
        polygamma, zeta
    ]
    @testset "$f" for f in binary_math_functions
        x = y = 2
        v = [f(x,y)]
        @test f.([x],y) == v
        @test f.(x,[y]) == v
        @test f.([x],[y]) == v
    end
end

@testset "MPFR" begin
    @testset "bessel functions" begin
        setprecision(53) do
            @test besselj(4, BigFloat(2)) ≈ besselj(4, 2.)
            @test besselj0(BigFloat(2)) ≈ besselj0(2.)
            @test besselj1(BigFloat(2)) ≈ besselj1(2.)
            @test bessely(4, BigFloat(2)) ≈ bessely(4, 2.)
            @test bessely0(BigFloat(2)) ≈ bessely0(2.)
            @test bessely1(BigFloat(2)) ≈ bessely1(2.)
        end
    end

    let err(z, x) = abs(z - x) / abs(x)
        @test 1e-60 > err(eta(parse(BigFloat,"1.005")), parse(BigFloat,"0.693945708117842473436705502427198307157819636785324430166786"))
        @test 1e-60 > err(exp(eta(big(1.0))), 2.0)
    end

    let a = parse(BigInt, "315135")
        @test typeof(erf(a)) == BigFloat
        @test typeof(erfc(a)) == BigFloat
    end

    # issue #101
    for i in 0:5
        @test gamma(big(i)) == gamma(i)
    end
end

@testset "Base Julia issue #17474" begin
    @test f64(complex(1f0,1f0)) === complex(1.0, 1.0)
    @test f64(1f0) === 1.0

    @test typeof(eta(big"2")) == BigFloat
    @test typeof(zeta(big"2")) == BigFloat
    @test typeof(digamma(big"2")) == BigFloat

    @test_throws MethodError trigamma(big"2")
    @test_throws MethodError trigamma(big"2.0")
    @test_throws MethodError invdigamma(big"2")
    @test_throws MethodError invdigamma(big"2.0")

    @test_throws MethodError eta(Complex(big"2"))
    @test_throws MethodError eta(Complex(big"2.0"))
    @test_throws MethodError zeta(Complex(big"2"))
    @test_throws MethodError zeta(Complex(big"2.0"))
    @test_throws MethodError zeta(1.0,big"2")
    @test_throws MethodError zeta(1.0,big"2.0")
    @test_throws MethodError zeta(big"1.0",2.0)
    @test_throws MethodError zeta(big"1",2.0)

    @test typeof(polygamma(3, 0x2)) == Float64
    @test typeof(polygamma(big"3", 2f0)) == Float32
    @test typeof(zeta(1, 2.0)) == Float64
    @test typeof(zeta(1, 2f0)) == Float64 # BitIntegers result in Float64 returns
    @test typeof(zeta(2f0, complex(2f0,0f0))) == Complex{Float32}
    @test typeof(zeta(complex(1,1), 2f0)) == Complex{Float64}
    @test typeof(zeta(complex(1), 2.0)) == Complex{Float64}
end

@test sprint(showerror, AmosException(1)) == "AmosException with id 1: input error."

@testset "missing data" begin
    for f in (digamma, erf, erfc, erfcinv, erfcx, erfi, erfinv, eta, gamma,
              invdigamma, logfactorial, trigamma)
        @test f(missing) === missing
    end
    @test beta(1.0, missing) === missing
    @test beta(missing, 1.0) === missing
    @test beta(missing, missing) === missing
    @test polygamma(4, missing) === missing
end


text:
#==============================================================================#
# XRay.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/xray-2016-04-12.normal.json
#==============================================================================#

__precompile__()

module XRay

using AWSCore

"""
    using AWSSDK.XRay.batch_get_traces
    batch_get_traces([::AWSConfig], arguments::Dict)
    batch_get_traces([::AWSConfig]; TraceIds=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/Traces", arguments::Dict)
    xray([::AWSConfig], "POST", "/Traces", TraceIds=, <keyword arguments>)

# BatchGetTraces Operation

Retrieves a list of traces specified by ID. Each trace is a collection of segment documents that originates from a single request. Use `GetTraceSummaries` to get a list of trace IDs.

# Arguments

## `TraceIds = [::String, ...]` -- *Required*
Specify the trace IDs of requests for which to retrieve segments.

## `NextToken = ::String`
Pagination token. Not used.

# Returns

`BatchGetTracesResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/BatchGetTraces)
"""
@inline batch_get_traces(aws::AWSConfig=default_aws_config(); args...) = batch_get_traces(aws, args)

@inline batch_get_traces(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/Traces", args)

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

"""
    using AWSSDK.XRay.get_encryption_config
    get_encryption_config([::AWSConfig], arguments::Dict)
    get_encryption_config([::AWSConfig]; )

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/EncryptionConfig", arguments::Dict)
    xray([::AWSConfig], "POST", "/EncryptionConfig", )

# GetEncryptionConfig Operation

Retrieves the current encryption configuration for X-Ray data.

# Arguments

# Returns

`GetEncryptionConfigResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/GetEncryptionConfig)
"""
@inline get_encryption_config(aws::AWSConfig=default_aws_config(); args...) = get_encryption_config(aws, args)

@inline get_encryption_config(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/EncryptionConfig", args)

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

"""
    using AWSSDK.XRay.get_service_graph
    get_service_graph([::AWSConfig], arguments::Dict)
    get_service_graph([::AWSConfig]; StartTime=, EndTime=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/ServiceGraph", arguments::Dict)
    xray([::AWSConfig], "POST", "/ServiceGraph", StartTime=, EndTime=, <keyword arguments>)

# GetServiceGraph Operation

Retrieves a document that describes services that process incoming requests, and downstream services that they call as a result. Root services process incoming requests and make calls to downstream services. Root services are applications that use the AWS X-Ray SDK. Downstream services can be other applications, AWS resources, HTTP web APIs, or SQL databases.

# Arguments

## `StartTime = timestamp` -- *Required*
The start of the time frame for which to generate a graph.

## `EndTime = timestamp` -- *Required*
The end of the time frame for which to generate a graph.

## `NextToken = ::String`
Pagination token. Not used.

# Returns

`GetServiceGraphResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/GetServiceGraph)
"""
@inline get_service_graph(aws::AWSConfig=default_aws_config(); args...) = get_service_graph(aws, args)

@inline get_service_graph(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/ServiceGraph", args)

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

"""
    using AWSSDK.XRay.get_trace_graph
    get_trace_graph([::AWSConfig], arguments::Dict)
    get_trace_graph([::AWSConfig]; TraceIds=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/TraceGraph", arguments::Dict)
    xray([::AWSConfig], "POST", "/TraceGraph", TraceIds=, <keyword arguments>)

# GetTraceGraph Operation

Retrieves a service graph for one or more specific trace IDs.

# Arguments

## `TraceIds = [::String, ...]` -- *Required*
Trace IDs of requests for which to generate a service graph.

## `NextToken = ::String`
Pagination token. Not used.

# Returns

`GetTraceGraphResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/GetTraceGraph)
"""
@inline get_trace_graph(aws::AWSConfig=default_aws_config(); args...) = get_trace_graph(aws, args)

@inline get_trace_graph(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/TraceGraph", args)

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

"""
    using AWSSDK.XRay.get_trace_summaries
    get_trace_summaries([::AWSConfig], arguments::Dict)
    get_trace_summaries([::AWSConfig]; StartTime=, EndTime=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/TraceSummaries", arguments::Dict)
    xray([::AWSConfig], "POST", "/TraceSummaries", StartTime=, EndTime=, <keyword arguments>)

# GetTraceSummaries Operation

Retrieves IDs and metadata for traces available for a specified time frame using an optional filter. To get the full traces, pass the trace IDs to `BatchGetTraces`.

A filter expression can target traced requests that hit specific service nodes or edges, have errors, or come from a known user. For example, the following filter expression targets traces that pass through `api.example.com`:

`service("api.example.com")`

This filter expression finds traces that have an annotation named `account` with the value `12345`:

`annotation.account = "12345"`

For a full list of indexed fields and keywords that you can use in filter expressions, see [Using Filter Expressions](http://docs.aws.amazon.com/xray/latest/devguide/xray-console-filters.html) in the *AWS X-Ray Developer Guide*.

# Arguments

## `StartTime = timestamp` -- *Required*
The start of the time frame for which to retrieve traces.

## `EndTime = timestamp` -- *Required*
The end of the time frame for which to retrieve traces.

## `Sampling = ::Bool`
Set to `true` to get summaries for only a subset of available traces.

## `FilterExpression = ::String`
Specify a filter expression to retrieve trace summaries for services or requests that meet certain requirements.

## `NextToken = ::String`
Specify the pagination token returned by a previous request to retrieve the next page of results.

# Returns

`GetTraceSummariesResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/GetTraceSummaries)
"""
@inline get_trace_summaries(aws::AWSConfig=default_aws_config(); args...) = get_trace_summaries(aws, args)

@inline get_trace_summaries(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/TraceSummaries", args)

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

"""
    using AWSSDK.XRay.put_encryption_config
    put_encryption_config([::AWSConfig], arguments::Dict)
    put_encryption_config([::AWSConfig]; Type=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/PutEncryptionConfig", arguments::Dict)
    xray([::AWSConfig], "POST", "/PutEncryptionConfig", Type=, <keyword arguments>)

# PutEncryptionConfig Operation

Updates the encryption configuration for X-Ray data.

# Arguments

## `KeyId = ::String`
An AWS KMS customer master key (CMK) in one of the following formats:

*   **Alias** - The name of the key. For example, `alias/MyKey`.

*   **Key ID** - The KMS key ID of the key. For example, `ae4aa6d49-a4d8-9df9-a475-4ff6d7898456`.

*   **ARN** - The full Amazon Resource Name of the key ID or alias. For example, `arn:aws:kms:us-east-2:123456789012:key/ae4aa6d49-a4d8-9df9-a475-4ff6d7898456`. Use this format to specify a key in a different account.

Omit this key if you set `Type` to `NONE`.

## `Type = "NONE" or "KMS"` -- *Required*
The type of encryption. Set to `KMS` to use your own key for encryption. Set to `NONE` for default encryption.

# Returns

`PutEncryptionConfigResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/PutEncryptionConfig)
"""
@inline put_encryption_config(aws::AWSConfig=default_aws_config(); args...) = put_encryption_config(aws, args)

@inline put_encryption_config(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/PutEncryptionConfig", args)

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

"""
    using AWSSDK.XRay.put_telemetry_records
    put_telemetry_records([::AWSConfig], arguments::Dict)
    put_telemetry_records([::AWSConfig]; TelemetryRecords=, <keyword arguments>)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/TelemetryRecords", arguments::Dict)
    xray([::AWSConfig], "POST", "/TelemetryRecords", TelemetryRecords=, <keyword arguments>)

# PutTelemetryRecords Operation

Used by the AWS X-Ray daemon to upload telemetry.

# Arguments

## `TelemetryRecords = [[ ... ], ...]` -- *Required*

```
 TelemetryRecords = [[
        "Timestamp" => <required> timestamp,
        "SegmentsReceivedCount" =>  ::Int,
        "SegmentsSentCount" =>  ::Int,
        "SegmentsSpilloverCount" =>  ::Int,
        "SegmentsRejectedCount" =>  ::Int,
        "BackendConnectionErrors" =>  [
            "TimeoutCount" =>  ::Int,
            "ConnectionRefusedCount" =>  ::Int,
            "HTTPCode4XXCount" =>  ::Int,
            "HTTPCode5XXCount" =>  ::Int,
            "UnknownHostCount" =>  ::Int,
            "OtherCount" =>  ::Int
        ]
    ], ...]
```

## `EC2InstanceId = ::String`

## `Hostname = ::String`

## `ResourceARN = ::String`

# Returns

`PutTelemetryRecordsResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/PutTelemetryRecords)
"""
@inline put_telemetry_records(aws::AWSConfig=default_aws_config(); args...) = put_telemetry_records(aws, args)

@inline put_telemetry_records(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/TelemetryRecords", args)

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

"""
    using AWSSDK.XRay.put_trace_segments
    put_trace_segments([::AWSConfig], arguments::Dict)
    put_trace_segments([::AWSConfig]; TraceSegmentDocuments=)

    using AWSCore.Services.xray
    xray([::AWSConfig], "POST", "/TraceSegments", arguments::Dict)
    xray([::AWSConfig], "POST", "/TraceSegments", TraceSegmentDocuments=)

# PutTraceSegments Operation

Uploads segment documents to AWS X-Ray. The X-Ray SDK generates segment documents and sends them to the X-Ray daemon, which uploads them in batches. A segment document can be a completed segment, an in-progress segment, or an array of subsegments.

Segments must include the following fields. For the full segment document schema, see [AWS X-Ray Segment Documents](https://docs.aws.amazon.com/xray/latest/devguide/xray-api-segmentdocuments.html) in the *AWS X-Ray Developer Guide*.

**Required Segment Document Fields**

*   `name` - The name of the service that handled the request.

*   `id` - A 64-bit identifier for the segment, unique among segments in the same trace, in 16 hexadecimal digits.

*   `trace_id` - A unique identifier that connects all segments and subsegments originating from a single client request.

*   `start_time` - Time the segment or subsegment was created, in floating point seconds in epoch time, accurate to milliseconds. For example, `1480615200.010` or `1.480615200010E9`.

*   `end_time` - Time the segment or subsegment was closed. For example, `1480615200.090` or `1.480615200090E9`. Specify either an `end_time` or `in_progress`.

*   `in_progress` - Set to `true` instead of specifying an `end_time` to record that a segment has been started, but is not complete. Send an in progress segment when your application receives a request that will take a long time to serve, to trace the fact that the request was received. When the response is sent, send the complete segment to overwrite the in-progress segment.

A `trace_id` consists of three numbers separated by hyphens. For example, 1-58406520-a006649127e371903a2de979. This includes:

**Trace ID Format**

*   The version number, i.e. `1`.

*   The time of the original request, in Unix epoch time, in 8 hexadecimal digits. For example, 10:00AM December 2nd, 2016 PST in epoch time is `1480615200` seconds, or `58406520` in hexadecimal.

*   A 96-bit identifier for the trace, globally unique, in 24 hexadecimal digits.

# Arguments

## `TraceSegmentDocuments = [::String, ...]` -- *Required*
A string containing a JSON document defining one or more segments or subsegments.

# Returns

`PutTraceSegmentsResult`

# Exceptions

`InvalidRequestException` or `ThrottledException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/xray-2016-04-12/PutTraceSegments)
"""
@inline put_trace_segments(aws::AWSConfig=default_aws_config(); args...) = put_trace_segments(aws, args)

@inline put_trace_segments(aws::AWSConfig, args) = AWSCore.Services.xray(aws, "POST", "/TraceSegments", args)

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

end # module XRay

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

ext: jl


text:
export Exact,
       Center,
       Subpixel

abstract type OverlapMethod end

function circular_overlap(xmin, xmax, ymin, ymax, nx, ny, r; method = :exact)
    out = fill(0.0, nx, ny)

    # width of each element
    dx = (xmax - xmin) / nx
    dy = (ymax - ymin) / ny

    # radius of one pixel
    pixel_radius = 0.5sqrt(dx^2 + dy^2)

    # bounding box
    bxmin = -r - 0.5dx
    bxmax = r + 0.5dx
    bymin = -r - 0.5dy
    bymax = r + 0.5dy

    for i in 1:nx
        # lower end of pixel
        pxmin = xmin + (i - 1) * dx
        pxcen = pxmin + 0.5dx
        # upper end of pixel
        pxmax = pxmin + dx

        if pxmax > bxmin && pxmin < bxmax
            for j in 1:ny
                pymin = ymin + (j - 1) * dy
                pycen = pymin + 0.5dy
                pymax = pymin + dy
                if pymax > bymin && pymin < bymax

                    # distance from circle to pixel
                    d = sqrt(pxcen^2 + pycen^2)

                    # fully within radius
                    if d < r - pixel_radius
                        @inbounds out[j, i] = 1
                    # partially within radius
                    elseif d < r + pixel_radius
                        if method === :exact
                            @inbounds out[j, i] = circular_overlap_single_exact(pxmin, pymin, pxmax, pymax, r) / (dx * dy)
                        elseif method === :center
                            @inbounds out[j, i] =  circular_overlap_single_subpixel(pxmin, pymin, pxmax, pymax, r, 1)
                        elseif method[1] === :subpixel
                            @inbounds out[j, i] =  circular_overlap_single_subpixel(pxmin, pymin, pxmax, pymax, r, method[2])
                        end
                    end
                end
            end
        end

    end
    return out
end

function circular_overlap_single_subpixel(xmin, ymin, xmax, ymax, r, subpixels)
    frac = 0
    dx = (xmax - xmin) / subpixels
    dy = (ymax - ymin) / subpixels
    r2 = r^2

    x = xmin - 0.5dx
    for i in 1:subpixels
        x += dx
        y = ymin - 0.5dy
        for j in 1:subpixels
            y += dy
            if x^2 + y^2 < r2
                frac += 1
            end
        end
    end
    return frac / subpixels^2
end

"""Area of overlap between a rectangle and a circle"""
function circular_overlap_single_exact(xmin, ymin, xmax, ymax, r)
    r ≤ 0 && return 0
    if 0 ≤ xmin
        0 ≤ ymin && return circular_overlap_core(xmin, ymin, xmax, ymax, r)
        0 ≥ ymax && return circular_overlap_core(-ymax, xmin, -ymin, xmax, r)
        return (circular_overlap_single_exact(xmin, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, xmax, ymax, r))
    elseif 0 ≥ xmax
        0 ≤ ymin && return circular_overlap_core(-xmax, ymin, -xmin, ymax, r)
        0 ≥ ymax && return circular_overlap_core(-xmax, -ymax, -xmin, -ymin, r)
        return (circular_overlap_single_exact(xmin, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, xmax, ymax, r))
    else
        0 ≤ ymin && return (circular_overlap_single_exact(xmin, ymin, 0, ymax, r) +
                            circular_overlap_single_exact(0, ymin, xmax, ymax, r))
        0 ≥ ymax && return (circular_overlap_single_exact(xmin, ymin, 0, ymax, r) +
                            circular_overlap_single_exact(0, ymin, xmax, ymax, r))
        return (circular_overlap_single_exact(xmin, ymin, 0, 0, r) +
                circular_overlap_single_exact(0, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, 0, ymax, r) +
                circular_overlap_single_exact(0, 0, xmax, ymax, r))
    end
end

"""Core of circular overlap routine"""
function circular_overlap_core(xmin, ymin, xmax, ymax, r)
    xmin^2 + ymin^2 > r^2 && return 0
    xmax^2 + ymax^2 < r^2 && return (xmax - xmin) * (ymax - ymin)

    d1 = sqrt(xmax^2 + ymin^2)
    d2 = sqrt(xmin^2 + ymax^2)
    if d1 < r && d2 < r
        x1, y1 = sqrt(r^2 - ymax^2), ymax
        x2, y2 = xmax, sqrt(r^2 - xmax^2)
        return ((xmax - xmin) * (ymax - ymin) -
                area_triangle(x1, y1, x2, y2, xmax, ymax) +
                area_arc(x1, y1, x2, y2, r))
    elseif d1 < r
        x1, y1 = xmin, sqrt(r^2 - xmin^2)
        x2, y2 = xmax, sqrt(r^2 - xmax^2)
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x1, ymin, xmax, ymin) +
                area_triangle(x1, y1, x2, ymin, x2, y2))
    elseif d2 < r
        x1, y1 = sqrt(r^2 - ymin^2), ymin
        x2, y2 = sqrt(r^2 - ymax^2), ymax
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, xmin, y1, xmin, ymax) +
                area_triangle(x1, y1, xmin, y2, x2, y2))
    else
        x1, y1 = sqrt(r^2 - ymin^2), ymin
        x2, y2 = xmin, sqrt(r^2 - xmin^2)
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x2, y2, xmin, ymin))
    end
end

####################################
# Some geometric helpers

"""Area of a triangle defined by three vertices"""
area_triangle(x0, y0, x1, y1, x2, y2) = 0.5abs(x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1))

"""
Area of a circular segment above a chord between two points with circle radius `r`

[Reference](http://mathworld.wolfram.com/CircularSegment.html)
"""
function area_arc(x0, y0, x1, y1, r)
    a = sqrt((x1 - x0)^2 + (y1 - y0)^2)
    θ = 2asin(0.5a / r)
    return 0.5r^2 * (θ - sin(θ))
end


text:
type ADAM
    decay1  :: Float64
    decay2  :: Float64
    epsilon :: Float64
    alpha   :: Float64
    weight1 :: Float64
    weight2 :: Float64
    moment1 :: Float64
    moment2 :: Float64
end

ADAM(;decay1=0.9, decay2=0.999, epsilon=1e-8, alpha=0.01) =
        ADAM(decay1, decay2, epsilon, alpha, decay1, decay2, 0.0, 0.0)

function gradient_step!(adam, gradient)
    adam.moment1 = adam.decay1 * adam.moment1 + (1 - adam.decay1) * gradient
    adam.moment2 = adam.decay2 * adam.moment2 + (1 - adam.decay2) * gradient^2

    corrected1 = adam.moment1/(1 - adam.weight1)
    corrected2 = adam.moment2/(1 - adam.weight2)

    adam.weight1 *= adam.decay1
    adam.weight2 *= adam.decay2

    delta  = corrected1 / ( sqrt(corrected2) + adam.epsilon)

    return adam.alpha * delta
end

max_stars_repo_name: adityam/periodic-sensor-scheduling

max_issues_repo_name: adityam/periodic-sensor-scheduling

max_forks_repo_name: adityam/periodic-sensor-scheduling


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using .DataFiles

# =========================================================================================
"""
    download_languages(
        languages::AbstractString = "eng";
        target::AbstractString = "tessdata",
        baseUrl::AbstractString = DATA_URL,
        force::Bool = false
    )::Bool

Download the data files for the specified languages.  Returns `false` if there is a problem
downloading the files.

__Arguments:__

| T | Name      | Default    | Description
|:--| :-------- | :--------- | :----------
| O | languages | `eng`      | The languages to download separated with "+".
| K | target    | `tessdata` | The directory to save the data files in.
| K | baseUrl   | ...        | The base URL to download the files from.
| K | force     | `false`    | Should the files be downloaded even if they exist?

__Details:__

By default the data files are saved in a "tessdata" directory under the current directory,
and the files are downloaded from https://github.com/tesseract-ocr/tessdata_best.

Normally if the file has already been downloaded then it is not downloaded again.  However
if `force` is true then the file is downloaded and the existing file is overwritten.

See also: [`update_languages`](@ref)
"""
function download_languages(
            languages::AbstractString = "eng";
            target::AbstractString = "tessdata",
            baseUrl::AbstractString = DATA_URL,
            force::Bool = false
        )::Bool
    local result = false
    local index = load_languages(target)

    try
        for lang in split(languages, "+")
            local filename = string(lang, DATA_FILE_EXT)
            local fullpath = joinpath(target, filename)
            local file     = get_language_file(index, lang)

            if force == true || ispath(fullpath) == false
                local url  = string(baseUrl, filename)
                local sha1 = download_file(url, fullpath)

                file.localsha1  = sha1
                file.timestamp  = mtime(fullpath)
            end
        end

        result = true
    catch ex
        if isa(ex, FileError) == false && isa(ex, NetworkError) == false
            rethrow(ex)
        end
    end
    save_languages(index, target)

    return result
end

# =================================================================================================
"""
    download_pdf_font(;
        target::AbstractString = "tessdata"
        baseUrl::AbstractString = PDF_FONT_URL,
        force::Bool = false
    )::Bool

Download the PDF font file.  Returns `false` if there is a problem downloading the files.

__Arguments:__

| T | Name    | Default    | Description
|:--| :------ | :--------- | :----------
| K | target  | `tessdata` | The directory to save the PDF font to.
| K | baseUrl | ...        | The base URL to download the file from.
| K | force   | `false`    | Should the file be downloaded even it it exists?

__Details:__

By default "pdf.ttf" is downloaded from https://github.com/tesseract-ocr/tessconfigs.
The client can specify a different URL to download the file from.  This file is needed when
asking Tesseract to generate a PDF.

Normally if the file has already been downloaded then it is not downloaded again.  However
if `force` is true then the file is downloaded and the existing file is overwritten.

See also: [`update_pdf_font`](@ref)
"""
function download_pdf_font(;
            target::AbstractString = "tessdata",
            baseUrl::AbstractString = PDF_FONT_URL,
            force::Bool = false
        )::Bool
    local result = false
    local index = load_pdf_font_info(target)

    try
        local fullpath = joinpath(target, "pdf.ttf")

        if force == true || ispath(fullpath) == false
            local fullurl = string(baseUrl, "pdf.ttf")
            local sha1    = download_file(fullurl, fullpath)

            index.font.localsha1  = sha1
            index.font.timestamp  = mtime(fullpath)
        end

        result = true
    catch ex
        if isa(ex, FileError) == false && isa(ex, NetworkError) == false
            rethrow(ex)
        end
    end

    save_pdf_font_info(index, target)

    return result
end


text:
# This file was generated, do not modify it. # hide
data2 = coerce(data, Count => Continuous)
schema(data2)


text:
# Functions for representing graphs with GraphPlot package
# https://github.com/JuliaGraphs/GraphPlot.jl

function plot_network(netw1::SparseMatrixCSC;
    layout = spring_layout,
    label = true,
    edgelinewidth = 1.0,
    color = Colors.parse(Colorant, "turquoise")
)

    directed = !is_symmetric(netw1)

    n = netw1.m
    swn = SimpleWeightedDiGraph(netw1')

    p = gplot(swn,
        layout = layout,
        NODESIZE = 0.20 / sqrt(n),
        nodefillc = color,
        nodelabel = 1:n,
        nodelabeldist = 0,
        nodelabelangleoffset = π / 4,
        NODELABELSIZE = (n < 50) ? 5 : 2,
        edgelabel = netw1.nzval,
        EDGELABELSIZE = label ? 3 : 0,
        edgestrokec = Colors.parse(Colorant, "gray74"),
        edgelinewidth = edgelinewidth,
        EDGELINEWIDTH = 1.7 / sqrt(n),
        arrowlengthfrac = (directed == false) ? 0.0 : maximum([0.028, 0.3 / length(netw1.nzval)]),
        arrowangleoffset = (π / 9),
    )
    return p
end

function plot_network(netw1::SimpleWeightedDiGraph;
    layout = spring_layout,
    label = true,
    edgelinewidth = 1.0,
    color = Colors.parse(Colorant, "turquoise"))
    return plot_network(netw1.weights; layout = layout, label = label,
        edgelinewidth = edgelinewidth, color = color)
end

function plot_network(netw1::SimpleWeightedGraph;
    layout = spring_layout,
    label = true,
    edgelinewidth = 1.0,
    color = Colors.parse(Colorant, "turquoise"))
    return plot_network(netw1.weights; layout = layout, label = label,
        edgelinewidth = edgelinewidth, color = color)
end

function plot_graph(g::SimpleGraph; color = false)
    # transform in a Simple Weighted Graph
    sparse_mat = adjacency_matrix(g)

    if color
        mat = SimpleWeightedGraph(sparse_mat)
        colors = palette(:default, sparse_mat.m)
        comp = dfs_find_components(mat)
        plot_network(sparse_mat; label = false, color = colors[comp])
    else
        plot_network(sparse_mat; label = false)
    end

end

function plot_graph(g::SimpleDiGraph; color = false)
    sparse_mat = adjacency_matrix(g)

    if color
        mat = SimpleWeightedDiGraph(sparse_mat')
        colors = palette(:default, sparse_mat.m)
        comp = dfs_find_components(mat)
        plot_network(sparse_mat; label = false, color = colors[comp])
    else
        plot_network(sparse_mat; label = false)
    end

end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots
using DiffEqSensitivity
using Zygote
using ForwardDiff
using LinearAlgebra
using Random
using Statistics
using ProgressBars, Printf
using Flux.Optimise: update!
using Flux.Losses: mae, mse
using BSON: @save, @load

is_restart  = false;
n_epoch     = 50000;
ntotal      = 50
n_plot      = 200;

batch_size  = 50;
opt         = ADAMW(0.005, (0.9, 0.999), 1.f-6);
lb          = 1.e-6;
ode_solver  = AutoTsit5(Rosenbrock23(autodiff=false));

function rober!(du, u, k, t)
    y1, y2, y3  = u
    k1, k2, k3  = k
    du[1]       = -k1 * y1 + k3 * y2 * y3
    du[2]       =  k1 * y1 - k3 * y2 * y3 - k2 * y2^2
    du[3]       =  k2 * y2^2
end

u0          = [1.0, 0, 0];
tspan       = (0.0, 1e5);
t_end       = tspan[2];
k           = [0.04, 3e7, 1e4];
tsteps      = 10 .^ (range(-5, log10(tspan[2]), length=ntotal));
prob_rober  = ODEProblem(rober!, u0, tspan, k);
sol_rober   = solve(prob_rober, ode_solver, saveat=tsteps, atol=1.f-8);
normdata    = Array(sol_rober)

yscale      = maximum(normdata, dims = 2)  #scale for each species

i_slow      = [1, 2, 3];
nslow       = length(i_slow);
node        = 5;
dudt2       = Chain(x -> x,
                  Dense(nslow, node, gelu),
                  Dense(node, node, gelu),
                  Dense(node, node, gelu),
                  Dense(node, node, gelu),
                  Dense(node, node, gelu),
                  Dense(node, node, gelu),
                  Dense(node, nslow))

p, re       = Flux.destructure(dudt2);
rep         = re(p)

yscale_     = yscale[:, 1]
function dudt!(du, u, p, t)
    du      .= rep(u) .* yscale_ /t_end
end

prob        = ODEProblem(dudt!, u0[i_slow], tspan)
sense       = BacksolveAdjoint(checkpointing=true; autojacvec=ZygoteVJP());
function predict_n_ode(p, sample)
    global rep  = re(p)
    _prob       = remake(prob, p=p, tspan=[0, tsteps[sample]])
    pred        = Array(solve(_prob, ode_solver, saveat=tsteps[1:sample], atol=lb, sensalg=sense))
end
pred        = predict_n_ode(p, ntotal)

function loss_n_ode(p, sample=ntotal)
    pred = predict_n_ode(p, sample)
    loss = mae(pred ./ yscale, normdata[i_slow, 1:size(pred)[2]] ./ yscale)
    return loss
end
loss_n_ode(p, ntotal)

list_loss = []
list_grad = []
iter = 1
cb = function (p, loss_mean, g_norm)
    global list_loss, list_grad, iter
    push!(list_loss, loss_mean)
    push!(list_grad, g_norm)

    if iter % n_plot == 0
        pred = predict_n_ode(p, ntotal)

        list_plt = []
        for i in 1:nslow
            j = i_slow[i]
            plt = scatter(tsteps[2:end], normdata[j,2:end], xscale=:log10, label="data")
            plot!(plt, tsteps[2:end], pred[i,2:end], lw=2, xscale=:log10, label="pred", framestyle=:box)
            plot!(plt, xtickfontsize=11, ytickfontsize=11, xguidefontsize=12, yguidefontsize=12)
            ylabel!(plt, "y$j")
            xlabel!(plt, "Time")
            if i == 1
                plot!(plt, legend=:best)
            else
                plot!(plt, legend=false)
            end
            push!(list_plt, plt)
        end
        plt_all = plot(list_plt..., layouts = (nslow, 1))
        png(plt_all, "figs/pred.png")

        plt_loss = plot(list_loss, xscale=:identity, yscale=:log10, label="loss", legend=:topright)
        plt_grad = plot(list_grad, xscale=:identity, yscale=:log10, label="grad_norm", legend=:bottomright)
        xlabel!(plt_loss, "Epoch")
        xlabel!(plt_grad, "Epoch")
        ylabel!(plt_loss, "Loss")
        ylabel!(plt_grad, "Gradient Norm")
        plt_all = plot([plt_loss, plt_grad]..., framestyle=:box, size = (1000, 400))
        plot!(plt_all, xtickfontsize=10, ytickfontsize=10, xguidefontsize=12, yguidefontsize=12)
        png(plt_all, "figs/loss_grad")

        @save "./checkpoint/mymodel.bson" p opt list_loss list_grad iter
    end
    iter += 1
    return false
end

if is_restart
    @load "./checkpoint/mymodel.bson" p opt list_loss list_grad iter
    iter += 1
end

epochs = ProgressBar(iter:n_epoch);
for epoch in epochs
    global p
    sample = rand(batch_size:ntotal)

    loss = loss_n_ode(p, sample)
    grad = ForwardDiff.gradient(x -> loss_n_ode(x, sample), p)

    grad_norm = norm(grad, 2)
    update!(opt, p, grad)

    set_description(epochs, string(@sprintf("Loss: %.4e grad: %.2e", loss, grad_norm)))
    cb(p, loss, grad_norm)
end

max_stars_repo_name: jiweiqi/StiffNeuralODE

max_issues_repo_name: jiweiqi/StiffNeuralODE

max_forks_repo_name: jiweiqi/StiffNeuralODE


text:
module WordTokenizers

using HTML_Entities
using StrTables
using Unicode

export poormans_tokenize, punctuation_space_tokenize,
       penn_tokenize, improved_penn_tokenize, nltk_word_tokenize,
       tweet_tokenize,
       tokenize,
       rulebased_split_sentences,
       split_sentences,
       set_tokenizer, set_sentence_splitter,
       rev_tokenize, rev_detokenize,
       toktok_tokenize

include("words/fast.jl")

include("words/simple.jl")
include("words/nltk_word.jl")
include("words/reversible_tokenize.jl")
include("words/sedbased.jl")
include("words/tweet_tokenizer.jl")
include("sentences/sentence_splitting.jl")
include("words/TokTok.jl")

include("set_method_api.jl")
include("split_api.jl")

end # module


text:
convert_name(x::Union{LibClang.CXCursor, LibClang.CXType}) = convert_name(string(x))

function convert_name(str::String)
	# TODO: some of these are actually safe to use
	issafe = !(str in (
		"begin", "let", "quote", "do", "end",
		"macro", "function", "return",
		"while", "for", "continue", "break",
		"if", "elseif", "else",
		"try", "catch", "finally",
		"local", "global", "const",
		"struct", "mutable", "abstract", "primitive", "type",
		"module", "baremodule",
		"using", "import", "export",
		"true", "false",
	))
	return issafe ? str : "var\"$(str)\""
end

function convert_name(typ::LibClang.CXType; nolongs::Bool = false)
	if typ.kind == LibClang.CXType_Bool
		return "𝐣𝐥.Cbool"
	elseif typ.kind == LibClang.CXType_Void
		return "𝐣𝐥.Cvoid"
	elseif typ.kind == LibClang.CXType_Char_S
		return "𝐣𝐥.Cchar"
	elseif typ.kind == LibClang.CXType_Char_U
		return "𝐣𝐥.Cuchar"
	elseif typ.kind == LibClang.CXType_SChar
		return "𝐣𝐥.Cchar"
	elseif typ.kind == LibClang.CXType_UChar
		return "𝐣𝐥.Cuchar"
	elseif typ.kind == LibClang.CXType_Short
		return "𝐣𝐥.Cshort"
	elseif typ.kind == LibClang.CXType_Int
		return "𝐣𝐥.Cint"
	elseif typ.kind == LibClang.CXType_Long
		return nolongs ? (sizeof(Clong) == sizeof(Cint) ? "𝐣𝐥.Cint" : "𝐣𝐥.Clonglong") : "𝐣𝐥.Clong"
	elseif typ.kind == LibClang.CXType_LongLong
		return "𝐣𝐥.Clonglong"
	elseif typ.kind == LibClang.CXType_UShort
		return "𝐣𝐥.Cushort"
	elseif typ.kind == LibClang.CXType_UInt
		return "𝐣𝐥.Cuint"
	elseif typ.kind == LibClang.CXType_ULong
		return nolongs ? (sizeof(Culong) == sizeof(Cuint) ? "𝐣𝐥.Cuint" : "𝐣𝐥.Culonglong") : "𝐣𝐥.Culong"
	elseif typ.kind == LibClang.CXType_ULongLong
		return "𝐣𝐥.Culonglong"
	elseif typ.kind == LibClang.CXType_Float
		return "𝐣𝐥.Cfloat"
	elseif typ.kind == LibClang.CXType_Double
		return "𝐣𝐥.Cdouble"
	elseif typ.kind == LibClang.CXType_LongDouble
		return "𝐣𝐥.Clongdouble"
	elseif typ.kind in (
		LibClang.CXType_Elaborated,
		LibClang.CXType_Enum,
	)
		decl = LibClang.clang_getTypeDeclaration(typ)
		if decl.kind == LibClang.CXCursor_EnumDecl
			kind = "𝐣𝐥.@cenum"
		elseif decl.kind == LibClang.CXCursor_StructDecl
			kind = "𝐣𝐥.@cstruct"
		elseif decl.kind == LibClang.CXCursor_UnionDecl
			kind = "𝐣𝐥.@cunion"
		else
			error("Unable to handle kind of type $(decl.kind)")
		end
		return "$(kind) $(convert_name(decl))"
	elseif typ.kind in (
		LibClang.CXType_Typedef,
		LibClang.CXType_Enum,
		LibClang.CXType_Record,
	)
		return convert_name(LibClang.clang_getTypeDeclaration(typ))
	end

	error("Unable to convert type name $(typ.kind), not yet implemented")
end


text:
using DEBTrait
using CSV, DataFrames, Statistics

n_polymers              = 0
n_monomers              = 6
n_microbes              = 2
n_enzymes               = 1
n_minerals              = 0

df_isolates             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/microbes2traits.csv", DataFrame, missingstring = "N/A")
df_monomers             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_class.csv", DataFrame, missingstring = "N/A")
df_p                    = filter(row -> row[:Rhizosphere_response] == "positive", df_isolates)
df_n                    = filter(row -> row[:Rhizosphere_response] == "negative", df_isolates)

L_DNA                   = [median(df_p.Genome_size), median(df_n.Genome_size)]
rRNA                    = [median(df_p.rRNA_genes), median(df_n.rRNA_genes)]
min_gentime             = [median(df_p.Min_gen_time), median(df_n.Min_gen_time)]

z_sugars_p              = median(convert(Array{Float64,1}, df_p.z_sugars))*ones(1,1)
z_organics_p            = median(convert(Array{Float64,1}, df_p.z_organic_acids))*ones(1,1)
z_aminos_p              = median(convert(Array{Float64,1}, df_p.z_amino_acids))*ones(1,1)
z_fattys_p              = median(convert(Array{Float64,1}, df_p.z_fatty_acids))*ones(1,1)
z_nucleos_p             = median(convert(Array{Float64,1}, df_p.z_nucleotides))*ones(1,1)
z_auxins_p              = median(convert(Array{Float64,1}, df_p.z_auxins))*ones(1,1)
genome_distr_p          = vcat(z_sugars_p, z_organics_p, z_aminos_p, z_fattys_p, z_nucleos_p, z_auxins_p)

z_sugars_n              = median(convert(Array{Float64,1}, df_n.z_sugars))*ones(1,1)
z_organics_n            = median(convert(Array{Float64,1}, df_n.z_organic_acids))*ones(1,1)
z_aminos_n              = median(convert(Array{Float64,1}, df_n.z_amino_acids))*ones(1,1)
z_fattys_n              = median(convert(Array{Float64,1}, df_n.z_fatty_acids))*ones(1,1)
z_nucleos_n             = median(convert(Array{Float64,1}, df_n.z_nucleotides))*ones(1,1)
z_auxins_n              = median(convert(Array{Float64,1}, df_n.z_auxins))*ones(1,1)
genome_distr_n          = vcat(z_sugars_n, z_organics_n, z_aminos_n, z_fattys_n, z_nucleos_n, z_auxins_n)

genome_distr            = hcat(genome_distr_p, genome_distr_n)

Isolate                 = BaseGenome(L_DNA, rRNA, min_gentime, genome_distr)
Monomer                 = BaseMonomer(df_monomers.D_S, df_monomers.N_C, df_monomers.N_N, df_monomers.y_DE)
Soil                    = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))

MetabolismC             = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
MetabolismC.y_EM        = ones(n_microbes)
MetabolismC.y_EX        = ones(n_microbes)
MetabolismC.α           = zeros(n_microbes)

AssimilationC           = PAssimilation(Isolate, Monomer, Soil, MetabolismC)

Setup                   = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
p                       = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)


text:
#config for build bert
const config = Dict(
  "attention_probs_dropout_prob" => 0.1,
  "initializer_range"            => 0.02,
  "hidden_act"                   => gelu,
  "num_hidden_layers"            => 12,
  "hidden_size"                  => 768,
  "max_position_embeddings"      => 512,
  "type_vocab_size"              => 2,
  "vocab_size"                   => 30522,
  "num_attention_heads"          => 12,
  "hidden_dropout_prob"          => 0.1,
  "intermediate_size"            => 3072,
)

#create bert model like pretrain struct
#you can define the bert model in the way you like and wrap it with TransformerModel
function create_bert()
  global config
  bert = Bert(
    config["hidden_size"],
    config["num_attention_heads"],
    config["intermediate_size"],
    config["num_hidden_layers"];
    act = config["hidden_act"],
    pdrop = config["hidden_dropout_prob"],
    attn_pdrop = config["attention_probs_dropout_prob"]
  )

  tok_emb = Embed(
    config["hidden_size"],
    config["vocab_size"]
  )

  seg_emb = Embed(
    config["hidden_size"],
    config["type_vocab_size"]
  )

  posi_emb = PositionEmbedding(
    config["hidden_size"],
    config["max_position_embeddings"];
    trainable = true
  )

  emb_post = Positionwise(
    LayerNorm(
      config["hidden_size"]
    ),
    Dropout(
      config["hidden_dropout_prob"]
    )
  )

  pooler = Dense(
    config["hidden_size"],
    config["hidden_size"],
    tanh
  )

  masklm = (
    transform = Chain(
      Dense(
        config["hidden_size"],
        config["hidden_size"],
        config["hidden_act"]
      ),
      LayerNorm(
        config["hidden_size"]
      )
    ),
    output_bias = param(randn(
      Float32,
      config["vocab_size"]
    ))
  )

  nextsentence = Chain(
    Dense(
      config["hidden_size"],
      2
    ),
    logsoftmax
  )

  emb = CompositeEmbedding(tok = tok_emb, pe = posi_emb, segment = seg_emb, postprocessor = emb_post)
  clf = (pooler = pooler, masklm = masklm, nextsentence = nextsentence)

  TransformerModel(emb, bert, clf)
end


text:
include(Pkg.dir("GenomeAnnotations","src","Config.jl"))
include(Pkg.dir("GenomeAnnotations","src","utils.jl"))

# db_id is normally associated with an organism, e.g. M2, H19
# it's tempting to call this organism_id or org_id or org_db or something
# but leave it general
function install_gencode_annotation( db_id )
    println("Installing annotations for $db_id")
    download_gencode(db_id)
end

function list_gencode_annotations()
   return ["GCM2","GCH19"]
end

# Can be M2 or H19
function download_gencode( db_id )
    LOCAL_STORE = LOCAL_STORE_PATH()
    local_store_dir = joinpath( LOCAL_STORE, "")
    path = ""
    if(db_id == "GCM2")
        ftp = "ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M2/gencode.vM2.annotation.gtf.gz"
        local_store_dir = joinpath( LOCAL_STORE, db_id)
        make_localstore(local_store_dir)
        basename_gencode = basename(ftp)
        path = joinpath(local_store_dir, basename_gencode)
        download_with_curl(ftp, path)
    elseif( db == "GCH19")
        error("Not implemented yet for Gencode $db")
    else
        error("Gencode db is not one of GCH19 or GCM2 (you gave $db)")
    end
end

# create some feature specific files for convenience, especially a gene coordinates file
# A non coding file?


text:
export generate_test_problems, generate_noisy_data!,
    generate_noisy_data, generate_clustered_noisy_data,
    generate_clustered_noisy_data!

"""

This dictionary represents the list of models used in the generation of random tests.
Return the tuple `(n, model, model_str)`, where

  - `n` is the number of parameters of the model
  - `model` is the model of the form `m(x, θ)`, where `x` are the
    variables and `θ` are the parameters
  - `model_str` is the string representing the model, used to build random generated problems

"""
const model_list = Dict(
    "linear" => (2, (x, θ) -> θ[1] * x[1] + θ[2],
                 "(x, θ) -> θ[1] * x[1] + θ[2]"),
    "cubic" => (4, (x, θ) -> θ[1] * x[1]^3 + θ[2] * x[1]^2 + θ[3] * x[1] + θ[4],
                "(x, θ) -> θ[1] * x[1]^3 + θ[2] * x[1]^2 + θ[3] * x[1] + θ[4]"),
    "expon" => (3, (x, θ) -> θ[1] + θ[2] * exp(- θ[3] * x[1]),
                "(x, θ) -> θ[1] + θ[2] * exp(- θ[3] * x[1])"),
    "logistic" => (4, (x, θ) -> θ[1] + θ[2] / (1.0 + exp(- θ[3] * x[1] + θ[4])),
                   "(x, θ) -> θ[1] + θ[2] / (1.0 + exp(- θ[3] * x[1] + θ[4]))"),
    "circle" => (3, (x, θ) -> (x[1] - θ[1])^2 + (x[2] - θ[2])^2 - θ[3]^2,
                 "(x, θ) -> (x[1] - θ[1])^2 + (x[2] - θ[2])^2 - θ[3]^2"),
    "ellipse" => (6, (x, θ) -> θ[1] * x[1]^2 + θ[2] * x[1] * x[2] + θ[3] * x[2]^2 +
                  θ[4] * x[1] + θ[5] * x[2] + θ[6],
                  "(x, θ) -> θ[1] * x[1]^2 + θ[2] * x[1] * x[2] + θ[3] * x[2]^2 + θ[4] * x[1] + θ[5] * x[2] + θ[6]")
)

"""

    interval_rand!(x::Vector{Float64},
        intervals::Vector{Tuple{Float64, Float64}})

Fill a vector `x` with uniformly distributed random numbers generated
in the interval given by `intervals`. It is assumed that `length(x) ==
length(intervals)`.

Throws an `ErrorException` if the dimension of `x` is smaller the
dimension of `intervals` or if the intervals are invalid.

"""
function interval_rand!(x::Vector{Float64},
    intervals::Vector{Tuple{Float64, Float64}})

    (length(x) < length(intervals)) &&
        error("Length of vector smaller than length of intervals.")

    for v in intervals

        (v[1] > v[2]) && error("Bad interval $(v[1])>$(v[2]).")

    end

    map!((v) -> v[1] + rand() * (v[2] - v[1]), x, intervals)

end

"""

    generate_test_problems(datFilename::String, solFilename::String,
        model::Function, modelStr::String, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(n), std::Float64=200.0,
        out_times::Float64=7.0)

    generate_test_problems(datFilename::String, solFilename::String,
        model::Function, modelStr::String, n::Int, np::Int, p::Int,
        cluster_interval::Tuple{Float64, Float64};
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(n), std::Float64=200.0,
        out_times::Float64=7.0)

Generate random data files for testing fitting problems.

  - `datFilename` and `solFilename` are strings with the name of the
    files for storing the random data and solution, respectively.
  - `model` is the model function and `modelStr` is a string
    representing this model function, e.g.

         model = (x, θ) -> θ[1] * x[1] + θ[2]
         modelStr = "(x, θ) -> θ[1] * x[1] + θ[2]"

    where vector `θ` represents the parameters (to be found) of the
    model and vector `x` are the variables of the model.
  - `n` is the number of parameters
  - `np` is the number of points to be generated.
  - `p` is the number of trusted points to be used in the LOVO
    approach.

If `cluster_interval` is provided, then generates outliers only in
this interval.

Additional parameters:

  - `xMin`, `xMax`: interval for generating points in one dimensional
    tests *Deprecated*
  - `x_interval`: interval for generating points in one dimensional
    tests
  - `θSol`: true solution, used for generating perturbed points
  - `std`: standard deviation
  - `out_times`: deviation for outliers will be `out_times * std`.

"""
function generate_test_problems(datFilename::String,
    solFilename::String, model::Function, modelStr::String, n::Int,
    np::Int, p::Int; gn_kwargs...)

    # Generate data file

    θSol = nothing

    open(datFilename, "w") do data

        vdata, θSol, outliers = generate_noisy_data(model, n, np, p;
                                                    gn_kwargs...)

        # Dimension of the domain of the function to fit
        @printf(data, "%d\n", 1)

        for k = 1:np

            @printf(data, "%20.15f %20.15f %1d\n",
                    vdata[k, 1], vdata[k, 2], Int(k in outliers))

        end

    end

    # Generate solution file

    open(solFilename, "w") do sol

        println(sol, n) # number of variables

        println(sol, θSol) # parameters

        println(sol, modelStr) # function expression

    end

end

function generate_test_problems(datFilename::String,
    solFilename::String, model::Function, modelStr::String, n::Int,
    np::Int, p::Int, x_interval::Tuple{Float64, Float64},
    cluster_interval::Tuple{Float64, Float64}; gn_kwargs...)

    # Generate data file

    θSol = nothing

    open(datFilename, "w") do data

        vdata, θSol, outliers = generate_clustered_noisy_data(model,
            n, np, p, x_interval, cluster_interval; gn_kwargs...)

        # Dimension of the domain of the function to fit
        @printf(data, "%d\n", 1)

        for k = 1:np

            @printf(data, "%20.15f %20.15f %1d\n",
                    vdata[k, 1], vdata[k, 2], Int(k in outliers))

        end

    end

    # Generate solution file

    open(solFilename, "w") do sol

        println(sol, n) # number of variables

        println(sol, θSol) # parameters

        println(sol, modelStr) # function expression

    end

end

"""

    get_unique_random_points(np::Int, npp::Int)

Choose exactly `npp` unique random points from a set containing `np`
points. This function is similar to `rand(vector)`, but does not allow
repetitions.

If `npp` < `np`, returns all the `np` points. Note that this function
is not very memory efficient, since the process of selecting unique
elements involves creating several temporary vectors.

Return a vector with the selected points.

"""
function get_unique_random_points(np::Int, npp::Int)

    # Check invalid arguments
    ((np <= 0) || (npp <=0)) && (return Vector{Int}())

    ntp = min(npp, np)

    v = Vector{Int}(undef, ntp)

    return get_unique_random_points!(v, np, npp)

end

"""

    get_unique_random_points!(v::Vector{Int}, np::Int, npp::Int)

Choose exactly `npp` unique random points from a set containing `np`
points. This function is similar to `rand(vector)`, but does not allow
repetitions.

If `npp` < `np`, returns all the `np` points. Note that this function
is not very memory efficient, since the process of selecting unique
elements involves creating several temporary vectors.

Return the vector `v` provided as argument filled with the selected
points.

"""
function get_unique_random_points!(v::Vector{Int}, np::Int, npp::Int)

    # Check invalid arguments
    ((np <= 0) || (npp <=0)) && return v

    ntp = min(np, npp)

    # Check invalid size
    (length(v) < ntp) && throw(DimensionMismatch("Incorrect size for vector."))

    # Check small np for efficiency
    if np == ntp

        for i = 1:np
            v[i] = i
        end

        return v

    end

    points = [1:np;]

    while ntp > 0

        u = rand(points, ntp)

        unique!(u)

        for i in u
            v[ntp] = i
            ntp -= 1
        end

        setdiff!(points, u)

    end

    return v

end

"""

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(Float64, n),
        std::Float64=200.0, out_times::Float64=7.0)

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
        x_interval::Tuple{Float64, Float64})

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
        θSol::Vector{Float64}, x_interval::Tuple{Float64, Float64})

Random generate a fitting one-dimensional data problem.

This function receives a `model(x, θ)` function, the number of parameters
`n`, the number of points `np` to be generated and the number of
trusted points `p`.

If the `n`-dimensional vector `θSol` is provided, then the exact
solution will not be random generated. The interval `[xMin, xMax]`
(*deprecated*) or `x_interval` for generating the values to evaluate
`model` can also be provided.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`.
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set

"""
function generate_noisy_data(model::Function, n::Int, np::Int,
    p::Int; gn_kwargs...)

    # Data matrix
    data = Array{Float64}(undef, np, 3)

    # Points selected to be outliers
    v = Vector{Int}(undef, np - p)

    return generate_noisy_data!(data, v, model, n, np, p;
                                gn_kwargs...)

end

# Deprecated

@deprecate(generate_noisy_data(model::Function, n, np, p,
    xMin::Float64, xMax::Float64), generate_noisy_data(model, n, np,
    p, (xMin, xMax)))

@deprecate(generate_noisy_data(model::Function, n::Int, np::Int,
    p::Int, θSol::Vector{Float64}, xMin::Float64, xMax::Float64),
    generate_noisy_data(model, n, np, p, θSol, (xMin, xMax)))

generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
    x_interval::Tuple{Float64, Float64}) = generate_noisy_data(model,
    n, np, p; x_interval=x_interval)

generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
    θSol::Vector{Float64}, x_interval::Tuple{Float64, Float64}) =
    generate_noisy_data(model, n, np, p; x_interval=x_interval,
    θSol=θSol)

"""

    generate_noisy_data!(data::AbstractArray{Float64, 2},
        v::Vector{Int}, model::Function, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(Float64, n),
        std::Float64=200.0, out_times::Float64=7.0)

Random generate a fitting one-dimensional data problem, storing the
data in matrix `data` and the outliers in vector `v`.

This function receives a `model(x, θ)` function, the number of parameters
`n`, the number of points `np` to be generated and the number of
trusted points `p`.

If the `n`-dimensional vector `θSol` is provided, then the exact
solution will not be random generated. The interval `[xMin, xMax]`
(*deprecated*) or `x_interval` for generating the values to evaluate
`model` can also be provided.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`.
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set

"""
function generate_noisy_data!(data::AbstractArray{Float64, 2},
    v::Vector{Int}, model::Function, n::Int, np::Int, p::Int;
    x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
    θSol::Vector{Float64}=10.0 * randn(Float64, n),
    std::Float64=200.0, out_times::Float64=7.0)

    @assert(x_interval[1] <= x_interval[2],
            "Invalid interval for random number generation.")
    @assert(size(data) == (np, 3),
            "Invalid size of data matrix. $(size(data)) != $((np, 3)).")
    @assert(length(v) >= np - p,
            "Invalid size for vector of outliers.")

    # Generate (x_i) where x_interval[1] <= x_i <= x_interval[2] (data)
    # Fix the problem of large interval with 1 element.
    x = (np == 1) ? sum(x_interval) / 2.0 : LinRange(x_interval[1], x_interval[2], np)

    # Points selected to be outliers
    get_unique_random_points!(v, np, np - p)

    # Add noise to some random points
    sgn = sign(randn())

    for k = 1:np

        y = model(x[k], θSol) + randn() * std

        noise = 0.0

        if k in v
            y = model(x[k], θSol)
            noise = (1.0 + 2 * rand()) * out_times * std * sgn
        end

        data[k, 1] = x[k]
        data[k, 2] = y + noise
        data[k, 3] = noise

    end

    return data, θSol, v

end

"""

    generate_clustered_noisy_data(model::Function, n::Int, np::Int,
        p::Int, x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

    generate_clustered_noisy_data(model::Function, n::Int,
        np::Int, p::Int, θSol::Vector{Float64},
        x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

Generate a test set with clustered outliers.

The arguments and optional arguments are the same for
[`generate_noisy_data!`](@ref), with exception of tuple
`cluster_interval` which is the interval to generate the clustered
outliers.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`. The same array given as argument
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set. The same vector given
    as argument.

"""
function generate_clustered_noisy_data(model::Function, n::Int,
    np::Int, p::Int, x_interval::Tuple{Float64,Float64},
    cluster_interval::Tuple{Float64, Float64}; kwargs...)

    data = Array{Float64, 2}(undef, np, 3)

    v = Vector{Int}(undef, np - p)

    return generate_clustered_noisy_data!(data, v, model, n, np, p,
               x_interval, cluster_interval; kwargs...)

end

generate_clustered_noisy_data(model::Function, n::Int, np::Int,
    p::Int, θSol::Vector{Float64}, x_interval::Tuple{Float64,Float64},
    cluster_interval::Tuple{Float64, Float64}; kwargs...) =
    generate_clustered_noisy_data(model, n, np, p, x_interval,
    cluster_interval, θSol=θSol; kwargs...)

"""

    generate_clustered_noisy_data!(data::Array{Float64, 2},
        v::Vector{Int}, model::Function, n::Int, np::Int, p::Int,
        x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

Generate a test set with clustered outliers. This version overwrites
the content of (`np` x `3`) matrix `data` and vector `v` with integer
indices to the position of outliers in `data`.

The arguments and optional arguments are the same for
[`generate_noisy_data!`](@ref), with exception of tuple
`cluster_interval` which is the interval to generate the clustered
outliers.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`. The same array given as argument
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set. The same vector given
    as argument.

"""
function generate_clustered_noisy_data!(data::Array{Float64, 2},
            v::Vector{Int}, model::Function, n::Int, np::Int,
            p::Int, x_interval::Tuple{Float64,Float64},
            cluster_interval::Tuple{Float64, Float64}; kwargs...)

    if (np - p > 0) &&
        !(x_interval[1] <= cluster_interval[1] <
          cluster_interval[2] <= x_interval[2])

        error("Bad interval for clustered data generation.")

    end

    interval_len = x_interval[2] - x_interval[1]

    fr1 = (cluster_interval[1] - x_interval[1]) / interval_len
    fr2 = (cluster_interval[2] - cluster_interval[1]) / interval_len
    fr3 = (x_interval[2] - cluster_interval[2]) / interval_len

    # Distribute data in a proportional way

    # Interval 2 will contain the clustered outliers
    np2 = max(np - p, Int(round(fr2 * np)))
    np1 = min(np - np2, Int(round(fr1 * np)))
    # Interval 3 will contain the remaining points
    np3 = max(0, np - np1 - np2)

    @debug("Clustered points: $(np1), $(np2), $(np3).")

    # Avoid repetition in the borders of the intervals

    δ_c = (cluster_interval[2] - cluster_interval[1]) / (np2 + 2)

    # Generate data

    tmpv = Vector{Int}()

    tmpd, θSol, tmpv = generate_noisy_data!(@view(data[1:np1, :]),
        tmpv, model, n, np1, np1; x_interval=(x_interval[1],
        cluster_interval[1]), kwargs...)

    generate_noisy_data!(@view(data[np1 + 1:np1 + np2, :]), v, model,
        n, np2, np2 - (np - p); x_interval=(cluster_interval[1] + δ_c,
        cluster_interval[2] - δ_c), θSol=θSol, kwargs...)

    # Update the outlier number with the correct number
    map!((x) -> x + np1, v, v)

    generate_noisy_data!(@view(data[np1 + np2 + 1:np, :]), tmpv,
        model, n, np3, np3; x_interval=(cluster_interval[2],
        x_interval[2]), θSol=θSol, kwargs...)

    return data, θSol, v

end


text:
# slab_rhs!
using ClimaCore

struct ThermalSlabParameters# <: CLIMAParameters.AbstractEarthParameterSet{F}
    h::FT
    ρ::FT
    c::FT
    T_init::FT
end

# domain
function ShellDomain(; radius = 6371e3, Nel = 8, Nq = 2)
    domain = ClimaCore.Domains.SphereDomain(radius)
    mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, Nel)
    topology = ClimaCore.Topologies.Topology2D(mesh)
    quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(topology, quad)
end

# init simulation
function slab_space_init(::Type{FT}, space, params) where {FT}

    coords = ClimaCore.Fields.coordinate_field(space)

    # initial condition
    T_sfc = map(coords) do coord
        T_sfc_0 = params.T_init
        anom_ampl = FT(2)
        radlat = coord.lat / FT(180) * pi
        lat_0 = FT(60) / FT(180) * pi
        lon_0 = FT(-90) / FT(180) * pi
        radlon = coord.long / FT(180) * pi
        stdev = FT(5) / FT(180) * pi
        anom = anom_ampl * exp(-((radlat - lat_0)^2 / 2stdev^2 + (radlon - lon_0)^2 / 2stdev^2))
        T_sfc = T_sfc_0 + anom
    end

    # prognostic variable
    Y = ClimaCore.Fields.FieldVector(T_sfc = T_sfc)

    return Y, space
end

# ode
function slab_rhs!(dY, Y, Ya, t)
    """
    Slab ocean:
    ∂_t T_sfc = F_sfc + G
    """
    p, F_sfc = Ya

    @. dY.T_sfc = (F_sfc) / (p.h * p.ρ * p.c)
end

struct SlabSimulation{P, Y, D, I}
    params::P
    Y_init::Y
    domain::D
    integrator::I
end

function slab_init(
    ::Type{FT},
    tspan;
    stepper = Euler(),
    nelements = 6,
    npolynomial = 4,
    dt = 0.02,
    saveat = 1.0e10,
    space = nothing,
) where {FT}

    params = ThermalSlabParameters(FT(0.5), FT(1500.0), FT(800.0), FT(281.0))

    Y, space = slab_space_init(FT, space, params)
    Ya = (params = params, F_sfc = ClimaCore.Fields.zeros(space)) #auxiliary
    problem = OrdinaryDiffEq.ODEProblem(slab_rhs!, Y, tspan, Ya)
    integrator = OrdinaryDiffEq.init(problem, stepper, dt = dt, saveat = saveat)

    SlabSimulation(params, Y, space, integrator)
end

get_slab_energy(slab_sim, T_sfc) = slab_sim.params.ρ .* slab_sim.params.c .* T_sfc .* slab_sim.params.h

max_stars_repo_name: CliMA/ClimaCoupler


text:
"""
    invariantManifolds(SYS, rv0, T; tf=1., nPts=20, α=1e-6, reltol=1e-12, integrator=TsitPap8())

Compute invariant manifolds of a periodic orbit in the CR3BP.
"""
function invariant_manifolds(sys::System, rv₀, T; tf=1., nPts=20, α=1e-6, reltol=1e-12, integrator=TsitPap8())

    Φ₀ = I(6)
    w₀ = vcat(rv₀,reshape(Φ₀,36,1))
    tspan = (0.,T)
    prob = ODEProblem(CR3BPstm!,w₀,tspan,sys)
    nomTraj = solve(prob, integrator, reltol=reltol)
    Φₜ = Matrix(reshape(nomTraj.u[end][7:42],6,6))
    rvₜ = nomTraj.u[end][1:6] # last time step
    Λ,V = eigen(Φₜ,sortby=isreal) # Λ is vector of eigenvalues, V is matrix of eigenvectors
    Yw = real(V[:,findall(isreal, Λ)]) # Eigenvectors corresponding to real eigenvalues
    Λ = Λ[findall(isreal, Λ)] # Purely real eigenvalues (have 0.0 imaginary component)
    Yws = Yw[:,findmin(real(Λ))[2]] # Eigenvector associated with stable eigenvalue λ < 1
    Ywu = Yw[:,findmax(real(Λ))[2]] # Eigenvector associated with unstable eigenvalue λ > 1

    tspan_forward = (0., tf)
    tspan_backward = (0., -tf)

    # Wsp = []
    # Wsn = []
    # Wup = []
    # Wun = []
    # for i = 1:length(nomTraj)
    #   Φᵢ = Matrix(reshape(nomTraj.u[i][7:42],6,6)) # do I need Matrix()?
    #   rv0sp = nomTraj[i][1:6] + α*Φᵢ*Yws/norm(Φᵢ*Yws);
    #   rv0sn = nomTraj[i][1:6] - α*Φᵢ*Yws/norm(Φᵢ*Yws);
    #   rv0up = nomTraj[i][1:6] + α*Φᵢ*Ywu/norm(Φᵢ*Ywu);
    #   rv0un = nomTraj[i][1:6] - α*Φᵢ*Ywu/norm(Φᵢ*Ywu);
    #
    #    prob_sp = ODEProblem(CR3BPdynamics!,rv0sp, tspan_backward, sys) # stable positive
    #    push!(Wsp, solve(prob_sp, integrator, reltol=reltol))
    #    prob_sn = ODEProblem(CR3BPdynamics!,rv0sn, tspan_backward, sys) # stable negative
    #    push!(Wsn, solve(prob_sn, integrator, reltol=reltol))
    #    prob_up = ODEProblem(CR3BPdynamics!,rv0up, tspan_forward, sys) # unstable positive
    #    push!(Wup, solve(prob_up, integrator, reltol=reltol))
    #    prob_un = ODEProblem(CR3BPdynamics!,rv0un, tspan_forward, sys) # unstable negative
    #    push!(Wun, solve(prob_un, integrator, reltol=reltol))
    # end

    Wsp = []
    Wsn = []
    Wup = []
    Wun = []
    for t = LinRange(0,T,nPts)
      Φᵢ = Matrix(reshape(nomTraj(t)[7:42],6,6)) # do I need Matrix()?
      rv0sp = nomTraj(t)[1:6] + α*Φᵢ*Yws/norm(Φᵢ*Yws);
      rv0sn = nomTraj(t)[1:6] - α*Φᵢ*Yws/norm(Φᵢ*Yws);
      rv0up = nomTraj(t)[1:6] + α*Φᵢ*Ywu/norm(Φᵢ*Ywu);
      rv0un = nomTraj(t)[1:6] - α*Φᵢ*Ywu/norm(Φᵢ*Ywu);

      prob_sp = ODEProblem(CR3BPdynamics!,rv0sp, tspan_backward, sys) # stable positive
      push!(Wsp, solve(prob_sp, integrator, reltol=reltol))
      prob_sn = ODEProblem(CR3BPdynamics!,rv0sn, tspan_backward, sys) # stable negative
      push!(Wsn, solve(prob_sn, integrator, reltol=reltol))
      prob_up = ODEProblem(CR3BPdynamics!,rv0up, tspan_forward, sys) # unstable positive
      push!(Wup, solve(prob_up, integrator, reltol=reltol))
      prob_un = ODEProblem(CR3BPdynamics!,rv0un, tspan_forward, sys) # unstable negative
      push!(Wun, solve(prob_un, integrator, reltol=reltol))
   end

    return Wsp, Wsn, Wup, Wun, Φₜ, nomTraj
end

function differential_corrector(sys::System, rv₀; myconst=3, iter=100, plot=false, t0=0., tf=1., dir=1, tol=1e-12)
   Φ₀ = I(6)

   tspan = (t0, tf)

   # event function
   condition(u, t, integrator) = u[2]
   affect!(integrator) = terminate!(integrator)
   cb = OrdinaryDiffEq.ContinuousCallback(condition, affect!)

   for k = 1:iter
      w₀ = vcat(rv₀, reshape(Φ₀,36,1))

      prob = ODEProblem(CR3BPstm!,w₀,tspan,sys)
      sol = solve(prob, TsitPap8(), reltol=tol, callback=cb)

      w = sol[end]
      rv = w[1:6]
      Φ = reshape(w[7:42],6,6)
      global T = 2*sol.t[end]

      ∂ẋ = 0 - rv[4]
      ∂ż = 0 - rv[6]

      if norm([∂ẋ, ∂ż]) < tol
         break
      end
      if k == iter
         @warn("Differential corrector did not converge")
      end

      ẋ,ẏ,ż,ẍ,ÿ,z̈ = CR3BPdynamics(rv, sys, 0)

      if myconst == 1
         ∂Φ = [Φ[4,3] Φ[4,5];
               Φ[6,3] Φ[6,5]]
         dyad = [ẍ;z̈]*[Φ[2,3] Φ[2,5]]

         ∂z, ∂ẏ = (∂Φ - dyad/ẏ)\[∂ẋ; ∂ż]
         rv₀[3] += ∂z
         rv₀[5] += ∂ẏ
      elseif myconst == 3
         ∂Φ = [Φ[4,1] Φ[4,5];
               Φ[6,1] Φ[6,5]]
         dyad = [ẍ;z̈]*[Φ[2,1] Φ[2,5]]

         ∂x, ∂ẏ = (∂Φ - dyad/ẏ)\[∂ẋ; ∂ż]
         rv₀[1] += ∂x
         rv₀[5] += ∂ẏ
      elseif myconst == 5
         ∂Φ = [Φ[4,1] Φ[4,3];
               Φ[6,1] Φ[6,3]]
         dyad = [ẍ;z̈]*[Φ[2,1] Φ[2,5]]

         ∂x, ∂z = (∂Φ - dyad/ẏ)\[∂ẋ; ∂ż]
         rv₀[1] += ∂x
         rv₀[3] += ∂ẏ
      else
         error("myconst should be 1, 3, or 5")
      end

   end

   return rv₀, T
end

""" condition(u, t, integrator) = u[2]
 affect!(integrator) = terminate!(integrator)
 cb = ContinousCallback(condition, affect!)
 sol = solve(prob, Tsit5(), calback=cb)

    rich3()
"""
function rich3(sys::System, Az, Lpt, NS, npts=10)
    μ = sys.μ
    γ = gammaL(sys, Lpt)

    if Lpt == 1;        won =  1;   primary = 1-μ;
    elseif Lpt == 2;    won = -1;   primary = 1-μ;
    elseif Lpt == 3;    won =  1;   primary = -μ;
    end

    c = zeros(4)

    if Lpt == 3
      for N = 2:4
        c[N]= (1/γ^3)*( 1-μ + (-primary*γ^(N+1))/((1+γ)^(N+1)) );
      end
    else
      for N = 2:4
        c[N]= (1/γ^3)*( (won^N)*μ + ((-1)^N)*((primary)*γ^(N+1))/((1+(-won)*γ)^(N+1)) );
      end
    end

    # polylambda = [1, 0, (c[2]-2), 0, -(c[2]-1)*(1+2*c[2])];
    polylambda = [-(c[2]-1)*(1+2*c[2]), 0, c[2]-2, 0, 1];
    lambdaroots = roots(polylambda); # lambda = frequency of orbit
 ## Stuck here 6/22/21 (need to select the right root)
    λ = real(sort(lambdaroots, by = x -> abs(imag(x))))[1]

    # if Lpt==3
    #    λ = lambdaroots[1];
    # else
    #    λ = lambdaroots[1];
    # end

    k   = 2*λ/(λ^2 + 1 - c[2]);

    del = λ^2 - c[2];

    d1  = ((3*λ^2)/k)*(k*( 6*λ^2 -1) - 2*λ);
    d2  = ((8*λ^2)/k)*(k*(11*λ^2 -1) - 2*λ);

    a21 = (3*c[3]*(k^2 - 2))/(4*(1 + 2*c[2]));
    a22 = 3*c[3]/(4*(1 + 2*c[2]));
    a23 = -(3*c[3]*λ/(4*k*d1))*( 3*(k^3)*λ - 6*k*(k-λ) + 4);
    a24 = -(3*c[3]*λ/(4*k*d1))*( 2 + 3*k*λ );

    b21 = -(3*c[3]*λ/(2*d1))*(3*k*λ - 4);
    b22 = 3*c[3]*λ/d1;
    d21 = -c[3]/(2*λ^2);

    a31 = -(9*λ/(4*d2))*(4*c[3]*(k*a23 - b21) + k*c[4]*(4 + k^2)) + ((9*λ^2 + 1 -c[2])/(2*d2))*(3*c[3]*(2*a23 - k*b21) + c[4]*(2 + 3*k^2));
    a32 = -(1/d2)*( (9*λ/4)*(4*c[3]*(k*a24 - b22) + k*c[4]) + 1.5*(9*λ^2 + 1 - c[2])*( c[3]*(k*b22 + d21 - 2*a24) - c[4]) );

    b31 = (.375/d2)*( 8*λ*(3*c[3]*(k*b21 - 2*a23) - c[4]*(2 + 3*k^2)) + (9*λ^2 + 1 + 2*c[2])*(4*c[3]*(k*a23 - b21) + k*c[4]*(4 + k^2)) );
    b32 = (1/d2)*( 9*λ*(c[3]*(k*b22 + d21 - 2*a24) - c[4]) + .375*(9*λ^2 + 1 + 2*c[2])*(4*c[3]*(k*a24 - b22) + k*c[4]) );

    d31 = (3/(64*λ^2))*(4*c[3]*a24 + c[4]);
    d32 = (3/(64*λ^2))*(4*c[3]*(a23- d21) + c[4]*(4 + k^2));

    s1  = (1/(2*λ*(λ*(1+k^2) - 2*k)))*( 1.5*c[3]*(2*a21*(k^2 - 2)-a23*(k^2 + 2) - 2*k*b21) - .375*c[4]*(3*k^4 - 8*k^2 + 8) );
    s2  = (1/(2*λ*(λ*(1+k^2) - 2*k)))*( 1.5*c[3]*(2*a22*(k^2 - 2)+a24*(k^2 + 2) + 2*k*b22 + 5*d21) + .375*c[4]*(12 - k^2) );

    a1  = -1.5*c[3]*(2*a21+ a23 + 5*d21) - .375*c[4]*(12-k^2);
    a2  =  1.5*c[3]*(a24-2*a22) + 1.125*c[4];

    l1 = a1 + 2*(λ^2)*s1;
    l2 = a2 + 2*(λ^2)*s2;

    # ADDITIONAL TERMS FROM GEOMETRY CENTER PAPER
    b33 = -k/(16*λ)*(12*c[3]*(b21-2*k*a21+k*a23)+3*c[4]*k*(3*k^2-4)+16*s1*λ*(λ*k-1));
    b34 = -k/(8*λ)*(-12*c[3]*k*a22+3*c[4]*k+8*s2*λ*(λ*k-1));
    b35 = -k/(16*λ)*(12*c[3]*(b22+k*a24)+3*c[4]*k);

    deltan = 2 - NS;

    Ax  = sqrt( (-del - l2*Az^2)/l1 );
    omg = 1+s1*Ax^2+s2*Az^2;
    freq=λ*omg;
    period=2*pi/freq;

    rvi  = zeros(npts,6);
    ss   = zeros(npts,1);
    if npts > 1
       dtau1= 2*pi/(npts-1);
    else
       dtau1= 2*pi;
    end
    tau1 = 0;
    for i=1:npts
       x = a21*Ax^2 + a22*Az^2 - Ax*cos(tau1) + (a23*Ax^2 - a24*Az^2)*cos(2*tau1) + (a31*Ax^3 - a32*Ax*Az^2)*cos(3*tau1);
       y = k*Ax*sin(tau1) + (b21*Ax^2 - b22*Az^2)*sin(2*tau1) + (b31*Ax^3 - b32*Ax*Az^2)*sin(3*tau1);
       z = deltan*Az*cos(tau1) + deltan*d21*Ax*Az*(cos(2*tau1) - 3) + deltan*(d32*Az*Ax^2 - d31*Az^3)*cos(3*tau1);
       y_plus = (b33*Ax^3 + b34*Ax*Az^2 - b35*Ax*Az^2)*sin(tau1);
       y = y + y_plus;     # ADD EXTRA TERMS FROM G.C. PAPER

       xdot = freq*Ax*sin(tau1) - 2*freq*(a23*Ax^2-a24*Az^2)*sin(2*tau1) - 3*freq*(a31*Ax^3 - a32*Ax*Az^2)*sin(3*tau1);
       ydot = freq*(k*Ax*cos(tau1) + 2*(b21*Ax^2 - b22*Az^2)*cos(2*tau1) + 3*(b31*Ax^3 - b32*Ax*Az^2)*cos(3*tau1));
       zdot = - freq*deltan*Az*sin(tau1) - 2*freq*deltan*d21*Ax*Az*sin(2*tau1) - 3*freq*deltan*(d32*Az*Ax^2 - d31*Az^3)*sin(3*tau1);
       ydot_plus = freq*(b33*Ax^3 + b34*Ax*Az^2 - b35*Ax*Az^2)*cos(tau1);
       ydot = ydot_plus + ydot; # ADD EXTRA TERMS FROM G.C. PAPER

       rvi[i,:]= γ*[ (primary+γ*(-won+x))/γ, y, z, xdot, ydot, zdot];
       ss[i]   = tau1/freq;
       tau1=tau1+dtau1;
    end

    return ss, rvi, period, Ax
end


text:
module CIFAR

# package code goes here
export traindata,
       testdata

include("get.jl")
include("load.jl")

end # module


text:
using DEMC
using Random
using LinearAlgebra
using Distributions
Random.seed!(31953150)
# set up target distribution: Multivariate Normal
ndim = 10 # Number of dimensions
μ = rand(ndim) # mean of each dimension
# log objective function
log_obj(mean) = -sum((mean.-μ).^2)

# set up of DEMCz chain
Npar = length(μ)
blockindex = [1:Npar] # parameter blocks: here choose all parameters to be updated simultaenously
Nblocks = length(blockindex)
eps_scale = 1e-5*ones(Npar) # scale of random error around DE update
γ = 2.38 # scale of DE update, 2.38 is the "optimal" number for a normal distribution
N = 5 # number of chains
K = 10 # every K steps add current N draws to Z
Z = randn((10*ndim, ndim)) # initial distribution (completely off to make a difficult test case)

# Number of iterations in Chain
Ngen = 5000
mc, Z = DEMC.demcz_anneal(log_obj, Z, N, K, Ngen, Nblocks, blockindex, eps_scale, γ; verbose=false, TN = 1e-4, T0 = 2)

# check best element
bestval = maximum(mc.log_obj)
bestel = findfirst(mc.log_obj.==bestval)
bestpar = mc.chain[bestel[1], :, bestel[2]]

@test abs(bestval)>-1e-1

#bestvals = [maximum(mc.log_obj[:, ig]) for ig = 1:Ngen]
#DEMC.plot(1:Ngen, bestvals)


text:
t(t0, tf, dt) = t0:dt:tf
t(dom::Domain) = t(dom.initial_year,dom.final_year,dom.dt)
t(params::ClimateModelParameters) = t(params.domain)
t(m::ClimateModel) = t(m.domain)

past_mask(t, present_year) = t .< present_year
past_mask(model) = past_mask(t(model), model.domain.present_year)


text:
abstract type TestRigidBody{B,R} <: RigidBody{R} end

############################################################################################
#                            OVERRIDE ALL ROTATION STUFF
############################################################################################

state_diff_size(model::TestRigidBody{false,<:UnitQuaternion}) = 13
state_diff_size(model::TestRigidBody{false}) = 12

function state_diff(model::TestRigidBody{false}, x::SVector, x0::SVector)
    x - x0
end

function state_diff_jacobian(model::TestRigidBody{false}, x::SVector)
    return I
end

function TrajectoryOptimization.∇²differential(model::TestRigidBody{false},
        x::SVector, dx::SVector)
    return I*0
end

function TrajectoryOptimization.inverse_map_jacobian(model::TestRigidBody{false}, x::SVector)
    return I
end

function TrajectoryOptimization.inverse_map_∇jacobian(model::TestRigidBody{false},
        x::SVector, b::SVector)
    return I*0
end


text:
using GeometricIntegrators
using GeometricExamples
using ImportMacros

@import GeometricProblems.LotkaVolterra2d as problem

include("lotka_volterra_2d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

equ = problem.lotka_volterra_2d_ode()

include("lotka_volterra_2d.jl")

run_tableau_list(tableau_list, equ, Δt, ntime)

generate_html(tableau_list, figures, figures_detail, prefix, fig_suff, fig_dir, out_dir, html_file;
              title=coll_title, params=get_params(Δt, ntime, problem))


text:
using XKPasswd
using Base.Test

srand(0)

@test XKPasswd.pwentropy(4, 2048, false) == 44.0
@test XKPasswd.generate(4) == ["simple curtain forgive society"]
@test XKPasswd.generate(4)[1] == "stay we besides someone"
@test XKPasswd.generate(8,npws=7,capitalize=true,delimstr="-")[3] == "Complete-Thread-Photograph-Everywhere-Steep-Paste-Over-Bag"
@test XKPasswd.generate(6, joinpath(dirname(@__FILE__), "..", "data", "jargon.txt"), append_digit=false) == ["welldesigned trigger tales decades historically listed"]
@test XKPasswd.generate(4, XKPasswd.google_10k_clean)[1] == "preventing feelings paying disorder"


text:
###########################################################################################
#                                    HNSW Algorithm                                       #
###########################################################################################

"""
    insert_point!(hnsw, q, l = get_random_level(hnsw.lgraph))
Insert index `query` referring to data point `data[q]` into the graph.
"""
function insert_point!(hnsw, query, l = get_random_level(hnsw.lgraph))
    add_vertex!(hnsw.lgraph, query, l)

    # Get enterpoint and highest level
    enter_point = get_enter_point(hnsw)
    L =  get_entry_level(hnsw)

    # Special Case for the very first entry
    if enter_point == 0
        set_enter_point!(hnsw, query)
        return nothing
    end

    ep = Neighbor(enter_point, distance(hnsw, enter_point, query))

    # Traverse through levels to l (assuming l < L)
    for level ∈ L:-1:l+1
        W = search_layer(hnsw, query, ep, 1,level)
        ep = nearest(W) #nearest element from q in W
    end

    # Insert query on all levels < min(L,l)
    for level ∈ min(L,l):-1:1
        W = search_layer(hnsw, query, ep, hnsw.efConstruction, level)
        add_connections!(hnsw, level, query, W)
        ep = nearest(W)
    end

    # Update enter point if inserted point has highest layer
    l > L && set_enter_point!(hnsw, query)
end

function search_layer(hnsw, query, enter_point, num_points, level)
    vl = get_list(hnsw.vlp) #Acquire VisitedList
    visit!(vl, enter_point)
    C = NeighborSet(enter_point) #set of candidates
    W = NeighborSet(enter_point) #dynamic list of found nearest neighbors
    while length(C) > 0
        c = pop_nearest!(C) # from query in C
        c.dist > furthest(W).dist && break #Stopping condition
        for e ∈ neighbors(hnsw.lgraph, level, c)
            if !isvisited(vl, e)
                visit!(vl, e)
                eN = Neighbor(e, distance(hnsw,query,e))
                if length(W) < num_points || eN.dist < furthest(W).dist
                    insert!(C, eN)
                    insert!(W, eN) #add optional maxlength feature?
                    length(W) > num_points && pop_furthest!(W)
                end
            end
        end
    end
    release_list(hnsw.vlp, vl)
    return W #num_points closest neighbors
end

function neighbor_heuristic(hnsw, level, candidates)
    M = max_connections(hnsw.lgraph, level)
    length(candidates) <= M  && return candidates

    chosen = typeof(candidates)()
    for e ∈ candidates
        length(chosen) < M  || break
        #Heuristic:
        good = true
        for r ∈ chosen
            if e.dist > distance(hnsw, e.idx, r.idx)
                good=false
                break
            end
        end
        good && insert!(chosen, e)
    end
    return chosen
end

###########################################################################################
#                                       KNN Search                                        #
###########################################################################################
function knn_search(hnsw, q, K)
    ef = max(K, hnsw.ef)
    @assert length(q)==length(hnsw.data[1])
    ep = get_enter_point(hnsw)
    epN = Neighbor(ep, distance(hnsw, q, ep))
    L = get_entry_level(hnsw) #layer of ep , top layer of hnsw
    for level ∈ L:-1:2 # Iterate from top to second lowest
        epN = search_layer(hnsw, q, epN, 1, level)[1]
    end
    W = search_layer(hnsw, q, epN, ef, 1)
    list = nearest(W, K)
    idx = map(x->x.idx, list)
    dist = map(x->x.dist, list)
    return idx, dist# K nearest elements to q
end

function knn_search(hnsw::HierarchicalNSW{T,F},
        q::AbstractVector{<:AbstractVector}, # query
        K) where {T,F}
    idxs = Vector{Vector{T}}(undef,length(q))
    dists = Vector{Vector{F}}(undef,length(q))
    @inbounds for n = 1:length(q)
        idxs[n], dists[n] = knn_search(hnsw, q[n], K)
    end
    idxs, dists
end


text:
import Base.==

mutable struct Basis{O, V, P} <: abstractBasis
    basis::O
    variables::V
    parameter::P
    f_
end

Base.show(io::IO, x::Basis) = print(io, "$(length(x.basis)) dimensional basis in ", "$(String.([v.op.name for v in x.variables]))")

@inline function Base.print(io::IO, x::Basis)
    show(io, x)
    println()
    if length(x.variables) == length(x.basis)
        for (i, bi) in enumerate(x.basis)
            println("d$(x.variables[i]) = $bi")
        end
    else
        for (i, bi) in enumerate(x.basis)
            println("f_$i = $bi")
        end
    end
end

is_independent(o::Operation) = isempty(o.args)

function Basis(basis::AbstractVector{Operation}, variables::AbstractVector{Operation};  parameters =  [])
    @assert all(is_independent.(variables)) "Please provide independent variables for base."

    bs = unique(basis)
    fix_single_vars_in_basis!(bs, variables)

    vs = [ModelingToolkit.Variable(Symbol(i)) for i in variables]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters]

    f_ = ModelingToolkit.build_function(bs, vs, ps, (), simplified_expr, Val{false})[1]
    return Basis(bs, variables, parameters, f_)
end

function update!(basis::Basis)
    vs = [ModelingToolkit.Variable(Symbol(i)) for i in variables(basis)]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters(basis)]

    basis.f_ = ModelingToolkit.build_function(basis.basis, vs, ps, (), simplified_expr, Val{false})[1]
    return
end

function Base.push!(b::Basis, ops::AbstractArray{Operation})
    @inbounds for o in ops
        push!(b.basis, o)
    end
    unique!(b.basis)
    update!(b)
    return
end

function Base.push!(b::Basis, op₀::Operation)
    op = simplify_constants(op₀)
    fix_single_vars_in_basis!(op, b.variables)
    push!(b.basis, op)
    # Check for uniqueness
    unique!(b)
    update!(b)
    return
end

function Base.deleteat!(b::Basis, inds)
    deleteat!(b.basis, inds)
    update!(b)
    return
end

function Base.merge(basis_a::Basis, basis_b::Basis)
    b =  unique(vcat(basis_a.basis, basis_b.basis))
    vs = unique(vcat(basis_a.variables, basis_b.variables))
    ps = unique(vcat(basis_a.parameter, basis_b.parameter))
    return Basis(b, vs, parameters = ps)
end

function Base.merge!(basis_a::Basis, basis_b::Basis)
    push!(basis_a, basis_b.basis)
    basis_a.variables = unique(vcat(basis_a.variables, basis_b.variables))
    basis_a.parameter = unique(vcat(basis_a.parameter, basis_b.parameter))
    update!(basis_a)
    return
end

Base.getindex(b::Basis, idx::Int64) = b.basis[idx]
Base.getindex(b::Basis, ids::UnitRange{Int64}) = b.basis[ids]
Base.getindex(b::Basis, ::Colon) = b.basis
Base.firstindex(b::Basis) = firstindex(b.basis)
Base.lastindex(b::Basis) = lastindex(b.basis)
Base.iterate(b::Basis) = iterate(b.basis)
Base.iterate(b::Basis, id::Int64) = iterate(b.basis, id)

function (==)(x::Basis, y::Basis)
    n = zeros(Bool, length(x.basis))
    @inbounds for (i, xi) in enumerate(x)
        n[i] = any(isequal.(xi, y.basis))
    end
    return all(n)
end

function count_operation(o::Expression, ops::AbstractArray)
    if isa(o, ModelingToolkit.Constant)
        return 0
    end
    k = o.op ∈ ops ? 1 : 0
    if !isempty(o.args)
        k += sum([count_operation(ai, ops) for ai in o.args])
    end
    return k
end

free_parameters(b::Basis; operations = [+]) = sum([count_operation(bi, operations) for bi in b.basis]) + length(b.basis)

(b::Basis)(x::AbstractArray{T, 1}; p::AbstractArray = []) where T <: Number = b.f_(x, isempty(p) ? parameters(b) : p)

function (b::Basis)(x::AbstractArray{T, 2}; p::AbstractArray = []) where T <: Number
    if (isempty(p) || eltype(p) <: Expression) && !isempty(parameters(b))
        pi = isempty(p) ? parameters(b) : p
        res = zeros(eltype(pi), length(b), size(x)[2])
    else
        pi = p
        res = zeros(eltype(x), length(b), size(x)[2])
    end
    @inbounds for i in 1:size(x)[2]
        res[:, i] .= b.f_(x[:, i], pi)
    end
    return res
end

Base.size(b::Basis) = size(b.basis)
Base.length(b::Basis) = length(b.basis)
ModelingToolkit.parameters(b::Basis) = b.parameter
variables(b::Basis) = b.variables

function jacobian(basis::Basis)

    vs = [ModelingToolkit.Variable(Symbol(i)) for i in variables(basis)]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters(basis)]

    j = calculate_jacobian(basis.basis, variables(basis))
    return ModelingToolkit.build_function(expand_derivatives.(j), vs, ps, (), simplified_expr, Val{false})[1]
end

function Base.unique!(b::Basis)
    N = length(b.basis)
    removes = Vector{Bool}()
    for i ∈ 1:N
        push!(removes, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    deleteat!(b, removes)
    update!(b)
end

function Base.unique(b::Basis)
    N = length(b.basis)
    returns = Vector{Bool}()
    for i ∈ 1:N
        push!(returns, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    returns = [!r for r in returns]
    return Basis(b.basis[returns], variables(b), parameters = parameters(b))
end

function Base.unique(b₀::AbstractVector{Operation})
    b = simplify_constants.(b₀)
    N = length(b)
    returns = Vector{Bool}()
    for i ∈ 1:N
        push!(returns, any([isequal(b[i], b[j]) for j in i+1:N]))
    end
    returns = [!r for r in returns]
    return b[returns]
end

function Base.unique!(b::AbstractArray{Operation})
    N = length(b)
    removes = Vector{Bool}()
    for i ∈ 1:N
        push!(removes, any([isequal(b[i], b[j]) for j in i+1:N]))
    end
    deleteat!(b, removes)
end

function fix_single_vars_in_basis!(basis,variables)
    for (ind, el) in enumerate(basis)
        for (ind_var, var) in enumerate(variables)
            if isequal(el,var)
                basis[ind] = 1var
            end
        end
    end
end

function dynamics(b::Basis)
    return b.f_
end

function ModelingToolkit.ODESystem(b::Basis)
    @assert length(b) == length(variables(b))
    # Define the time
    @parameters t
    @derivatives D'~t

    vs = similar(b.variables)
    dvs = similar(b.variables)
    for (i, vi) in enumerate(b.variables)
        vs[i] = ModelingToolkit.Operation(vi.op, [t])
        dvs[i] = D(vs[i])
    end
    eqs = dvs .~ b(vs, p = b.parameter)
    return ODESystem(eqs)
end

function ModelingToolkit.ODESystem(b::Basis, independent_variable::Operation)
    @assert length(b) == length(variables(b))-1
    @derivatives D'~independent_variable

    vars = [vi for vi in variables(b) if ! isequal(vi, independent_variable)]

    vs = similar(vars)
    dvs = similar(vars)

    for (i, vi) in enumerate(vars)
        vs[i] = ModelingToolkit.Operation(vi.op, [independent_variable])
        dvs[i] = D(vs[i])
    end

    eqs = dvs .~ b([vs..., independent_variable], p = b.parameter)
    return ODESystem(eqs, independent_variable, vs, b.parameter)
end


text:
# Implements the WENO method used in Dumbser et al (DOI 10.1016/j.cma.2013.09.022)
module Weno

export coefficient_matrices, oscillation_indicator, weno

import Polynomials: polyint, polyder

import Basis: basis
include("array_manipulation.jl")

function coefficient_matrices(N)
  # Generate linear systems governing the coefficients of the basis polynomials
  fHalfN = floor(N/2)
  cHalfN = ceil(N/2)
  M1 = zeros(N+1, N+1); M2 = zeros(N+1, N+1); M3 = zeros(N+1, N+1); M4 = zeros(N+1, N+1);
  ψ = basis(N)

  for p=1:N+1, a=1:N+1
    ψintp = polyint(ψ[p])
    M1[a,p] = ψintp(a-N)      - ψintp(a-N-1)
    M2[a,p] = ψintp(a)        - ψintp(a-1)
    M3[a,p] = ψintp(a-fHalfN) - ψintp(a-fHalfN-1)
    if Bool(N%2)
      M4[a,p] = ψintp(a-cHalfN) - ψintp(a-cHalfN-1)
    end
  end
  return factorize(M1), factorize(M2), factorize(M3), factorize(M4)
end

function oscillation_indicator(N)
  # Generate the oscillation indicator matrix
  Σ = zeros(N+1, N+1)
  ψ = basis(N)

  for a=1:N
    ψdera = [polyder(ψp, a) for ψp in ψ]
    for p=1:N+1, m=1:N+1
      antiderivative = polyint(ψdera[p] * ψdera[m])
      Σ[p,m] += antiderivative(1) - antiderivative(0)
    end
  end
  return Σ
end

@inline function indicator(chΣT, w, n)
  # The indicator is w' * Σ * w, ie |w' * chol(Σ)'|^2
  # NB The Cholesky decomposition always exists as Σ is SPD
  ret = 0.
  for j=1:n
    tmp = 0.
    for i=j:n
      @inbounds tmp += w[i] * chΣT[i,j]
    end
    ret += tmp^2
  end
  return ret
end

function coeffs3(ret, N, w1,w2,w3, M1,M2,M3, λs,λc, chΣT, ε, r)
  # Calculate coefficients of basis polynomials and weights
  A_ldiv_B!(M1, w1)
  A_ldiv_B!(M2, w2)
  A_ldiv_B!(M3, w3)
  o1 = λs / (indicator(chΣT, w1, N+1) + ε)^r
  o2 = λs / (indicator(chΣT, w2, N+1) + ε)^r
  o3 = λc / (indicator(chΣT, w3, N+1) + ε)^r
  oSum = o1 + o2 + o3
  for i = 1:N+1
    @inbounds ret[i] = (o1*w1[i] + o2*w2[i] + o3*w3[i]) / oSum
  end
end

function coeffs4(ret, N, w1,w2,w3,w4, M1,M2,M3,M4, λs,λc, chΣT, ε, r)
  # Calculate coefficients of basis polynomials and weights
  A_ldiv_B!(M1, w1)
  A_ldiv_B!(M2, w2)
  A_ldiv_B!(M3, w3)
  A_ldiv_B!(M4, w4)
  o1 = λs / (indicator(chΣT, w1, N+1) + ε)^r
  o2 = λs / (indicator(chΣT, w2, N+1) + ε)^r
  o3 = λc / (indicator(chΣT, w3, N+1) + ε)^r
  o4 = λc / (indicator(chΣT, w4, N+1) + ε)^r
  oSum = o1 + o2 + o3 + o4
  for i = 1:N+1
    @inbounds ret[i] = (o1*w1[i] + o2*w2[i] + o3*w3[i] + o4*w4[i]) / oSum
  end
end

function weno(u, N, M1=nothing, M2=nothing, M3=nothing, M4=nothing, chΣT=nothing,
              λc=1e5, λs=1, r=8, ε=1e-14)
  # Find reconstruction coefficients of u to order N+1

  if M1==nothing
    M1,M2,M3,M4 = coefficient_matrices(N)
  end
  if chΣT==nothing
    Σ = oscillation_indicator(N)
    chΣT = chol(Σ)'
  end

  nx, ny, nz, nvar = size(u)
  fHalfN = Int(floor(N/2))
  cHalfN = Int(ceil(N/2))
  oddN = Bool(N%2)

  w1 = Array{Float64}(N+1)
  w2 = Array{Float64}(N+1)
  w3 = Array{Float64}(N+1)
  w4 = Array{Float64}(N+1)

  Wx = zeros(nx, ny, nz, N+1, nvar)
  Wx0 = extend(u, N, nx, ny, nz, nvar, 1)

  for v=1:nvar, k=1:nz, j=1:ny
    w0 = view(Wx0, :, j, k, v)
    for i=1:nx
      for ind = 0:N
        @inbounds w1[ind+1] = w0[i+ind]
        @inbounds w2[ind+1] = w0[i+N+ind]
        @inbounds w3[ind+1] = w0[i+cHalfN+ind]
      end
      if oddN
        for ind = 0:N
          @inbounds w4[ind+1] = w0[i+fHalfN+ind]
        end
        coeffs4(view(Wx,i,j,k,:,v), N, w1,w2,w3,w4, M1,M2,M3,M4, λs,λc, chΣT,ε,r)
      else
        coeffs3(view(Wx,i,j,k,:,v), N, w1,w2,w3, M1,M2,M3, λs,λc, chΣT,ε,r)
      end
    end
  end
  if ny==1 && nz==1
    return Wx
  end

  Wxy = zeros(nx, ny, nz, N+1, N+1, nvar)
  Wxy0 = extend(Wx, N, nx, ny, nz, nvar, 2)

  for v=1:nvar, a=1:N+1, k=1:nz, i=1:nx
    w0 = view(Wxy0, i, :, k, a, v)
    for j=1:ny
      for ind = 0:N
        @inbounds w1[ind+1] = w0[j+ind]
        @inbounds w2[ind+1] = w0[j+N+ind]
        @inbounds w3[ind+1] = w0[j+cHalfN+ind]
      end
      if oddN
        for ind = 0:N
          @inbounds w4[ind+1] = w0[j+fHalfN+ind]
        end
        coeffs4(view(Wxy,i,j,k,a,:,v), N, w1,w2,w3,w4, M1,M2,M3,M4, λs,λc, chΣT,ε,r)
      else
        coeffs3(view(Wxy,i,j,k,a,:,v), N, w1,w2,w3, M1,M2,M3, λs,λc, chΣT,ε,r)
      end
    end
  end
  if nz==1
    return Wxy
  end

  Wxyz = zeros(nx, ny, nz, N+1, N+1, N+1, nvar)
  Wxyz0 = extend(Wxy, N, nx, ny, nz, nvar, 3)

  for v=1:nvar, b=1:N+1, a=1:N+1, j=1:ny, i=1:nx
    w0 = view(Wxyz0, i, j, :, a, b, v)
    for k=1:nz
      for ind = 0:N
        @inbounds w1[ind+1] = w0[k+ind]
        @inbounds w2[ind+1] = w0[k+N+ind]
        @inbounds w3[ind+1] = w0[k+cHalfN+ind]
      end
      if oddN
        for ind = 0:N
          @inbounds w4[ind+1] = w0[k+fHalfN+ind]
        end
        coeffs4(view(Wxyz,i,j,k,a,b,:,v), N, w1,w2,w3,w4, M1,M2,M3,M4, λs,λc, chΣT,ε,r)
      else
        coeffs3(view(Wxyz,i,j,k,a,b,:,v), N, w1,w2,w3, M1,M2,M3, λs,λc, chΣT,ε,r)
      end
    end
  end
  return Wxyz

end

end   # module Weno

max_stars_repo_name: haranjackson/Julia-WENO


text:
using ClimateMachine, JLD2, GLMakie, MPI
ClimateMachine.init()
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Ocean
using ClimateMachine.Ocean.Domains
using ClimateMachine.Ocean.Fields
using ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
using ClimateMachine.GenericCallbacks: EveryXSimulationTime
using ClimateMachine.GenericCallbacks: EveryXSimulationSteps
using ClimateMachine.Ocean: current_step, Δt, current_time
using ClimateMachine.Ocean: JLD2Writer, OutputTimeSeries, write!


text:
using WaterLily
using LinearAlgebra: norm2
include("TwoD_plots.jl")
include("TwoD_circle.jl")
using SmoothLivePlot

function v_plot(t,v,i=length(t))
    sleep(0.001)
    plot(t[1:i],v[1:i],xlims=(first(t),last(t)),legend=false)
    scatter!(t[i:i],v[i:i])
    plot!(xaxis=("time"),yaxis=("cross-stream velocity"))
end
"""
    sim_measure!(sim, I; duration=1, step=0.1)

Example function to run a simulation `sim` over time `duration`,
measuring the velocity `u[I]` every `step` and visualize `using
SmoothLivePlot`. The function modifies `sim` and returns the time
and velocity history `t,v`.

# Examples
```jldoctest
julia> sim, I = circle(128,64,Re=250), CartesianIndex(60,40,2);
julia> t, v = sim_measure!(sim, I, duration=100);
```
"""
function sim_measure!(sim,I;duration=1,step=0.1)
    t₀ = round(sim_time(sim))
    t = range(t₀,t₀+duration;step)
    v = Vector{Float64}(undef,length(t))
    plt = @makeLivePlot v_plot(t,v,1)
    for i ∈ 1:length(t)
        sim_step!(sim,t[i])
        v[i] = sim.flow.u[I]
        modifyPlotObject!(plt,arg2=v,arg3=i)
    end
    t,v
end


text:
using PyPlot
using HDF5, JLD
using Statistics

using SubmodularMaximization

pygui(false)
name = "compare_solvers"
fig_path = "./fig/$name"
data_path = "./data/$name"
mkpath(fig_path)
mkpath(data_path)

num_trials = 100

num_agents = 50
num_sensors = 10
nominal_area = 2.0

sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 2 * sensor_radius

agent_specification = CircleAgentSpecification(sensor_radius, station_radius,
                                         num_sensors)

f(x) = mean_area_coverage(x, 100)

solvers = Any[]

#push!(solvers, (solve_optimal, "Optimal"))
#push!(solvers, (solve_worst, "Worst-case"))
push!(solvers, (solve_random, "Random"))
push!(solvers, (solve_myopic, "Myopic"))

partitions = [2, 4, 8]
for num_partitions in partitions
  solve_n(p) = solve_n_partitions(num_partitions, p)
  push!(solvers, (solve_n, "Partition-$num_partitions"))
end

push!(solvers, (solve_sequential, "Sequential"))

results = zeros(num_trials, length(solvers))

problems = map(1:num_trials) do unused
  agents = generate_agents(agent_specification, num_agents)
  problem = ExplicitPartitionProblem(f, agents)
end
@save "$data_path/partition_matroids" map(x->x.partition_matroid, problems)

for trial_num in 1:length(problems)
  println("Trial: $trial_num")

  for solver_num in 1:length(solvers)
    solver = solvers[solver_num][1]
    name = solvers[solver_num][2]

    solution = solver(problems[trial_num])

    results[trial_num, solver_num] = solution.value
  end
end
@save "$data_path/results" results

# now analyze weights
weight_matrices = map(problems, 1:length(problems)) do problem, ii
  println("Weight matrix $ii")
  compute_weight_matrix(problem)
end

@save "$data_path/weights" weight_matrices

total_weights = map(weight_matrices, 1:length(problems)) do weight_matrix, ii
  println("Total weights $ii")
  total_weight(weight_matrix)
end

edge_sets = map(weight_matrices, 1:length(problems)) do weight_matrix, ii
  println("Triangle $ii")
  extract_triangle(weight_matrix)
end

figure()
boxplot(results, notch=false, vert=false)
yticks(1:length(solvers), map(x->x[2], solvers))
xlabel("Area coverage")
tight_layout()

save_fig(fig_path, "results")

# plot histograms of edge weights

figure()
PyPlot.plt.hist(total_weights, 20)
tt = "Total Graph Weight Frequency"
ylabel("Frequency")
xlabel("Redundancy graph weight")
tight_layout()

save_fig(fig_path, tt)
title(tt)

figure()
PyPlot.plt.hist(vcat(edge_sets...), 20)
tt = "Edge Weight Frequency"
ylabel("Frequency")
xlabel("Edge weight (redundancy)")
tight_layout()

save_fig(fig_path, tt)
title(tt)

# sequential solutions and edge weights
sequential_mean = mean(results[:,end][:])

partition_values = mean(results[:,3:end-1]; dims=1)[:]

mean_weight = mean(total_weights)
println("Mean weight: $mean_weight")
bounds = map(x->mean_weight/x, partitions)

#figure()
#plot([partitions[1], partitions[end]], [sequential_mean, sequential_mean])
#plot(partitions, partition_values)
#plot(partitions, partition_values + bounds)
#legend(["Sequential", "Partition-n"])

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing