text:
using Random
using BenchmarkTools
using LightOSM
using OpenStreetMapX
using Graphs
using DataStructures
using JSON

"""
Setup
"""

Random.seed!(1234) # Set seed so experiment is reproducible

point = GeoLocation(-37.8142176, 144.9631608) # Melbourne, Australia
radius = 5 # km
data_file = joinpath(@__DIR__, "benchmark_map.osm")

g_losm = LightOSM.graph_from_download(:point,
                                      point=point,
                                      radius=radius,
                                      weight_type=:distance,
                                      save_to_file_location=data_file)

# g_losm = LightOSM.graph_from_file(data_file, weight_type=:distance)
@time g_losm_precompute = LightOSM.graph_from_file(data_file, weight_type=:distance, precompute_dijkstra_states=true)

g_osmx = OpenStreetMapX.get_map_data(data_file,
                                     use_cache=false,
                                     only_intersections=false,
                                     trim_to_connected_graph=true)

"""
Define benchmark functions
"""

function losm_shortest_path(g::LightOSM.OSMGraph, o_d_nodes, algorithm)
    for (o, d) in o_d_nodes
        try
            LightOSM.shortest_path(g, o, d, algorithm=algorithm)
        catch
            # Error exception will be thrown if path does not exist from origin to destination node
        end
    end
end

function osmx_shortest_path(g::OpenStreetMapX.MapData, o_d_nodes, algorithm)
    for (o, d) in o_d_nodes
        try
            OpenStreetMapX.find_route(g, o, d, g.w, routing=algorithm, heuristic=(u, v) -> OpenStreetMapX.get_distance(u, v, g_osmx.nodes, g_osmx.n), get_distance=false, get_time=false)
        catch
            # Error exception will be thrown if path does not exist from origin to destination node
        end
    end
end

function lg_shortest_path(g::LightOSM.OSMGraph, o_d_indices, algorithm)
    if algorithm == :astar
        for (o, d) in o_d_indices
            try
                Graphs.a_star(g.graph, o, d, g.weights)
            catch
                # Error exception will be thrown if path does not exist from origin to destination node
            end
        end
    elseif algorithm == :dijkstra
        for (o, d) in o_d_indices
            try
                state = Graphs.dijkstra_shortest_paths(g.graph, o, g.weights)
                Graphs.enumerate_paths(state, d)
            catch
                # Error exception will be thrown if path does not exist from origin to destination node
            end
        end
    end
end

function extract(bmark_trial)
    io = IOBuffer()
    show(io, "text/plain", bmark_trial)

    s = String(take!(io))
    s = split.((split(s, "\n")), ":")

    result = Dict{String,String}()

    for item in s
        if length(item) >= 2
            result[strip(item[1])] = strip(item[2])
        end
    end
    @info "Extracted benchmark result: $result"
    return result
end

"""
Define experiment
"""

n_paths = [1, 10, 100, 1000, 10000]
o_d_indices = OrderedDict()
o_d_nodes = OrderedDict()

losm_nodes = collect(keys(g_losm.nodes))
osmx_nodes = collect(keys(g_osmx.nodes))
common_nodes = intersect(losm_nodes, osmx_nodes)

for n in n_paths
    rand_o_d_indices = rand(1:length(common_nodes), n, 2)
    o_d_indices[n] = [[o, d] for (o, d) in eachrow(rand_o_d_indices) if o != d]
    o_d_nodes[n] = [[common_nodes[o], common_nodes[d]] for (o, d) in eachrow(rand_o_d_indices) if o != d]
end

"""
Run experiment
"""
results = Dict(:dijkstra => DefaultDict(Dict), :astar => DefaultDict(Dict))

results[:dijkstra][1][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[1], :dijkstra))
results[:dijkstra][1][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[1], :dijkstra))
results[:dijkstra][1][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[1], :dijkstra))
results[:dijkstra][1][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[1], :dijkstra))

results[:dijkstra][10][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[10], :dijkstra))
results[:dijkstra][10][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[10], :dijkstra))
results[:dijkstra][10][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[10], :dijkstra))
results[:dijkstra][10][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[10], :dijkstra))

results[:dijkstra][100][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[100], :dijkstra))
results[:dijkstra][100][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[100], :dijkstra))
results[:dijkstra][100][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[100], :dijkstra))
results[:dijkstra][100][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[100], :dijkstra))

results[:dijkstra][1000][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[1000], :dijkstra))
results[:dijkstra][1000][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[1000], :dijkstra))
results[:dijkstra][1000][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[1000], :dijkstra))
results[:dijkstra][1000][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[1000], :dijkstra))

results[:dijkstra][10000][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[10000], :dijkstra))
results[:dijkstra][10000][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[10000], :dijkstra))
results[:dijkstra][10000][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[10000], :dijkstra))
results[:dijkstra][10000][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[10000], :dijkstra))

results[:astar][1][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[1], :astar))
results[:astar][1][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[1], :astar))
results[:astar][1][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[1], :astar))
results[:astar][1][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[1], :astar))

results[:astar][10][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[10], :astar))
results[:astar][10][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[10], :astar))
results[:astar][10][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[10], :astar))
results[:astar][10][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[10], :astar))

results[:astar][100][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[100], :astar))
results[:astar][100][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[100], :astar))
results[:astar][100][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[100], :astar))
results[:astar][100][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[100], :astar))

results[:astar][1000][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[1000], :astar))
results[:astar][1000][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[1000], :astar))
results[:astar][1000][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[1000], :astar))
results[:astar][1000][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[1000], :astar))

results[:astar][10000][:losm] = extract(@benchmark losm_shortest_path(g_losm, o_d_nodes[10000], :astar))
results[:astar][10000][:losm_precompute] = extract(@benchmark losm_shortest_path(g_losm_precompute, o_d_nodes[10000], :astar))
results[:astar][10000][:osmx] = extract(@benchmark osmx_shortest_path(g_osmx, o_d_nodes[10000], :astar))
results[:astar][10000][:lg] = extract(@benchmark lg_shortest_path(g_losm, o_d_indices[10000], :astar))

"""
Export Results
"""

open(joinpath(@__DIR__, "benchmark_results.json"), "w") do io
    write(io, json(results))
end


text:
#==============================================================================#
# Polly.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/polly-2016-06-10.normal.json
#==============================================================================#

__precompile__()

module Polly

using AWSCore

"""
    using AWSSDK.Polly.delete_lexicon
    delete_lexicon([::AWSConfig], arguments::Dict)
    delete_lexicon([::AWSConfig]; LexiconName=)

    using AWSCore.Services.polly
    polly([::AWSConfig], "DELETE", "/v1/lexicons/{LexiconName}", arguments::Dict)
    polly([::AWSConfig], "DELETE", "/v1/lexicons/{LexiconName}", LexiconName=)

# DeleteLexicon Operation

Deletes the specified pronunciation lexicon stored in an AWS Region. A lexicon which has been deleted is not available for speech synthesis, nor is it possible to retrieve it using either the `GetLexicon` or `ListLexicon` APIs.

For more information, see [Managing Lexicons](http://docs.aws.amazon.com/polly/latest/dg/managing-lexicons.html).

# Arguments

## `LexiconName = ::String` -- *Required*
The name of the lexicon to delete. Must be an existing lexicon in the region.

# Returns

`DeleteLexiconOutput`

# Exceptions

`LexiconNotFoundException` or `ServiceFailureException`.

# Example: To delete a lexicon

Deletes a specified pronunciation lexicon stored in an AWS Region.

Input:
```
[
    "Name" => "example"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/DeleteLexicon)
"""
@inline delete_lexicon(aws::AWSConfig=default_aws_config(); args...) = delete_lexicon(aws, args)

@inline delete_lexicon(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "DELETE", "/v1/lexicons/{LexiconName}", args)

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

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

    using AWSCore.Services.polly
    polly([::AWSConfig], "GET", "/v1/voices", arguments::Dict)
    polly([::AWSConfig], "GET", "/v1/voices", <keyword arguments>)

# DescribeVoices Operation

Returns the list of voices that are available for use when requesting speech synthesis. Each voice speaks a specified language, is either male or female, and is identified by an ID, which is the ASCII version of the voice name.

When synthesizing speech ( `SynthesizeSpeech` ), you provide the voice ID for the voice you want from the list of voices returned by `DescribeVoices`.

For example, you want your news reader application to read news in a specific language, but giving a user the option to choose the voice. Using the `DescribeVoices` operation you can provide the user with a list of available voices to select from.

You can optionally specify a language code to filter the available voices. For example, if you specify `en-US`, the operation returns a list of all available US English voices.

This operation requires permissions to perform the `polly:DescribeVoices` action.

# Arguments

## `LanguageCode = "cy-GB", "da-DK", "de-DE", "en-AU", "en-GB", "en-GB-WLS", "en-IN", "en-US", "es-ES", "es-US", "fr-CA", "fr-FR", "is-IS", "it-IT", "ko-KR", "ja-JP", "nb-NO", "nl-NL", "pl-PL", "pt-BR", "pt-PT", "ro-RO", "ru-RU", "sv-SE" or "tr-TR"`
The language identification tag (ISO 639 code for the language name-ISO 3166 country code) for filtering the list of voices returned. If you don't specify this optional parameter, all available voices are returned.

## `NextToken = ::String`
An opaque pagination token returned from the previous `DescribeVoices` operation. If present, this indicates where to continue the listing.

# Returns

`DescribeVoicesOutput`

# Exceptions

`InvalidNextTokenException` or `ServiceFailureException`.

# Example: To describe available voices

Returns the list of voices that are available for use when requesting speech synthesis. Displayed languages are those within the specified language code. If no language code is specified, voices for all available languages are displayed.

Input:
```
[
    "LanguageCode" => "en-GB"
]
```

Output:
```
Dict(
    "Voices" => [
        Dict(
            "Gender" => "Female",
            "Id" => "Emma",
            "LanguageCode" => "en-GB",
            "LanguageName" => "British English",
            "Name" => "Emma"
        ),
        Dict(
            "Gender" => "Male",
            "Id" => "Brian",
            "LanguageCode" => "en-GB",
            "LanguageName" => "British English",
            "Name" => "Brian"
        ),
        Dict(
            "Gender" => "Female",
            "Id" => "Amy",
            "LanguageCode" => "en-GB",
            "LanguageName" => "British English",
            "Name" => "Amy"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/DescribeVoices)
"""
@inline describe_voices(aws::AWSConfig=default_aws_config(); args...) = describe_voices(aws, args)

@inline describe_voices(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "GET", "/v1/voices", args)

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

"""
    using AWSSDK.Polly.get_lexicon
    get_lexicon([::AWSConfig], arguments::Dict)
    get_lexicon([::AWSConfig]; LexiconName=)

    using AWSCore.Services.polly
    polly([::AWSConfig], "GET", "/v1/lexicons/{LexiconName}", arguments::Dict)
    polly([::AWSConfig], "GET", "/v1/lexicons/{LexiconName}", LexiconName=)

# GetLexicon Operation

Returns the content of the specified pronunciation lexicon stored in an AWS Region. For more information, see [Managing Lexicons](http://docs.aws.amazon.com/polly/latest/dg/managing-lexicons.html).

# Arguments

## `LexiconName = ::String` -- *Required*
Name of the lexicon.

# Returns

`GetLexiconOutput`

# Exceptions

`LexiconNotFoundException` or `ServiceFailureException`.

# Example: To retrieve a lexicon

Returns the content of the specified pronunciation lexicon stored in an AWS Region.

Input:
```
[
    "Name" => ""
]
```

Output:
```
Dict(
    "Lexicon" => Dict(
        "Content" => "<?xml version="1.0" encoding="UTF-8"?>
<lexicon version="1.0"
      xmlns="http://www.w3.org/2005/01/pronunciation-lexicon"
      xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
      xsi:schemaLocation="http://www.w3.org/2005/01/pronunciation-lexicon
        http://www.w3.org/TR/2007/CR-pronunciation-lexicon-20071212/pls.xsd"
      alphabet="ipa"
      xml:lang="en-US">
  <lexeme>
    <grapheme>W3C</grapheme>
    <alias>World Wide Web Consortium</alias>
  </lexeme>
</lexicon>",
        "Name" => "example"
    ),
    "LexiconAttributes" => Dict(
        "Alphabet" => "ipa",
        "LanguageCode" => "en-US",
        "LastModified" => 1.478542980117e9,
        "LexemesCount" => 1,
        "LexiconArn" => "arn:aws:polly:us-east-1:123456789012:lexicon/example",
        "Size" => 503
    )
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/GetLexicon)
"""
@inline get_lexicon(aws::AWSConfig=default_aws_config(); args...) = get_lexicon(aws, args)

@inline get_lexicon(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "GET", "/v1/lexicons/{LexiconName}", args)

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

"""
    using AWSSDK.Polly.get_speech_synthesis_task
    get_speech_synthesis_task([::AWSConfig], arguments::Dict)
    get_speech_synthesis_task([::AWSConfig]; TaskId=)

    using AWSCore.Services.polly
    polly([::AWSConfig], "GET", "/v1/synthesisTasks/{TaskId}", arguments::Dict)
    polly([::AWSConfig], "GET", "/v1/synthesisTasks/{TaskId}", TaskId=)

