text:
#==============================================================================#
# AWSLambda.jl
#
# AWS Lambda API.
#
# See http://docs.aws.amazon.com/lambda/latest/dg/API_Reference.html
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

#=
TODO

 - Test on windows.

TODO later:

 - can bindeps be deployed easily from BinaryBuilder/Linux??

 - rename functions to shorter e.g. AWSLambda.create
=#

#FIXME https://github.com/awslabs/serverless-application-model
#FIXME https://github.com/awslabs/aws-sam-local
#FIXME DLQ
#FIXME tagging
#FIXME environment variables

__precompile__()

module AWSLambda

using AWSCore
using AWSIAM
using JSON
using InfoZIP
using Retry
using SymDict
using DataStructures
using Glob
using FNVHash
using Base.Pkg
using HTTP
import HTTP: @require, precondition_error

const jl_version = "JL_$(replace(string(VERSION), ".", "_"))"
const aws_lamabda_jl_version = "0.3.0"

#-------------------------------------------------------------------------------
# AWS Lambda REST API.
#-------------------------------------------------------------------------------

function lambda(aws::AWSConfig, verb; path="", query=[], headers=Dict())

    aws = copy(aws)
    aws[:ordered_json_dict] = false

    resource = HTTP.escapepath("/2015-03-31/functions/$path")

    query = @SymDict(headers, query...)

    r = AWSCore.Services.lambda(aws, verb, resource, query)

    if isa(r, Dict)
        r = symboldict(r)
    end

    return r
end

list_lambdas(aws::AWSConfig) =
    [symboldict(f) for f in lambda(aws, "GET")[:Functions]]

list_lambdas() = list_lambdas(default_aws_config())

function lambda_configuration(aws::AWSConfig, name)

    @protected try

        return lambda(aws, "GET", path="$name/configuration")

    catch e
        @ignore if ecode(e) == "404" end
    end

    return nothing
end

lambda_configuration(name) = lambda_configuration(default_aws_config(), name)

function lambda_update_configuration(aws::AWSConfig, name, options)

    lambda(aws, "PUT", path="$name/configuration", query=options)
end

lambda_update_configuration(name, options) =
    lambda_update_configuration(default_aws_config(), name, options)

lambda_exists(aws::AWSConfig, name) = lambda_configuration(aws, name) != nothing

lambda_exists(name) = lambda_exists(default_aws_config(), name)

function create_lambda(aws::AWSConfig, name;
                       ZipFile=nothing,
                       S3Key="$name.zip",
                       S3Bucket=get(aws, :lambda_bucket, nothing),
                       Handler="lambda_function.lambda_handler",
                       Role=role_arn(aws, "jl_lambda_eval_lambda_role"),
                       Runtime="python2.7",
                       MemorySize = 1536,
                       Timeout=300,
                       args...)

    if ZipFile != nothing
        ZipFile = base64encode(ZipFile)
        Code = @SymDict(ZipFile)
    else
        Code = @SymDict(S3Key, S3Bucket)
    end

    query = @SymDict(FunctionName = name,
                     Code,
                     Handler,
                     Role,
                     Runtime,
                     MemorySize,
                     Timeout,
                     args...)

    @repeat 5 try

        lambda(aws, "POST", query=query)

    catch e
        # Retry in case Role was just created and is not yet active...
        @delay_retry if ecode(e) == "400" end
    end
end

create_lambda(name; kw...) = create_lambda(default_aws_config(), name; kw...)

function update_lambda(aws::AWSConfig, name;
                       ZipFile=nothing,
                       S3Key="$name.zip",
                       S3Bucket=get(aws, :lambda_bucket, nothing),
                       args...)

    if ZipFile != nothing
        ZipFile = base64encode(ZipFile)
        query = @SymDict(ZipFile)
    else
        query = @SymDict(S3Key, S3Bucket)
    end

    lambda(aws, "PUT", path="$name/code", query=query)

    if !isempty(args)
        lambda(aws, "PUT", path="$name/configuration", query=@SymDict(args...))
    end
end

update_lambda(name; kw...) = update_lambda(default_aws_config(), name; kw...)

function lambda_publish_version(aws::AWSConfig, name, alias)

    r = lambda(aws, "POST", path="$name/versions")
    @protected try
        lambda_create_alias(aws, name, alias, FunctionVersion=r[:Version])
    catch e
        @ignore if ecode(e) == "409"
            lambda_update_alias(aws, name, alias, FunctionVersion=r[:Version])
        end
    end
end

lambda_publish_version(name, alias) =
    lambda_publish_version(default_aws_config(), name, alias)

function lambda_create_alias(aws::AWSConfig, name, alias;
                             FunctionVersion="\$LATEST")

    lambda(aws, "POST", path="$name/aliases",
                        query=@SymDict(FunctionVersion, Name=alias))
end

lambda_create_alias(name, alias; kw...) =
    lambda_create_alias(default_aws_config(), name, alias; kw...)

function lambda_update_alias(aws::AWSConfig, name, alias;
                             FunctionVersion="\$LATEST")

    lambda(aws, "PUT", path="$name/aliases/$alias",
                       query=@SymDict(FunctionVersion))
end

lambda_update_alias(name, alias; kw...) =
    lambda_update_alias(default_aws_config(), name, alias; kw...)

function lambda_add_permission(aws::AWSConfig, name, permission)

    lambda(aws, "POST", path="$name/policy", query=permission)
end

lambda_add_permission(name, permission) =
    lambda_add_permission(default_aws_config(), name, permission)

function lambda_delete_permission(aws::AWSConfig, name, id)

    lambda(aws, "DELETE", path="$name/policy/$id")
end

lambda_delete_permission(name, id) =
    lambda_delete_permission(default_aws_config(), name, id)

function lambda_delete_permissions(aws::AWSConfig, name)

    for p in lambda_get_permissions(aws, name)
        lambda_delete_permission(aws, name, p["Sid"])
    end
end

lambda_delete_permissions(name) =
    lambda_delete_permissions(default_aws_config(), name)

function lambda_get_permissions(aws::AWSConfig, name)

    @protected try
        r = lambda(aws, "GET", path="$name/policy")
        return JSON.parse(r[:Policy])["Statement"]
    catch e
        @ignore if ecode(e) == "404"
            return Dict[]
        end
    end
end

lambda_get_permissions(name) =
    lambda_get_permissions(default_aws_config(), name)

function delete_lambda(aws::AWSConfig, name)

    @protected try
        lambda(aws, "DELETE", path=name)
    catch e
        @ignore if ecode(e) == "404" end
    end

    return
end

delete_lambda(name) = delete_lambda(default_aws_config(), name)

export AWSLambdaException

type AWSLambdaException <: Exception
    name::String
    message::String
end

function Base.show(io::IO, e::AWSLambdaException)

    info = try
        JSON.parse(e.message)
    catch
        Dict("message" => e.message)
    end
    println(io, string("AWSLambdaException \"", e.name, "\":\n",
                       info["message"], "\n"))
    for (k, v) in info
        if k != "message"
            println(io, "$k: $v")
        end
    end
end

function invoke_lambda(aws::AWSConfig, name, args::Dict; async=false)

    @protected try

        r = lambda(aws, "POST",
                        path="$name/invocations",
                        headers=Dict("X-Amz-Invocation-Type" =>
                                     async ? "Event" : "RequestResponse"),
#FIXME async not working?
                        query=args)

        if isa(r, Dict) && haskey(r, :errorMessage)
            throw(AWSLambdaException(string(name), r[:errorMessage]))
        end

        return r

    catch e
        if ecode(e) == "429"
            e.message = e.message *
                " See http://docs.aws.amazon.com/lambda/latest/dg/limits.html"
        end
    end

    @assert false # Unreachable
end

invoke_lambda(name, args::Dict; kw...) = invoke_lambda(default_aws_config(),
                                                       name, args; kw...)

invoke_lambda(aws::AWSConfig, name; args...) =
    invoke_lambda(aws, name, symboldict(args))

invoke_lambda(name; args...) =
    invoke_lambda(default_aws_config(), name; args...)

async_lambda(aws::AWSConfig, name, args) =
    invoke_lambda(aws, name, args; async=true)

async_lambda(name, args) = async_lambda(default_aws_config(), name, args)

function create_lambda_role(aws::AWSConfig, name, policy="")

    name = "$(name)_lambda_role"

    @protected try

        r = AWSIAM.iam(aws, Action = "CreateRole",
                        Path = "/",
                        RoleName = name,
                        AssumeRolePolicyDocument = """{
                            "Version": "2012-10-17",
                            "Statement": [ {
                                "Effect": "Allow",
                                "Principal": {
                                    "Service": "lambda.amazonaws.com"
                                },
                                "Action": "sts:AssumeRole"
                            } ]
                        }""")

    catch e
        @ignore if ecode(e) == "EntityAlreadyExists" end
    end

    AWSIAM.iam(aws, Action = "AttachRolePolicy",
                    RoleName = name,
                    PolicyArn = "arn:aws:iam::aws:policy" *
                                "/service-role/AWSLambdaBasicExecutionRole")

    if policy != ""
        AWSIAM.iam(aws, Action = "PutRolePolicy",
                        RoleName = name,
                        PolicyName = name,
                        PolicyDocument = policy)
    end

    return role_arn(aws, name)
end

#-------------------------------------------------------------------------------
# Python Code Support.
#-------------------------------------------------------------------------------

function create_py_lambda(aws::AWSConfig, name, py_code;
                          Role = create_lambda_role(aws, name))

    options = @SymDict(Role, ZipFile = create_zip("lambda_function.py" => py_code))

    old_config = lambda_configuration(aws, name)

    if old_config == nothing
        create_lambda(aws, name; options...)
    else
        update_lambda(aws, name; options...)
    end
end

create_py_lambda(name, py_code; kw...) =
    create_py_lambda(default_aws_config(), name, py_code; kw...)

#-------------------------------------------------------------------------------
# Julia Modules.
#-------------------------------------------------------------------------------

function local_module_cache()

    # List of modules in local ".ji" cache.
    r = [Symbol(splitext(f)[1]) for f in
            [[readdir(p) for p in
                filter(isdir, Base.LOAD_CACHE_PATH)]...;]]

    # List of modules compiled in to sys image.
    append!(r, filter(x->Main.eval(:(try isa($x, Module)
                                     catch ex
                                         if typeof(ex) != UndefVarError
                                            rethrow(ex)
                                         end
                                         false
                                     end)), names(Main)))

    return unique(r)
end

function __init__()
    global _lambda_module_cache = Symbol[]
end

global _lambda_module_cache = Symbol[]

global default_lambda_module_cache = [
    :AWSCore,
    :AWSIAM,
    :AWSLambda,
    :AWSS3,
    :Base,
    :Compat,
    :Core,
    :DataStructures,
    :FNVHash,
    :Glob,
    :HTTP,
    :InfoZIP,
    :IniFile,
    :IterTools,
    :JSON,
    :LightXML,
    :Main,
    :MbedTLS,
    :Nullables,
    :Retry,
    :SymDict,
    :XMLDict]

# List of modules in the Lambda sandbox ".ji" cache.