# GetSpeechSynthesisTask Operation

Retrieves a specific SpeechSynthesisTask object based on its TaskID. This object contains information about the given speech synthesis task, including the status of the task, and a link to the S3 bucket containing the output of the task.

# Arguments

## `TaskId = ::String` -- *Required*
The Amazon Polly generated identifier for a speech synthesis task.

# Returns

`GetSpeechSynthesisTaskOutput`

# Exceptions

`InvalidTaskIdException`, `ServiceFailureException` or `SynthesisTaskNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/GetSpeechSynthesisTask)
"""
@inline get_speech_synthesis_task(aws::AWSConfig=default_aws_config(); args...) = get_speech_synthesis_task(aws, args)

@inline get_speech_synthesis_task(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "GET", "/v1/synthesisTasks/{TaskId}", args)

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

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

    using AWSCore.Services.polly
    polly([::AWSConfig], "GET", "/v1/lexicons", arguments::Dict)
    polly([::AWSConfig], "GET", "/v1/lexicons", <keyword arguments>)

# ListLexicons Operation

Returns a list of pronunciation lexicons stored in an AWS Region. For more information, see [Managing Lexicons](http://docs.aws.amazon.com/polly/latest/dg/managing-lexicons.html).

# Arguments

## `NextToken = ::String`
An opaque pagination token returned from previous `ListLexicons` operation. If present, indicates where to continue the list of lexicons.

# Returns

`ListLexiconsOutput`

# Exceptions

`InvalidNextTokenException` or `ServiceFailureException`.

# Example: To list all lexicons in a region

Returns a list of pronunciation lexicons stored in an AWS Region.

Input:
```
[

]
```

Output:
```
Dict(
    "Lexicons" => [
        Dict(
            "Attributes" => Dict(
                "Alphabet" => "ipa",
                "LanguageCode" => "en-US",
                "LastModified" => 1.478542980117e9,
                "LexemesCount" => 1,
                "LexiconArn" => "arn:aws:polly:us-east-1:123456789012:lexicon/example",
                "Size" => 503
            ),
            "Name" => "example"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/ListLexicons)
"""
@inline list_lexicons(aws::AWSConfig=default_aws_config(); args...) = list_lexicons(aws, args)

@inline list_lexicons(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "GET", "/v1/lexicons", args)

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

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

    using AWSCore.Services.polly
    polly([::AWSConfig], "GET", "/v1/synthesisTasks", arguments::Dict)
    polly([::AWSConfig], "GET", "/v1/synthesisTasks", <keyword arguments>)

# ListSpeechSynthesisTasks Operation

Returns a list of SpeechSynthesisTask objects ordered by their creation date. This operation can filter the tasks by their status, for example, allowing users to list only tasks that are completed.

# Arguments

## `MaxResults = ::Int`
Maximum number of speech synthesis tasks returned in a List operation.

## `NextToken = ::String`
The pagination token to use in the next request to continue the listing of speech synthesis tasks.

## `Status = "scheduled", "inProgress", "completed" or "failed"`
Status of the speech synthesis tasks returned in a List operation

# Returns

`ListSpeechSynthesisTasksOutput`

# Exceptions

`InvalidNextTokenException` or `ServiceFailureException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/ListSpeechSynthesisTasks)
"""
@inline list_speech_synthesis_tasks(aws::AWSConfig=default_aws_config(); args...) = list_speech_synthesis_tasks(aws, args)

@inline list_speech_synthesis_tasks(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "GET", "/v1/synthesisTasks", args)

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

"""
    using AWSSDK.Polly.put_lexicon
    put_lexicon([::AWSConfig], arguments::Dict)
    put_lexicon([::AWSConfig]; LexiconName=, Content=)

    using AWSCore.Services.polly
    polly([::AWSConfig], "PUT", "/v1/lexicons/{LexiconName}", arguments::Dict)
    polly([::AWSConfig], "PUT", "/v1/lexicons/{LexiconName}", LexiconName=, Content=)

# PutLexicon Operation

Stores a pronunciation lexicon in an AWS Region. If a lexicon with the same name already exists in the region, it is overwritten by the new lexicon. Lexicon operations have eventual consistency, therefore, it might take some time before the lexicon is available to the SynthesizeSpeech operation.

For more information, see [Managing Lexicons](http://docs.aws.amazon.com/polly/latest/dg/managing-lexicons.html).

# Arguments

## `LexiconName = ::String` -- *Required*
Name of the lexicon. The name must follow the regular express format [0-9A-Za-z]{1,20}. That is, the name is a case-sensitive alphanumeric string up to 20 characters long.

## `Content = ::String` -- *Required*
Content of the PLS lexicon as string data.

# Returns

`PutLexiconOutput`

# Exceptions

`InvalidLexiconException`, `UnsupportedPlsAlphabetException`, `UnsupportedPlsLanguageException`, `LexiconSizeExceededException`, `MaxLexemeLengthExceededException`, `MaxLexiconsNumberExceededException` or `ServiceFailureException`.

# Example: To save a lexicon

Stores a pronunciation lexicon in an AWS Region.

Input:
```
[
    "Content" => "file://example.pls",
    "Name" => "W3C"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/PutLexicon)
"""
@inline put_lexicon(aws::AWSConfig=default_aws_config(); args...) = put_lexicon(aws, args)

@inline put_lexicon(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "PUT", "/v1/lexicons/{LexiconName}", args)

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

"""
    using AWSSDK.Polly.start_speech_synthesis_task
    start_speech_synthesis_task([::AWSConfig], arguments::Dict)
    start_speech_synthesis_task([::AWSConfig]; OutputFormat=, OutputS3BucketName=, Text=, VoiceId=, <keyword arguments>)

    using AWSCore.Services.polly
    polly([::AWSConfig], "POST", "/v1/synthesisTasks", arguments::Dict)
    polly([::AWSConfig], "POST", "/v1/synthesisTasks", OutputFormat=, OutputS3BucketName=, Text=, VoiceId=, <keyword arguments>)

# StartSpeechSynthesisTask Operation

Allows the creation of an asynchronous synthesis task, by starting a new `SpeechSynthesisTask`. This operation requires all the standard information needed for speech synthesis, plus the name of an Amazon S3 bucket for the service to store the output of the synthesis task and two optional parameters (OutputS3KeyPrefix and SnsTopicArn). Once the synthesis task is created, this operation will return a SpeechSynthesisTask object, which will include an identifier of this task as well as the current status.

# Arguments

## `LexiconNames = [::String, ...]`
List of one or more pronunciation lexicon names you want the service to apply during synthesis. Lexicons are applied only if the language of the lexicon is the same as the language of the voice.

## `OutputFormat = "json", "mp3", "ogg_vorbis" or "pcm"` -- *Required*
The format in which the returned output will be encoded. For audio stream, this will be mp3, ogg_vorbis, or pcm. For speech marks, this will be json.

## `OutputS3BucketName = ::String` -- *Required*
Amazon S3 bucket name to which the output file will be saved.

## `OutputS3KeyPrefix = ::String`
The Amazon S3 Key prefix for the output speech file.

## `SampleRate = ::String`
The audio frequency specified in Hz.

The valid values for mp3 and ogg_vorbis are "8000", "16000", and "22050". The default value is "22050".

Valid values for pcm are "8000" and "16000" The default value is "16000".

## `SnsTopicArn = ::String`
ARN for the SNS topic optionally used for providing status notification for a speech synthesis task.

## `SpeechMarkTypes = ["sentence", "ssml", "viseme" or "word", ...]`
The type of speech marks returned for the input text.

## `Text = ::String` -- *Required*
The input text to synthesize. If you specify ssml as the TextType, follow the SSML format for the input text.

## `TextType = "ssml" or "text"`
Specifies whether the input text is plain text or SSML. The default value is plain text.

## `VoiceId = "Geraint", "Gwyneth", "Mads", "Naja", "Hans", "Marlene", "Nicole", "Russell", "Amy", "Brian", "Emma", "Raveena", "Ivy", "Joanna", "Joey", "Justin", "Kendra", "Kimberly", "Matthew", "Salli", "Conchita", "Enrique", "Miguel", "Penelope", "Chantal", "Celine", "Lea", "Mathieu", "Dora", "Karl", "Carla", "Giorgio", "Mizuki", "Liv", "Lotte", "Ruben", "Ewa", "Jacek", "Jan", "Maja", "Ricardo", "Vitoria", "Cristiano", "Ines", "Carmen", "Maxim", "Tatyana", "Astrid", "Filiz", "Vicki", "Takumi", "Seoyeon" or "Aditi"` -- *Required*
Voice ID to use for the synthesis.

# Returns

`StartSpeechSynthesisTaskOutput`

# Exceptions

`TextLengthExceededException`, `InvalidS3BucketException`, `InvalidS3KeyException`, `InvalidSampleRateException`, `InvalidSnsTopicArnException`, `InvalidSsmlException`, `LexiconNotFoundException`, `ServiceFailureException`, `MarksNotSupportedForFormatException` or `SsmlMarksNotSupportedForTextTypeException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/StartSpeechSynthesisTask)
"""
@inline start_speech_synthesis_task(aws::AWSConfig=default_aws_config(); args...) = start_speech_synthesis_task(aws, args)

@inline start_speech_synthesis_task(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "POST", "/v1/synthesisTasks", args)

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

"""
    using AWSSDK.Polly.synthesize_speech
    synthesize_speech([::AWSConfig], arguments::Dict)
    synthesize_speech([::AWSConfig]; OutputFormat=, Text=, VoiceId=, <keyword arguments>)

    using AWSCore.Services.polly
    polly([::AWSConfig], "POST", "/v1/speech", arguments::Dict)
    polly([::AWSConfig], "POST", "/v1/speech", OutputFormat=, Text=, VoiceId=, <keyword arguments>)

# SynthesizeSpeech Operation

Synthesizes UTF-8 input, plain text or SSML, to a stream of bytes. SSML input must be valid, well-formed SSML. Some alphabets might not be available with all the voices (for example, Cyrillic might not be read at all by English voices) unless phoneme mapping is used. For more information, see [How it Works](http://docs.aws.amazon.com/polly/latest/dg/how-text-to-speech-works.html).

# Arguments

## `LexiconNames = [::String, ...]`
List of one or more pronunciation lexicon names you want the service to apply during synthesis. Lexicons are applied only if the language of the lexicon is the same as the language of the voice. For information about storing lexicons, see [PutLexicon](http://docs.aws.amazon.com/polly/latest/dg/API_PutLexicon.html).

## `OutputFormat = "json", "mp3", "ogg_vorbis" or "pcm"` -- *Required*
The format in which the returned output will be encoded. For audio stream, this will be mp3, ogg_vorbis, or pcm. For speech marks, this will be json.

## `SampleRate = ::String`
The audio frequency specified in Hz.

The valid values for `mp3` and `ogg_vorbis` are "8000", "16000", and "22050". The default value is "22050".

Valid values for `pcm` are "8000" and "16000" The default value is "16000".

## `SpeechMarkTypes = ["sentence", "ssml", "viseme" or "word", ...]`
The type of speech marks returned for the input text.

## `Text = ::String` -- *Required*
Input text to synthesize. If you specify `ssml` as the `TextType`, follow the SSML format for the input text.

## `TextType = "ssml" or "text"`
Specifies whether the input text is plain text or SSML. The default value is plain text. For more information, see [Using SSML](http://docs.aws.amazon.com/polly/latest/dg/ssml.html).

## `VoiceId = "Geraint", "Gwyneth", "Mads", "Naja", "Hans", "Marlene", "Nicole", "Russell", "Amy", "Brian", "Emma", "Raveena", "Ivy", "Joanna", "Joey", "Justin", "Kendra", "Kimberly", "Matthew", "Salli", "Conchita", "Enrique", "Miguel", "Penelope", "Chantal", "Celine", "Lea", "Mathieu", "Dora", "Karl", "Carla", "Giorgio", "Mizuki", "Liv", "Lotte", "Ruben", "Ewa", "Jacek", "Jan", "Maja", "Ricardo", "Vitoria", "Cristiano", "Ines", "Carmen", "Maxim", "Tatyana", "Astrid", "Filiz", "Vicki", "Takumi", "Seoyeon" or "Aditi"` -- *Required*
Voice ID to use for the synthesis. You can get a list of available voice IDs by calling the [DescribeVoices](http://docs.aws.amazon.com/polly/latest/dg/API_DescribeVoices.html) operation.

# Returns

`SynthesizeSpeechOutput`

# Exceptions

`TextLengthExceededException`, `InvalidSampleRateException`, `InvalidSsmlException`, `LexiconNotFoundException`, `ServiceFailureException`, `MarksNotSupportedForFormatException` or `SsmlMarksNotSupportedForTextTypeException`.

# Example: To synthesize speech

Synthesizes plain text or SSML into a file of human-like speech.

Input:
```
[
    "LexiconNames" => [
        "example"
    ],
    "OutputFormat" => "mp3",
    "SampleRate" => "8000",
    "Text" => "All Gaul is divided into three parts",
    "TextType" => "text",
    "VoiceId" => "Joanna"
]
```

Output:
```
Dict(
    "AudioStream" => "TEXT",
    "ContentType" => "audio/mpeg",
    "RequestCharacters" => 37
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/polly-2016-06-10/SynthesizeSpeech)
"""
@inline synthesize_speech(aws::AWSConfig=default_aws_config(); args...) = synthesize_speech(aws, args)

@inline synthesize_speech(aws::AWSConfig, args) = AWSCore.Services.polly(aws, "POST", "/v1/speech", args)

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

end # module Polly

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

ext: jl


text:
export SARTS
export SART
export RLDataSet

abstract type RLDataSet end

const SARTS = (:state, :action, :reward, :terminal, :next_state)
const SART = (:state, :action, :reward, :terminal)


text:
# LP programming using MATLAB-like construction
#
#   maximize 1000 x + 350 y
#
#       s.t. x >= 30, y >= 0
#            x - 1.5y >= 0  (i.e. -x + 1.5 y <= 0)
#            12 x + 8 y <= 1000
#            1000 x + 300 y <= 70000
#
#   solution: (59.0909, 36.3636)
#   objv: 71818.1818
#

using MathProgBase
using Gurobi

env = Gurobi.Env()

model = gurobi_model(env;
	name="lp_02",
	sense=:maximize,
	f = [1000., 350.],
	A = [-1. 1.5; 12. 8.; 1000. 300.],
	b = [0., 1000., 70000.],
	lb = [0., 30.])

println(model)

optimize(model)

println()
println("soln = $(get_solution(model))")
println("objv = $(get_objval(model))")


text:
###### customizable show ########

const DECIMALS = 2
const SHOWINT = false


text:
# Average fields along the field directions (i.e., Fw along the w-direction).  See the
# description of create_m̂() below.

# create_mean() assumes the space dimension and field dimension are the same.  In other
# words, when the space coordinate indices are (i,j,k), then the field has three vector
# components.  Therefore, for the input field array F[i,j,k,w], we assume w = 1:3.  This
# assumption is reasonable, because line-averaging the w-component of the field along the
# w-direction requires the interpolation points to be a half-grid point away from the
# original points along the w-direction, which means that the domain has a nonzero thickness
# along the w-direction.  Therefore, if the field has the w-component, the domain should
# have the w-dimension as well.
#
# There are cases where the averaging needs to be done in the direction orthogonal to the
# field direciton.  For example, in the TM equation, we have Ez at pixel corners, and we may
# want to know Ez at the pixel edges (at Hx- and Hy-points).  In this case the field
# dimension is 1 and different from the shape dimension = 2.  However, it doesn't make sense
# to have Ez at the Hx- and Hy-points simultaneously, because then the density of Ez points
# is doubled.  We usually need to know Ez at either the Hx-points or the Hy-points, but not
# both.  Then the needed operator can be easily generated using create_m̂().
export create_m̂, create_mean

# Assumption: we don't calculate averages for interpolated fields.  In other words, we
# calculate averages for the original fields defined on Yee's grid.  (There are exceptions;
# see below.)  This assumption simplifies implementation of the symmetry boundary condition
# in the averaging operators greatly.
#
# For example, suppose we use the primal grid planes as the E-field planes.  Then, for the
# symmetry boundary condition (which is the PEC boundary condition in this case), the fields
# defined on the domain boundaries are the tangential  E-fields and normal H-fields.  Both
# fields are zero on the symmetry boundaries.  Therefore, if we assume that we calculate
# averages only for those fields defined on Yee's grid, we do not have to construct
# separate averaging operators for the tangential E-field and normal H-field.
#
# See the comments in differential.jl as well.

# Note for SALT users
# The averaging scheme in SALT will be mostly used with the Bloch boundary condition, so the
# implementation of the impact of boundary conditions in the averaging scheme is actually
# not that important, because field averaging at the boundaries is trivial for the Bloch
# boundary condition.

# The forward and (-1) × backward averaging operators are not the transpose of each other
# when the symmetry boundary condition is used: one has element with the factor of 2
# multiplied while the other with the factor of zero.  Still, the symmetry boundary enforces
# a special structure on the subpixel smoothed material parameter tensor on the symmetry
# boundary.  This tensor of a special structure has zero entries where the non-transpose
# relationship between the forward and backward averaging matters, and eventually makes the
# effect of this non-transpose relationship does not propagate up to the symmetry of the
# final operator (i.e., the final operator is still symmetric even though the forward and
# (-1) × backward averaging operators are not the transpose of each other). See my notes
# entitled [Beginning of the part added on Aug/14/2018] in RN - Subpixel Smoothing.nb.

# Wrapper to create the arithmetic averaging operator by default
# This corresponds to the wrappers to create discrete versions of differential operators,
# but the nonzero enries in the matrix are 0.5 rather than 1.0.
#
# Note that when the positional arguments are SVector's, this function may call the concrete
# implementation directly instead of the next wrapper function taking ∆l and ∆l′⁻¹ even if
# the keyword arguments are not SVector's, because methods are dispatched based only on the
# types of positional arguments.  Therefore, we convert the keyword arguments into SVector's.
create_mean(isfwd::AbsVecBool,  # isfwd[w] = true|false: m̂_w is forward|backward averaging
            N::AbsVecInteger,  # size of grid
            isbloch::AbsVecBool=fill(true,length(isfwd)),  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(length(isfwd));  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
            order_cmpfirst::Bool=true) =  # true to use Cartesian-component-major ordering for more tightly banded matrix
    (K = length(isfwd); create_mean(isfwd, fill.(1.0,(N...,)), fill.(1.0,(N...,)), isbloch, e⁻ⁱᵏᴸ; order_cmpfirst))

# Wrapper to convert AbstractVector's to SVec's.
create_mean(isfwd::AbsVecBool,  # isfwd[w] = true|false: m̂_w is forward|backward averaging
            ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with; vectors of length N
            ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by; vectors of length N
            isbloch::AbsVecBool=fill(true,K),  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
            order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
            ) where {K} =
    create_mean(SBool{K}(isfwd), ∆l, ∆l′⁻¹, SBool{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); order_cmpfirst)

# Creates the field-averaging operator for all three Cartegian components.
#
# The function used to support setting nonzero matrix blocks at off-diagonal blocks, but the
# need for doing it is absorbed in MaxwellFDM.jl/src/param.jl/create_param_matrix() by
# putting the material parameters at off-diagonal blocks.  Therefore, the field-averaging
# operators always the diagonal blocks as nonzero blocks, as they perform averaging on the
# already generated input and output fields.  See RN - Subpixel Smoothing > [Update (May/13/2018)].
function create_mean(isfwd::SBool{K},  # isfwd[w] = true|false: m̂_w is forward|backward averaging
                     ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with; vectors of length N
                     ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by; vectors of length N
                     isbloch::SBool{K},  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
                     e⁻ⁱᵏᴸ::SNumber{K};  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
                     order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
                     ) where {K}
    T = promote_type(eltype.(∆l)..., eltype.(∆l′⁻¹)..., eltype(e⁻ⁱᵏᴸ))  # eltype(eltype(∆l)) can be Any if ∆l is inhomogeneous
    N = SVec(length.(∆l))
    M = prod(N)
    KM = K * M

    Itot = VecInt(undef, 2KM)
    Jtot = VecInt(undef, 2KM)
    Vtot = Vector{T}(undef, 2KM)

    indblk = 0  # index of matrix block
    for nw = 1:K  # direction of averaging
        I, J, V = create_minfo(nw, isfwd[nw], N, ∆l[nw], ∆l′⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw])  # averaging along nv-direction

        nv = nw  # component of output field generated by w-directional averaging operator m̂_w (row index of matrix block)
        istr, ioff = order_cmpfirst ? (K, nv-K) : (1, M*(nv-1))  # (row stride, row offset)

        nu = nw  # component of input field to feed to w-directional averaging operator m̂_w (column index of matrix block)
        jstr, joff = order_cmpfirst ? (K, nu-K) : (1, M*(nu-1))  # (column stride, column offset)

        @. I = istr * I + ioff
        @. J = jstr * J + joff

        # For some reason, using .= below is slower because it uses 1 allocatiotn.  On the
        # other hand, using = does not use allocation and therefore faster.
        indₛ, indₑ = indblk*2M + 1, (indblk+1)*2M  # each I, J, V is length-2M
        Itot[indₛ:indₑ] = I
        Jtot[indₛ:indₑ] = J
        Vtot[indₛ:indₑ] = V
        indblk += 1
    end

    return dropzeros!(sparse(Itot, Jtot, Vtot, KM, KM))  # KM×KM matrix with 2 entries per row (so 2KM entries for I, J, V)
end

## Field-averaging operators "m̂" (as used in de Moerloose and de Zutter)
#
# This creates the averaging operator for a single Cartesian component.  For the operator
# for all three Cartesian components, use create_mean.
#
# Construction of these operators are similar to that of difference operators.  However,
# unlike the difference operators that are primarilly used for curl and hence differentiates
# the fields along the direction normal to the fields, the field-averaging operators average
# fields along the field direction.  As a result, backward (rather than forward) averaging
# is for primal fields.
create_m̂(nw::Integer,  # for length(N) = 3, 1|2|3 for averaging along x|y|z
         isfwd::Bool,  # true|false for forward|backward averaging
         N::AbsVecInteger,  # size of grid
         isbloch::Bool=true,  # boundary condition in w-direction
         e⁻ⁱᵏᴸ::Number=1.0) =  # Bloch phase factor
    (K = length(N); ∆w = ones(N[nw]); ∆w′⁻¹ = ones(N[nw]); create_m̂(nw, isfwd, SInt{K}(N), ∆w, ∆w′⁻¹, isbloch, e⁻ⁱᵏᴸ))

create_m̂(nw::Integer,  # for length(N) = 3, 1|2|3 for averaging along x|y|z
         isfwd::Bool,  # true|false for forward|backward averaging
         N::AbsVecInteger,  # size of grid
         ∆w::AbsVecNumber,  # line segments to multiply with; vector of length N[nw]
         ∆w′⁻¹::AbsVecNumber,  # inverse of line segments to divide by; vector of length N[nw]
         isbloch::Bool=true,  # boundary condition in w-direction
         e⁻ⁱᵏᴸ::Number=1.0) =  # Bloch phase factor
    (K = length(N); M = prod(N); dropzeros!(sparse(create_minfo(nw, isfwd, SInt{K}(N), ∆w, ∆w′⁻¹, isbloch, e⁻ⁱᵏᴸ)..., M, M)))

function create_minfo(nw::Integer,  # for K = 3, 1|2|3 for averaging along x|y|z
                      isfwd::Bool,  # true|false for forward|backward averaging
                      N::SInt{K},  # size of grid
                      ∆w::AbsVecNumber,  # line segments to multiply with; vector of length N[nw]
                      ∆w′⁻¹::AbsVecNumber,  # inverse of line segments to divide by; vector of length N[nw]
                      isbloch::Bool,  # boundary condition in w-direction
                      e⁻ⁱᵏᴸ::Number  # Bloch phase factor
                      ) where {K}
    M = prod(N)
    Nw = N[nw]
    ŵ = SVec(ntuple(identity,Val(K))) .== nw  # [0,true,0] for w == y
    ns = isfwd ? 1.0 : -1.0  # number for sign

    # Below, when constructing I, J's, V's, note that a tuple of array subscripts (i,j,k)
    # indicates a row index of the final matrix.  In other words, the entries with the same
    # subscripts correspond to the same *output* field.  So, when setting I[i,j,k], J₀[i,j,k],
    # Jₛ[i,j,k], V₀[i,j,k], Vₛ[i,j,k], we must know that they are for the entries in the
    # same row index (or the same output field).  Specifically,
    #
    # - I is the row index itself, so it is the identity map: I[i,j,k] = LinearIndices(N)[i,j,k].
    #
    # - J₀[i,j,k] is the column index of the diagonal entry in the row subscripted (i,j,k).
    # Because the row and column indices are the same for the diagonal entries, J₀ = I and
    # not created separately.
    #
    # - Jₛ[i,j,k] is the column index of the off-diagonal entry in the row subscripted
    # (i,j,k).
    #
    # - V₀[i,j,k] is the value of the diagonal entry in the row subscripted (i,j,k).
    #
    # - Vₛ[i,j,k] is the value fo the off-diagonal entry in the row subscripted (i,j,k).

    I = reshape(collect(1:M), N.data)

    Jₛ = reshape(collect(1:M), N.data)
    shifts = -ns * ŵ  # shift vector for sign; [0,-1,0] for w == y and ns = +1
    Jₛ = circshift(Jₛ, shifts.data)

    # Align ∆w and ∆w′⁻¹ in the w-direction.
    vec1 =  @SVector ones(Int,K)
    sizew = @. !ŵ * vec1 + ŵ * N  # [1,Ny,1] for w == y
    ∆W = reshape(∆w, sizew.data)
    ∆W′⁻¹ = reshape(∆w′⁻¹, sizew.data)

    # Construct the values of the diagonal and off-diagonal nonzero entries of the matrix
    # for the periodic (≠ Bloch) boundary condition.  (These values will be updated for
    # specific boundary conditions later.)
    #
    # To figure out the entries of V₀ and Vₛ, it is most convenient to figure out the entire
    # 3D arrays V₀ and Vₛ, instead of entry by entry.  In doing so, note that the entries
    # with the same subscripts correspond to the same *output* field (= same row index).
    T = promote_type(eltype(∆w), eltype(∆w′⁻¹), eltype(e⁻ⁱᵏᴸ))
    V₀ = fill(T(0.5), N.data) .* ∆W .* ∆W′⁻¹  # values of diagonal entries

    Vₛ = fill(T(0.5), N.data) .* ∆W  # values of off-diagonal entries (division later)
    Vₛ = circshift(Vₛ, shifts.data)
    Vₛ .*= ∆W′⁻¹

    # Modify V₀ and Vₛ according to the boundary condition.  What we do is basically to take
    # the operator for the periodic boundary condition (constructed above) as a template and
    # then modify it for the actual boundary condition, either by applying the Bloch phase
    # factors for the Bloch boundary condition or by multiplying 0 or 2 to some diagonal and
    # off-diagonal entries for the symmetry boundary condition.  Because this modification
    # is due to the boundary condition, the modification needs to be made only at the array
    # entries with the first or last subscript of the Cartesion direction in which averaging
    # is taken.
    #
    # The outline of the procedure is as follows:
    #
    # 1. Detertmine the rows in which the modification should be made.  This is done by
    # figuring out the array subscripts (i,j,k) of the output average array entry to which
    # the input boundary fields contribute.
    #
    # 2. In each of these rows, determine the column in which the modification should be
    # made.  Because each row has only two nonzero entries (one diagoal and one off-diagonal),
    # this step is to determine whether the modification should be made at the diagonal
    # entry or the off-diagonal entry in the row currently concerned.  If the (i,j,k) of the
    # input fields contributing to the currently concerned output averages are the same as
    # the (i,j,k) of the output averages, the modification is made in V₀ (diagonal).  If
    # they are different, the modification is made in Vₛ (off-diagonal).
    #
    # 3. In the determined rows, modify either V₀ or Vₛ depending on the above step.  This
    # is to modify V₀[i,j,k] or Vₛ[i,j,k], where (i,j,k) is the array subscripts specifying
    # the determined row.
    #
    # For the final sparsity patterns of the operators, see my notes entitled [Beginning of
    # the part added on Aug/14/2018] in RN - Subpixel Smoothing.nb.
    #
    # Below, V[Base.setindex(axes(V), iw:iw, nw)...] mimics the implementation of slicedim
    # and means V[:,iw:iw,:] for w = y.  (axes(V) = (1:Nx,1:Ny,1:Nz) and setindex(t,v,i)
    # does t[i] = v, so Base.setindex(axes(V),iw:iw,nw) creates (1:Nx,iw:iw,1:Nz) for w = y.)
    # The use of iw:iw instead of iw is to support 1D arrays.  If V is 1D, then V[iw] is a
    # scalar and the dot equal on V[iw] fails.
    if isbloch
        # Application of the aformentioned procedure:
        #
        # 0. Goal
        # Our goal is to create the ghost fields on one boundary by multiplying the Bloch
        # phase factors to the nonghost fields on the opposite boundary.  Note that w is the
        # direction of averaging.
        #
        # 1. Determination of the indices of the rows to modify
        # For isfwd = true, we are taking forward averaging between primal grid planes.  The
        # ghost input fields are at the positive-end boundary, and they contritube to the
        # averages at the last dual grid points, so the rows to modify are subscripted (i,j,k)
        # with iw ≡ (i,j,k)[w] = Nw.
        #
        # For isfwd = false, we are taking backward averaging between dual grid planes.  The
        # ghost input fields are before the negative-end boundary, and they contribute to
        # the averages at the first primal grid points, so the rows to modify are
        # subscripted (i,j,k) with iw ≡ (i,j,k)[w] = 1.
        #
        # 2. Determination of the column in each row to modify
        # For isfwd = true, the w-subscript of the ghost input field is greater by 1 than
        # that of the output average.
        #
        # For isfwd = false, the w-subscript of the ghost input field is less by 1 than that
        # of the output average.
        #
        # The nonghost input fields are brought to these ghost fields (from the opposite
        # boundary), an they are different from the input fields subscripted the same as the
        # output averages.  Therefore, the matrix entries to modify by multiplying the Bloch
        # phase factors with are off-diagonal entries in Vₛ.
        #
        # 3. Modification of the values in V
        # For isfwd = true, the nonghost input fields have a smaller w-subscript than the
        # ghost input fields that they are brought into.  Therefore, e⁻ⁱᵏᴸ must be
        # multiplied to the nonghost input fields to create the ghost input fields.
        #
        # For isfwd = false, the nonghost input fields have a larger w-subscript than the
        # ghost input fields that they are brought into.  Therefore, e⁺ⁱᵏᴸ must be
        # multiplied to the nonghost input fields to create the ghost input fields.
        iw = isfwd ? Nw : 1
        Vₛ[Base.setindex(axes(Vₛ), iw:iw, nw)...] .*= e⁻ⁱᵏᴸ^ns
    else  # symmetry bounndary
        # A. isfwd = true (forward averaging)
        #
        # A-0. Goal
        # Forward averaging is taken between primal grid planes.  Because the domain
        # boundaries are always primal grid planes, some input fields (either tangential or
        # normal) are on the domain boundaries.  Those input fields are zero for the
        # symmetry boundary condition and do not contribute to the output averages.  Our
        # goal is to make sure that the forward averaging operator properly ignore those
        # fields that should be zero.
        #
        # A-1. Determination of the indices of the the rows to modify
        # The negative-end-boundary input fields contribute to the w-directional averages
        # evaluated at the first (nonghost) dual grid points in the w-direction.  The
        # positive-end-boundary (ghost) input fields contribute to the w-directional
        # averages evaluated at the last dual grid points in the w-direction.  Therefore,
        # the array subscripts (i,j,k) we need to consider are the ones with iw ≡ (i,j,k)[w] = 1, Nw.
        #
        # A-2. Determination of the column in each row to modify
        # The negative-end-boundary input fields have subscripts the same as the averages
        # they contribute to, so for iw = 1 we modify the diagonal entries stored in V₀.
        # The positive-end-boundary input fields have w-subscript greater by 1 than the
        # averages they contribute to, so for iw = Nw we modify the off-diagonal entries
        # stored in Vₛ.
        #
        # A-3. Modification of the values in V
        # The negative-end-boundary input fields must be zeroed, so for (i,j,k) with iw = 1
        # we set V₀[i,j,k] = 0.
        # The positive-end-boundary input fields must be also zeroed, so for (i,j,k) with iw = Nw
        # we set Vₛ[i,j,k] = 0.  Note that the positive-end-boundary input fields are ghost
        # fields that are brought from the negative-end-boundary (nonghost) input fields.
        #
        #
        # B. isfwd = false (backward averaging)
        #
        # B-0. Goal
        # Backward averaging is taken between dual grid planes.  The domain boundaries are
        # always primal grid planes.  Therefore, no input fields (either tangential or
        # normal) are on the domain boundaries, and thus no input fields needs to be zeroed.
        #
        # However, the field at dual grid points are symmetric (rather than anti-symmetric)
        # around the symmetry boundaries (so that their interpolated values at the symmetry
        # boundary are nonzero). Therefore, the average between the two is double the
        # contribution from one input field, leading to multiplication of a factor of 2 at
        # the symmetry boundaries.
        #
        # On the positive-end boundary, the average of the dual field is not evaluated,
        # so we don't have to do anything to realize the positivte-end symmetry boundary.
        # On the other hand, on the negative-end boundary the average of the dual field is
        # evaluated.  We have to make sure that the resulting average is double the
        # contribution from one input field there.
        #
        # B-1. Determination of the indices of the the rows to modify
        # The output averages at the negative-end boundary are subscripted (i,j,k) with
        # iw ≡ (i,j,k)[w] = 1.
        #
        # B-2. Determination of the column in each row to modify
        # The output averages at the negative-end boundary is taken between the first
        # (nonghost) dual grid point and the ghost point before that point.  We can realize
        # the symmetry boundary by doubling the effect of the input field at the first
        # dual grid point and zeroing the effect of the input field at the ghost point.
        #
        # The input field at the first dual grid point is subscripted the same as the output
        # average, so its effect is described in a diagonal entry in V₀.  The input field at
        # the ghost point is subscripted 1 less in the w-direction, so its effect is
        # described in an off-diagonal entry in Vₛ.
        #
        # B-3. Modification of the values in V
        # We double V₀[i,j,k] and zero Vₛ[i,j,k] for (i,j,k) with iw = 1.
        #
        # The above procedure A and B are written compactly below without separating the
        # isfwd = true and false cases.

        # Replace some diagonal entries (V₀) with 0 or 2.
        iw = 1
        val = isfwd ? 0 : 2
        V₀[Base.setindex(axes(V₀), iw:iw, nw)...] .*= val

        # Replace some off-diagonal entries (Vₛ) with 0.
        iw = isfwd ? Nw : 1
        Vₛ[Base.setindex(axes(Vₛ), iw:iw, nw)...] .= 0

        # Note that the forward and backward averaging operators are not the transpose of
        # each other because they have different diagonals.  This may seem to make the final
        # material parameter matrix (with the input and output field averaging operators
        # multiplied) nonsymmetric.  However, zeros in the material parameter tensor
        # enforced by the symmetry boundary condition nullifies the possibility of
        # nonsymmetry.  See my notes entitled [Beginning of the part added on Aug/14/2018]
        # in RN - Subpixel Smoothing.nb.  See also kottke_input_accurate in smoothing.jl.
    end

    i = reshape(I, M)
    jₛ = reshape(Jₛ, M)
    v₀ = reshape(V₀, M)
    vₛ = reshape(Vₛ, M)

    Is = [i; i]  # row indices of [diagonal; off-diagonal] entries
    Js = [i; jₛ]  # column indices of [diagonal; off-diagonal] entries
    Vs = [v₀; vₛ]  # matrix entries of [diagonal, off-diagonal] entris

    return Is, Js, Vs
end


text:
import FLOWUnsteady
uns = FLOWUnsteady
vlm = uns.vlm
import GeometricTools
gt = GeometricTools
using PyPlot
using LinearAlgebra: norm, dot, cross
using LaTeXStrings                      #Need for setting up plots

# velocities = [];       # Total Velocity over wing

function main(;disp_plot = true)

    #Define aspects of wing. Single, symmetric wing.

    #Forward Swept Wing
    span = 1.0;              #wing span
    aspectratio = 10.0;      #wing aspect ratio
    taperratio = 0.5;        #wing taper ratio
    wingtwist = 4.0;         #wing twist
    wingsweep = -5.0;        #wing sweep in degrees
    wingdihedral = 6.0;      #wing dihedral in degrees

    #Glider
    # span = 1.0;
    # aspectratio = 15.0;
    # taperratio = 0.9;
    # wingtwist = 1.0;
    # wingsweep = 1.0;
    # wingdihedral = 1.0;

    #This bulids the wing using the variables we gave it.
    mainwing = vlm.simpleWing(span,aspectratio,taperratio,wingtwist,wingsweep,wingdihedral);

    #Now we create a wing system
    system = vlm.WingSystem();       #creates an empty system

    vlm.addwing(system, "mainwing", mainwing);  #Adds wing to wystem

    #We will now save the system to a .vtk file. MAKE SURE YOU ARE IN A NON-IMPORANT LOCATION TO SAVE THIS.
    #I added an empty folder for this.

    #Prep file
    Vinf(X,t) = [1,0,0];        #non-dimensional function defining free stream velocity
    vlm.setVinf(system, Vinf)   #set freestream velocity for the system
    "Freestream Velocity set"
    #"$Vinf"

    run_name = "tutorial"           #define identifier at beginning of file names
    save_path = "./simplewing/"     #define directory where files will be saved

    run(`rm -rf $save_path`)        #clear out directory where files will be saved
    run(`mkdir $save_path`)         #re-create directory fresh

    vlm.save(system, run_name; path=save_path)  #save geometry in a .vtk file format

    #now we will show the wing using paraview

    # run(`paraview --data="$(save_path)/$(run_name)_mainwing_vlm.vtk"`)

    #We now have the wing set up and running in ParaView.

    #We will now add a rotor
    rotor_file = "apc10x7.csv"          # Rotor geometry
    data_path = uns.def_data_path       # Path to rotor database

    #this creates a rotor, and generates a heck-a-ton of output, which we suppress.
    R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]] #get the radius for later

    Vinf(X,t) = [1,0,0];
    # Vinf = [1,0,0];

    rotor = uns.generate_rotor(rotor_file; pitch=0.0,
                                            n=10, CW=true, ReD=1.5e6,
                                            verbose=true, xfoil=false,
                                            data_path=data_path,
                                            plot_disc=false);

    #Create a rotor object bsed on parameters above.
    rotors = vlm.Rotor[rotor];
    #We now set the location of the vehicle to be in line with the rotor.

    #First we set the vehicle at the origin.
    vehicleorigin = [0.0; 0.0; 0.0];
    vehicleaxis = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0];

    #Now we set the rotor origin in line with the vehicle
    rotororigin = [-0.1; 0.0; 0.0];
    RPM_Val = 1000

    #And now line the rest of the rotor up, part by part.
    for rotor in rotors
        vlm.setcoordsystem(rotor, rotororigin, vehicleaxis; user=true);
        vlm.setRPM(rotor, RPM_Val)
    end

    #Now we will store the rotor in a tuple
    rotor_systems = (rotors,);
    # vlm.setRPM(rotors, RPM);       #Maybe this.

    #finally, we add the rotor to the wing system. Again, part by part.
    for rotor in rotors
        vlm.addwing(system, run_name, rotor);
    end

    #Now we run it to see the visualization
    run(`rm -rf $save_path`)
    run(`mkdir $save_path`)
    #
    vlm.setVinf(system, Vinf)

    # vlm.setRPM(system, RPM);

    vlm.save(system, run_name; path=save_path)
    # run(`paraview --data="$(save_path)/tutorial_mainwing_vlm.vtk;tutorial_tutorial_Blade1_vlm.vtk;tutorial_tutorial_Blade2_vlm.vtk;tutorial_tutorial_Blade1_loft.vtk;tutorial_tutorial_Blade2_loft.vtk;"`)

    #Starting the new code block. There is still problems in the above code.

    ######################### Unrun code. #######################################
    #Should work, as it is basically line for line from online

    #create a vlm system to add our main wing to.
    vlm_system = vlm.WingSystem();
    vlm.addwing(vlm_system, "mainwing", mainwing);

    #Create a wake system. Add VLM system and rotor to this wake system.
    wake_system = vlm.WingSystem();
    # vlm.addwing(wake_system, "SolveVlm", vlm_system)

    # for rotor in rotors
    #     vlm.addwing(wake_system, run_name, rotor);
    # end

    #We need to define tilting objects (as they are part of the code) but ours will be empty

    tilting_systems = ();   #Empty tuple

    #The vehicle, rotor, and freestream are all set up now.

    #We will now define a maneuver.

    #Start with Vehicle velocity function
    Vvehicle(t) = [-1.0,0.0,0.0];

    #Define Vehicle angle
    anglevehicle(t) = zeros(3); #Flat angle

    #define angle of tilting system. We have no tilting system, so make it empty
    angle = ();

    #Now, lastly, we define the rotation rate function for the rotors
    RPM_fun(t) = 10.0;   #1.0 is constant rotation
    RPM = (RPM_fun,);

    #using each individual element, we can create a maneuver now.
    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle);

    #We can plot the maneuver now, which is good.
    uns.plot_maneuver(maneuver);

    #This is cool, but lets put it in paraview.
    vlm.setVinf(system, Vinf);
    vehicle = uns.QVLMVehicle( system;
                              tilting_systems = tilting_systems,
                              rotor_systems = rotor_systems,
                              vlm_system = vlm_system,
                              wake_system = wake_system,
                            );
    #Now we need some reference information: referenc velocity, total time, and initial conditions
    Forces_wing = [];      # Total force over wing.        Done, line 70
    forces_per_span = [];  # Force per unit span on wing.  Done, line 82
    cL_wing = [];          # Coefficient of lift of wing.  Done, line 87
    L_wing = [];           # Lift over wing.               Done, line 76
    cD_wing = [];          # Coefficint of drag over wing. Done Line 107
    D_wing = [];           # Effective drage over wing.    Done Line 106
    dt = 0;
    Time = 0;
    nsteps = 300;

    Flag = false;           #testing for the forces output

    Vref = 10.0         #define a reference velocity for the vehicle
    ttot = 1.0          #define a total simulation time, in seconds
    nsteps = 300        #define the number of steps the simulation will take
    prev_wing = deepcopy(mainwing)
    rhoinf = 1.23
    b = span
    qinf = 0.5*Vref^2 * rhoinf
    ar = aspectratio
    wake_coupled = false
    function monitor(sim, PFIELD, T, DT; figname="monitor_plot", nsteps_plot=1, disp_plot = disp_plot, figsize_factor = 1)
        #We need to write the monitor function, as seen in the examples, such the heavingwing example.

        #WRITE YOUR OWN MONITOR FUNCTION, as it is specific to your simulation.

        #To get the forces/moments/lift/lift coefficitnets, either declare variables outside of the function that
        #   you add to with the push!() command each iteration, or print all the numbers to a text file which you
        #   can use later. This is how you get these values for the wing.

        #Example of monitor function starts on line 199 of heavingwing.jl

        #There is a lot going on in here, especially a lot of plotting. Do you need this?
        #Maybe deleate anything you don't need to find the Forces and coefficients.

        aux = PFIELD.nt/nsteps
            clr = (1-aux, 0, aux)

            if PFIELD.nt==0 && disp_plot
                figure(figname, figsize=[7*2, 5*2]*figsize_factor)
                subplot(221)
                xlim([0,1])
                xlabel(L"$\frac{2y}{b}$")
                ylabel(L"$\frac{Cl}{CL}$")
                title("Spanwise lift distribution")

                subplot(222)
                xlim([0,1])
                xlabel(L"$\frac{2y}{b}$")
                ylabel(L"$\frac{Cd}{CD}$")
                title("Spanwise drag distribution")

                subplot(223)
                xlabel("Simulation time (s)")
                ylabel(L"Lift Coefficient $C_L$")
                ylim([0,1]);

                subplot(224)
                xlabel("Simulation time (s)")
                ylabel(L"Drag Coefficient $C_D$")
                ylim([0,0.1]);

                figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
                subplot(121)
                xlabel(L"$\frac{2y}{b}$")
                ylabel(L"Circulation $\Gamma$")
                subplot(122)
                xlabel(L"$\frac{2y}{b}$")
                ylabel(L"Effective velocity $V_\infty$")

                figure(figname*"_3", figsize=[7*2, 5*1]*figsize_factor)
                subplot(121)
                xlabel("Simulation time")
                ylabel("Velocity")
                subplot(122)
                xlabel("Simulation time")
                ylabel("Angular velocity")

                figure(figname*"_4", figsize=[7*2, 5*1]*figsize_factor)
                subplot(121)
                xlabel("Simulation time")
                ylabel("Total Force")
                subplot(122)
                xlabel("Simulation time")
                ylabel("Force per unit span")

            end

            if PFIELD.nt!=0 && PFIELD.nt%nsteps_plot==0 && disp_plot
                figure(figname)

                # Force at each VLM element
                Ftot = uns.calc_aerodynamicforce(mainwing, prev_wing, PFIELD, Vinf, DT,
                                                                rhoinf; t=PFIELD.t)
                # Add Forces to forces vector
                resultant_Force = sqrt(Ftot[1]^2 + Ftot[2]^2 + Ftot[3]^2);  #Resultnat vector
                push!(Forces_wing, restultant_Force);
                L, D, S = uns.decompose(Ftot, [0,0,1], [-1,0,0]);
                vlm._addsolution(mainwing, "L", L)
                vlm._addsolution(mainwing, "D", D)
                vlm._addsolution(mainwing, "S", S)

                # #Add Lift to lift vector
                push!(L_wing, L);

                # Force per unit span at each VLM element
                ftot = uns.calc_aerodynamicforce(mainwing, prev_wing, PFIELD, Vinf, DT,
                                            rhoinf; t=PFIELD.t, per_unit_span=true);
                l, d, s = uns.decompose(ftot, [0,0,1], [-1,0,0]);

                # Add Force/span to f vector
                resultant_force = sqrt(ftot[1]^2 + ftot[2]^2 + ftot[3]^2);
                push!(forces_per_span, resultant_force);
                # Lift of the wing
                Lwing = norm(sum(L))
                CLwing = Lwing/(qinf*b^2/ar)
                ClCL = norm.(l) / (Lwing/b)

                # Add cl to vector
                push!(cL_wing, CLwing);

                # Drag of the wing
                Dwing = norm(sum(D))
                CDwing = Dwing/(qinf*b^2/ar)
                CdCD = [sign(dot(this_d, [1,0,0])) for this_d in d].*norm.(d) / (Dwing/b) # Preserves the sign of drag

                # Add Dwing and CDwing to vectors
                push!(D_wing, Dwing)
                push!(cD_wing, CDwing);

                vlm._addsolution(mainwing, "Cl/CL", ClCL)
                vlm._addsolution(mainwing, "Cd/CD", CdCD)

                y2b = 2*mainwing._ym/b

                subplot(221)
                # plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
                plot(y2b, ClCL, "-", label="FLOWVLM", alpha=0.5, color=clr)

                subplot(222)
                # plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
                plot(y2b, CdCD, "-", label="FLOWVLM", alpha=0.5, color=clr)

                subplot(223)
                # plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
                plot([T], [CLwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

                subplot(224)
                # plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
                plot([T], [CDwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

                figure(figname*"_2")
                subplot(121)
                plot(y2b, mainwing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
                if wake_coupled && PFIELD.nt!=0
                    subplot(122)
                    plot(y2b, norm.(mainwing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
                    if VehicleType==uns.VLMVehicle
                        plot(y2b, norm.(mainwing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr)
                    end
                    plot(y2b, [norm(Vinf(vlm.getControlPoint(mainwing, i), T)) for i in 1:vlm.get_m(mainwing)],
                                                                "-k", label="FLOWVLM", alpha=0.5)
                end

                figure(figname*"_3")
                subplot(121)
                plot([sim.t], [sim.vehicle.V[1]], ".r", label=L"V_x", alpha=0.5)
                plot([sim.t], [sim.vehicle.V[2]], ".g", label=L"V_y", alpha=0.5)
                plot([sim.t], [sim.vehicle.V[3]], ".b", label=L"V_z", alpha=0.5)
                if PFIELD.nt==1; legend(loc="best", frameon=false); end;
                subplot(122)
                plot([sim.t], [sim.vehicle.W[1]], ".r", label=L"\Omega_x", alpha=0.5)
                plot([sim.t], [sim.vehicle.W[2]], ".g", label=L"\Omega_y", alpha=0.5)
                plot([sim.t], [sim.vehicle.W[3]], ".b", label=L"\Omega_z", alpha=0.5)
                if PFIELD.nt==1; legend(loc="best", frameon=false); end;

                #ADDED CODE HERE.
                #The forces aren't coming through, and I don't know why.

                figure(figname*"_4")
                subplot(121)
                plot([sim.t], resltant_Force, label="Total forces on wing", alpha = 0.5)
                subplot(122)
                plot([sim.t], ftot[1],label = "Force per unit span", alpha = 0.5)

            end

            prev_wing = deepcopy(mainwing)

            return false

    end
    #initial conditions
    tinit = 0.0                                  #initial time
    Vinit = Vref*maneuver.Vvehicle(tinit/ttot)   #initial linear velocity
    Winit = zeros(3)                             #initial angular velocity
    RPMreference = RPM_Val

    #Everything is now defined! Let's create the simulation

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMreference, ttot; Vinit=Vinit, Winit=Winit, t=tinit);

    #Visulaize the simulation in paraview
    # files = uns.visualize_kinematics(   simulation, nsteps, save_path;
    #                                     run_name=run_name,
    #                                     prompt=false,
    #                                     paraview=false
    #                                 )

    # run(`paraview --data="$save_path/$files"`)

    #The simulation also has some outputs we can recieve. The following code does this.

    # nullfunc(args...) = false
    pfield = uns.run_simulation(simulation, nsteps;
                                        surf_sigma=R/10,
                                        Vinf=Vinf,
                                        save_path=save_path,
                                        run_name=run_name,
                                        prompt=false,
                                        verbose=true,
                                        extra_runtime_function=monitor     #THIS LINE IS WHERE YOU CALL THE FORCES
                                )
                                #See line 325 of heavingwing.jl

    #run(`paraview --data="$save_path/$(files);tutorial_pfield...vtk"`)

    #THIS WORKS

    #Now, we go into the code and pull the lift/drag/and moments from the wing.
    #Plot those over time.

    #Look into the examples on FLOWUnsteady for the function call that returns ALL of those values.

    #We need to write the monitor function, as seen in the examples, such the heavingwing example.

    #WRITE YOUR OWN MONITOR FUNCTION, as it is specific to your simulation.

    #To get the forces/moments/lift/lift coefficitnets, either declare variables outside of the function that
    #   you add to with the push!() command each iteration, or print all the numbers to a text file which you
    #   can use later. This is how you get these values for the wing.

    #Example of monitor function starts on line 199 of heavingwing.jl

    dt = ttot/nsteps;
    time_vec = 0:dt:ttot;

    return cD_wing,cL_wing,D_wing,L_wing,time_vec,Forces_wing, forces_per_span

end

cD_wing,cL_wing,D_wing,L_wing,time_vec, Forces_wing, forces_per_span = main()

#
println(time_vec)
println(Forces_wing)    #This is not being added to. Is it a problem with it as a global variable?
                        #Can Julia push into this array in the function, then print it here?

FirstForce = zeros();            #Empty Tuple

Forces_wing = insert!(Forces_wing,1,FirstForce);     #add azero to the start

println(Forces_wing(1));
#FROM WHAT I UNDERSTNAD, THESE PLOTS ARE DONE, WITHOUT LABELS
#Plots of Forces
#Forces_wing is a very convoluted vector. It is a tuple.
# figure("Forces")
# subplot(211)
# plot(time_vec, Forces_wing[:,1])
# xlabel("Time(s)");
# ylabel("Force")
# subplot(212)
# plot(time_vec, forces_per_span[:,1]);
# xlabel("Time(s)")
# ylabel("Force")
#
# #Plots of lift
# figure("Lift")
# subplot(211)
# plot(time_vec, cl_wing);
# xlabel("time (s)");
# ylabel("cL")
# subplot(212)
# plot(time_vec, L_wing, label = cl);
# xlabel("time(s)")
# ylabel("L")

max_stars_repo_name: brigoste/FlowVLMForces

max_issues_repo_name: brigoste/FlowVLMForces

max_forks_repo_name: brigoste/FlowVLMForces


text:
function testing_frame()
    frame = Frame()
    resize!(frame, 4)

    topology = Topology()
    add_atom!(topology, Atom("H"))
    add_atom!(topology, Atom("O"))
    add_atom!(topology, Atom("O"))
    add_atom!(topology, Atom("H"))

    add_bond!(topology, 0, 1)
    add_bond!(topology, 1, 2)
    add_bond!(topology, 2, 3)

    set_topology!(frame, topology)
    return frame
end

@testset "Selection" begin
    frame = testing_frame()
    selection = Selection("name O")
    @test size(selection) == 1

    matches = evaluate(selection, frame)
    @test size(matches, 1) == 2
    @test matches[1] == 1
    @test matches[2] == 2

    selection = Selection("bonds: all")
    @test size(selection) == 2

    matches = evaluate(selection, frame)
    @test size(matches, 1) == 3
    @test ((0, 1) in matches) == true
    @test ((1, 2) in matches) == true
    @test ((2, 3) in matches) == true

    selection = Selection("angles: all")
    @test size(selection) == 3

    matches = evaluate(selection, frame)
    @test size(matches, 1) == 2
    @test ((0, 1, 2) in matches) == true
    @test ((1, 2, 3) in matches) == true

    selection = Selection("dihedrals: all")
    @test size(selection) == 4

    matches = evaluate(selection, frame)
    @test size(matches, 1) == 1
    @test ((0, 1, 2, 3) in matches) == true

    @test selection_string(selection) == "dihedrals: all"

    copy = deepcopy(selection)
    @test selection_string(copy) == "dihedrals: all"
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "GAP"
upstream_version = v"4.12.0-dev"
version = v"400.1191.001"

julia_versions = [v"1.6.0", v"1.7.0", v"1.8.0"]

# Collection of sources required to complete build
sources = [
    # snapshot of GAP master branch leading up to GAP 4.12:
    GitSource("https://github.com/gap-system/gap.git", "401c797476b787e748a3890be4ce95ae4e5d52ae"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(upstream_version)/gap-$(upstream_version)-core.tar.gz",
#                  "2b6e2ed90fcae4deb347284136427105361123ac96d30d699db7e97d094685ce"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done

# must run autogen.sh if compiling from git snapshot and/or if configure was patched;
# it doesn't hurt otherwise, too, so just always do it
./autogen.sh

# compile a native version of GAP to generate c_oper1.c and c_type1.c
mkdir native-build
cd native-build
rm ${host_libdir}/*.la  # delete *.la, they hardcode libdir='/workspace/destdir/lib'
../configure --build=${MACHTYPE} --host=${MACHTYPE} \
    --with-gmp=${host_prefix} \
    --without-readline \
    --with-zlib=${host_prefix} \
    CC=${CC_BUILD} CXX=${CXX_BUILD}
make -j${nproc}
cp build/c_*.c ../src/
cd ..

# remove the native build, it has done its job
rm -rf native-build

# compile GAP
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gmp=${prefix} \
    --with-readline=${prefix} \
    --with-zlib=${prefix} \
    --with-gc=julia \
    --with-julia
make -j${nproc}

# install GAP binaries
make install-bin install-headers install-libgap

# also install config.h
cp build/config.h ${prefix}/include/gap

# the license
install_license LICENSE

# get rid of *.la files, they just cause trouble
rm ${prefix}/lib/*.la

# get rid of the wrapper shell script, which is useless for us
mv ${prefix}/bin/gap.real ${prefix}/bin/gap

# install gac and sysinfo.gap
mkdir -p ${prefix}/share/gap/
cp gac sysinfo.gap ${prefix}/share/gap/

# We deliberately do NOT install the GAP library, documentation, etc. because
# they are identical across all platforms; instead, we use another platform
# independent artifact to ship them to the user.
"""

include("../../L/libjulia/common.jl")
platforms = vcat(libjulia_platforms.(julia_versions)...)
filter!(!Sys.iswindows, platforms)

# we only care about 64bit builds
filter!(p -> nbits(p) == 64, platforms)

# Windows is not supported
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("gap", :gap),
    LibraryProduct("libgap", :libgap),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # for the "native" build that generates c_oper1.c and c_type1.c
    HostBuildDependency("GMP_jll"),
    HostBuildDependency("Zlib_jll"),

    Dependency("GMP_jll"),
    Dependency("Readline_jll"),
    Dependency("Zlib_jll"),
    BuildDependency("libjulia_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7", julia_compat="1.6", init_block="""

    sym = dlsym(libgap_handle, :GAP_InitJuliaMemoryInterface)
    ccall(sym, Nothing, (Any, Ptr{Nothing}), @__MODULE__, C_NULL)
""")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
module MyModule

export f, g, h

f(x) = x^3
g(x) = f(f(x))
h(x) = 1/f(x)

end


text:
module TestCompute

using LibCloud.Compute
using Base.Test
using PyCall

function ec2(key, secret)
    nd = NodeDriver(ComputeProvider.EC2, key, secret)
    println("Regions:")
    for reg in list_regions(nd)
        println("    - ", reg)
    end

    println("Sizes:")
    for sz in list_sizes(nd)
        println("    - ", sz)
    end

    println("Locations:")
    for loc in list_locations(nd)
        println("    - ", loc)
    end

    println("Volumes:")
    for vol in list_volumes(nd)
        println(vol)
        for snapshot in list_volume_snapshots(nd, vol)
            println("    - ", snapshot)
        end
    end
end

function dummy()
    driver = NodeDriver(ComputeProvider.DUMMY, 0)

    driver_features = features(driver)
    @test "create_node" in keys(driver_features)
    @test isempty(driver_features["create_node"])

    println("Creating a node...")
    node = Node(create_node(driver))
    @test !isempty(get_uuid(node))
    @test state(node) == NodeState.RUNNING

    println("Listing nodes...")
    nodes = list_nodes(driver)
    @test !isempty(nodes)
    for n in nodes
        @test isa(n, Node)
    end

    println("Listing sizes...")
    sizes = list_sizes(driver)
    @test !isempty(sizes)
    for s in sizes
        @test isa(s, NodeSize)
    end

    println("Listing images...")
    images = list_images(driver)
    @test !isempty(images)
    for i in images
        @test isa(i, NodeImage)
    end

    println("Listing locations...")
    locs = list_locations(driver)
    @test !isempty(locs)
    for l in locs
        @test isa(l, NodeLocation)
    end

    println("Rebooting node...")
    @test reboot(node) == true
    @test state(node) == NodeState.REBOOTING
    println("Destroying node...")
    @test destroy(node) == true
    @test state(node) == NodeState.TERMINATED

    println("Creating key...")
    k = NodeAuthSSHKey("my pub key")
    @test k.o[:pubkey] == "my pub key"
end

function gce(service_account_id, key, project)
    driver = NodeDriver(ComputeProvider.GCE, service_account_id, key, project=project)

    driver_features = features(driver)
    @test "create_node" in keys(driver_features)
    @test isempty(driver_features["create_node"])

    println("Listing nodes...")
    nodes = list_nodes(driver)
    @test !isempty(nodes)
    for n in nodes
        @test isa(n, Node)
    end

    println("Listing sizes...")
    sizes = list_sizes(driver)
    @test !isempty(sizes)
    for s in sizes
        @test isa(s, NodeSize)
    end

    println("Listing images...")
    images = list_images(driver)
    @test !isempty(images)
    for i in images
        @test isa(i, NodeImage)
    end

    println("Listing locations...")
    locs = list_locations(driver)
    @test !isempty(locs)
    for l in locs
        @test isa(l, NodeLocation)
    end
end

end # module TestCompute


text:
function get_robust_path(start_node::Array, end_node::Array, c::Array, d::Array, Gamma, origin, destination)
	@assert length(start_node)==length(end_node)
	@assert length(start_node)==length(c)
	@assert length(start_node)==length(d)

	no_arc = length(start_node)

	dset = [sort(d, rev=true); 0]
	best_obj = Inf
	best_path = []
	best_x = []

	# Constructing Lset as in Lee and Kwon (2014) http://dx.doi.org/10.1007/s10288-014-0270-7
	k_max = round(Int, ceil((no_arc - Gamma)/2))

	Lset = Array{Int}(k_max+1)
	for k=1:k_max
		Lset[k] = Gamma + 2*k-1
	end
	Lset[k_max+1]=no_arc + 1

	# This for-loop examines the smallest theta first.
	for i=length(Lset):-1:1
		theta = dset[Lset[i]]
		link_cost = c + max(d - theta, 0)

		# println("$i : $theta : $(Gamma*theta) : $best_obj")

		if Gamma*theta > best_obj
			# If this condition is true, next theta values will always have a bigger objective value.
			# We can stop examining theta values.
			break
		end

		(pp, xx, ll) = get_shortest_path(start_node, end_node, link_cost, origin, destination)
		# current_obj = Gamma * theta + dot(link_cost, xx)
		current_obj = Gamma * theta + ll

		if current_obj < best_obj
			best_obj = current_obj
			best_path = pp
			best_x = xx
		end

	end

	return best_path, best_x, best_obj
end


text:
###################################################
## routing/virtualtimes.jl
## Create virtual routing times and paths in a network
###################################################

"""
    `roadTypePaths`
    - simple version of RoutingPaths
    - returns a timing object with a constant speed for each road type
"""
function roadTypeRouting(n::Network)
    # Timings are 60% of maximum speed
    times = maxSpeedTimes(n) / 0.6
    return RoutingPaths(n, times)
end

"""
    `randomTimeRouting`
    - returns a timing object with every edge having a randomly selected speed between 0 and 130 km/h
"""
function randomTimeRouting(n::Network)
    times = randomTimes(n)
    return RoutingPaths(n, times)
end

"""
    `fixedSpeedTimes` : given a speeds for each road-type, returns times for the whole network
    speeds are kph
"""
function fixedSpeedTimes(n::Network, speeds::Vector{Float64})
    g = n.graph
    times = spzeros(nv(g),nv(g))
    for ((o,d),r) in n.roads
        times[o,d] = 3.6*r.distance/speeds[r.roadType]
    end
    return times
end

"""
    `maxSpeedTimes`
    - returns road times corresponding to maximum allowed speed (sparse array)
    - `maxspeed`: km/h maximum speeds for each road-type
"""
maxSpeedTimes(n::Network; maxspeed::Vector{Float64} = [130.,110.,90.,50.,50.,20.,0.,0.]) =
fixedSpeedTimes(n, maxspeed)

"""
    `uniformTimes`
    - returns road times corresponding to constant speed (sparse array)
    - `speed` in km/h
"""
function uniformTimes(n::Network, speed::Real=90.)
    g = n.graph
    times = spzeros(nv(g),nv(g))
    for ((o,d),r) in n.roads
        times[o,d] = 3.6*r.distance/speed
    end
    return times
end

"""
    `randomTimes`
    - returns road times corresponding to random speed (sparse array) between 0 and 130 km/h
"""
function randomTimes(n::Network)
    srand()
    g = n.graph
    times = spzeros(nv(g), nv(g))
    for ((o,d), r) in n.roads
        times[o,d] = 3.6 * r.distance / (rand() * 130)
    end
    return times
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


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

mutable struct TrustedIdProviderProperties <: SwaggerModel
    idProvider::Any # spec type: Union{ Nothing, String } # spec name: idProvider

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

const _property_map_TrustedIdProviderProperties = Dict{Symbol,Symbol}(Symbol("idProvider")=>Symbol("idProvider"))
const _property_types_TrustedIdProviderProperties = Dict{Symbol,String}(Symbol("idProvider")=>"String")
Base.propertynames(::Type{ TrustedIdProviderProperties }) = collect(keys(_property_map_TrustedIdProviderProperties))
Swagger.property_type(::Type{ TrustedIdProviderProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_TrustedIdProviderProperties[name]))}
Swagger.field_name(::Type{ TrustedIdProviderProperties }, property_name::Symbol) =  _property_map_TrustedIdProviderProperties[property_name]

function check_required(o::TrustedIdProviderProperties)
    (getproperty(o, Symbol("idProvider")) === nothing) && (return false)
    true
end

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


text:
module AKNS
using Base, ApproxFun, ApproxFunBase, SpecialFunctions, ApproxFunFourier, Reexport, LinearAlgebra,
    ApproxFunRational, BandedMatrices, SparseArrays
#, MacroTools#, Reexport, AbstractFFTs, FFTW, InfiniteArrays, FillArrays, FastTransforms, IntervalSets,

export AKNSscattering, ODEFT, AKNSscattering_rat, discreteRHP, Residue, ODEFT_Gaussian, ODEFT_solve, grideval

u1 = (x,w) -> w ≈ 0 ? 1im*pi + 2im*atan(x) : exp(1im*x*w)*( exp(-w + loggamma(0,1im*x*w-w)) - exp(w + loggamma(0,1im*x*w+w)) + exp(w)*( x >= 0 ? 2im*pi : 0.0))
ϕrat = (x,w) -> (w <= 0 ? u1(x,w) : -u1(x,-w) |> conj)/(2im)
ϕrat∞ = (w) -> pi*exp(-abs(w))
function Φrat(x,w)
    return ϕrat(x,w)/ϕrat∞(w)
end
rat = x -> 1/(1 + x^2)

σ = 1.;
ϕGl = (x,w) -> (0.5*sqrt(pi/σ))*exp.(-w^2/(4*σ))*(1 .+ erf.(x*sqrt(σ) .+ 1im*w/(2*sqrt(σ))))
ϕG_best = (x,w) -> (sqrt(pi/σ))*exp.(-w^2/(4*σ) + 1im*w*x) .- (0.5*sqrt(pi/σ))*exp(-σ*x^2)*erfcx.(x*sqrt(σ) .+ 1im*w/(2*sqrt(σ)))
ϕG = (x,w) -> ( x < -8 ? 0.0im : ϕG_best(x,w))

ϕG∞ = (w) -> (sqrt(pi/σ))*exp.(-w^2/(4*σ))
function ΦG(x,w)
    exp.(1im*x*w).*(0.5*(1 .+ erf.(x*sqrt(σ) .+ 1im*w/(2*sqrt(σ)))))
end
Ga = x -> exp(-σ*x^2)

ϕsech∞ = (w) -> 2*atan(exp(w))

struct ODEFT
    D::SparseMatrixCSC
    T::SparseMatrixCSC
    G::Fun
    ω
end

struct ODEFTsol
    uu::Fun
    c0::ComplexF64
    ϕ::ComplexF64
    k::Float64
end

function ODEFT(g,n::Integer,m::Integer,L=1.0)
    sp1 = OscRational(0.,L);
    G = Fun(zai(g),sp1,m+1)
    D = Derivative(sp1)
    T = Array(I,n,n)
    T = hcat([1;1;zeros(n-2)],T)
    T = vcat(T,zeros(m-n,n+1))
    Dn = hcat(zeros(m),D[1:m,1:n])
    ω = k -> (4*pi*G.space.domain.L*exp(-G.space.domain.L*abs(k)))
    ODEFT(sparse(Dn),sparse(T),G,ω)
end

function ODEFT_Gaussian(g,n::Integer,m::Integer,L=1.0)
    sp1 = OscRational(0.,L);
    G = Fun(zai(g),sp1,m+1)
    D = Derivative(sp1)
    T = Array(I,n,n)
    FF = Fun(zai(Ga),sp1,m+1)
    T = vcat(T,zeros(m-n,n))
    T = hcat(-FF.coefficients[1:m],T)
    Dn = hcat(zeros(m),D[1:m,1:n])
    ODEFT(sparse(Dn),sparse(T),G,ϕG∞)
end

function  (F::ODEFT)(k)
    #ω = (4*pi*F.G.space.domain.L*exp(-F.G.space.domain.L*abs(k)))
    S = copy(F.D)
    S[:,2:end] -= 1im*k*F.T[:,2:end]
    S[:,1] -= F.T[:,1]
    m = size(S)[1]
    S = SparseMatrixCSC(S)
    u = S\pad(F.G.coefficients,m)
    (u[1]*F.ω(k), norm(S*u-pad(F.G.coefficients,m)))
end

function ODEFT_solve(F::ODEFT,k)
    S = copy(F.D)
    S[:,2:end] -= 1im*k*F.T[:,2:end]
    S[:,1] -= F.T[:,1]
    m = size(S)[1]
    u = (SparseMatrixCSC(S)\pad(F.G.coefficients,m))
    c0 = u[1]
    sp = F.G.space
    uu = Fun(sp |> ApproxFunRational.OscRational_to_OscLaurent, u[2:end] |> ApproxFunRational.coefs_to_OscLaurent)
    ODEFTsol(uu,c0,F.ω(k),k)
end

function (F::ODEFTsol)(x)
    pts = points(F.uu)
    F.c0*ϕ*ΦG(x, F.k) + F.uu(x)
end

function grideval(F::ODEFTsol)
    pts = points(F.uu)
    (F.c0*map(x -> ϕG(x, F.k), pts) + values(F.uu),pts)
end

function normgrideval(F::ODEFTsol)
    pts = points(F.uu)
    (map(x -> ΦG(x, F.k), pts) + values(F.uu)/F.c0,pts)
end

struct AKNSscattering
    qΦG
    D::SparseMatrixCSC
    T::SparseMatrixCSC
    GaFun::Fun
    Mqmat::SparseMatrixCSC
    Mrmat::SparseMatrixCSC
    Q::Fun
    R::Fun
    ϕG∞
end

function rat_qΦG(q,n,w,sp)
    Fhat = ODEFT_Gaussian(rat,n+1,n+101)
    y = ODEFT_solve(Fhat,w) |> grideval
    u = transform(sp |> ApproxFunRational.OscRational_to_OscLaurent,y[1].*map(q,y[2]))
    Fun(sp, u |> ApproxFunRational.coefs_to_OscRational)
end

function sech_qΦG(q,n,w,sp)
    Fhat = ODEFT_Gaussian(sech,n+1,n+101)
    y = ODEFT_solve(Fhat,w) |> grideval
    u = transform(sp |> ApproxFunRational.OscRational_to_OscLaurent,y[1].*map(q,y[2]))
    Fun(sp, u |> ApproxFunRational.coefs_to_OscRational)
end

function AKNSscattering_rat(f1,f2,n,m)
    L = 1.0
    sp1 = OscRational(0.,L)
    GaFun = Fun(zai(rat),sp1,m)
    Q = Fun(zai(f1),sp1,n+2)
    R = Fun(zai(f2),sp1,n+2)
    #qΦG = (q,n,w) -> Fun(x -> ϕrat(x,w)*q(x),sp1,n+2)
    qΦG = (q,n,w) -> rat_qΦG(q,n,w,sp1)
    #return (qΦG1,qΦG)
    D = Derivative(sp1)
    Dn = D[1:m,1:n+1] |> sparse
    #Dn = hcat(pad(GaFun.coefficients[1:end]/ϕrat∞(0.0),m),Dn)
    Dn = hcat(pad(GaFun.coefficients[1:end],m),Dn)
    Mr = Multiplication(R,sp1)
    Mq = Multiplication(Q,sp1)
    T = zeros(Complex{Float64},m,n+2)
    T[1:n+1,2:n+2] += Array(I,n+1,n+1)
    T = T |> SparseMatrixCSC
    Mqmat = Mq[1:m,1:n+1] |> sparse
    Mrmat = Mr[1:m,1:n+1] |> sparse
    AKNSscattering(qΦG,Dn,T,GaFun,Mqmat,Mrmat,Q,R,ϕrat∞)
end

function AKNSscattering_sech(f1,f2,n,m)
    L = 1.0
    sp1 = OscRational(0.,L)
    GaFun = Fun(zai(sech),sp1,m)
    Q = Fun(zai(f1),sp1,n+2)
    R = Fun(zai(f2),sp1,n+2)
    #qΦG = (q,n,w) -> Fun(x -> ϕrat(x,w)*q(x),sp1,n+2)
    qΦG = (q,n,w) -> sech_qΦG(q,n,w,sp1)
    #return (qΦG1,qΦG)
    D = Derivative(sp1)
    Dn = D[1:m,1:n+1] |> sparse
    #Dn = hcat(pad(GaFun.coefficients[1:end]/ϕrat∞(0.0),m),Dn)
    Dn = hcat(pad(GaFun.coefficients[1:end],m),Dn)
    Mr = Multiplication(R,sp1)
    Mq = Multiplication(Q,sp1)
    T = zeros(Complex{Float64},m,n+2)
    T[1:n+1,2:n+2] += Array(I,n+1,n+1)
    T = T |> SparseMatrixCSC
    Mqmat = Mq[1:m,1:n+1] |> sparse
    Mrmat = Mr[1:m,1:n+1] |> sparse
    AKNSscattering(qΦG,Dn,T,GaFun,Mqmat,Mrmat,Q,R,ϕsech∞)
end

function AKNSscattering(f1,f2,n,m,L = 1.0)
    sp1 = OscRational(0.,L)
    GaFun = Fun(zai(Ga),sp1,m)
    Q = Fun(zai(f1),sp1,n+2)
    R = Fun(zai(f2),sp1,n+2)
    qΦG = (q,n,w) -> Fun(x -> ϕG(x,w)*q(x),sp1,n+2)
    D = Derivative(sp1)
    Dn = D[1:m,1:n+1] |> sparse
    Dn = hcat(pad(GaFun.coefficients[1:end],m),Dn)
    Mr = Multiplication(R,sp1)
    Mq = Multiplication(Q,sp1)
    T = zeros(Complex{Float64},m,n+2)
    T[1:n+1,2:n+2] += Array(I,n+1,n+1)
    T = T |> SparseMatrixCSC
    Mqmat = Mq[1:m,1:n+1] |> sparse
    Mrmat = Mr[1:m,1:n+1] |> sparse
    AKNSscattering(qΦG,Dn,T,GaFun,Mqmat,Mrmat,Q,R,ϕG∞)
end

function _evalAKNS(A,k,bool)
    if abs(k) > 1e15
        return Array(I,2,2)*(1+0.0im)
    end
    m = size(A.D)[1]
    n = size(A.D)[2] -1
    σ = 1
    osc = A.qΦG(A.Q,m,2k*σ)
    non_osc = A.qΦG(A.R,m,0.0)
    D_osc = copy(A.D)
    #D_osc[:,1] = D_osc[:,1]*A.ϕG∞(0.0)
    D_osc -= 2im*k*σ*A.T
    Mn1 = -hcat(pad(osc.coefficients[1:end],m),A.Mqmat)
    Mn2 = -hcat(pad(non_osc.coefficients[1:end],m),A.Mrmat)

    S1 = vcat(hcat(A.D[:,1:end-1],Mn1[:,1:end-1]),hcat(Mn2[:,1:end-1],D_osc[:,1:end-1]))
    #S1 = vcat(hcat(D_osc[:,1:end-1],Mn2[:,1:end-1]),hcat(Mn1[:,1:end-1],A.D[:,1:end-1]))
    #S1[:,n+1] *= A.ϕG∞(2k)
    b = [zeros(m,1);pad(A.R.coefficients[1:end],m)]
    u1 = (S1\b)
    c1 = u1[[1;n+1]]

    σ = -1
    osc = A.qΦG(A.R,n,2k*σ)
    non_osc = A.qΦG(A.Q,n,0.0)
    D_osc -= 4im*k*σ*A.T
    Mn1 = -hcat(pad(osc.coefficients[1:end],m),A.Mrmat)
    Mn2 = -hcat(pad(non_osc.coefficients[1:end],m),A.Mqmat)

    S1 = vcat(hcat(A.D[:,1:end-1],Mn1[:,1:end-1]),hcat(Mn2[:,1:end-1],D_osc[:,1:end-1]))
    #S1[:,n+1] *= A.ϕG∞(2k)
    b = [zeros(m,1);pad(A.Q.coefficients[1:end],m)]
    u2 = (S1\b)
    c2 = u2[1:n:end]
    J = hcat(c1 ,c2 |>reverse)
    S = copy(J)
    J[1,2] *= A.ϕG∞(-2k)
    J[2,1] *= A.ϕG∞(2k)
    J[1,1] *= A.ϕG∞(0)
    J[2,2] *= A.ϕG∞(0)
    #(J + I,u1,u2)
    if bool
        return (S,u1,u2,non_osc.space)
    else
        return J + I
    end
end

function _evalAKNS1(A,k)
    if abs(k) > 1e15
        return Array(I,2,2)*(1+0.0im)
    end
    m = size(A.D)[1]
    n = size(A.D)[2] -1
    σ = 1
    osc = A.qΦG(A.Q,m,2k*σ)
    non_osc = A.qΦG(A.R,m,0.0)
    D_osc = copy(A.D)
    #D_osc[:,1] = D_osc[:,1]*A.ϕG∞(0.0)
    D_osc -= 2im*k*σ*A.T
    Mn1 = -hcat(pad(osc.coefficients[1:end],m),A.Mqmat)
    Mn2 = -hcat(pad(non_osc.coefficients[1:end],m),A.Mrmat)

    S1 = vcat(hcat(A.D[:,1:end-1],Mn1[:,1:end-1]),hcat(Mn2[:,1:end-1],D_osc[:,1:end-1]))
    #S1 = vcat(hcat(D_osc[:,1:end-1],Mn2[:,1:end-1]),hcat(Mn1[:,1:end-1],A.D[:,1:end-1]))
    #S1[:,n+1] *= A.ϕG∞(2k)
    b = [zeros(m,1);pad(A.R.coefficients[1:end],m)]
    u1 = (S1\b)
    c1 = u1[[1;n+1]]

    c1[2] *= A.ϕG∞(2k)
    c1[1] *= A.ϕG∞(0)
    c1[1] += 1.0
    return c1
end

function _evalAKNS2(A,k)
    if abs(k) > 1e15
        return Array(I,2,2)*(1+0.0im)
    end
    m = size(A.D)[1]
    n = size(A.D)[2] -1
    # σ = 1
    # osc = A.qΦG(A.Q,m,2k*σ)
    # non_osc = A.qΦG(A.R,m,0.0)

    # #D_osc[:,1] = D_osc[:,1]*A.ϕG∞(0.0)
    # D_osc -= 2im*k*σ*A.T
    # Mn1 = -hcat(pad(osc.coefficients[1:end],m),A.Mqmat)
    # Mn2 = -hcat(pad(non_osc.coefficients[1:end],m),A.Mrmat)
    #
    # S1 = vcat(hcat(A.D[:,1:end-1],Mn1[:,1:end-1]),hcat(Mn2[:,1:end-1],D_osc[:,1:end-1]))
    # #S1 = vcat(hcat(D_osc[:,1:end-1],Mn2[:,1:end-1]),hcat(Mn1[:,1:end-1],A.D[:,1:end-1]))
    # #S1[:,n+1] *= A.ϕG∞(2k)
    # b = [zeros(m,1);pad(A.R.coefficients[1:end],m)]
    # u1 = (S1\b)
    # c1 = u1[[1;n+1]]

    σ = -1
    D_osc = copy(A.D)
    osc = A.qΦG(A.R,n,2k*σ)
    non_osc = A.qΦG(A.Q,n,0.0)
    D_osc -= 2im*k*σ*A.T
    Mn1 = -hcat(pad(osc.coefficients[1:end],m),A.Mrmat)
    Mn2 = -hcat(pad(non_osc.coefficients[1:end],m),A.Mqmat)

    S1 = vcat(hcat(A.D[:,1:end-1],Mn1[:,1:end-1]),hcat(Mn2[:,1:end-1],D_osc[:,1:end-1]))
    #S1[:,n+1] *= A.ϕG∞(2k)
    b = [zeros(m,1);pad(A.Q.coefficients[1:end],m)]
    u2 = (S1\b)
    c2 = u2[1:n:end] |> reverse
    c2[1] *= A.ϕG∞(-2k)
    #J[2,1] *= A.ϕG∞(2k)
    #J[1,1] *= A.ϕG∞(0)
    c2[2] *= A.ϕG∞(0)
    c2[2] += 1.0
    #(J + I,u1,u2)
    return c2
end

function (A::AKNSscattering)(k)
    _evalAKNS(A,k,false)
end

function (A::AKNSscattering)(k,j)
    if j == 1
        _evalAKNS1(A,k)
    else
        _evalAKNS2(A,k)
    end
end

struct Jost
    k1::Number
    k2::Number
    c1::Number
    c2::Number
    u1::Fun
    u2::Fun
end

function (J::Jost)(x)
    [J.c1*ϕG(x,J.k1) + J.u1(x); J.c2*ϕG(x,J.k2) + J.u2(x)]
end

function Jost(A::AKNSscattering,k)
    (S,u1,u2,sp) = _evalAKNS(A,k,true)
    n = (length(u1) ÷ 2)-1

    (c1,c2) = S[:,1]
    f1 = Fun(sp,u1[2:n])
    f2 = Fun(sp,u1[n+3:end])
    J1 = Jost(0.0,2k,c1,c2,f1,f2)

    (c1,c2) = S[:,2]
    f2 = Fun(sp,u2[2:n])
    f1 = Fun(sp,u2[n+3:end])
    J2 = Jost(-2k,0.0,c1,c2,f1,f2)
    return (J1,J2)
end

function λ(A::AKNSscattering)
    n = size(A.D)[2]-1
    S1 = vcat(hcat(1im*A.D[1:n,2:n+1],-1im*A.Mqmat[1:n,1:n]),hcat(1im*A.Mrmat[1:n,1:n],-1im*A.D[1:n,2:n+1]))
    lam = eigvals(S1 |> Array)
    S = abs.(imag(lam)) .> 1e-1
    lam[S]
end

struct discreteRHP
    zp::Vector
    zm::Vector
    U::Array
end

function discreteRHP(zp::Vector,zm::Vector,cp::Vector,cm::Vector)
    ZP = repeat(zp,1,length(zm))
    ZM = repeat(zm,1,length(zp)) |> transpose

    M12 = -1.0./(ZP - ZM)
    M21 = -1.0./(ZM - ZP) |> transpose
    M12 = cp.*M12
    M21 = cm.*M21
    sol = vcat(hcat(I,M12),hcat(M21,I))\vcat(hcat(0*cp,cp),hcat(cm,0*cm))
    discreteRHP(zp,zm,sol)
end

function (M::discreteRHP)(z)
    r1 = [1+sum(M.U[1:length(M.zp),1]./(z .- M.zp));sum(M.U[length(M.zp)+1:end,1]./(z .- M.zm))]
    r2 = [sum(M.U[1:length(M.zp),2]./(z .- M.zp));1+sum(M.U[length(M.zp)+1:end,2]./(z .- M.zm))]
    hcat(r1,r2) |> transpose
end

function Residue(M::discreteRHP)
    r1 = [sum(M.U[1:length(M.zp),1]);sum(M.U[length(M.zp)+1:end,1])]
    r2 = [sum(M.U[1:length(M.zp),2]);sum(M.U[length(M.zp)+1:end,2])]
    hcat(r1,r2) |> transpose
end

end#module


text:
# Effective waves are related to the ansatz u ~ amps.*exp(im*k_eff*x), where k_eff is the effective wavenumber with Im(k_eff) > 0.

export wavematrix2D, wavematrix3D, wavematrix3DPlane # supplies matrices used for disperision equation and effective eignvectors
export far_field_pattern, pair_field_pattern, diff_far_field_pattern
export wavenumbers_path, wavenumbers_mesh, reduce_kvecs
export wavenumber_challis, one_species_low_wavenumber,
        two_species_approx_wavenumber

export gray_square!, gray_square

export reflection_coefficient_low_volfrac, wavenumber_low_volfrac, wavenumber_very_low_volfrac

export wienerhopf_reflection_coefficient, wienerhopf_wavevectors

include("far_fields.jl")

include("effective_waves.jl")
include("effective_wavevectors.jl")
include("reflection_effective.jl")
include("dispersion.jl")
include("wavenumber_effective.jl")
include("wavenumber_path.jl")
include("wavenumber_mesh.jl")
include("utils.jl")

include("low_frequency.jl")
include("low_volfrac.jl")
include("alternative_wavenumbers.jl")
include("two_species_approximate.jl")