function lambda_module_cache(aws::AWSConfig = default_aws_config())

    global _lambda_module_cache
    if isempty(_lambda_module_cache)

        @protected try
            _lambda_module_cache = lambda_eval(aws, :(filter(x->Main.eval(:(
                try
                    isa($x, Module) && !isfile(string("/tmp/julia/", $x, ".ji"))
                catch ex
                    if typeof(ex) != UndefVarError
                        rethrow(ex)
                    end
                    false
                end
            )), names(Main))))
        catch e
            if ecode(e) == "404"
                return default_lambda_module_cache
            end
        end

    end
    return _lambda_module_cache
end

# List of source files required by "modules".

function precompiled_module_files(aws::AWSConfig, modules::Vector{Symbol})

    exclude = lambda_module_cache(aws)

    modules = collect(filter(m->!(m in exclude), modules))

    return unique([[_precompiled_module_files(i, exclude) for i in modules]...;])
end

function _precompiled_module_files(m, exclude_modules)

    # Module must be precompiled...
    if !(m in local_module_cache())
        error("$m not precompiled. Do \"using $m\".")
    end

    r = Dict()
    for p in Base.find_all_in_cache_path(m)
        modules, files = Base.cache_dependencies(p)
        for (f,t) in files
            r[f] = nothing
        end
        for (m,t) in modules
            if m in exclude_modules
                continue
            end
            for f in _precompiled_module_files(m, exclude_modules)
                r[f] = nothing
            end
        end
    end
    return collect(keys(r))
end

# Split paths from common prefix...

function path_prefix_split(paths::Vector)

    if length(paths) == 0
        return "", []
    end

    # Find longest common prefix...
    i = 1
    for i = 1:length(paths[1])
        if any(f -> length(f) <= i || f[1:i] != paths[1][1:i], paths)
            break
        end
    end

    # Ensure prefix is a dir path...
    while i > 0 && any(f->!isdirpath(f[1:i]), paths)
        i -= 1
    end

    return paths[1][1:i], [p[i+1:end] for p in paths]
end

# Find load path and source files for "modules"...

function module_files(aws::AWSConfig, modules::Vector{Symbol})

    if length(modules) == 0
        return [], OrderedDict()
    end

    pkgd = realpath(Pkg.dir())

    # Build a dict of files for each load path location...
    d = Dict()
    add_file = (p,f) -> (if !haskey(d, p); d[p] = [] end; push!(d[p], f))
    for file in precompiled_module_files(aws, modules)
        file = realpath(file)

        if startswith(file, pkgd)
            add_file(pkgd, file[length(pkgd)+2:end])
            continue
        end

        for p in filter(isdir, LOAD_PATH)
            p = realpath(abspath(p))
            if joinpath(p, basename(file)) == file
                add_file(p, basename(file))
                break
            end
        end
    end

    if isempty(d)
        return [], OrderedDict()
    end

    # Remove common prefix from load path locations...
    load_path = collect(Base.Iterators.filter(p->p != pkgd, keys(d)))
    prefix, short_load_path = path_prefix_split(load_path)
    push!(short_load_path, basename(pkgd))

    # Build archive of file content...
    archive = OrderedDict()
    for (p, s) in zip([load_path...; pkgd],
                      [short_load_path...; basename(pkgd)])
        for f in get(d,p,[])
            path = joinpath("julia", s, f)
            path = replace(path, "\\", "/")
            archive[path] = read(joinpath(p, f))
        end
    end

    return short_load_path, archive
end

#-------------------------------------------------------------------------------
# Julia Code Support.
#-------------------------------------------------------------------------------

# Base64 representation of Julia objects...

function serialize64(x)

    buf = IOBuffer()
    b64 = Base64EncodePipe(buf)
    serialize(b64, x)
    close(b64)
    String(take!(buf))
end

# Invoke a Julia AWS Lambda function.
# Serialise "args" and deserialise result.

function invoke_jl_lambda(aws::AWSConfig, name::String, args...)

    r = invoke_lambda(aws, name, jl_data = serialize64(args))
    try
        println(r[:stdout])
    end
    try
        return deserialize(Base64DecodePipe(IOBuffer(r[:jl_data])))
    end
    return r
end

invoke_jl_lambda(name::String, args...) =
    invoke_jl_lambda(default_aws_config(), name, args...)

# Returns an anonymous function that calls "func" in the Lambda sandbox.
# e.g. lambda_function(println)("Hello")
#      lambda_function(x -> x*x)(4)
#      lambda_function(s -> JSON.parse(s))("{}")

lambda_function(aws, f) =
    (a...) -> invoke_jl_lambda(aws, "jl_lambda_eval:$jl_version",
                                    eval(Main,:(()->$f($a...))))

lambda_function(f) = lambda_function(default_aws_config(), f)

# Prepend `body` with expression to install `modules` under "/tmp"...
function embed_modules_in_body(aws, body, modules)

    load_path, mod_files = module_files(aws, modules)
    quote
        for (path, code) in $([(k, v) for (k,v) in mod_files])
            path = "/tmp/$path"
            println("$path, $(length(code)) bytes")
            mkpath(dirname(path))
            write(path, code)
        end
        if !("/tmp/julia" in Base.LOAD_CACHE_PATH)
            insert!(Base.LOAD_CACHE_PATH, 1, "/tmp/julia")
        end
        for dir in ["/tmp/julia/$p" for p in $load_path]
            if !(dir in Base.LOAD_PATH)
                push!(Base.LOAD_PATH, dir)
            end
        end
        for u in $modules
            eval(Main, :(using $u))
        end
        Base.invokelatest(()->$body)
    end
end

"""
    @lambda [using ...] function name(args...) body end [aws_config]

Define local `function` that executes `body` on AWS Lambda.
"""
macro lambda(args...)

    @require 1 <= length(args) <= 3
    @require all(x-> x isa Expr, args[1:end-1])
    @require args[1].head in (:using, :toplevel, :function)
    @require args[1].head != :toplevel || all(x->x.head == :using, args[1].args)
    @require args[1].head == :function || length(args) > 1 &&
                                          args[2].head == :function

    # Optional last argument is AWSConfig.
    if !(args[end] isa Expr) || args[end].head != :function
        aws = args[end]
    else
        aws = :(Main.AWSCore.default_aws_config())
    end

    if args[1].head == :function
        f = args[1]
        body = f.args[2]
        modules = Symbol[]
    else
        # Optional first argument is `using Foo, Bar, ...`
        f = args[2]
        modules = args[1]
        if modules.head == :using
            modules = Expr(:toplevel, modules)
        end
        @assert modules.head == :toplevel
        modules = [x.args[1] for x in modules.args]
        body = f.args[2]
        body = embed_modules_in_body(eval(aws), body, modules)
    end

    # Replace function body with Lambda invocation.
    args = f.args[1].args[2:end]
    arg_names = [isa(a, Expr) ? a.args[1] : a for a in args]
    l = Expr(:quote, Expr(:->, Expr(:tuple, args...), body))
    f.args[2] = :(AWSLambda.lambda_function($aws, eval(Main, $l))($(arg_names...)))

    return esc(f)
end

# Evaluate "expr" in the Lambda sandbox.

lambda_eval(aws, expr) =
    invoke_jl_lambda(aws, "jl_lambda_eval:$jl_version", expr)

lambda_eval(expr) = lambda_eval(default_aws_config(), expr)

macro lambda_eval(expr)
    if expr.head == :block
        expr = [e for e in expr.args]
    else
        expr = QuoteNode(expr)
    end
    :(lambda_eval($expr))
end

# Evaluate "code" in the Lambda sandbox.

lambda_include_string(aws, code) = lambda_function(aws, include_string)(code)

lambda_include_string(code) = lambda_include_string(default_aws_config(), code)

# Evaluate "filename" in the Lambda sandbox.

function lambda_include(aws::AWSConfig, filename)
    code = readstring(filename)
    lambda_function(aws, include_string)(code, filename)
end

lambda_include(filename) = lambda_include(default_aws_config(), filename)

                                                           # For build_sysimg.jl
                                       @static if isdefined(Base, :uv_eventloop)

# Create an AWS Lambda to run "jl_code".

function create_jl_lambda(name, jl_code, modules=Symbol[], options=SymbolDict())

    options = Dict{Symbol,Any}(copy(options))

    if !haskey(options, :aws)
        aws = default_aws_config()
    else
        aws = options[:aws]
        delete!(options, :aws)
    end

    # Find files and load path for required modules...
    load_path, mod_files = module_files(aws, modules)
    full_load_path = join([":/var/task/julia/$p" for p in load_path])
    py_config = ["os.environ['JULIA_LOAD_PATH'] += '$full_load_path'\n"]

    if AWSCore.debug_level > 0
        println("create_jl_lambda($name)")
        for f in keys(mod_files)
            println("    $f")
        end
    end

    # Get env vars from "aws" or "options"...
    env = get(options, :env,
          get(aws, :lambda_env, Dict()))
    options[:env] = env
    for (n,v) in env
        push!(py_config, "os.environ['$n'] = \"\"\"$v\"\"\"\n")
    end

    # Start with ZipFile from options...
    if !haskey(options, :ZipFile)
        options[:ZipFile] = UInt8[]
    end

    if AWSCore.debug_level > 1
        @show py_config
    end

    # Add lambda source and module files to zip...
    open_zip(options[:ZipFile]) do z
        merge!(z, mod_files)
        z["lambda_config.py"] = join(py_config)
        z["module_$name.jl"] = jl_code
    end

    # FNV hash of deployed Julia code is stored in the Description field.
    old_config = lambda_configuration(aws, name)
    new_code_hash = options[:ZipFile] |> open_zip |> Dict |>
                    serialize64 |> fnv32 |> hex
    old_code_hash = old_config == nothing ? nothing :
                    get(old_config, :Description, nothing)

    # Don't create a new lambda if one already exists with same code...
    if new_code_hash == old_code_hash && !get(aws, :lambda_force_update, false)
        return
    end

    options[:Description] = new_code_hash

    deploy_lambda(aws, name, load_path, options, old_config == nothing)
end

@lambda using AWSLambda, AWSS3, InfoZIP function deploy_lambda(
                                          aws, name, load_path, options, is_new)

    println("Foo")

    AWSCore.set_debug_level(1)

    # Unzip lambda source and modules files to /tmp...
    mktempdir() do tmpdir

        for (n, v) in options[:env]
            ENV[n] = replace(v, r"^/var/task", tmpdir)
        end

        if haskey(options, :ZipURL)
            InfoZIP.unzip(Requests.get(options[:ZipURL]).data, tmpdir)
        end

        InfoZIP.unzip(options[:ZipFile], tmpdir)

        # Create module precompilation directory under /tmp...
        v = "v$(VERSION.major).$(VERSION.minor)"
        ji_path = "$tmpdir/julia/lib/$v"
        mkpath(ji_path)

        # Run precompilation...
        cmd = "push!(LOAD_PATH, \"$tmpdir\")\n"
        for p in load_path
            cmd *= "push!(LOAD_PATH, \"$tmpdir/julia/$p\")\n"
        end
        cmd *= "insert!(Base.LOAD_CACHE_PATH, 1, \"$ji_path\")\n"
        cmd *= "using module_$name\n"
        println(cmd)
        run(`$JULIA_HOME/julia -e $cmd`)
        run(`chmod -R a+r $ji_path`)

        # Create new ZIP combining base image and new files from /tmp...
        run(`rm -f /tmp/lambda.zip`)
        run(Cmd(`zip -q --symlinks -r -9 /tmp/lambda.zip .`,
                dir="/var/task"))
        run(Cmd(`zip -q --symlinks -r -9 /tmp/lambda.zip .`,
                dir=tmpdir))
        options[:ZipFile] = read("/tmp/lambda.zip")

        run(`rm -f /tmp/lambda.zip`)
    end

    # Deploy the lambda to AWS...
    if is_new
        r = AWSLambda.create_lambda(aws, name; options...)
    else
        r = AWSLambda.update_lambda(aws, name; options...)
    end

end

# Create an AWS Lambda function.
#
# e.g.
#
#   AWSLambda.@deploy function hello(a, b)
#
#       message = "Hello $a$b"
#
#       println("Log: $message")
#
#       return message
#   end aws
#
#   hello("World", "!") # FIXME invoke lambda
#   Hello World!
#
# @deploy deploys an AWS Lambda that contains the body of the Julia function.

"""
    AWSLambda.@deploy [(option=value, ...)] [using ...] function name(args...)
        function_body
    end

Deploy `function` to AWS Lambda.
"""
macro deploy(options_ex, using_ex, function_ex)
    deploy_ex(options_ex, using_ex, function_ex)
end

macro deploy(ex, function_ex)
    if ex.head in (:toplevel, :using)
        deploy_ex(:(Dict{Symbol,Any}()), ex, function_ex)
    else
        deploy_ex(ex, Expr(:toplevel), function_ex)
    end
end

macro deploy(function_ex)
    deploy_ex(:(Dict{Symbol,Any}()), Expr(:toplevel), function_ex)
end

function deploy_ex(options_ex, using_ex, function_ex)

    # Wrap single-module using expression with :toplevel
    if using_ex.head == :using
        using_ex = Expr(:toplevel, using_ex)
    end

    @require using_ex.head == :toplevel
    @require all(x->x.head == :using, using_ex.args)
    @require function_ex.head == :function

    modules = Symbol[x.args[1] for x in using_ex.args]

    # Rewrite function name to be :lambda_function...
    call, body = function_ex.args
    name = call.args[1]
    call.args[1] = :lambda_function
    args = function_ex.args[1].args[2:end]

    # Generate code to extract args from event Dict...
    arg_names = [isa(a, Expr) ? a.args[1] : a for a in args]
    get_args = Expr(:tuple, [:(event[$(string(a))]) for a in arg_names]...)

    jl_code = """
        __precompile__()

        module $(Symbol("module_$name"))

        $(join(["using $m" for m in modules], "; "))

        $function_ex

        function lambda_function_with_event(event::Dict{String,Any})
            lambda_function$get_args
        end

        end
        """

    :(create_jl_lambda($(string(name)), $jl_code, $modules, $(esc(options_ex))))
end
                                            end # isdefined(Base, :uv_eventloop)
                                                           # For build_sysimg.jl

#-------------------------------------------------------------------------------
# Deploy pre-cooked Julia Base Lambda
#-------------------------------------------------------------------------------

function deploy_jl_lambda_eval(aws::AWSConfig = default_aws_config();
    bucket = "octech.com.au.$(aws[:region]).awslambda.jl.deploy",
    base_zip = "jl_lambda_eval_$(VERSION)_$(aws_lamabda_jl_version).zip")

    old_config = lambda_configuration(aws, "jl_lambda_eval")

    options = @SymDict(S3Key=base_zip,
                       S3Bucket=bucket,
                       Role = create_lambda_role(aws, "jl_lambda_eval"))

    if old_config == nothing
        create_lambda(aws, "jl_lambda_eval"; options...)
    else
        update_lambda(aws, "jl_lambda_eval"; options...)
    end

    lambda_publish_version(aws, "jl_lambda_eval", jl_version)
end

end # module AWSLambda

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


text:
struct PIR <: FileFormat end

# PIR Parser
# ==========

# Description of PIR/NBRF format in :
# - https://salilab.org/modeller/9v7/manual/node445.html
# - http://iubio.bio.indiana.edu/soft/molbio/readseq/classic/src/Formats
# - http://www.bioinformatics.nl/tools/crab_pir.html
# - http://emboss.sourceforge.net/docs/themes/seqformats/NbrfFormat.html

# They could be annotations/comments after the sequence, we are ignoring that at the moment.

function _pre_readpir(io::Union{IO, AbstractString})
    IDS = String[]
    SEQS = String[]
    GS = Dict{Tuple{String,String},String}()

    # Sequence types from http://iubio.bio.indiana.edu/soft/molbio/readseq/classic/src/Formats
    nucleic_types = Set(["DL",  # DNA, linear
                         "DC",  # DNA, circular
                         "RL",  # RNA, linear
                         "RC",  # RNA, circular
                         "N1",  # functional RNA, other than tRNA
                         "N3"]) # tRNA

    nonres = r"[^A-Za-z-\.]+" # large negated character class becuase seqs are free text

    finished = false
    line_number = 0
    seq_number = 0
    for line::String in lineiterator(io)
        line_number += 1
        if startswith(line, '>')
            line_number = 1
            seq_number += 1
            finished = false
            m = match(r"^>([A-Z][A-Z0-9]);(\S+)", line) # e.g. >P1;5fd1
            if m !== nothing
                seq_type = m[1]
                seq_id   = m[2]
                push!(IDS, seq_id)
                push!(GS, (seq_id, "Type") => seq_type)
                if seq_type in nucleic_types
                    warn("Type $seq_type indicates that $seq_id is a nucleic acid.")
                end
            else
                throw(ErrorException("Identifier doesn't match the PIR format: $line"))
            end
        elseif line_number == 2 && !isempty(IDS)
            push!(GS, (IDS[seq_number], "Title") => line)
        elseif line_number == 3
            if endswith(line, '*')
                finished = true
            end
            push!(SEQS, replace(line, nonres => ""))
        elseif !finished && length(SEQS) == seq_number && seq_number != 0
            if endswith(line, '*')
                finished = true
            end
            @inbounds SEQS[seq_number] = SEQS[seq_number] * replace(line, nonres => "")
        end
    end

    (IDS, SEQS, GS)
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{PIR},
                   output::Type{AnnotatedMultipleSequenceAlignment};
                   generatemapping::Bool=false,
                   useidcoordinates::Bool=false,
                   deletefullgaps::Bool=true,
                   keepinserts::Bool=false)::AnnotatedMultipleSequenceAlignment
    IDS, SEQS, GS = _pre_readpir(io)
    annot = Annotations()
    annot.sequences = GS
    _generate_annotated_msa(annot, IDS, SEQS, keepinserts, generatemapping,
                            useidcoordinates, deletefullgaps)
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{PIR},
                   output::Type{NamedResidueMatrix{Array{Residue,2}}};
                   deletefullgaps::Bool=true)::NamedResidueMatrix{Array{Residue,2}}
    IDS, SEQS, _ = _pre_readpir(io)
    msa = _generate_named_array(SEQS, IDS)
    if deletefullgaps
        return deletefullgapcolumns(msa)
    end
    msa
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{PIR},
                   output::Type{MultipleSequenceAlignment};
                   deletefullgaps::Bool=true)::MultipleSequenceAlignment
    msa = parse(io,
                format,
                NamedResidueMatrix{Array{Residue,2}},
                deletefullgaps=deletefullgaps)
    MultipleSequenceAlignment(msa)
end

function Base.parse(io::Union{IO,AbstractString},
                   format::Type{PIR},
                   output::Type{Matrix{Residue}};
                   deletefullgaps::Bool=true)::Matrix{Residue}
    IDS, SEQS, _ = _pre_readpir(io)
    _strings_to_matrix_residue_unsafe(SEQS, deletefullgaps)
end

function Base.parse(io, format::Type{PIR};
                    generatemapping::Bool=false,
                    useidcoordinates::Bool=false,
                    deletefullgaps::Bool=true,
                    keepinserts::Bool=false)::AnnotatedMultipleSequenceAlignment
    parse(io, PIR, AnnotatedMultipleSequenceAlignment,
          generatemapping=generatemapping,
          useidcoordinates=useidcoordinates,
          deletefullgaps=deletefullgaps,
          keepinserts=keepinserts)
end

# Print Pfam
# ==========

function _print_pir_seq(io::IO, seq_type, seq_id, seq_title, seq)
    println(io, ">", seq_type, ';', seq_id)
    println(io, seq_title)
    char_counter = 0
    for char in seq
        print(io, char)
        char_counter += 1
        if char_counter % 80 == 0
            println(io)
        end
    end
    println(io, '*')
end

function Base.print(io::IO, msa::AnnotatedMultipleSequenceAlignment,
                    format::Type{PIR})
    seqann = getannotsequence(msa)
    seqnames = sequencenames(msa)
    for i in 1:nsequences(msa)
        seq_id = seqnames[i]
        if haskey(seqann, (seq_id, "Type"))
            seq_type = seqann[(seq_id, "Type")]
        else
            warn("There is not sequence Type annotation for $seq_id, using XX (Unknown).")
            seq_type = "XX"
        end
        if haskey(seqann, (seq_id, "Title"))
            seq_title = seqann[(seq_id, "Title")]
        else
            warn("There is not sequence Title annotation for $seq_id.")
            seq_title = ""
        end
        seq = stringsequence(msa, i)
        _print_pir_seq(io, seq_type, seq_id, seq_title, seq)
    end
end

function Base.print(io::IO, msa::AbstractMatrix{Residue}, format::Type{PIR})
    seqnames = sequencenames(msa)
    for i in 1:nsequences(msa)
        _print_pir_seq(io, "XX", seqnames[i], "", stringsequence(msa, i))
    end
end

ext: jl


text:
module TestDiversity
using Test

using Diversity

@testset "Deprecations" begin
    # Only run deprecation checks on macs
end

end


text:
module WhittakerTests

using FactCheck, Whittaker

include("whittaker1.jl")
include("whittaker2.jl")

isinteractive() || FactCheck.exitstatus()

end


text:
struct JC69abs <: JC69
    λ::Float64
    function JC69abs(λ::Float64)
        if λ <= 0.
            error("JC69 parameter λ must be positive")
        end
        new(λ)
    end
end

function show(io::IO, object::JC69abs)
  print(io, "\r\e[0m\e[1mJ\e[0mukes and \e[1mC\e[0mantor 19\e[1m69\e[0m model (absolute rate form)
λ = $(object.λ)")
end

JC69(λ) = JC69abs(λ)

function Q(mod::JC69abs)
  λ = mod.λ
  Q₁ =  0.25 * λ
  Q₂ = -(Q₁ * 3)
  return Qmatrix(Q₂, Q₁, Q₁, Q₁,
                 Q₁, Q₂, Q₁, Q₁,
                 Q₁, Q₁, Q₂, Q₁,
                 Q₁, Q₁, Q₁, Q₂)
end

function P(mod::JC69abs, t::Float64)
  if t < 0.0
    error("t must be positive")
  end
  λ = mod.λ
  ω = exp(-t * λ)
  P₁ = 0.25 + 0.75 * ω
  P₂ = 0.25 + 0.25 * ω
  return Pmatrix(P₁, P₂, P₂, P₂,
                 P₂, P₁, P₂, P₂,
                 P₂, P₂, P₁, P₂,
                 P₂, P₂, P₂, P₁)
end


text:
using Test, SimpleNL, MadNLP, Ipopt, ForwardDiff, SpecialFunctions
import Random: seed!
const eps = 1e-6
seed!(0)

compare(p1,p2) = maximum(abs.(p1 .- p2)) <= eps

function compare(f1,f2,xpoints,ppoints)
    for i=1:length(xpoints)
        compare(f1(xpoints[i],ppoints[i]),f2(xpoints[i],ppoints[i])) || return false
    end
    return true
end

function compare(f1,f2,y1,y2,xpoints,ppoints)
    for i=1:length(xpoints)
        f1(y1,xpoints[i],ppoints[i])
        f2(y2,xpoints[i],ppoints[i])
        compare(y1,y2) || return false
    end
    return true
end

@testset "SimpleNLModels test" begin
    @testset "Function Test" begin
        include("function_test.jl")
    end
    @testset "Print Test" begin
        include("print_test.jl")
    end
    @testset "Model Test" begin
        include("model_test.jl")
    end
end


text:
function sizeforfft(len::Int, T::DataType=Float64)
    len2 = nextpow(2,len)
    zeros(Complex{T}, len2)
end

mutable struct PaddedFFT{T} <: Function
    U::Vector{Complex{T}}
    PaddedFFT(n::Int, D=Float64) = n > 0 ? new{D}(sizeforfft(n, D)) : error("n=$n must be than 0.")
    PaddedFFT{T}(v::Vector{Complex{T}}) where {T <: Union{Float32, Float64}} = new{T}(sizeforfft(length(v), T))
end

function (pfft::PaddedFFT{T})(u2::Vector{Complex{T}}) where {T <: Real}
    len = length(u2)
    fill!(pfft.U, 0+0*im)
    @inbounds pfft.U[1:len] = u2
    fft!(pfft.U)
    nothing
end

# using BenchmarkTools
# u2 = randn(2^20 - 20000) .+ 0im;
# ffft64 = PaddedFFT{Float64}(u2)
# size(ffft64.U)
# @btime ffft64(u2)
# ffft64(u2)
# u32 = randn(Float32, 2^20 - 20000) .+ 0im;
# ffft32 = PaddedFFT{Float32}(length(u32),Float32)
# @btime ffft32(u32)


text:
abstract LevelMethods

type CenterLevels <: LevelMethods
  nlevels :: Int64
  rangeCalculator :: Function
end

function CenterLevels(; nlevels :: Int64 = 5, center_by :: Function = x-> mean(x))
  rangeCalculator = function ( x :: Array{Float64})
    center = center_by(x)
    r1 = 3*std(x)
    return (center - r1, center + r1)
  end
  return CenterLevels(nlevels, rangeCalculator)
end
export CenterLevels

type PDFlevels <: LevelMethods
  nlevels :: Int64
  low_cutoff :: Float64
  high_cutoff :: Float64
end

function PDFlevels(; nlevels :: Int64 = 5, low_cutoff :: Float64 = 0.15)
  return PDFlevels( nlevels, low_cutoff, 1-low_cutoff)
end

function levels( xarr :: Array{Float64} , m :: PDFlevels)
  x1 = sort(xarr[:])
  lx = length(xarr)
  low = x1[floor(Int64,m.low_cutoff * lx)]
  high = x1[floor(Int64,m.high_cutoff * lx)]
  rr = collect( low:(high-low)/(m.nlevels):high )
  return ((:levels, rr), (:extend, "both"))
end

function levels( xarr :: Array{Float64} , m :: CenterLevels)
  (low, high) = m.rangeCalculator(xarr)
  rr = collect( low:(high-low)/(m.nlevels):high )
  return ((:levels, rr), (:extend, "both"))
end

export levels
export PDFlevels
export CenterLevels


text:
# Compute network force components
# Alex Tam, 12/10/2020

"Compute the contractile/expansive forces in the network"
function network_force(s, s_old, af, xl, mm, parN, parA, parM, Lxx, Lxy, Lyx, Lyy)
    x::Vector{Real} = build_dof(s); # Build DOF vector
    # Compute the network force (+ve expansive, -ve contractile)
    fxx = -ForwardDiff.derivative(a -> energy_functional(x, s_old, af, mm, parN, parA, parM, a, Lxy, Lyx, Lyy), Lxx::Real);
    fxy = -ForwardDiff.derivative(a -> energy_functional(x, s_old, af, mm, parN, parA, parM, Lxx, a, Lyx, Lyy), Lxy::Real);
    fyx = -ForwardDiff.derivative(a -> energy_functional(x, s_old, af, mm, parN, parA, parM, Lxx, Lxy, a, Lyy), Lyx::Real);
    fyy = -ForwardDiff.derivative(a -> energy_functional(x, s_old, af, mm, parN, parA, parM, Lxx, Lxy, Lyx, a), Lyy::Real);
    return [fxx, fxy, fyx, fyy]
end


text:
"""
This struct contains data for the energy market product.
    name: Energy.
    capacity_factors: Dictionary of expected capacity factors for each scenario.
    marginal_cost: Marginal cost per unit of energy provided.
"""
mutable struct Energy <: OperatingProduct
    name::Symbol
    capacity_factors::Dict{String, Array{Float64, 2}}
    marginal_cost::Float64
end

# Capacity factors only returned when product is of type Energy
function get_capacity_factors(prod::T) where T<:Product
    return
end

get_capacity_factors(prod::Energy) = prod.capacity_factors

# Capacity factors only set when product is of type Energy
function set_capacity_factors!(product::T,
                              scenario_name::String,
                              capacity_factors::Array{Float64, 2}) where T <: Product
    return
end

function set_capacity_factors!(product::Energy,
                              scenario_name::String,
                              capacity_factors::Array{Float64, 2})
    product.capacity_factors[scenario_name] = capacity_factors
    return
end


text:
"""
This `OCIM2` module is used to load the OCIM2 matrices and grid for use in AIBECS.

!!! tip
    To load the default OCIM2 matrix and grid, do
    ```
    julia> grd, T = OCIM2.load()
    ```
    But you can also load the other matrices by specifying which version you want, e.g.,
    ```
    julia> grd, T = OCIM2.load(version="OCIM2_KiHIGH_noHe")
    ```
    See *DeVries and Holzer* (2019) for more details

!!! note
    The files, that are downloaded from a public and persistant URL in FigShare,
    were created with the code available at https://github.com/briochemc/OceanCirculations.
"""
module OCIM2

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using BSON                  # For saving circulation as BSON format
using Unitful               # for units
using Reexport
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

# OCIM2 URLs
function url(; version="CTL_He")
    url_start = "https://files.figshare.com/"
    version == "CTL_He"             ? "$(url_start)22475534/OCIM2_CTL_He.bson" :
    version == "CTL_noHe"           ? "$(url_start)23442083/OCIM2_CTL_noHe.bson" :
    version == "KiHIGH_He"          ? "$(url_start)23419721/OCIM2_KiHIGH_He.bson" :
    version == "KiHIGH_noHe"        ? "$(url_start)23418215/OCIM2_KiHIGH_noHe.bson" :
    version == "KiLOW_He"           ? "$(url_start)23416925/OCIM2_KiLOW_He.bson" :
    version == "KiLOW_noHe"         ? "$(url_start)23417861/OCIM2_KiLOW_noHe.bson" :
    version == "KvHIGH_He"          ? "$(url_start)23416583/OCIM2_KvHIGH_He.bson" :
    version == "KvHIGH_KiHIGH_noHe" ? "$(url_start)23417639/OCIM2_KvHIGH_KiHIGH_noHe.bson" :
    version == "KvHIGH_KiLOW_He"    ? "$(url_start)23416181/OCIM2_KvHIGH_KiLOW_He.bson" :
    version == "KvHIGH_KiLOW_noHe"  ? "$(url_start)23416055/OCIM2_KvHIGH_KiLOW_noHe.bson" :
    version == "KvHIGH_noHe"        ? "$(url_start)23417396/OCIM2_KvHIGH_noHe.bson" :
    OCIM2versionerror(version)
end

OCIM2versionerror(version) = error("""`$version` is not a valid OCIM2 version name.

                                   Valid versions are `CTL_He`, `CTL_noHe`, `KiHIGH_He`, `KiHIGH_noHe`, `KiLOW_He`, `KiLOW_noHe`, `KvHIGH_He`, `KvHIGH_KiHIGH_noHe`, `KvHIGH_KiLOW_He`, `KvHIGH_KiLOW_noHe`, and `KvHIGH_noHe`.

                                   See *DeVries and Holzer* (2019) for more details on OCIM2 configurations.""")

# OCIM2 Hashes
function sha(; version="CTL_He")
    version == "CTL_He"             ? "2786a56e89315feb44102629470ec565bd67ee765cb0a49364a51e6206a15490" :
    version == "CTL_noHe"           ? "bdfbdc67b224ea8a588450454feb1f440cfb3745c191ff3929ac66897a036e02" :
    version == "KiHIGH_He"          ? "bf89d657b690bbfc8055c2182c04d50f5a408b7c610b4dc8daad52fba2482bf4" :
    version == "KiHIGH_noHe"        ? "9ec3f7237e541185bbf577bc7bee8c6308f2e4a2e0457a7460e73a962902cd8c" :
    version == "KiLOW_He"           ? "2a7094a5fc0bf5aa6d5d4e0f7ad4ed89901a7431a8d039df74835d84d4259a91" :
    version == "KiLOW_noHe"         ? "a732a382365235446091f77c2cf2961af641b10662a653edec05295349151e8f" :
    version == "KvHIGH_He"          ? "bf89d657b690bbfc8055c2182c04d50f5a408b7c610b4dc8daad52fba2482bf4" :
    version == "KvHIGH_KiHIGH_noHe" ? "8ee547b2fd4ddc0b61796f0c4bc0e57622615c1e492939377b44a899f86a3122" :
    version == "KvHIGH_KiLOW_He"    ? "c207e88e88a155972f51a82e1ec216efb9bd28073464074a7b023460ff9805b6" :
    version == "KvHIGH_KiLOW_noHe"  ? "22c3de72f6d630901d929c769d1c88c09ba15264f489c34c1b0d1069fc2538c5" :
    version == "KvHIGH_noHe"        ? "419cc207a0e104427677b816c645e8fdbca83d1201204dd4dd10a95200fbc070" :
    OCIM2versionerror(version)
end

# Create registry entry for OCIM in BSON format
function register_OCIM2(; version="CTL_He")
    register(
        DataDep(
            "AIBECS-OCIM2_$version",
            """
            References:
            - $(citation())
            """,
            url(version=version),
            sha(version=version),
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load

Returns the grid, the transport matrix, and the He fluxes (in that order).

!!! tip
    To load the default OCIM2 matrix and grid, do
    ```
    julia> grd, T = OCIM2.load()
    ```
    But you can also load the other matrices by specifying which version you want, e.g.,
    ```
    julia> grd, T = OCIM2.load(version="KiHIGH_noHe")
    ```
    See *DeVries and Holzer* (2019) for more details
"""
function load(; version="CTL_He")
    register_OCIM2(version=version)
    bson_file = @datadep_str string("AIBECS-OCIM2_$version/", "OCIM2_$version.bson")
    BSON.@load bson_file grid T He3Flux He4Flux
    @info """You are about to use the OCIM2_$version model.
          If you use it for research, please cite:

          - $(citation())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "DeVries_Holzer_2019" key.)
          """
    return grid, ustrip.(T), He3Flux, He4Flux
end

citation() = "DeVries, T., & Holzer, M. (2019). Radiocarbon and helium isotope constraints on deep ocean ventilation and mantle‐³He sources. Journal of Geophysical Research: Oceans, 124, 3036–3057. https://doi.org/10.1029/2018JC014716"

versions() = ["CTL_He", "CTL_noHe", "KiHIGH_He", "KiHIGH_noHe", "KiLOW_He", "KiLOW_noHe", "KvHIGH_He", "KvHIGH_KiHIGH_noHe", "KvHIGH_KiLOW_He", "KvHIGH_KiLOW_noHe", "KvHIGH_noHe"]

end # end module

export OCIM2


text:
ENV["GKSwstype"] = "100" #magic fix tests on github headless machines with GR as plotting backend

using ThreeBodyTB
using Test

include("crystal_testing.jl")
include("test_ewald.jl")
include("sym_test.jl")
include("test_laguerre.jl")
include("test_U.jl")
include("test_scf.jl")
include("test_forces.jl")
include("test_proto.jl")
include("test_makescf.jl")
include("test_atomicproj.jl")

include("test_examples.jl")

delete!(ENV, "GKSwstype") #undo magic to avoid side effects

Nothing


text:
##### Diffusive fluxes
@inline diffusive_flux_x(i, j, k, grid, κˣ::Number, c) = κˣ * Ax(i, j, k, grid) * ∂x⁰(i, j, k, grid, c)
@inline diffusive_flux_y(i, j, k, grid, κʸ::Number, c) = κʸ * Ay(i, j, k, grid) * ∂y⁰(i, j, k, grid, c)
@inline diffusive_flux_z(i, j, k, grid, κᶻ::Number, c) = κᶻ * Az(i, j, k, grid) * ∂z⁰(i, j, k, grid, c)

##### Laplacian diffusion operator
@inline function κ∇²(i, j, k, grid, κˣ, κʸ, κᶻ, c)
    return 1/volume(i, j, k, grid) * (δx⁺(i, j, k, grid, diffusive_flux_x, κˣ, c) +
                       δy⁺(i, j, k, grid, diffusive_flux_y, κʸ, c) +
                       δz⁺(i, j, k, grid, diffusive_flux_z, κᶻ, c))
end

##### calculate the tendency by diffusion for tracer c
@kernel function calc_diffusion!(Gc, grid, κˣ, κʸ, κᶻ, c, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + grid.Hx
    jj = j + grid.Hy
    kk = k + grid.Hz
    @inbounds Gc[ii, jj, kk] = Gc[ii, jj, kk] + κ∇²(ii, jj, kk, grid, κˣ, κʸ, κᶻ, c) * ΔT
end

function nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κˣ, κʸ, κᶻ, ΔT)
    calc_diffusion_kernel! = calc_diffusion!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))

    events=[]
    for name in nut_names
        event = calc_diffusion_kernel!(Gcs[name].data, g, κˣ, κʸ, κᶻ, nutrients[name].data, ΔT, dependencies=barrier)
        push!(events, event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end
nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κ, ΔT) =
    nut_diffusion!(Gcs, arch::Architecture, g, nutrients, κ, κ, κ, ΔT)


text:
"""
    WarpAffine(σ = 0.1) <: ProjectiveTransform

A three-point affine warp calculated by randomly moving 3 corners
of an item. Similar to a random translation, shear and rotation.
"""
struct WarpAffine <: ProjectiveTransform
    σ
end

getrandstate(::WarpAffine) = abs(rand(Int))

function getprojection(
        tfm::WarpAffine,
        bounds::Bounds{2};
        randstate = getrandstate(tfm))
    T = Float32
    rng = Random.seed!(Random.MersenneTwister(), randstate)
    scale = sqrt(prod(length.(bounds.rs)))

    srcps = shuffle(rng, SVector{2, T}.(corners(bounds)))[1:3]
    offsets = rand(rng, SVector{2, T}, 3)
	offsets = map(v -> (v .* (2one(T)) .- one(T)) .* convert(T, scale * tfm.σ), offsets)

	return threepointwarpaffine(srcps, srcps .+ offsets)
end

# Adapted from
"""
    threepointwarpaffine(srcps, dstps)

Calculate an affine [`CoordinateTransformations.LinearMap`](#)
from 3 source points to 3 destination points.

Adapted from  [CoordinateTransformations.jl#30](https://github.com/JuliaGeometry/CoordinateTransformations.jl/issues/30#issuecomment-610337378).
"""
function threepointwarpaffine(
        srcps::T, dstps::T) where {V, T<:AbstractArray{<:SVector{2,V}}}
	X = vcat(hcat(dstps...), ones(1,3))'
    Y = hcat(srcps...)'
    c = (X \ Y)'
    A = SMatrix{2, 2, V}(c[:, 1:2])
    b = SVector{2, V}(c[:, 3])
    return AffineMap(A, b)
end

function corners(bounds::Bounds)
    (Tuple(x) for x in ImageTransformations.CornerIterator(CartesianIndices(bounds.rs)))
end


text:
#using JLD

function generate_theta(plan::abc_pmc_plan_type, sampler::Distribution, ss_true, epsilon::Float64; num_max_attempt = plan.num_max_attempt)
      @assert(epsilon>=0.0)
      dist_best = Inf
      local theta_best
      #summary_stats_log = abc_log_type()
      #push!(summary_stats_log.generation_starts_at, 1)
      accept_log = abc_log_type()
      reject_log = abc_log_type()
      push!(accept_log.generation_starts_at, 1)
      push!(reject_log.generation_starts_at, 1)
      attempts = 0
      all_attempts = num_max_attempt
      for a in 1:num_max_attempt
         theta_star = rand(sampler)
         if (typeof(theta_star) <: Real)   # in case return a scalar, make into array
            theta_star = fill(theta_star, length(plan.prior))  ### TODO: Univariate uniform prior returns length 1 -> need to generalize for multiple bins.
            #theta_star = [theta_star]
         end
         plan.normalize(theta_star)
         if(!plan.is_valid(theta_star)) continue end
         attempts += 1
         data_star = plan.gen_data(theta_star)
         ss_star = plan.calc_summary_stats(data_star)

         accept_prob = 0.
         local dist_star
         for d in 1:plan.num_dist_per_obs
           dist_star = plan.calc_dist(ss_true,ss_star)
         #=
         if plan.save_params
           push!(summary_stats_log.theta, theta_star)
         end
         if plan.save_summary_stats
           push!(summary_stats_log.ss, ss_star)
         end
         if plan.save_distances
           push!(summary_stats_log.dist, dist_star)
         end
         =#
           if dist_star < dist_best
              dist_best = dist_star
              theta_best = copy(theta_star)
           end
           if(dist_star < epsilon)
             accept_prob += 1.0
           end
         end
         accept_prob /= plan.num_dist_per_obs

         accept = rand() < accept_prob

         if(!accept)
              push_to_abc_log!(reject_log,plan,theta_star,ss_star,dist_star)
         else
              push_to_abc_log!(accept_log,plan,theta_star,ss_star,dist_star)
              all_attempts = a
              break
         end
      end
      if dist_best == Inf
        error("# Failed to generate any acceptable thetas.")
      end
      # println("gen_theta: d= ",dist_best, " num_valid_attempts= ",attempts, " num_all_attempts= ", all_attempts, " theta= ", theta_best)
      #return (theta_best, dist_best, attempts, summary_stats_log)
      return (theta_best, dist_best, attempts, accept_log, reject_log)
end

function generate_abc_sample(plan::abc_pmc_plan_type, sampler::Distribution, ss_true, epsilon::Real; n::Integer = 100)
  if in_parallel
   nw = nworkers()
   @assert (nw > 1)
   return generate_abc_sample_distributed_map(plan,ss_true)
  else
   return generate_abc_sample_serial(plan,ss_true)
  end
end

# Generate initial abc population from prior, aiming for d(ss,ss_true)<epsilon
function init_abc(plan::abc_pmc_plan_type, ss_true; in_parallel::Bool = plan.in_parallel)
  if in_parallel
    nw = nworkers()
    @assert (nw > 1)
    #@assert (plan.num_part > 2*nw)  # Not really required, but seems more likely to be a mistake
    #return init_abc_parallel_map(plan,ss_true)
    return init_abc_distributed_map(plan,ss_true)
  else
    return init_abc_serial(plan,ss_true)
  end
end

# run the ABC algorithm matching to summary statistics ss_true, starting from an initial population (e.g., output of previous call)
function run_abc(plan::abc_pmc_plan_type, ss_true, pop::abc_population_type; verbose::Bool = false, print_every::Integer=1, in_parallel::Bool = plan.in_parallel )
  attempts = zeros(Int64,plan.num_part)
  # Set initial epsilon tolerance based on current population
  epsilon = quantile(pop.dist,plan.init_epsilon_quantile)
  eps_diff_count = 0

  # Uncomment for generation info output to log file
  #
  if verbose
     f_log = open("generation_log.txt", "w")
  end
  #

  # Uncomment for history output at end of run to terminal
  #
  eps_arr = []
  mean_arr = []
  std_arr = []
  #

  for t in 1:plan.num_max_times
    local new_pop
    sampler = plan.make_proposal_dist(pop, plan.tau_factor)
    if plan.adaptive_quantiles
      min_quantile = t==1 ? 4.0/size(pop.theta,2) : 1.0/sqrt(size(pop.theta,2))
      epsilon = choose_epsilon_adaptive(pop, sampler, min_quantile=min_quantile)
    end
    if in_parallel
      #new_pop = update_abc_pop_parallel_darray(plan, ss_true, pop, epsilon, attempts=attempts)
      new_pop = update_abc_pop_parallel_distributed_map(plan, ss_true, pop, sampler, epsilon, attempts=attempts)
    else
      #new_pop = update_abc_pop_serial(plan, ss_true, pop, epsilon, attempts=attempts)
      new_pop = update_abc_pop_serial(plan, ss_true, pop, sampler, epsilon, attempts=attempts)
    end
    pop = new_pop
    if verbose && (t%print_every == 0)
       println("# t= ",t, " eps= ",epsilon, " med(d)= ",median(pop.dist), " attempts= ",median(attempts), " ",maximum(attempts), " reps= ", sum(pop.repeats), " ess= ",ess(pop.weights,pop.repeats)) #," mean(theta)= ",mean(pop.theta,dims=2) )#) #, " tau= ",diag(tau) ) #
       println("Mean(theta)= ", mean(pop.theta, dims=2), " Stand. Dev.(theta)= ", std(pop.theta, dims=2))

       # Uncomment for generation info output to log file
       #
       println(f_log, "# t= ",t, " eps= ",epsilon, " med(d)= ",median(pop.dist), " attempts= ",median(attempts), " ",maximum(attempts), " reps= ", sum(pop.repeats), " ess= ",ess(pop.weights,pop.repeats))
       println(f_log, "Mean(theta)= ", mean(pop.theta, dims=2), " Stand. Dev.(theta)= ", std(pop.theta, dims=2))
       flush(f_log)
       #save(string("gen-",t,".jld"), "pop_out", pop, "ss_true", ss_true)
       #

    end
    # Uncomment for history output at end of run to terminal
    #
    push!(eps_arr, epsilon)
    push!(mean_arr, mean(pop.theta,dims=2)[1])
    push!(std_arr, std(pop.theta,dims=2)[1])
    #
    #if epsilon < plan.target_epsilon  # stop once acheive goal
    if maximum(pop.dist) < plan.target_epsilon  # stop once acheive goal
       println("# Reached ",epsilon," after ", t, " generations.")
       break
    end
    if median(attempts)>0.2*plan.num_max_attempt
      println("# Halting due to ", median(attempts), " median number of valid attempts.")
      break
    end
    if sum(pop.repeats)>plan.num_part
      println("# Halting due to ", sum(pop.repeats), " repeats.")
      break
    end
    eps_old = epsilon
    if maximum(attempts)<0.75*plan.num_max_attempt
      epsilon = minimum([maximum(pop.dist),epsilon * plan.epsilon_reduction_factor])
    end
    if ((abs(eps_old-epsilon)/epsilon) < 1.0e-5)
      eps_diff_count += 1
    else
      eps_diff_count = 0
    end
    if eps_diff_count > 1
      println("# Halting due to epsilon goal remaining the same for 3 consecutive generations.")
      break
    end  end # t / num_times
    #println("mean(theta) = ",[ sum(pop.theta[i,:])/size(pop.theta,2) for i in 1:size(pop.theta,1) ])

  # Uncomment for history output at end of run to terminal
  #
  if verbose
     println("Epsilon history = ", eps_arr)
     println("Mean history = ", mean_arr)
     println("Std Dev. history = ", std_arr)
  end
  #

  # Uncomment for generation info output to log file
  #
  if verbose
     close(f_log)
  end
  #

  return pop
end

# run the ABC algorithm matching to summary statistics ss_true
function run_abc(plan::abc_pmc_plan_type, ss_true; verbose::Bool = false, print_every::Integer=1, in_parallel::Bool =  plan.in_parallel )                         # Initialize population, drawing from prior
  pop::abc_population_type = init_abc(plan,ss_true, in_parallel=in_parallel)
  if verbose
      println("pop_init: ",pop)
  end
  run_abc(plan, ss_true, pop; verbose=verbose, print_every=print_every, in_parallel=in_parallel )
end

function choose_epsilon_adaptive(pop::abc_population_type, sampler::Distribution; min_quantile::Real = 1.0/sqrt(size(pop.theta,2)) )
  sampler_logpdf = logpdf(sampler, pop.theta)
  target_quantile_this_itteration = min(1.0, exp(-maximum(sampler_logpdf .- pop.logpdf)) )
  if target_quantile_this_itteration > 1.0
     target_quantile_this_itteration = 1.0
  end
  optimal_target_quantile = target_quantile_this_itteration
  if target_quantile_this_itteration < min_quantile
     target_quantile_this_itteration = min_quantile
  end
  epsilon = quantile(pop.dist,target_quantile_this_itteration)
  println("# Estimated target quantile is ", optimal_target_quantile, " using ",target_quantile_this_itteration, " resulting in epsilon = ", epsilon)
  return epsilon
end


text:
module Scale

using Color
using Compose
using DataFrames
using Gadfly

import Gadfly.element_aesthetics, Gadfly.isconcrete, Gadfly.concrete_length,
       Gadfly.nonzero_length

include("color.jl")

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#   datas: Zero or more data objects. (Yes, I know "datas" is not a real word.)
#
# Returns:
#   nothing
#
function apply_scales(scales,
                      aess::Vector{Gadfly.Aesthetics},
                      datas::Gadfly.Data...)
    for scale in scales
        apply_scale(scale, aess, datas...)
    end
end

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   datas: Zero or more data objects.
#
# Returns:
#   A vector of Aesthetics of the same length as datas containing scaled data.
#
function apply_scales(scales, datas::Gadfly.Data...)
    aess = [Gadfly.Aesthetics() for _ in datas]
    apply_scales(scales, aess, datas...)
    aess
end

# Transformations on continuous scales
type ContinuousScaleTransform
    f::Function     # transform function
    finv::Function  # f's inverse
    label::Function # produce a string given some value f(x)
end

const identity_transform =
    ContinuousScaleTransform(identity, identity, Gadfly.fmt_float)
const log10_transform =
    ContinuousScaleTransform(log10, x -> 10^x,
                             x -> @sprintf("10%s", Gadfly.fmt_float(x)))
const log2_transform =
    ContinuousScaleTransform(log10, x -> 2^x,
                             x -> @sprintf("2%s", Gadfly.fmt_float(x)))
const ln_transform =
    ContinuousScaleTransform(log, exp,
                             x -> @sprintf("e%s", Gadfly.fmt_float(x)))
const asinh_transform =
    ContinuousScaleTransform(asinh, sinh, x -> Gadfly.fmt_float(sinh(x)))
const sqrt_transform =
    ContinuousScaleTransform(sqrt, x -> x^2,
                             x -> @sprintf("√%s", Gadfly.fmt_float(x)))

# Continuous scale maps data on a continuous scale simple by calling
# `convert(Float64, ...)`.
type ContinuousScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
    trans::ContinuousScaleTransform
end

const x_vars = [:x, :x_min, :x_max]
const y_vars = [:y, :y_min, :y_max]

# Commonly used scales.
const x_continuous = ContinuousScale(x_vars, identity_transform)
const y_continuous = ContinuousScale(y_vars, identity_transform)
const x_log10      = ContinuousScale(x_vars, log10_transform)
const y_log10      = ContinuousScale(y_vars, log10_transform)
const x_log2       = ContinuousScale(x_vars, log2_transform)
const y_log2       = ContinuousScale(y_vars, log2_transform)
const x_log        = ContinuousScale(x_vars, ln_transform)
const y_log        = ContinuousScale(y_vars, ln_transform)
const x_asinh      = ContinuousScale(x_vars, asinh_transform)
const y_asinh      = ContinuousScale(y_vars, asinh_transform)
const x_sqrt       = ContinuousScale(x_vars, sqrt_transform)
const y_sqrt       = ContinuousScale(y_vars, sqrt_transform)

function element_aesthetics(scale::ContinuousScale)
    return scale.vars
end

# Apply a continuous scale.
#
# Args:
#   scale: A continuos scale.
#   datas: Zero or more data objects.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#
# Return:
#   nothing
#
function apply_scale(scale::ContinuousScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            label_var = symbol(@sprintf("%s_label", string(var)))
            if getfield(data, var) === nothing
                continue
            end

            ds = Array(Float64, concrete_length(getfield(data, var)))
            i = 1
            for d in getfield(data, var)
                if isconcrete(d)
                    ds[i] = scale.trans.f(convert(Float64, d))
                    i += 1
                end
            end

            setfield(aes, var, ds)
            if contains(Set(names(aes)...), label_var)
                setfield(aes, label_var, scale.trans.label)
            end
        end
    end
end

discretize(values::Vector) = PooledDataArray(values)
discretize(values::DataArray) = PooledDataArray(values)
discretize(values::PooledDataArray) = values

type DiscreteScaleTransform
    f::Function
end

type DiscreteScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
end

element_aesthetics(scale::DiscreteScale) = scale.vars

const x_discrete = DiscreteScale(x_vars)
const y_discrete = DiscreteScale(y_vars)

function apply_scale(scale::DiscreteScale, aess::Vector{Gadfly.Aesthetics},
                     datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            label_var = symbol(@sprintf("%s_label", string(var)))

            if getfield(data, var) === nothing
                continue
            end

            disc_data = discretize(getfield(data, var))
            setfield(aes, var, Int64[r for r in disc_data.refs])

            # The labeler for discrete scales is a closure over the discretized data.
            function labeler(i)
                if 0 < i <= length(levels(disc_data))
                    string(levels(disc_data)[int(i)])
                else
                    ""
                end
            end

            if contains(Set(names(aes)...), label_var)
                setfield(aes, label_var, labeler)
            end
        end
    end
end

type DiscreteColorScale <: Gadfly.ScaleElement
    f::Function # A function f(n) that produces a vector of n colors.
end

function element_aesthetics(scale::DiscreteColorScale)
    [:color]
end

# Common discrete color scales
const color_hue = DiscreteColorScale(
    h -> convert(Vector{ColorValue},
        distinguishable_colors(h, c -> deuteranopic(c, 0.8),
                               LCHab(70, 60, 240),
                               Float64[65, 70, 75, 80, 85],
                               Float64[0, 50, 60],
                               Float64[h for h in 0:30:360])))

function apply_scale(scale::DiscreteColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        if data.color === nothing
            continue
        end
        ds = discretize(data.color)
        colors = convert(Vector{ColorValue}, scale.f(length(levels(ds))))
        colorvals = Array(ColorValue, nonzero_length(ds.refs))
        i = 1
        for k in ds.refs
            if k != 0
                colorvals[i] = colors[k]
                i += 1
            end
        end

        colored_ds = PooledDataArray(colorvals, colors)
        aes.color = colored_ds

        color_map = {color => string(label)
                     for (label, color) in zip(levels(ds), colors)}
        aes.color_label = c -> color_map[c]
        aes.color_key_colors = colors
    end
end

type ContinuousColorScale <: Gadfly.ScaleElement
    # A function of the form f(p) where 0 <= p <= 1, that returns a color.
    f::Function
end

element_aesthetics(::ContinuousColorScale) = [:color]

# Common continuous color scales
# TODO: find a good color combo
const color_gradient = ContinuousColorScale(
        lab_gradient(LCHab(20, 44, 262), LCHab(100, 44, 262)))

function apply_scale(scale::ContinuousColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    cmin = Inf
    cmax = -Inf
    for data in datas
        if data.color === nothing
            continue
        end

        for c in data.color
            if c === NA
                continue
            end

            c = convert(Float64, c)
            if c < cmin
                cmin = c
            end

            if c > cmax
                cmax = c
            end
        end
    end

    # It's a little weird to be doing this here. If a more elegant solution
    # arises, it's worth reorganizing.

    if cmin == Inf || cmax == -Inf
        return nothing
    end

    ticks = Gadfly.optimize_ticks(cmin, cmax)
    if ticks[1] == 0 && cmin >= 1
        ticks[1] = 1
    end

    cmin = ticks[1]
    cmax = ticks[end]
    cspan = cmax != cmin ? cmax - cmin : 1

    for (aes, data) in zip(aess, datas)
        if data.color === nothing
            continue
        end

        nas = [c === NA for c in data.color]
        cs = Array(ColorValue, length(data.color))
        for (i, c) in enumerate(data.color)
            if c === NA
                continue
            end
            cs[i] = scale.f((convert(Float64, c) - cmin) / cspan)
        end

        aes.color = PooledDataArray(cs, nas)

        color_labels = Dict{ColorValue, String}()
        for tick in ticks
            r = (tick - cmin) / cspan
            color_labels[scale.f(r)] = Gadfly.fmt_float(tick)
        end

        # Add a gradient of steps between labeled colors.
        num_steps = 1
        for (i, j) in zip(ticks, ticks[2:end])
            span = j - i
            for step in 1:num_steps
                k = i + span * (step / (1 + num_steps))
                r = (k - cmin) / cspan
                color_labels[scale.f(r)] = ""
            end
        end

        aes.color_label = c -> color_labels[c]
        aes.color_key_colors = [k for k in keys(color_labels)]
        sort!(aes.color_key_colors, order=Sort.Reverse)
        aes.color_key_continuous = true
    end
end

# Label scale is always discrete, hence we call it 'label' rather
# 'label_discrete'.
type LabelScale <: Gadfly.ScaleElement
end

function apply_scale(scale::LabelScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        if data.label === nothing
            continue
        end

        aes.label = discretize(data.label)
    end
end

element_aesthetics(::LabelScale) = [:scale]

const label = LabelScale()

end # module Scale


text:
function DolittleLU(A::Array; pivot=true)

  #A = the matrix of coefficients 'A' must be squared
  #B = The Vector 'B' right hand side vector Vertical
    m,n =size(A)

    if m != n
        error("Matrix must be squared")
    end

    # System of equations
    U = A
    P=Matrix{Float64}(I,m,m)
    p=zeros(m)
    L=zeros(m,m)
    #Forward Elimination
    for j=1:m-1
        if pivot==true
        k=argmax(abs.(U[j:m,j]))   #Find the Max value in the colum k to change the order
        Ind=j+k-1                  #Indicates de index of the row to change
        U[[j,Ind],:]=U[[Ind,j],:]  #Change the order of rows
        P[[j,Ind],:]=P[[Ind,j],:]
        L[[j,Ind],:]=L[[Ind,j],:]
        end
       for i=j+1:m

            f=U[i,j]/U[j,j]
            L[i,j]=f
            U[i,:]=U[i,:].-f.*U[j,:]
        end
    end
    for i=1:m
        L[i,i]=1
        p[i]=argmax(P[i,:])
    end
    U=UpperTriangular(U)
    L=LowerTriangular(L)

    return L,U, P, p
end


text:
################################################################################
#
#  Map/AbGrp.jl : Types for maps with domains of type GrpAbFinGen
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015 - 2018 Tommy Hofmann, Carlo Sircana
#
################################################################################

################################################################################
#
#  Morphisms between finitely generated abelian groups
#
################################################################################

mutable struct GrpAbFinGenMap <: Map{GrpAbFinGen, GrpAbFinGen,
                                     HeckeMap, GrpAbFinGenMap}
  header::MapHeader{GrpAbFinGen, GrpAbFinGen}

  map::fmpz_mat
  imap::fmpz_mat
  im::GrpAbFinGen  # if set
  ke::GrpAbFinGen  # if set

  function GrpAbFinGenMap(G::GrpAbFinGen)
    r = new()
    r.header = MapHeader(G, G)
    r.map = identity_matrix(FlintZZ, ngens(G))
    r.imap = identity_matrix(FlintZZ, ngens(G))
    return r
  end

  function GrpAbFinGenMap(From::GrpAbFinGen, To::GrpAbFinGen, M::fmpz_mat)
    r = new()
    r.header = MapHeader(From, To)
    r.map = M
    return r
  end

  function GrpAbFinGenMap(From::GrpAbFinGen, To::GrpAbFinGen, M::fmpz_mat, Mi::fmpz_mat)
    r = new()
    r.header = MapHeader(From, To)
    r.map = M
    r.imap = Mi
    return r
  end
end

function image(f::Map(GrpAbFinGenMap), a::GrpAbFinGenElem)
  return GrpAbFinGenElem(codomain(f), a.coeff * f.map)
end

function preimage(f::Map(GrpAbFinGenMap), a::GrpAbFinGenElem)
  return GrpAbFinGenElem(domain(f), a.coeff * f.imap)
end

################################################################################
#
#  Isomorphisms from abelian groups onto multliplicative groups of residue
#  rings of maximal orders
#
################################################################################

mutable struct GrpAbFinGenToNfOrdQuoRingMultMap <: Map{GrpAbFinGen, NfOrdQuoRing, HeckeMap, GrpAbFinGenToNfOrdQuoRingMultMap}
  header::MapHeader{GrpAbFinGen, NfOrdQuoRing}
  generators::Vector{NfOrdQuoRingElem}
  discrete_logarithm::Function

  # The multiplicative group, tame part
  tame::Dict{NfOrdIdl,Tuple{NfOrdElem,fmpz,Function}}

  # Multiplicative group, wild part
  wild::Dict{NfOrdIdl,Tuple{Array{NfOrdElem,1},Array{fmpz,1},Function}}

  function GrpAbFinGenToNfOrdQuoRingMultMap(Q::NfOrdQuoRing,
                              generators::Vector{NfOrdQuoRingElem},
                              snf_structure::Vector{fmpz},
                              disc_log::Function)
    @assert length(generators) == length(snf_structure)
    @hassert :NfOrdQuoRing 1 all(g->parent(g)==Q,generators)

    G = DiagonalGroup(snf_structure)
    @assert isa(G,GrpAbFinGen)
    @assert issnf(G)

    function _image(a::GrpAbFinGenElem)
      @assert parent(a) == G
      y = one(Q)
      for i in 1:length(generators)
        a[i] == 0 && continue
        y *= generators[i]^a[i]
      end
      return y
    end

    function _preimage(a::NfOrdQuoRingElem)
      @assert parent(a) == Q
      return G(disc_log(a))
    end

    z = new()
    z.header = MapHeader(G, Q, _image, _preimage)
    z.generators = generators
    z.discrete_logarithm = disc_log
    return z
  end

  function GrpAbFinGenToNfOrdQuoRingMultMap(G::GrpAbFinGen, Q::NfOrdQuoRing, exp::Function, disc_log::Function)
    z = new()
    z.header = MapHeader(G, Q, exp, disc_log)
    return z
  end
end

################################################################################
#
#  Morphisms from finite abelian groups onto units of orders
#
################################################################################

mutable struct AbToNfOrdMultGrp <: Map{GrpAbFinGen, NfOrd, SetMap, AbToNfOrdMultGrp}
  domain::GrpAbFinGen
  codomain::NfOrd
  generator::NfOrdElem

  function AbToNfOrdMultGrp(O::NfOrd, order::Int, generator::NfOrdElem)
    G = DiagonalGroup([order])
    z = new()
    z.domain = G
    z.codomain = O
    z.generator = generator
    return z
  end
end

function (f::AbToNfOrdMultGrp)(a::GrpAbFinGenElem)
  return f.generator^a[1]
end


text:
import JSON, JSON3
using Test

function learn_JSON()
    s = "{\"a_number\" : 5.0, \"an_array\" : [\"string\", 9]}"
    j = JSON.parse(s)
    JSON.print(stdout, j, 4)
    # out = IOBuffer()
    # JSON.print(out, j)
    # @show String(take!(out))
    @show JSON.json(j)
end

learn_JSON()

nothing


text:
using DiffEqBase, Markdown

#DAE Example Problems
export prob_dae_resrob

### DAE Problems

f = function (r, yp, y, p, tres)
    r[1]  = -0.04*y[1] + 1.0e4*y[2]*y[3]
    r[2]  = -r[1] - 3.0e7*y[2]*y[2] - yp[2]
    r[1] -=  yp[1]
    r[3]  =  y[1] + y[2] + y[3] - 1.0
end
u0 = [1.0, 0, 0]
du0 = [-0.04, 0.04, 0.0]

@doc doc"""
The Robertson biochemical reactions in DAE form

```math
\frac{dy₁}{dt} = -k₁y₁+k₃y₂y₃
```
```math
\frac{dy₂}{dt} =  k₁y₁-k₂y₂^2-k₃y₂y₃
```
```math
1 = y₁ + y₂ + y₃
```
where ``k₁=0.04``, ``k₂=3\times10^7``, ``k₃=10^4``. For details, see:
Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Problems Page 129
Usually solved on ``[0,1e11]``
"""
prob_dae_resrob = DAEProblem(f,du0,u0,(0.0,100000.0))


text:
using AccurateArithmetic
using Test
using Random: seed!

using Base:IEEEFloat

seed!(1602)
include("bigfloats.jl")
include("randfloats.jl")

const nrands = 16
const flbig = randbf(nrands)
const fl64s = map(Float64, flbig)
const fl32s = map(Float32, flbig)

hi64, lo64 = Floats2(Float64, flbig[1])
hi32, lo32 = Floats2(Float32, flbig[2])

@test add_(lo64, hi64) == (hi64, lo64)
@test add_hilo_(hi32, lo32) == (hi32, lo32)

@test sub_(lo64, hi64) == Floats2(Float64, BigFloat(lo64) - BigFloat(hi64))
@test sub_hilo_(hi32, lo32) == Floats2(Float32, BigFloat(hi32) - BigFloat(lo32))

hi64, lo64 = Floats2(Float64, flbig[3])
hi32, lo32 = Floats2(Float32, flbig[4])

@test powr2_(hi64) == Floats2(Float64, BigFloat(hi64)^2)

@test mul_(hi64, lo64) == Floats2(Float64, BigFloat(hi64)*BigFloat(lo64))
@test dvi_(hi32, lo32) == Floats2(Float32, BigFloat(hi32)/BigFloat(lo32))

@test root2_(hi64) == Floats2(Float64, sqrt(BigFloat(hi64)))
@test inv_(hi32)  == Floats2(Float32, inv(BigFloat(hi32)))

hi64, md64, lo64 = Floats3(Float64, flbig[5])
hi32, md32, lo32 = Floats3(Float32, flbig[5])

@test mul_(hi64, lo64, md64) == Floats3(Float64, BigFloat(hi64) * BigFloat(lo64) * BigFloat(md64))
@test mul_(hi64, lo64, md64) == Floats3(Float64, BigFloat(hi64) * BigFloat(lo64) * BigFloat(md64))

@test fma_(hi64, lo64, md64) == Floats3(Float64, BigFloat(hi64) * BigFloat(lo64) + BigFloat(md64))
@test fms_(md32, lo32, hi32) == Floats3(Float32, BigFloat(md32) * BigFloat(lo32) - BigFloat(hi32))

flbig64s = map(BigFloat, fl64s)
@test sum_(fl64s) == Floats2(Float64, sum(flbig64s))


text:
# This file was generated, do not modify it. # hide
RANDOM_SEED = 42;

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using DataFrames, CSV, HTTP, Random, LinearAlgebra, Statistics, Optim, DataFramesMeta, GLM

include("create_grids.jl")

function PS6()
#Question 1

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)

df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#question 2

θ=glm(@formula(Y ~ Odometer * Odometer^2 * RouteUsage * RouteUsage^2 * Branded * time * time^2), df_long, Binomial(), LogitLink())

#Question 3

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df2 = CSV.read(HTTP.get(url).body)
Y = Matrix(df2[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df2[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df2[:,:RouteUsage])
B = Vector(df2[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df2[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df2[:,:Zst])
zval,zbin,xval,xbin,xtran = create_grids()

log_df=DataFrame()
log_df.Odometer=kron(ones(zbin),xval) #20301
log_df.RouteUsage =kron(ones(xbin),zval) #20301
log_df.Branded=zeros((size(log_df,1)))
log_df.time=zeros(size(log_df,1))

function values(din,Zstate,Xstate,xtran,zbin,xbin,xval)
    FV=zeros(size(xtran,1),2,T+1)
    for t=2:T
        for b=0:1
            din.time[t]=t
            din.Branded[t]=b
            p0 = predict(θ, din)

            FV[:, b+1, t] = - .9 .*log.(p0)
        end
    end
    FVT1 = zeros(size(df,1), T)
    for i=1:size(df,1)
        row0 = (Zstate[i]-1)*xbin+1
        for t=1:T
            row1  = Xstate[i,t] + (Zstate[i]-1)*xbin
            FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,df.Branded[i]+1,t+1]
        end
    end

    return FVT1'[:]
end

fvt1=values(log_df,Zstate,Xstate,xtran,zbin,xbin,xval)
df_long = @transform(df_long,fv = fvt1)

theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),df_long, Binomial(), LogitLink(),offset=df_long.fv)

end

@time PS6()


text:
module Day8

using Advent: save_input

struct Instruction
    op::Symbol
    arg::Int
    Instruction(op::Union{AbstractString,Symbol}, arg::Int) = new(Symbol(op), arg)
end

function parse_instructions(filename)
    insts = Instruction[]
    open(filename) do io
        while !eof(io)
            op, arg = readuntil(io, ' '), parse(Int, readline(io))
            push!(insts, Instruction(op, arg))
        end
    end
    insts
end

function simulate!(insts)
    pos = 1
    acc = 0
    seen = BitSet()
    while checkbounds(Bool, insts, pos) && !(pos in seen)
        push!(seen, pos)
        inst = insts[pos]
        inst.op == :acc && (acc += inst.arg; pos += 1)
        inst.op == :jmp && (pos += inst.arg)
        inst.op == :nop && (pos += 1)
    end
    (pos == length(insts) + 1, acc)
end

function solve_a(filename)
    _, res = filename |> parse_instructions |> simulate!
    res
end

function solve_b(filename)
    insts = parse_instructions(filename)
    for (i, inst) in enumerate(insts)
        if inst.op in (:jmp, :nop)
            save = insts[i]
            insts[i] = Instruction((insts[i].op == :jmp ? :nop : :jmp), insts[i].arg)
            success, res = simulate!(insts)
            success && return res
            insts[i] = save
        end
    end
    -1
end

solve(filename = save_input(2020, 8)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
"""
    RiemannMatrix

Type containing information about the Riemann Matrix, for use in computing the Riemann theta function.
"""
struct RiemannMatrix
    τ::Array{<:Number} # matrix in Siegel upper-half space
    g::Integer # genus
    X::Array{<:Real} # real part of τ
    Y::Array{<:Real} # imaginary part of τ
    T::Array{<:Real} # upper triangular matrix in Cholesky factorization of Y
    ellipsoid::Array{} # points to sum over for computing theta functions

    """
        RiemannMatrix(τ, siegel=false, ϵ=1.0e-12, nderivs=4)

    Construct a RiemannMatrix type using a matrix τ in the Siegel upper-half space, for use in computing the Riemann theta function.

    # Arguments
    - `τ::Array{<:Number}`: 2-dimensional array of complex numbers
    - `siegel::Bool=false`: do a siegel transformation on τ if true, and use the input matrix otherwise
    - `ϵ::Real=1.0e-12`: absolute error in value of theta function or its derivatives
    - `nderivs::Integer=4`: highest order of the derivatives of the theta function

    # Examples
    ```julia
    julia> RiemannMatrix([1+im -1; -1 1+im], siegel=true, nderivs=2)
    ```
    """
    function RiemannMatrix(τ::Array{<:Number}; siegel::Bool=false, ϵ::Real=1.0e-12, nderivs::Integer=4)
        τ = (siegel ? siegel_transform(τ)[2] : τ);
        X = real(τ);
        Y = imag(τ);
        g = size(Y)[1];
        T = (g > 1 ? convert(Array{Float64}, cholesky(Symmetric(Y)).U) : sqrt.(Y));
        ρ = svp(T)[2]*sqrt(π); # length of shortest vector in lattice generated by √π*T
        radius = [radius_ellipsoid(ϵ, T, Y, ρ, i) for i=0:nderivs]; # radius of ellipsoids for computing theta function, and up to nderivs derivatives.
        ellipsoid = [ellipsoid_uniform(T, r) for r in radius];
        new(τ, g, X, Y, T, ellipsoid)
    end
end

"""
    random_siegel(g)

Sample a random matrix in the Siegel upper half space with genus g.

# Arguments
- `g::Integer`: genus.

# Examples
```julia
julia> random_siegel(5)
```
"""
function random_siegel(g::Integer)
    Mx = 2*rand(g, g) - ones(g,g); # entries are random between [-1,1)
    My = 2*rand(g, g) - ones(g,g); # entries are random between [-1,1)
    τ = Symmetric((1/2)*(transpose(Mx) + Mx)) + Symmetric(transpose(My)*My)*im; # symmetric real part and psd imaginary part
    return convert(Array, τ);
end


text:
struct stSubMX{T} <: subArray{T}
    nID::Int64
    ranges::Vector{Tuple{UnitRange{Int64},UnitRange{Int64}}}
    MXfun::Vector{T} # [R^k][i,j] # T<:AbstractMatrix{<:AbstractVector{<:Real}}
end
struct stSubV{T} <: subArray{T}
    nID::Int64
    Vfun::T # T<:AbstractVector{<:Real}
end
stSubV() = stSubV(Vector{Float64}(undef,0),Vector{Vector{SubV}}(undef,0))

struct Result{d,lddepT<:AbstractLDDEProblem{d},mT,submxT,stsubmxT,subvT,stsubvT,tT,AavgT} <: AbstractResult{d}
    LDDEP::lddepT
    method::mT
    itoisometrymethod::ItoIsometryMethod
    subMXs::Vector{Vector{submxT}} # [[A(t1),A(t2),...],[B1(t1),B1(t2),...],[B2(t1),B2(t2),...],...]
    stsubMXs::Vector{Vector{stsubmxT}} # [[A(t1),A(t2),...],[B1(t1),B1(t2),...],[B2(t1),B2(t2),...],...]
    #
    subVs::Vector{subvT} # [c(t1),c(t2),...]
    stsubVs::Vector{Vector{stsubvT}} # [[σ1(t1),σ1(t2),...], [σ2(t1),σ2(t2),...],...]
    #
    A_avgs::AavgT
    ts::Vector{tT} # [0,t1,t2,...]
    n_steps::Int64 # number of time steps
    n::Int64 # Large discretisation matrix size
    calculate_additive::Bool
end

function Result(LDDEP::LDDEProblem{d,AT,BT, cT}, method::DiscretizationMethod{fT}, DiscretizationLength::Real; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false,im::ItoIsometryMethod{K} = Trapezoidal(20,method)) where {d,AT,BT,cT,N,fT,K}
    # DiscretizationLength discretisated time interval length
    # n_steps: how many mapping matrix to calculate
    ts =  collect(fT,take(n_steps + 1, iterated(x -> method.Δt + x, zero(method.Δt))))
    n = (rOfDelay(DiscretizationLength, method) + 1) * d
    A_avgs = calculate_Aavgs(LDDEP.A, ts, method.Δt)

    subMXs = [Vector{SubMX{eltype(A_avgs)}}(undef, n_steps) for i in 1:(length(LDDEP.Bs) + 1)] # []
    stsubMXs = [Vector{stSubMX{SizedArray{Tuple{d,d},SArray{Tuple{K},eltype(eltype(A_avgs)),1,K},2,2}}}(undef,n_steps) for i in 1:(length(LDDEP.αs)+length(LDDEP.βs))]

    if calculate_additive
        subVs = Vector{SubV{SVector{d,eltype(eltype(A_avgs))}}}(undef, n_steps) # []
        stsubVs = Vector{Vector{stSubV{SizedArray{Tuple{d},SArray{Tuple{K},Float64,1,K},1,1}}}}(undef, LDDEP.w);
        # stSubV(LDDEP.dnoise.delayMX.τ.τ.ws, [Vector{SubV}(undef,n_steps) for i in LDDEP.dnoise.delayMX.τ.τ.ws]) # []
    else
        subVs = Vector{SubV{SVector{d,eltype(eltype(A_avgs))}}}(undef, 0)
        stsubVs = Vector{Vector{stSubV{SizedArray{Tuple{d},SArray{Tuple{K},Float64,1,K},1,1}}}}(undef, 0)
    end
    Result(LDDEP, method, im, subMXs, stsubMXs, subVs, stsubVs, A_avgs, ts, n_steps, n,calculate_additive)
end

struct stDiscreteMapping{tT,dmxT,stmxT,dvT,stvT}
    ts::Vector{tT}
    detMXs::Vector{dmxT} # F [time]
    detVs::Vector{dvT} # f [time]
    stMXs::Vector{Vector{stmxT}} # G [time][noiseID]
    stVs::Vector{Vector{stvT}} # g [time][noiseID]
end

struct DiscreteMapping_M1{tT,mxT,vT}
    ts::Vector{tT}
    M1_MXs::Vector{mxT}
    M1_Vs::Vector{vT}
end

struct DiscreteMapping_M2{tT,mx1T,mx2T,mx12T,v1T,v2T}#,mx12T,v2T
    ts::Vector{tT}
    M1_MXs::Vector{mx1T} # F [time]
    M1_Vs::Vector{v1T} # c_1 [time]
    M2_MXs::Vector{mx2T} # H [time]
    M1toM2_MXs::Vector{mx12T} # C [time]
    M2_Vs::Vector{v2T} # c_2 [time]
end


text:
include("graph.jl")
include("DFS.jl")
M = [0 1 0 1 0 0;
     0 0 0 0 1 0;
     0 0 0 0 1 1;
     0 1 0 0 0 0;
     0 0 0 1 0 0;
     0 0 0 0 0 1]
G = nonEmptyGraph(M)
DFS!(G)
printGraph(G)


text:
# Types
include("QuadForm/Types.jl")
include("QuadForm/Quad/Types.jl")
include("QuadForm/Herm/Types.jl")

# Binary quadratic forms
include("QuadForm/QuadBin.jl")

# Functionality for spaces
include("QuadForm/Spaces.jl")

# Functionality for lattices
include("QuadForm/Lattices.jl")

# Misc things that need to be sorted
include("QuadForm/Misc.jl")

# Functionality for pseudo matrices
include("QuadForm/PseudoMatrices.jl")

# Functionality for IO with Hecke/Magma
include("QuadForm/IO.jl")

# Quadratic
include("QuadForm/Quad/Spaces.jl")
include("QuadForm/Quad/Lattices.jl")
include("QuadForm/Quad/ZLattices.jl")
include("QuadForm/Quad/Genus.jl")
include("QuadForm/Quad/GenusRep.jl")
include("QuadForm/Quad/Legacy.jl")

# Hermitian
include("QuadForm/Herm/Spaces.jl")
include("QuadForm/Herm/Lattices.jl")
include("QuadForm/Herm/Genus.jl")
include("QuadForm/Herm/GenusRep.jl")
include("QuadForm/Herm/Mass.jl")
include("QuadForm/Herm/Legacy.jl")
include("QuadForm/Herm/LocallyIsometricSublattice.jl")

include("QuadForm/Morphism.jl")
include("QuadForm/Database.jl")
include("QuadForm/Enumeration.jl")
include("QuadForm/LineOrbits.jl")
include("QuadForm/MassQuad.jl")

# Torsion
include("QuadForm/Torsion.jl")