text:
export Connection, readrequestline!, readheaders!, readbody!

const DEFAULT_MAX_HEADER = UInt64(80 * 1024)
const DEFAULT_MAX_TRAILER = UInt64(20 * 1024)
const DEFAULT_MAX_BODY = typemax(UInt64)
const DEFAULT_MAX_URI = UInt64(8000)
const DEFAULT_MAX_METHOD = UInt64(16)
const DEFAULT_MAX_VERSION = UInt64(10)

const OCTETS = Set(UInt8(0):UInt8(255))
const CHAR = Set(UInt8(0):UInt8(127))
const VCHAR = Set(UInt8(0x21):UInt8(0x7E))  # visible characters
const UCALPHA = Set(UInt8('A'):UInt8('Z'))
const LCALPHA = Set(UInt8('a'):UInt8('z'))
const ALPHA = union(UCALPHA, LCALPHA)
const DIGIT = Set(UInt8('0'):UInt8('9'))
const CTL = union(Set(UInt8(0):UInt8(39)), Set((UInt8(127))))

const CR = UInt8('\r')
const LF = UInt8('\n')
const SP = UInt8(' ')
const HT = UInt8('\t')
const QUOTE = UInt8('"')
const BACKSLASH = UInt8('\\')   # used for escapes in strings and comments
const LEFTPAREN = UInt8('(')           # used for beginning of comments in some headers
const RIGHTPAREN = UInt8(')')           # used for end of comments in some headers
const COLON = UInt8(':')
const COMMA = UInt8(',')
const EQUAL = UInt8('=')
const LWS = Set([SP,HT])
const TEXT = union(setdiff(OCTETS, CTL), LWS)
const HEX = union(DIGIT, Set(Vector{UInt8}("abcdefABCDEF")))
const SEPARATOR = Set(Vector{UInt8}("()<>@,;:\\\"/[]?={} \t"))
# const TOKEN = setdiff(setdiff(CHAR, SEPARATOR), CTL)
const TOKEN = union(Set(Vector{UInt8}("!#\$%&'*+-.^_`|~")), ALPHA, DIGIT)

const URI_GENDELIMS = Set(Vector{UInt8}(":/?#[]@"))
const URI_SUBDELIMS = Set(Vector{UInt8}("!\$&'()*+,;="))
const URI_RESERVED = union(URI_GENDELIMS, URI_SUBDELIMS)
const URI_UNRESERVED = union(ALPHA, DIGIT, Set(Vector{UInt8}("-._~")))
const URI_CHARS = union(URI_RESERVED, URI_UNRESERVED)

const CRLF = [UInt8('\r'), UInt8('\n')]

const V1_0 = Vector{UInt8}("1.0")
const V1_1 = Vector{UInt8}("1.1")
const V2_0 = Vector{UInt8}("2.0")

const GENERAL_HEADER_FIELDS = [
    "cache-control",
    "connection",
    "date",
    "pragma",
    "trailer",
    "transfer-encoding",
    "upgrade",
    "via",
    "warning"
]

const COMMENTABLE_HEADER_FIELDS = Set([
    "user-agent",
    "server",
    "via"
])

mutable struct Connection
    socket::IO
    method::String
    uri::String
    httpversion::String
    rawheaders::Vector{Pair{String,String}}
    headers::Headers
    bodystream::IO
    rawtrailers::Vector{Pair{String,String}}
    trailers::Headers

    requestlineread::Bool
    headersread::Bool
    bodystartedread::Bool
    bodyfullyread::Bool
    responsesent::Bool
end

function Connection(socket::IO)
    Connection(socket,
        "",
        "",
        "",
        Vector{Pair{String,String}}(),
        Headers(),
        IOBuffer(""),
        Vector{Pair{String,String}}(),
        Headers(),
        false,
        false,
        false,
        false,
        false)
end

function readrequestline!(connection::Connection; maxuri=DEFAULT_MAX_URI, maxmethod=DEFAULT_MAX_METHOD, maxversion=DEFAULT_MAX_VERSION, strict=false)
    if !isequal(false, connection.requestlineread)
        throw(HTTPError(500; message="Cannot read http connection request line when it has already been read", shouldclose=true))
    end

    connection.requestlineread = true
    method = Vector{UInt8}()
    uri = Vector{UInt8}()
    version = Vector{UInt8}()
    socket = connection.socket
    foundend = false

    # parse the method
    ind = 1
    while isopen(socket) && !eof(socket)
        b = read(socket, UInt8)
        if b == CR && ind == 1 && !strict
            # https://tools.ietf.org/html/rfc7230#section-3.5 (robustness)
            # this is technically wrong, but expected bad behavior
            ind += 1
            b = read(socket, UInt8)
            if b != LF
                throw(HTTPError(400; message="Invalid character in http method", shouldclose=true))
            end
            ind += 1
        end
        if b == SP
            # method is done
            if length(method) == 0
                throw(HTTPError(400; message="Missing http method", shouldclose=true))
            end
            break
        end
        if ind > maxmethod
            throw(HTTPError(405; message="Unwilling to accept methods longer than $maxheader bytes", shouldclose=true))
        end
        if !in(b, TOKEN)
            throw(HTTPError(400; message="Invalid character in http method", shouldclose=true))
        end
        push!(method, b)
        ind += 1
    end

    ind = 1
    while isopen(socket)
        b = read(socket, UInt8)
        if b == SP
            # method is done
            break
        end
        if ind > maxuri
            throw(HTTPError(414; message="Unwilling to accept a uri longer than $maxuri bytes", shouldclose=true))
        end
        if !in(b, URI_CHARS)
            throw(HTTPError(400; message="Invalid character in uri", shouldclose=true))
        end
        push!(uri, b)
        ind += 1
    end

    ind = 1
    founddot = false
    while isopen(socket)
        b = read(socket, UInt8)
        if (ind < 6)
            if ind == 1 && b != UInt8('H')
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            elseif (ind == 2 || ind == 3) && b != UInt8('T')
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            elseif ind == 4 && b != UInt8('P')
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            elseif ind == 5 && b != UInt8('/')
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            end
        elseif b == UInt8('.')
            if founddot
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            end

            founddot = true
            push!(version, b)
        elseif b == CR
            if !founddot
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            end

            b = read(socket, UInt8)
            if b != LF
                throw(HTTPError(400; message="Invalid request line (line terminator)", shouldclose=true))
            end

            foundend = true

            break
        elseif b == LF
            if strict
                throw(HTTPError(400; message="Invalid request line (line terminator)", shouldclose=true))
            end

            if !founddot
                throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
            end

            foundend = true
            break

        elseif !in(b, DIGIT)
            throw(HTTPError(400; message="Invalid http version text", shouldclose=true))
        else
            push!(version, b)
        end

        if (ind > maxversion)
            throw(HTTPError(505; message="Unwilling to accept a version string longer than $maxversion bytes long", shouldclose=true))
        end
        ind += 1
    end

    if !foundend
        throw(HTTPError(400; message="Unexpected end of file", shouldclose=true))
    end
    connection.method = String(method)
    connection.uri = String(uri)
    # copy here is necessary because `String` method seems to be consuming contents of the byte array
    # this has to be a big in Julia
    connection.httpversion = String(copy(version))

    if version != V1_0 && version != V1_1
        # TODO: support HTTP/2.0
        throw(HTTPError(505, message="HTTP Version not supported", shouldclose=true))
    end
    return connection
end

function processheaders!(connection::Connection; maxheader=DEFAULT_MAX_HEADER, keepcomments=true, strict=false)
    headers = Vector{Pair{String,String}}()
    headername = Vector{UInt8}()
    headervalue = Vector{UInt8}()
    foundend = false
    foundcolon = false
    instring = false
    commentlevel = 0
    socket = connection.socket
    ind = 1
    while isopen(socket)

        b = read(socket, UInt8)

        if foundcolon && length(headervalue) == 0 && commentlevel == 0 && !instring
            # trim leading white space for value
            while in(b, LWS)
                b = read(socket, UInt8)
                if ind > maxheader
                    throw(HTTPError(431; message="Unwilling to process headers with total size greater than $maxheaders bytes", shouldclose=true))
                end
                ind += 1
            end
        end

        if b == CR
            b = read(socket, UInt8)
            if b != LF
                throw(HTTPError(400; message="Invalid header line (line terminator)", shouldclose=true))
            elseif instring
                throw(HTTPError(400; message="Invalid header line (unterminated quoted string)", shouldclose=true))
            elseif commentlevel != 0
                throw(HTTPError(400; message="Invalid header line (unterminated comment)", shouldclose=true))
            elseif !foundcolon
                if length(headername) == 0
                    # last header found
                    foundend = true
                    break
                else
                    throw(HTTPError(400; message="Invalid header line (missing name separator)", shouldclose=true))
                end
            end

            push!(headers, String(headername) => String(headervalue))

            headername = Vector{UInt8}()
            headervalue = Vector{UInt8}()
            foundcolon = false
        elseif b == LF
            if strict
                throw(HTTPError(400; message="Invalid header line (line terminator)", shouldclose=true))
            elseif instring
                throw(HTTPError(400; message="Invalid header line (unterminated quoted string)", shouldclose=true))
            elseif commentlevel != 0
                throw(HTTPError(400; message="Invalid header line (unterminated comment)", shouldclose=true))
            elseif !foundcolon
                if length(headername) == 0
                    # last header found
                    foundend = true
                    break
                else
                    throw(HTTPError(400; message="Invalid header line (missing name separator)", shouldclose=true))
                end
            end

            push!(headers, String(headername) => String(headervalue))

            headername = Vector{UInt8}()
            headervalue = Vector{UInt8}()
            foundcolon = false
        elseif !foundcolon
            if b == COLON
                foundcolon = true
            elseif !in(b, TOKEN)
                throw(HTTPError(400; message="Invalid character in header name (charcode $b)", shouldclose=true))
            else
                push!(headername, b)
            end
        else
            if !instring && commentlevel == 0 && in(b, LWS)
                # TODO: LWS is only valid after or before a separator in a non-text header field
                #       or as part of a text header field.  If it is part of a text header
                #       field, should it be condensed?  Really don't want to backtrack to
                #       condense this after the fact.  Also, want to be consistent. So,
                #       probably always condense, or never condense.  Maybe allow this as
                #       an option.
                while in(b, LWS)
                    b = read(socket, UInt8)
                    if ind > maxheader
                        throw(HTTPError(431; message="Unwilling to process headers with size greater than $maxheaders bytes", shouldclose=true))
                    end
                    ind += 1
                end

                if b == CR
                    #  Trailing whitespace may be ignored
                    b = read(socket, UInt8)
                    if b != LF
                        throw(HTTPError(400; message="Invalid header line (line terminator)", shouldclose=true))
                    end
                    push!(headers, String(headername) => String(headervalue))

                    headername = Vector{UInt8}()
                    headervalue = Vector{UInt8}()
                    foundcolon = false
                    continue
                elseif b == LF
                    if strict
                        throw(HTTPError(400; message="Invalid header line (line terminator)", shouldclose=true))
                    end
                    push!(headers, String(headername) => String(headervalue))
                    headername = Vector{UInt8}()
                    headervalue = Vector{UInt8}()
                    foundcolon = false
                    continue
                else
                    # Per rfc2616 section 4.2, multiple whitespace LWS characters may
                    # be replaced by a single space character without altering its meaning
                    push!(headervalue, SP)
                end
            end

            if instring
                if b == QUOTE
                    push!(headervalue, b)
                    instring = false
                elseif b == BACKSLASH
                    # quoted pair
                    ind += 1
                    push!(headervalue, b)
                    b = read(socket, UInt8)
                    if !in(b, CHAR)
                        throw(HTTPError(400; message="Invalid quoted pair when parsing header value quoted string (charcode $b)", shouldclose=true))
                    end
                    push!(headervalue, b)
                elseif !in(b, TEXT)
                    throw(HTTPError(400; message="Invalid character when parsing header value quoted string (charcode $b)", shouldclose=true))
                else
                    push!(headervalue, b)
                end
            elseif commentlevel != 0
                if b == RIGHTPAREN
                    if keepcomments
                        push!(headervalue, b)
                    end
                    commentlevel -= 1
                elseif b == LEFTPAREN
                    if keepcomments
                        push!(headervalue, b)
                    end
                    commentlevel += 1
                elseif b == BACKSLASH
                    # quoted pair
                    if keepcomments
                        push!(headervalue, b)
                    end
                    ind += 1
                    b = read(socket, UInt8)
                    if !in(b, CHAR)
                        throw(HTTPError(400; message="Invalid quoted pair when parsing header value comment (charcode $b)", shouldclose=true))
                    end
                    if keepcomments
                        push!(headervalue, b)
                    end
                elseif !in(b, TEXT)
                    throw(HTTPError(400; message="Invalid character when parsing header value comment (charcode $b)", shouldclose=true))
                else
                    if keepcomments
                        push!(headervalue, b)
                    end
                end
            elseif b == QUOTE
                # TODO: This is not valid in the middle of a text header field
                push!(headervalue, b)
                instring = true
            elseif b == LEFTPAREN && in(lowercase(String(copy(headername))), COMMENTABLE_HEADER_FIELDS)
                # Dirty dirty bug consumes my
                if keepcomments
                    push!(headervalue, b)
                end
                commentlevel += 1
            elseif !in(b, TEXT) && !in(b, CTL)
                throw(HTTPError(400; message="Invalid character when parsing header value (charcode $b)", shouldclose=true))
            else
                # Per the spec the value may be *TEXT or combinations
                # of token, separators, and quoted-string.
                # We aren't really handling these two divergent possibilities here.
                # Currently, we are only validating that it MAY be consistent with one of these
                # two possibilities based on the characters given.  Enforcing a strict mode
                # where only one possibility is allowed is a TODO
                push!(headervalue, b)
            end
        end

        if ind > maxheader
            throw(HTTPError(431; message="Unwilling to process headers with size greater than $maxheaders bytes", shouldclose=true))
        end
        ind += 1
    end

    if !foundend
        throw(HTTPError(400; message="Unexpected end of stream", shouldclose=true))
    end

    headers
end

function readheaders!(connection::Connection; maxheader=DEFAULT_MAX_HEADER, keepcomments=true, strict=false)
    if !isequal(false, connection.headersread)
        throw(HTTPError(500, message="Cannot read http connection headers when they have already been read", shouldclose=true))
    end

    connection.headersread = true

    if isequal(false, connection.requestlineread)
        throw(HTTPError(500, message="Cannot read http connection headers before the requestline has been read", shouldclose=true))
    end

    headers = processheaders!(connection; maxheader=maxheader, keepcomments=keepcomments, strict=strict)
    connection.rawheaders = headers
    connection.headers = Headers(headers...)

    connection
end

function readtrailers!(connection::Connection; maxtrailer=DEFAULT_MAX_TRAILER, keepcomments=true, strict=false)
    if !connection.bodyfullyread
        throw(HTTPError(500, message="Cannot read http connection trailers when the body has not been processed", shouldclose=true))
    end

    trailers = processheaders!(connection; maxheader=maxtrailer, keepcomments=keepcomments, strict=strict)
    connection.rawtrailers = trailers
    connection.trailers = Headers(trailers...)

    connection
end

function createbodystream!(connection::Connection; maxbody=DEFAULT_MAX_BODY, maxtrailer=DEFAULT_MAX_TRAILER)
    if isequal(false, connection.headersread)
        throw(HTTPError(500; message="Cannot read body before headers have been parsed", shouldclose=true))
    end

    if isequal(true, connection.bodystartedread)
        throw(HTTPError(500; message="Cannot read body twice", shouldclose=true))
    end

    headers = connection.headers

    function onstart()
        connection.bodystartedread = true
    end

    function makesizedstream()
        lengths = []
        try
            lengths = parse_header_value(headers["Content-Length"]; parameterized=false, permittedchars=DIGIT)
        catch e
            throw(HTTPError(400; message="Invalid Content-Length header (invalid syntax)", shouldclose=true))
        end

        function iterator(res, pair)
            (v,) = pair
            try
                v = parse(UInt64, v, base=10)
            catch e
                throw(HTTPError(400; message="Invalid Content-Length header (field value must be a decimal number)", shouldclose=true))
            end
            if res == nothing || res == v
                v
            else
                throw(HTTPError(400; message="Invalid Content-Length header (multiple values do not match)", shouldclose=true))
            end
        end
        length = reduce(iterator, lengths; init=nothing)
        if length > maxbody
            throw(HTTPError(413; message="Content-Length too large", shouldclose=true))
        end
        function ondone()
            connection.bodyfullyread = true
        end
        connection.bodystream = SizedStream(connection.socket, length; ondone=ondone, onstart=onstart)
    end

    if haskey(headers, "Transfer-Encoding")
        encodings = nothing
        try
            # can extension encodings be parameterized?
            encodings = parse_header_value(headers["Transfer-Encoding"])
        catch e
            throw(HTTPError(400; message="Invalid Transfer-Encoding header (invalid syntax)", shouldclose=true))
        end
        hastransferencoding = reduce((res, encoding) -> res || (first(encoding) == "chunked"), encodings; init = false)
        if !hastransferencoding && haskey(headers, "Content-Length")
            makesizedstream()
        elseif first(encodings[end]) != "chunked"
            throw(HTTPError(400; message="Invalid Transfer-Encoding header (chunked must be the last value)", shouldclose=true))
        else
            function ondone()
                connection.bodyfullyread = true
                readtrailers!(connection; maxtrailer=maxtrailer)
            end
            connection.bodystream = ChunkedTransferDecodeStream(connection.socket; maxsize=maxbody, ondone=ondone, onstart=onstart)
        end
    elseif haskey(headers, "Content-Length")
        makesizedstream()
    else
        # no body
        connection.bodyfullyread = true
        connection.bodystartedread = true
    end

    connection
end

function processrequest!(connection::Connection;
    maxheader=DEFAULT_MAX_HEADER,
    maxuri=DEFAULT_MAX_URI,
    maxmethod=DEFAULT_MAX_METHOD,
    maxversion=DEFAULT_MAX_VERSION,
    maxbody=DEFAULT_MAX_BODY,
    maxtrailer=DEFAULT_MAX_TRAILER,
    keepcomments=true,
    strict=false)

    readrequestline!(connection; maxuri=maxuri, maxmethod=maxmethod, maxversion=maxversion, strict=strict)
    readheaders!(connection; maxheader=maxheader, strict=strict)
    createbodystream!(connection; maxbody=maxbody, maxtrailer=maxtrailer)

    connection
end

function writeresponse!(connection::Connection, statusCode::Integer, statusReason::AbstractString, headers::Headers, body::IO)
    if connection.responsesent
        warn("Response has already been sent")
        return connection
    end
    statusCode = UInt16(statusCode)
    socket = connection.socket
    write(socket, "HTTP/1.1")
    write(socket, SP)
    write(socket, string(statusCode))
    write(socket, SP)
    write(socket, string(statusReason))
    write(socket, CRLF)

    for (name,value) in headers
        write(socket, name)
        write(socket, COLON)
        write(socket, SP)
        write(socket, value)
        write(socket, CRLF)
    end

    write(socket, CRLF)

    while !eof(body)
        write(socket, read(body, UInt8))
    end

    connection.responsesent = true

    connection
end


text:
using ITensors
using LinearAlgebra
using Printf
using Random
using TBLIS

Random.seed!(1234)

let
  nthreads = 4

  N = 100
  sites = siteinds("S=1", N)
  ampo = AutoMPO()
  for j in 1:(N - 1)
    ampo += 0.5, "S+", j, "S-", j + 1
    ampo += 0.5, "S-", j, "S+", j + 1
    ampo += "Sz", j, "Sz", j + 1
  end
  H = MPO(ampo, sites)
  psi0 = randomMPS(sites, 10)

  sweeps_compile = Sweeps(2)
  maxdim!(sweeps_compile, 10)
  cutoff!(sweeps_compile, 1E-15)

  sweeps = Sweeps(6)
  maxdim!(sweeps, 20, 100, 200, 300)
  cutoff!(sweeps, 0.0)
  @show sweeps

  #
  # Using BLAS backend
  #

  ITensors.disable_tblis()
  BLAS.set_num_threads(nthreads)

  # Compile
  dmrg(H, psi0, sweeps_compile; outputlevel=0)

  println("Using BLAS with $nthreads threads\n")
  energy, psi = @time dmrg(H, psi0, sweeps)
  @printf("Final energy = %.12f\n", energy)
  println()

  #
  # Using TBLIS backend
  #

  ITensors.enable_tblis()
  BLAS.set_num_threads(1)
  TBLIS.set_num_threads(nthreads)

  # Compile
  dmrg(H, psi0, sweeps_compile; outputlevel=0)

  println("Using TBLIS with $(TBLIS.get_num_threads()) threads (and 1 BLAS thread)\n")
  energy, psi = @time dmrg(H, psi0, sweeps)
  @printf("Final energy = %.12f\n", energy)
end

ext: jl


text:
"""
    apple_alert!(;title="Julia", subtitle="", sound="")

Use a nicer appearing backend for MacOS systems.

# Arguments
- `title::AbstractString="Julia"`: Notification title
- `subtitle::AbstractString=""`: Notification subtitle
- `sound::AbstractString=""`: Notification sound
- `clear::Bool`: if true, clears all apple_alert settings and reverts to default
  backend

# Examples
```jldoctest
julia> Alert.apple_alert!(title="𝒥𝓊𝓁𝒾𝒶", subtitle="Alert.jl", sound="Crystal");

julia> alert("From a fancy backend!")
```
"""
function apple_alert!(;clear=false,kwds...)
    if clear
        set_alert_backend!()
    else
        Alert.set_alert_backend!(apple_backend(;kwds...))
    end
end

"""

    apple_backend(title="Julia", subtitle="", sound="", test=false)

Helper function to define a Script Editor backend for MacOS systems.

# Arguments
- `title::AbstractString="Julia"`: Notification title
- `subtitle::AbstractString=""`: Notification subtitle
- `sound::AbstractString=""`: Notification sound
- `test::Bool=false`: Display a sample notification when initializing

# Examples
```jldoctest
julia> mybackend = Alert.apple_backend(title="𝒥𝓊𝓁𝒾𝒶", subtitle="Alert.jl", sound="Crystal");

julia> Alert.set_alert_backend!(mybackend);

julia> alert("From a fancy backend!")
```
"""
function apple_backend(;
    title::AbstractString="Julia",
    subtitle::AbstractString="",
    sound::AbstractString="",
    test::Bool=false
)
    @static if !Sys.isapple()
        @warn "You are attempting to set a MacOS backend on a different system!"
    end

    options = ""
    options = length(title) > 0 ? options * " with title \"$title\"" : options
    options = length(subtitle) > 0 ? options * " subtitle \"$subtitle\"" : options
    options = length(sound) > 0 ? options * " sound name \"$sound\"" : options

    @debug "Current command" cmd=`osascript -e 'display notification "'message'"'$options''`

    if test
        testmsg = "Your notifications will now look like this!"
        run(`osascript -e 'display notification "'$testmsg'"'$options''`)
    end

    return message -> run(`osascript -e 'display notification "'$message'"'$options''`)
end


text:
#= MIT License

Copyright (c) 2020, 2021 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

module Utils

# data structures for the flight state and the flight log
# functions for creating a demo flight state, demo flight log, loading and saving flight logs
# function se() for reading the settings
# in addition helper functions for working with rotations

using Rotations, StaticArrays, StructArrays, RecursiveArrayTools, Arrow, YAML
export SysState, ExtSysState, SysLog, MyFloat

export demo_state, demo_syslog, demo_log, load_log, syslog2extlog, save_log, rot3d, se

const MyFloat = Float32               # type to use for postions
const DATA_PATH = "./data"            # path for log files and other data

mutable struct Settings
    log_file::String
    model::String
    segments::Int64       # number of tether segments
    sample_freq::Int64
    time_lapse::Float64
    zoom::Float64
    fixed_font::String
end
const SETTINGS = [Settings("","",0,0,0,0,"")]

# getter function for the Settings struct
function se()
    if SETTINGS[1].segments == 0
        # load settings from YAML
        dict = YAML.load_file(joinpath(DATA_PATH, "settings.yaml"))
        SETTINGS[1].log_file=dict["system"]["log_file"]
        SETTINGS[1].model=dict["kite"]["model"]
        SETTINGS[1].segments=dict["system"]["segments"]
        SETTINGS[1].sample_freq=dict["system"]["sample_freq"]
        SETTINGS[1].time_lapse=dict["system"]["time_lapse"]
        SETTINGS[1].zoom=dict["system"]["zoom"]
        SETTINGS[1].fixed_font=dict["system"]["fixed_font"]
    end
    return SETTINGS[1]
end

# basic system state; one of these will be saved per time step
struct SysState
    time::Float64                          # time since start of simulation in seconds
    orient::MVector{4, Float32}            # orientation of the kite (quaternion)
    elevation::MyFloat                     # elevation angle in radians
    azimuth::MyFloat                       # azimuth angle in radians
    l_tether::MyFloat                      # tether length [m]
    v_reelout::MyFloat                     # reel out velocity [m/s]
    force::MyFloat                         # tether force [N]
    depower::MyFloat                       # depower settings
    v_app::MyFloat                         # apparent wind speed [m/s]
    X::MVector{se().segments+1, MyFloat}   # vector of particle positions in x
    Y::MVector{se().segments+1, MyFloat}   # vector of particle positions in y
    Z::MVector{se().segments+1, MyFloat}   # vector of particle positions in z
end

# extended SysState containing derived values for plotting
struct ExtSysState
    time::Float64                          # time since launch in seconds
    orient::UnitQuaternion{Float32}        # orientation of the kite
    X::MVector{se().segments+1, MyFloat}   # vector of particle positions in x
    Y::MVector{se().segments+1, MyFloat}   # vector of particle positions in y
    Z::MVector{se().segments+1, MyFloat}   # vector of particle positions in z
    x::MyFloat                             # kite position in x
    y::MyFloat                             # kite position in y
    z::MyFloat                             # kite position in z
end

# flight log, containing the basic data as struct of arrays
# and in addition an extended view on the data that includes derived/ calculated values for plotting
# finally meta data like the file name of the log file is included
struct SysLog
    name::String
    syslog::StructArray{SysState}    # struct of vectors
    extlog::StructArray{ExtSysState} # struct of vectors, containing derived values
end

# functions
"""
Calculate the rotation of reference frame (ax, ay, az) so that it matches the reference frame (bx, by, bz).
All parameters must be 3-element vectors. Both refrence frames must be orthogonal,
all vectors must already be normalized.
Source: http://en.wikipedia.org/wiki/User:Snietfeld/TRIAD_Algorithm
"""
function rot3d(ax, ay, az, bx, by, bz)
    R_ai = hcat(ax, az, ay)
    R_bi = hcat(bx, bz, by)
    return R_bi * R_ai'
end

# create a demo state with a given height and time
function demo_state(height=6.0, time=0.0)
    a = 10
    X = range(0, stop=10, length=(se().segments)+1)
    Y = zeros(length(X))
    Z = (a .* cosh.(X./a) .- a) * height/ 5.430806
    r_xyz = RotXYZ(pi/2, -pi/2, 0)
    q = UnitQuaternion(r_xyz)
    orient = MVector{4, Float32}(q.w, q.x, q.y, q.z)
    elevation = asin(Z[end]/X[end])
    return SysState(time, orient, elevation,0.,0.,0.,0.,0.,0.,X, Y, Z)
end

# create a demo flight log with given name [String] and duration [s]
function demo_syslog(name="Test flight"; duration=10)
    max_height = 6.03
    steps   = Int(duration * se().sample_freq) + 1
    time_vec = Vector{Float64}(undef, steps)
    myzeros = zeros(MyFloat, steps)
    elevation = Vector{Float64}(undef, steps)
    orient_vec = Vector{MVector{4, Float32}}(undef, steps)
    X_vec = Vector{MVector{se().segments+1, MyFloat}}(undef, steps)
    Y_vec = Vector{MVector{se().segments+1, MyFloat}}(undef, steps)
    Z_vec = Vector{MVector{se().segments+1, MyFloat}}(undef, steps)
    for i in range(0, length=steps)
        state = demo_state(max_height * i/steps, i/se().sample_freq)
        time_vec[i+1] = state.time
        orient_vec[i+1] = state.orient
        elevation[i+1] = asin(state.Z[end]/state.X[end])
        X_vec[i+1] = state.X
        Y_vec[i+1] = state.Y
        Z_vec[i+1] = state.Z
    end
    return StructArray{SysState}((time_vec, orient_vec, elevation, myzeros,myzeros,myzeros,myzeros,myzeros,myzeros, X_vec, Y_vec, Z_vec))
end

# extend a flight systom log with the fieds x, y, and z (kite positions) and convert the orientation to the type UnitQuaternion
function syslog2extlog(syslog)
    x_vec = @view VectorOfArray(syslog.X)[end,:]
    y_vec = @view VectorOfArray(syslog.Y)[end,:]
    z_vec = @view VectorOfArray(syslog.Z)[end,:]
    orient_vec = Vector{UnitQuaternion{Float32}}(undef, length(syslog.time))
    for i in range(1, length=length(syslog.time))
        orient_vec[i] = UnitQuaternion(syslog.orient[i])
    end
    return StructArray{ExtSysState}((syslog.time, orient_vec, syslog.X, syslog.Y, syslog.Z, x_vec, y_vec, z_vec))
end

# create an artifical log file for demonstration purposes
function demo_log(name="Test_flight"; duration=10)
    syslog = demo_syslog(name, duration=duration)
    return SysLog(name, syslog, syslog2extlog(syslog))
end

function save_log(flight_log::SysLog)
    Arrow.ArrowTypes.registertype!(SysState, SysState)
    filename=joinpath(DATA_PATH, flight_log.name) * ".arrow"
    Arrow.write(filename, flight_log.syslog, compress=:lz4)
end

function load_log(filename::String)

    Arrow.ArrowTypes.registertype!(SysState, SysState)
    Arrow.ArrowTypes.registertype!(MVector{4, Float32}, MVector{4, Float32})
    if isnothing(findlast(isequal('.'), filename))
        fullname = joinpath(DATA_PATH, filename) * ".arrow"
    else
        fullname = joinpath(DATA_PATH, filename)
    end
    table = Arrow.Table(fullname)
    myzeros = zeros(MyFloat, length(table.time))
    syslog = StructArray{SysState}((table.time, table.orient, table.elevation, table.azimuth, table.l_tether, table.v_reelout, table.force, table.depower, table.v_app, table.X, table.Y, table.Z))
    return SysLog(basename(fullname[1:end-6]), syslog, syslog2extlog(syslog))
end

function test(save=false)
    if save
        log_to_save=demo_log()
        save_log(log_to_save)
    end
    return(load_log("Test_flight.arrow"))
end

end


text:
@testset "Eps" begin
    @testset for T in (Posit8, Posit16, Posit32,
                    # pX1_mul is not supported yet by C library
                    # Posit8_1, Posit16_1, Posit24_1,
                    Posit8_2, Posit16_2, Posit24_2)

        @test eps(T) == eps(one(T))
        @test floatmin(T) == eps(zero(T))
        if ~(T == Posit8)
            # Excluding Posit(8,0) as this test is only true for exponent bits > 0
            @test floatmax(T) == eps(floatmax(T))
        end
        @test notareal(T) == eps(notareal(T))
    end
end


text:
module DiscFESpaces

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery

export DiscFESpace

import Gridap: num_free_dofs
import Gridap: num_diri_dofs
import Gridap: diri_tags
import Gridap: apply_constraints
import Gridap: apply_constraints_rows
import Gridap: apply_constraints_cols
import Gridap: celldofids
import Gridap: interpolate_values
import Gridap: interpolate_diri_values
import Gridap: CellField
import Gridap: CellBasis
import Gridap: Triangulation

"""
Concrete implementation of a discontinuous FE space
with no boundary conditions and for arbitrary
reference elements. Even in the case that the used reference element
has dofs in low dimensional n-faces, the resulting FE space
is always discontinuous.
"""
struct DiscFESpace{D,Z,T} <: FESpace{D,Z,T}
  reffe::RefFE{Z,T}
  nfree::Int
  celldofs::IndexCellVector{Int}
  cellbasis::CellBasis{Z,T}
  trian::Triangulation{D,Z}
end

function DiscFESpace(
  reffe::RefFE{Z,T}, trian::Triangulation{D,Z}) where {D,Z,T}
  cellbasis = _setup_cellbasis(reffe,trian)
  celldofs, nfree = _setup_celldofs(reffe,ncells(trian))
  DiscFESpace{D,Z,T}(reffe,nfree,celldofs,cellbasis,trian)
end

function DiscFESpace(
  reffe::RefFE{D,T}, model::DiscreteModel{D}) where {D,T}

  trian = Triangulation(model)
  DiscFESpace(reffe,trian)
end

num_free_dofs(this::DiscFESpace) = this.nfree

num_diri_dofs(::DiscFESpace) = 0

diri_tags(::DiscFESpace) = Int[]

apply_constraints(
  this::DiscFESpace, cellvec::CellVector, cellids::CellNumber) = cellvec

apply_constraints_rows(
  this::DiscFESpace, cellmat::CellMatrix, cellids::CellNumber) = cellmat

apply_constraints_cols(
  this::DiscFESpace, cellmat::CellMatrix, cellids::CellNumber) = cellmat

function celldofids(this::DiscFESpace)
  this.celldofs
end

function interpolate_values(this::DiscFESpace,fun::Function)
  _setup_interpolate_values(this, fun)
end

function interpolate_values(this::DiscFESpace{D,Z,T},val::T) where {D,Z,T}
  fun(x) = val
  _setup_interpolate_values(this, fun)
end

function interpolate_diri_values(this::DiscFESpace, funs::Vector{<:Function})
  _interpolate_diri_values(this::DiscFESpace, funs)
end

function interpolate_diri_values(this::DiscFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
  _interpolate_diri_values(this::DiscFESpace, vals)
end

function _interpolate_diri_values(this::DiscFESpace, funs)
  @assert length(funs) == 0
  T = value_type(this)
  E = dof_type(T)
  zeros(E,0)
end

function CellField(
  this::DiscFESpace{D,Z,T},
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E}) where {D,Z,T,E}

  @assert length(diri_dofs) == 0

  _setup_cellfield(this,free_dofs)

end

function CellBasis(this::DiscFESpace{D,Z,T}) where {D,Z,T}
  this.cellbasis
end

function Triangulation(this::DiscFESpace{D,Z}) where {Z,D}
  this.trian
end

# Helpers

function _setup_cellbasis(reffe,trian)

  trian_reffes = CellRefFEs(trian)
  @notimplementedif !isa(trian_reffes,ConstantCellValue)
  @assert polytope(reffe).extrusion == polytope(trian_reffes.value).extrusion

  shb = ConstantCellValue(shfbasis(reffe), ncells(trian))
  phi = CellGeomap(trian)
  attachgeomap(shb,phi)
end

function _setup_celldofs(reffe,ncells)
  ncelldofs = length(shfbasis(reffe))
  nfree = ncelldofs*ncells
  celldofs_data = collect(1:nfree)
  celldofs = CellVectorFromDataAndStride(celldofs_data,ncelldofs)
  celldofs, nfree
end

function _setup_cellfield(this,free_dofs)

  ncell = ncells(this.trian)
  ncelldofs = length(shfbasis(this.reffe))
  @assert length(free_dofs) == ncell*ncelldofs

  celldofsvals = CellVectorFromDataAndStride(free_dofs,ncelldofs)

  lincomb(this.cellbasis, celldofsvals)

end

function  _setup_interpolate_values(fesp, fun)

  reffe = fesp.reffe
  dofb = dofbasis(reffe)
  trian = fesp.trian
  phi = CellGeomap(trian)
  uphys = fun ∘ phi

  T = value_type(fesp)
  E = dof_type(T)
  free_dofs = zeros(E, num_free_dofs(fesp))
  diri_dofs = zeros(E, num_diri_dofs(fesp))
  aux = zeros(E, length(dofb))

  _interpolate_kernel!(free_dofs,dofb,uphys,aux)

  (free_dofs, diri_dofs)

end

function  _interpolate_kernel!(free_dofs,dofb,uphys,aux)
  i = 1
  for u in uphys
    evaluate!(dofb,u,aux)
    for a in aux
      free_dofs[i] = a
      i += 1
    end
  end
end

end # module


text:
println("=============================================================================")
println("Beginning execution of script GMHExamples.jl/target/sincos/scripts/SinCos2.jl")
println("=============================================================================")

import GeneralizedMetropolisHastings #pre-include to avoid race condition when running in parallel
import GMHExamples #pre-include to avoid race conditions

println("Number of parallel processes running: ",nprocs())
println("Use addprocs() in the REPL if you want to run on more processes")

@everywhere using GeneralizedMetropolisHastings
@everywhere using GMHExamples
using PyPlot

###Same as SinCos1.jl, but with an in-place target function
nproposals = 30

#MCMC iteration specifications
nburnin = 400
niterations = 1000
ntunerperiod = 20

#Time points to simulate the sine-cosine model
timepoints = linspace(0.0,10.0,500)

###Values of the model parameters (sine and cosine coefficients)
a = 1/2
b = 2/3
lows = [0.2,0.3]
highs = [0.8,1.0]

###The variance of the normal noise on the input data
variance = [0.09,0.09]

println("==========================================")
println("Simulation parameters defined successfully")
println("==========================================")

###Create a Spring-Mass model with measurement data and ODE function
m = sincosmodel!(timepoints,[a,b],variance,lows,highs)

###Show the model
println("==========================")
println("Model defined successfully")
println("==========================")
show(m)

###Plot the measurement data (simmulated data + noise)
figure(string("SinCos2-",nprocs())) ; clf()
subplot(221)
plot(dataindex(m),measurements(m)[:,1];label="sine")
plot(dataindex(m),measurements(m)[:,2];label="cosine")
xlabel("Time")
ylabel("Amplitude")
title("Sine-Cosine data")
grid("on")
legend(loc="upper right",fancybox="true")

###Create a Metropolis sampler with a Normal proposal density
s = sampler(:mh,:normal,0.1,ones(2))
println("============================")
println("Sampler defined successfully")
println("============================")
show(s)

###Create a tuner that scales the proposal density
t = tuner(:scale,ntunerperiod,0.5,:erf)
println("==========================")
println("Tuner defined successfully")
println("==========================")
show(t)

###Create a Generalized Metropolis-Hastings runner (which will default to Standard MH when nproposals=1)
p = policy(:mh,nproposals;initialize=:prior)
r = runner(p,niterations,nproposals;numburnin = nburnin)
println("===========================")
println("Runner defined successfully")
println("===========================")
show(r)

###Run the MCMC (can take quite a bit of time)
println("=======================")
println("Run the MCMC simulation")
println("=======================")
@time c = run!(r,m,s,t)
println("=========================")
println("Completed MCMC simulation")
println("=========================")

###Show the result of the simulations
show(c)

nparas = numparas(m)
meanparamvals = mean(samples(c),2)
stdparamvals = std(samples(c),2)

println("Results of the MCMC simulation:")
println(" mean a: ",meanparamvals[1])
println(" mean b: ",meanparamvals[2])

println("================")
println("Plotting results")
println("================")

###Plot the loglikelihood values across samples
###After an initial few low values, this should remain relatively high
subplot(222)
plot(1:numsamples(c),logposterior(c,m))
title("Log-Posterior values across samples")
xlabel("Samples")
ylabel("Log-Posterior")

###Plot the histograms of a,b values
for i=1:nparas
  subplot(222 + i)
    h = PyPlot.plt[:hist](vec(getindex(samples(c),i,:)),20)
  grid("on")
  title("Parameter $(parameters(m)[i].key)")
  xlabel("Values")
  ylabel("Samples")
end

###Finally, plot the average model results in the data window
subplot(221)
modeldata = evaluate!(m,vec(meanparamvals))
plot(dataindex(m),modeldata[:,1];label="model sine")
plot(dataindex(m),modeldata[:,2];label="model cosine")
legend(loc="upper right",fancybox="true")


text:
@recipe function f(tm::TorusMantle; n_arc = 40, subn = 10)
    seriestype --> :mesh3d
    if haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :mesh3d
        @series begin
            label --> "Torus Mantle"
            mesh(tm, n_arc = n_arc)
        end
    else
        ls = lines(tm)
        linecolor --> :black
        @series begin
            label --> "Torus Mantle"
            ls[1]
        end
        for i in 2:length(ls)
            @series begin
                label := nothing
                ls[i]
            end
        end
    end
    if (haskey(plotattributes, :show_normal) && plotattributes[:show_normal]) &&
            tm.φ === nothing && tm.θ === nothing
        @series begin
            label := nothing
            seriestype := :vector
            pts = _plt_points_for_normals(tm)
            ns = broadcast(p -> normal(tm, p) / 5, pts)
            [(pts[i], ns[i]) for i in eachindex(pts)]
        end
    end
end

function _plt_points_for_normals(tm::TorusMantle{T}) where {T}
    pts = [ CartesianPoint{T}( tm.r_torus+tm.r_tube, zero(T), zero(T)),
            CartesianPoint{T}( tm.r_torus-tm.r_tube, zero(T), zero(T)),
            CartesianPoint{T}( tm.r_torus, zero(T), tm.r_tube ),
            CartesianPoint{T}( tm.r_torus, zero(T),-tm.r_tube ),
            CartesianPoint{T}(-tm.r_torus+tm.r_tube, zero(T), zero(T)),
            CartesianPoint{T}(-tm.r_torus-tm.r_tube, zero(T), zero(T)),
            CartesianPoint{T}(-tm.r_torus, zero(T), tm.r_tube ),
            CartesianPoint{T}(-tm.r_torus, zero(T),-tm.r_tube ),
            CartesianPoint{T}( zero(T), tm.r_torus+tm.r_tube, zero(T)),
            CartesianPoint{T}( zero(T), tm.r_torus-tm.r_tube, zero(T)),
            CartesianPoint{T}( zero(T), tm.r_torus, tm.r_tube ),
            CartesianPoint{T}( zero(T), tm.r_torus,-tm.r_tube ),
            CartesianPoint{T}( zero(T),-tm.r_torus+tm.r_tube, zero(T)),
            CartesianPoint{T}( zero(T),-tm.r_torus-tm.r_tube, zero(T)),
            CartesianPoint{T}( zero(T),-tm.r_torus, tm.r_tube ),
            CartesianPoint{T}( zero(T),-tm.r_torus,-tm.r_tube ) ]
    _transform_into_global_coordinate_system(pts, tm)
end


text:
isincident(v::VRepElement, h::HRepElement) = v in hyperplane(h)
isincident(h::HRepElement, v::VRepElement) = v in hyperplane(h)

abstract type Incident{N, T, ElemT<:Element{N, T}, PT<:Polyhedron{N, T}, IdxT<:Index{N, T}} end

function Base.get(p::Polyhedron{N, T}, inc::Incident{N, T, ElemT}) where {N, T, ElemT}
    el = get(p, inc.idx)
    incs = _inctype(inc)[]
    for idx in Indices{N, T, ElemT}(p)
        eli = get(p, idx)
        if isincident(el, eli)
            push!(incs, _incel(inc, idx, eli))
        end
    end
    incs
end

struct IncidentElements{N, T, ElemT, PT, IdxT} <: Incident{N, T, ElemT, PT, IdxT}
    p::PT
    idx::IdxT
end
IncidentElements{N, T, ElemT}(p, idx) where {N, T, ElemT<:Element{N, T}} = IncidentElements{N, T, ElemT, typeof(p), typeof(idx)}(p, idx)
_inctype(inc::IncidentElements{N, T, ElemT}) where {N, T, ElemT} = ElemT
_incel(inc::IncidentElements, idx, eli) = eli

struct IncidentIndices{N, T, ElemT, PT, IdxT} <: Incident{N, T, ElemT, PT, IdxT}
    p::PT
    idx::IdxT
end
IncidentIndices{N, T, ElemT}(p, idx) where {N, T, ElemT<:Element{N, T}} = IncidentIndices{N, T, ElemT, typeof(p), typeof(idx)}(p, idx)
_inctype(inc::IncidentIndices{N, T, ElemT}) where {N, T, ElemT} = Index{N, T, ElemT}
_incel(inc::IncidentIndices, idx, eli) = idx


text:
## regular expressions to use
whiteRe = r"\s*"
spaceRe = r"\s+"
nonSpaceRe = r"\S"
eqRe = r"\s*="
curlyRe = r"\s*\}"

# # section
# ^ inversion
# / close section
# > partials
# { dont' escape
# & unescape a variable
# = set delimiters {{=<% %>=}} will set delimiters to <% %>
# ! comments
# | lamda "section" with *evaluated* value
tagRe = r"^[#^/<>{&=!|]"

function asRegex(txt)
    for i in ("[","]")
        txt = replace(txt, Regex("\\$i") => "\\$i")
    end
    for i in ("(", ")", "{","}", "|")
        txt = replace(txt, Regex("[$i]") => "[$i]")
    end
    Regex(txt)
end

isWhitespace(x) = occursin(whiteRe, x)
function stripWhitespace(x)
    y = replace(x, r"^\s+" => "")
    replace(y, r"\s+$" => "")
end

## this is for falsy value
## Falsy is true if x is false, 0 length, "", ...
falsy(x::Bool) = !x
falsy(x::Array) = isempty(x) || all(falsy, x)
falsy(x::AbstractString) = x == ""
falsy(x::Nothing) = true
falsy(x::Missing) = true
#falsy(x) = (x == nothing) || false                #  default
function falsy(x)
    Tables.istable(x) && isempty(Tables.rows(x)) && return true
    x == nothing && return true
    false
end

## escape_html with entities

entityMap = [("&", "&"),
             ("<", "<"),
             (">", ">"),
             ("'", "'"),
             ("\"", """),
             ("/", "/")]

function escape_html(x)
    y = string(x)
    for (k,v) in entityMap
        y = replace(y, k => v)
    end
    y
end

## Make these work
function escapeRe(string)
    replace(string, r"[\-\[\]{}()*+?.,\\\^$|#\s]" => "\\\$&");
end

function escapeTags(tags)
   [Regex(escapeRe(tags[1]) * "\\s*"),
    Regex("\\s*" * escapeRe(tags[2]))]
end

# key may be string or a ":symbol"
function normalize(key)
    if occursin(r"^:", key)
        key = key[2:end]
        key = Symbol(key)
    end
    return key
end

## hueristic to avoid loading DataFrames
## Once `Tables.jl` support for DataFrames is available, this can be dropped
is_dataframe(x) = !isa(x, Dict) && !isa(x, Module) &&!isa(x, Array) && occursin(r"DataFrame", string(typeof(x)))


text:
"""
    ODESolvers

Ordinary differential equation solvers
"""
module ODESolvers

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays
using ..SystemSolvers
using ..MPIStateArrays: array_device, realview
using ..GenericCallbacks

export solve!, updatedt!, gettime, getsteps

abstract type AbstractODESolver end

"""
    gettime(solver::AbstractODESolver)

Returns the current simulation time of the ODE solver `solver`
"""
gettime(solver::AbstractODESolver) = solver.t

"""
    getdt(solver::AbstractODESolver)

Returns the current simulation time step of the ODE solver `solver`
"""
getdt(solver::AbstractODESolver) = solver.dt

"""
    getsteps(solver::AbstractODESolver)

Returns the number of completed time steps of the ODE solver `solver`
"""
getsteps(solver::AbstractODESolver) = solver.steps

"""
    ODESolvers.general_dostep!(Q, solver::AbstractODESolver, p,
                               timeend::Real, adjustfinalstep::Bool)

Use the solver to step `Q` forward in time from the current time, to the time
`timeend`. If `adjustfinalstep == true` then `dt` is adjusted so that the step
does not take the solution beyond the `timeend`.
"""
function general_dostep!(
    Q,
    solver::AbstractODESolver,
    p,
    timeend::Real;
    adjustfinalstep::Bool,
)
    time, dt = gettime(solver), getdt(solver)
    final_step = false
    if adjustfinalstep && time + dt > timeend
        orig_dt = dt
        dt = timeend - time
        updatedt!(solver, dt)
        final_step = true
    end
    @assert dt > 0

    dostep!(Q, solver, p, time)

    if !final_step
        updatetime!(solver, time + dt)
    else
        updatedt!(solver, orig_dt)
        updatetime!(solver, timeend)
    end
end

"""
    updatetime!(solver::AbstractODESolver, time)

Change the current time to `time` for the ODE solver `solver`.
"""
updatetime!(solver::AbstractODESolver, time) = (solver.t = time)

"""
    updatedt!(solver::AbstractODESolver, dt)

Change the time step size to `dt` for the ODE solver `solver`.
"""
updatedt!(solver::AbstractODESolver, dt) = (solver.dt = dt)

"""
    updatesteps!(solver::AbstractODESolver, dt)

Set the number of elapsed time steps for the ODE solver `solver`.
"""
updatesteps!(solver::AbstractODESolver, steps) = (solver.steps = steps)

isadjustable(solver::AbstractODESolver) = true

# {{{ run!
"""
    solve!(Q, solver::AbstractODESolver; timeend,
           stopaftertimeend=true, numberofsteps, callbacks)

Solves an ODE using the `solver` starting from a state `Q`. The state `Q` is
updated inplace. The final time `timeend` or `numberofsteps` must be specified.

A series of optional callback functions can be specified using the tuple
`callbacks`; see the `GenericCallbacks` module.
"""
function solve!(
    Q,
    solver::AbstractODESolver,
    param = nothing;
    timeend::Real = Inf,
    adjustfinalstep = true,
    numberofsteps::Integer = 0,
    callbacks = (),
)

    @assert isfinite(timeend) || numberofsteps > 0
    if adjustfinalstep && !isadjustable(solver)
        error("$solver does not support time step adjustments. Can only be used with `adjustfinalstep=false`.")
    end
    t0 = gettime(solver)

    # Loop through an initialize callbacks (if they need it)
    GenericCallbacks.init!(callbacks, solver, Q, param, t0)

    step = 0
    time = t0
    while time < timeend
        step += 1
        updatesteps!(solver, step)

        time = general_dostep!(
            Q,
            solver,
            param,
            timeend;
            adjustfinalstep = adjustfinalstep,
        )

        val = GenericCallbacks.call!(callbacks, solver, Q, param, time)
        if val !== nothing && val > 0
            return gettime(solver)
        end

        # Figure out if we should stop
        if step == numberofsteps
            break
        end
    end

    # Loop through to fini callbacks
    GenericCallbacks.fini!(callbacks, solver, Q, param, time)

    return gettime(solver)
end
# }}}

include("BackwardEulerSolvers.jl")
include("MultirateInfinitesimalGARKExplicit.jl")
include("MultirateInfinitesimalGARKDecoupledImplicit.jl")
include("LowStorageRungeKuttaMethod.jl")
include("LowStorageRungeKutta3NMethod.jl")
include("StrongStabilityPreservingRungeKuttaMethod.jl")
include("AdditiveRungeKuttaMethod.jl")
include("MultirateInfinitesimalStepMethod.jl")
include("MultirateRungeKuttaMethod.jl")
include("SplitExplicitMethod.jl")
include("DifferentialEquations.jl")

end # module

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
"""
  Utils.jl File

Machine Learning shared utility functions (Module BetaML.Utils)

`?BetaML.Utils` for documentation

- Part of [BetaML](https://github.com/sylvaticus/BetaML.jl)
- New to Julia? [A concise Julia tutorial](https://github.com/sylvaticus/juliatutorial) - [Julia Quick Syntax Reference book](https://julia-book.com)

"""

"""
    Utils module

Provide shared utility functions for various machine learning algorithms.

For the complete list of functions provided see below. The main ones are:

## Helper functions for logging
- Most BetAML functions accept a parameter `verbosity` that expect one of the element in the `Verbosity` enoum (`NONE`, `LOW`, `STD`, `HIGH` or `FULL`)
- Writing complex code and need to find where something is executed ? Use the macro [`@codeLocation`](@ref)

## Stochasticity management
- Utils provide [`FIXEDSEED`], [`FIXEDRNG`] and [`generateParallelRngs`](@ref). All stochastic functions accept a `rng` paraemter. See the "Getting started" section in the tutorial for details.

## Data processing
- Various small and large utilities for helping processing the data, expecially before running a ML algorithm
- Includes [`getPermutations`](@ref), [`oneHotEncoder`](@ref), [`integerEncoder`](@ref) (and [`integerDecoder`](@ref)), [`partition`](@ref), [`scale`](@ref) (and [`getScaleFactors`](@ref)), [`pca`](@ref), [`crossValidation`](@ref)

## Samplers
- Utilities to sample from data (e.g. for neural network training or for cross-validation)
- Include the "generic" type [`SamplerWithData`](@ref), together with the sampler implementation [`KFold`](@ref) and the function [`batch`](@ref)

## Transformers
- Funtions that "transform" a single input (that can be also a vector or a matrix)
- Includes varios NN "activation" functions ([`relu`](@ref), [`celu`](@ref), [`sigmoid`](@ref), [`softmax`](@ref), [`pool1d`](@ref)) and their derivatives (`d[FunctionName]`), but also [`gini`](@ref), [`entropy`](@ref), [`variance`](@ref), [`BIC`](@ref bic), [`AIC`](@ref aic)

## Measures
- Several functions of a pair of parameters (often `y` and `ŷ`) to measure the goodness of `ŷ`, the distance between the two elements of the pair, ...
- Includes "classical" distance functions ([`l1_distance`](@ref), [`l2_distance`](@ref), [`l2²_distance`](@ref) [`cosine_distance`](@ref)), "cost" functions for continuous variables ([`squaredCost`](@ref), [`meanRelError`](@ref)) and comparision functions for multui-class variables ([`crossEntropy`](@ref), [`accuracy`](@ref), [`ConfusionMatrix`](@ref)).

"""
module Utils

using LinearAlgebra, Printf, Random, Statistics, Combinatorics, Zygote, CategoricalArrays, StableRNGs

using ForceImport
@force using ..Api

export Verbosity, NONE, LOW, STD, HIGH, FULL,
       FIXEDSEED, FIXEDRNG, @codeLocation, generateParallelRngs,
       reshape, makeColVector, makeRowVector, makeMatrix, issortable, getPermutations,
       oneHotEncoder, integerEncoder, integerDecoder, colsWithMissing, getScaleFactors, scale, scale!, batch, partition, shuffle, pca,
       didentity, relu, drelu, elu, delu, celu, dcelu, plu, dplu,  #identity and rectify units
       dtanh, sigmoid, dsigmoid, softmax, dsoftmax, pool1d, softplus, dsoftplus, mish, dmish, # exp/trig based functions
       bic, aic,
       autoJacobian,
       squaredCost, dSquaredCost, mse, crossEntropy, dCrossEntropy, classCounts, classCountsWithLabels, meanDicts, mode, gini, entropy, variance,
       error, accuracy, meanRelError, ConfusionMatrix,
       crossValidation, AbstractDataSampler, SamplerWithData, KFold,
       l1_distance,l2_distance, l2²_distance, cosine_distance, lse, sterling,
       #normalFixedSd, logNormalFixedSd,
       radialKernel, polynomialKernel

# Various functions that we add a method to
import Base.print, Base.println, Base.findfirst, Base.findall, Base.error, Random.shuffle

#include("Miscelanneous.jl")
include("Logging_utils.jl")
include("Processing.jl")
include("Stochasticity.jl")
include("Samplers.jl")
include("Transformers.jl")
include("Measures.jl")

end # end module


text:
function backward_pass!(policy::PolicyData, problem::ProblemData;
    mode=:nominal)

    # horizon
    H = length(problem.states)

    # gradients
    fx = problem.model.jacobian_state
    fu = problem.model.jacobian_action
    gx = problem.objective.gradient_state
    gu = problem.objective.gradient_action

    # Hessians
    gxx = problem.objective.hessian_state_state
    guu = problem.objective.hessian_action_action
    gux = problem.objective.hessian_action_state

    # policy
    if mode == :nominal
        K = policy.K
        k = policy.k
    else
        K = policy.K_candidate
        k = policy.k_cand
    end

    # value function approximation
    P = policy.value.hessian
    p = policy.value.gradient

    # state-action value function approximation
    Qx = policy.action_value.gradient_state
    Qu = policy.action_value.gradient_action
    Qxx = policy.action_value.hessian_state_state
    Quu = policy.action_value.hessian_action_action
    Qux = policy.action_value.hessian_action_state

    # terminal value function
    P[H] .= gxx[H]
    p[H] .= gx[H]

    for t = H-1:-1:1
        # Qx[t] .= gx[t] + fx[t]' * p[t+1]
        mul!(Qx[t], transpose(fx[t]), p[t+1])
        Qx[t] .+= gx[t]

        # Qu[t] .= gu[t] + fu[t]' * p[t+1]
        mul!(Qu[t], transpose(fu[t]), p[t+1])
        Qu[t] .+= gu[t]

        # Qxx[t] .= gxx[t] + fx[t]' * P[t+1] * fx[t]
        mul!(policy.xx̂_tmp[t], transpose(fx[t]), P[t+1])
        mul!(Qxx[t], policy.xx̂_tmp[t], fx[t])
        Qxx[t] .+= gxx[t]

        # Quu[t] .= guu[t] + fu[t]' * P[t+1] * fu[t]
        mul!(policy.ux̂_tmp[t], transpose(fu[t]), P[t+1])
        mul!(Quu[t], policy.ux̂_tmp[t], fu[t])
        Quu[t] .+= guu[t]

        # Qux[t] .= gux[t] + fu[t]' * P[t+1] * fx[t]
        mul!(policy.ux̂_tmp[t], transpose(fu[t]), P[t+1])
        mul!(Qux[t], policy.ux̂_tmp[t], fx[t])
        Qux[t] .+= gux[t]

        # K[t] .= -1.0 * Quu[t] \ Qux[t]
        # k[t] .= -1.0 * Quu[t] \ Qu[t]
		policy.uu_tmp[t] .= Quu[t]
        LAPACK.potrf!('U', policy.uu_tmp[t])
        K[t] .= Qux[t]
        k[t] .= Qu[t]
        LAPACK.potrs!('U', policy.uu_tmp[t], K[t])
		LAPACK.potrs!('U', policy.uu_tmp[t], k[t])
		K[t] .*= -1.0
		k[t] .*= -1.0

        # P[t] .=  Qxx[t] + K[t]' * Quu[t] * K[t] + K[t]' * Qux[t] + Qux[t]' * K[t]
        # p[t] .=  Qx[t] + K[t]' * Quu[t] * k[t] + K[t]' * Qu[t] + Qux[t]' * k[t]
		mul!(policy.ux_tmp[t], Quu[t], K[t])

		mul!(P[t], transpose(K[t]), policy.ux_tmp[t])
		mul!(P[t], transpose(K[t]), Qux[t], 1.0, 1.0)
		mul!(P[t], transpose(Qux[t]), K[t], 1.0, 1.0)
		P[t] .+= Qxx[t]

		mul!(p[t], transpose(policy.ux_tmp[t]), k[t])
		mul!(p[t], transpose(K[t]), Qu[t], 1.0, 1.0)
		mul!(p[t], transpose(Qux[t]), k[t], 1.0, 1.0)
		p[t] .+= Qx[t]
    end
end


text:
using Compat

if VERSION.minor == 7
  @eval using LinearAlgebra, SparseMatrices
end

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

A = Float64[10 1 -5;-20 3 20;5 3 5]
b = Float64[1; 2; 6]

(Al, Au, Ap) = lu(A)

@show y = Al \ b[Ap]
println()

@show x = Au \ y
println()

@show As = sparse(A)
println()

@show As \ y
println()

@show F = lufact(As)
println()

@show ys = F[:L] \ (F[:Rs] .* b)
println()

@show x = F[:U] \ ys
println()

@show F[:L]*F[:U] == (F[:Rs] .* A)[F[:p], F[:q]]
println()

@show F[:L]*F[:U]
println()

@show sparse((F[:Rs] .* A)[F[:p], F[:q]])
println()

@show lufact(A)
println()

F = lufact(A)

if VERSION.minor == 7
  F.L * F.U == A[p, :]
else
  F[:L]* F[:U] == A[F[:p], :]
end

end


text:
export TV_cuda

f_inds(rs, b) = ntuple(i -> i == b ? rs[i] .+ 1 : rs[i], length(rs))

"""
    TV_cuda(; num_dims=2)
This function returns a function to calculate the Total Variation regularizer
of a 2 or 3 dimensional array.
`num_dims` can be either `2` or `3`.

```julia-repl
julia> using CUDA

julia> reg = TV_cuda(num_dims=2);

julia> reg(CuArray([1 2 3; 4 5 6; 7 8 9]))
12.649111f0
```
"""
function TV_cuda(; num_dims=2, weights=ones(Float32, num_dims), ϵ=1f-8)
    if num_dims == 3
        return arr -> TV_3D_view(arr, weights, ϵ)
    elseif num_dims == 2
        return arr -> TV_2D_view(arr, weights, ϵ)
    else
        throw(ArgumentError("num_dims must be 2 or 3"))
    end

    return reg_TV
end

function TV_2D_view(arr::AbstractArray{T, N}, weights, ϵ=1f-8) where {T, N}
    as = ntuple(i -> axes(arr, i), Val(N))
    rs = map(x -> first(x):last(x)-1, as)
    arr0 = view(arr, f_inds(rs, 0)...)
    arr1 = view(arr, f_inds(rs, 1)...)
    arr2 = view(arr, f_inds(rs, 2)...)
    return @fastmath sum(sqrt.(ϵ .+ weights[1] .* (arr1 .- arr0).^2 .+ weights[2] .* (arr0 .- arr2).^2))
end

function TV_3D_view(arr::AbstractArray{T, N}, weights, ϵ=1f-8) where {T, N}
    as = ntuple(i -> axes(arr, i), Val(N))
    rs = map(x -> first(x):last(x)-1, as)
    arr0 = view(arr, f_inds(rs, 0)...)
    arr1 = view(arr, f_inds(rs, 1)...)
    arr2 = view(arr, f_inds(rs, 2)...)
    arr3 = view(arr, f_inds(rs, 3)...)
    return @fastmath sum(sqrt.(ϵ .+ weights[1] .* (arr1 .- arr0).^2 .+
                               weights[2] .* (arr2 .- arr0).^2 .+  weights[3] .* (arr3 .- arr0).^2 ))
end


text:
"""
    AbstractPowerRepresentation

An abstract representation type of points and tangent vectors on a power manifold.
"""
abstract type AbstractPowerRepresentation end

"""
    ArrayPowerRepresentation

Representation of points and tangent vectors on a power manifold using multidimensional
arrays where first dimensions are equal to [`representation_size`](@ref) of the
wrapped manifold and the following ones are equal to the number of elements in each
direction.

[`Torus`](@ref) uses this representation.
"""
struct ArrayPowerRepresentation <: AbstractPowerRepresentation end

"""
    NestedPowerRepresentation

Representation of points and tangent vectors on a power manifold using arrays
of size equal to `TSize` of a [`PowerManifold`](@ref).
Each element of such array stores a single point or tangent vector.

[`GraphManifold`](@ref) uses this representation.
"""
struct NestedPowerRepresentation <: AbstractPowerRepresentation end

@doc raw"""
    AbstractPowerManifold{M,TPR} <: Manifold

An abstract [`Manifold`](@ref) to represent manifolds that are build as powers
of another [`Manifold`](@ref) `M` with representation type `TPR`, a subtype of
[`AbstractPowerRepresentation`](@ref).
"""
abstract type AbstractPowerManifold{M<:Manifold,TPR<:AbstractPowerRepresentation} <:
              Manifold end

@doc raw"""
    PowerManifold{TM<:Manifold, TSize<:Tuple, TPR<:AbstractPowerRepresentation} <: AbstractPowerManifold{TM}

The power manifold $\mathcal M^{n_1× n_2 × … × n_d}$ with power geometry
 `TSize` statically defines the number of elements along each axis.

For example, a manifold-valued time series would be represented by a power manifold with
$d$ equal to 1 and $n_1$ equal to the number of samples. A manifold-valued image
(for example in diffusion tensor imaging) would be represented by a two-axis power
manifold ($d=2$) with $n_1$ and $n_2$ equal to width and height of the image.

While the size of the manifold is static, points on the power manifold
would not be represented by statically-sized arrays. Operations on small
power manifolds might be faster if they are represented as [`ProductManifold`](@ref).

# Constructor

    PowerManifold(M, N_1, N_2, ..., N_n)
    PowerManifold(M, NestedPowerRepresentation(), N_1, N_2, ..., N_n)
    M^(N_1,N_2, ..., N_n)

Generate the power manifold $M^{N_1 × N_2 × … × N_n}$.
By default, the [`ArrayPowerRepresentation`](@ref) of points
and tangent vectors is used, although a different one, for example
[`NestedPowerRepresentation`](@ref), can be given as the second argument to the
constructor.
"""
struct PowerManifold{TM<:Manifold,TSize,TPR<:AbstractPowerRepresentation} <:
       AbstractPowerManifold{TM,TPR}
    manifold::TM
end

function PowerManifold(M::Manifold, size::Int...)
    return PowerManifold{typeof(M),Tuple{size...},ArrayPowerRepresentation}(M)
end
function PowerManifold(
    M::Manifold,
    ::TPR,
    size::Int...,
) where {TPR<:AbstractPowerRepresentation}
    PowerManifold{typeof(M),Tuple{size...},TPR}(M)
end

@doc raw"""
    PowerMetric <: Metric

Represent the [`Metric`](@ref) on an [`AbstractPowerManifold`](@ref), i.e. the inner
product on the tangent space is the sum of the inner product of each elements
tangent space of the power manifold.
"""
struct PowerMetric <: Metric end

"""
    PowerRetraction(retraction::AbstractRetractionMethod)

Power retraction based on `retraction`. Works on [`AbstractPowerManifold`](@ref)s.
"""
struct PowerRetraction{TR<:AbstractRetractionMethod} <: AbstractRetractionMethod
    retraction::TR
end

"""
    InversePowerRetraction(inverse_retractions::AbstractInverseRetractionMethod...)

Power inverse retraction of `inverse_retractions`. Works on [`AbstractPowerManifold`](@ref)s.
"""
struct InversePowerRetraction{TR<:AbstractInverseRetractionMethod} <:
       AbstractInverseRetractionMethod
    inverse_retraction::TR
end

"""
    PowerPointDistribution(M::AbstractPowerManifold, distribution)

Power distribution on manifold `M`, based on `distribution`.
"""
struct PowerPointDistribution{TM<:AbstractPowerManifold,TD<:MPointDistribution,TX} <:
       MPointDistribution{TM}
    manifold::TM
    distribution::TD
    point::TX
end

"""
    PowerFVectorDistribution([type::VectorBundleFibers], [x], distr)

Generates a random vector at a `point` from vector space (a fiber of a tangent
bundle) of type `type` using the power distribution of `distr`.

Vector space type and `point` can be automatically inferred from distribution `distr`.
"""
struct PowerFVectorDistribution{
    TSpace<:VectorBundleFibers{<:VectorSpaceType,<:AbstractPowerManifold},
    TD<:FVectorDistribution,
    TX,
} <: FVectorDistribution{TSpace,TX}
    type::TSpace
    point::TX
    distribution::TD
end

"""
    PowerBasisData{TB<:AbstractArray}

Data storage for an array of basis data.
"""
struct PowerBasisData{TB<:AbstractArray}
    bases::TB
end

const POWER_BASIS_LIST_CACHED = [
    CachedBasis{<:AbstractBasis{ℝ},<:PowerBasisData},
    CachedBasis{<:ManifoldsBase.AbstractOrthogonalBasis{ℝ},<:PowerBasisData},
    CachedBasis{<:ManifoldsBase.AbstractOrthonormalBasis{ℝ},<:PowerBasisData},
    CachedBasis{<:AbstractBasis{ℂ},<:PowerBasisData},
]

const PowerManifoldMultidimensional =
    AbstractPowerManifold{<:Manifold,ArrayPowerRepresentation} where {TSize}
const PowerManifoldNested =
    AbstractPowerManifold{<:Manifold,NestedPowerRepresentation} where {TSize}

_access_nested(x, i::Int) = x[i]
_access_nested(x, i::Tuple) = x[i...]

function allocate_result(M::PowerManifoldNested, f, x...)
    return [allocate_result(M.manifold, f, map(y -> _access_nested(y, i), x)...) for i in get_iterator(M)]
end
function allocate_result(M::PowerManifoldNested, f::typeof(flat), w::TFVector, x)
    alloc = [allocate(_access_nested(w.data, i)) for i in get_iterator(M)]
    return FVector(CotangentSpace, alloc)
end
function allocate_result(M::PowerManifoldNested, f::typeof(get_vector), p, X)
    return [allocate_result(M.manifold, f, _access_nested(p, i)) for i in get_iterator(M)]
end
function allocate_result(M::PowerManifoldNested, f::typeof(sharp), w::CoTFVector, x)
    alloc = [allocate(_access_nested(w.data, i)) for i in get_iterator(M)]
    return FVector(TangentSpace, alloc)
end
function allocate_result(M::PowerManifoldNested, f::typeof(get_coordinates), p, X, B)
    return invoke(
        allocate_result,
        Tuple{Manifold,typeof(get_coordinates),Any,Any,typeof(B)},
        M,
        f,
        p,
        X,
        B,
    )
end

function allocation_promotion_function(M::AbstractPowerManifold, f, args::Tuple)
    return allocation_promotion_function(M.manifold, f, args)
end

^(M::Manifold, n) = PowerManifold(M, n...)

"""
    check_manifold_point(M::AbstractProductManifold, p; kwargs...)

Check whether `p` is a valid point on an [`AbstractPowerManifold`](@ref) `M`,
i.e. each element of `p` has to be a valid point on the base manifold.

The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::AbstractPowerManifold, p; kwargs...)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        imp = check_manifold_point(M.manifold, _read(M, rep_size, p, i); kwargs...)
        imp === nothing || return imp
    end
    return nothing
end

"""
    check_tangent_vector(M::AbstractPowerManifold, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to `p` an the [`AbstractPowerManifold`](@ref)
`M`, i.e. atfer [`check_manifold_point`](@ref)`(M, p)`, and all projections to
base manifolds must be respective tangent vectors.
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(
    M::AbstractPowerManifold,
    p,
    X;
    check_base_point = true,
    kwargs...,
)
    if check_base_point
        mpe = check_manifold_point(M, p)
        mpe === nothing || return mpe
    end
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        imp = check_tangent_vector(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i);
            kwargs...,
        )
        imp === nothing || return imp
    end
    return nothing
end

function det_local_metric(
    M::MetricManifold{<:AbstractPowerManifold,PowerMetric},
    p::AbstractArray,
)
    result = one(number_eltype(p))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result *= det_local_metric(M.manifold, _read(M, rep_size, p, i))
    end
    return result
end

@doc raw"""
    distance(M::AbstractPowerManifold, p, q)

Compute the distance between `q` and `p` on an [`AbstractPowerManifold`](@ref),
i.e. from the element wise distances the Forbenius norm is computed.
"""
function distance(M::AbstractPowerManifold, p, q)
    sum_squares = zero(number_eltype(p))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sum_squares +=
            distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))^2
    end
    return sqrt(sum_squares)
end

@doc raw"""
    exp(M::AbstractPowerManifold, p, X)

Compute the exponential map from `p` in direction `X` on the [`AbstractPowerManifold`](@ref) `M`,
which can be computed using the base manifolds exponential map elementwise.
"""
exp(::AbstractPowerManifold, ::Any...)

function exp!(M::AbstractPowerManifold, q, p, X)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        exp!(
            M.manifold,
            _write(M, rep_size, q, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
        )
    end
    return q
end

@doc raw"""
    flat(M::AbstractPowerManifold, p, X::FVector{TangentSpaceType})

use the musical isomorphism to transform the tangent vector `X` from the tangent space at
`p` on an [`AbstractPowerManifold`](@ref) `M` to a cotangent vector.
This can be done elementwise for each entry of `X` (and `p`).
"""
flat(::AbstractPowerManifold, ::Any...)

function flat!(M::AbstractPowerManifold, ξ::CoTFVector, p, X::TFVector)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        flat!(
            M.manifold,
            FVector(CotangentSpace, _write(M, rep_size, ξ.data, i)),
            _read(M, rep_size, p, i),
            FVector(TangentSpace, _read(M, rep_size, X.data, i)),
        )
    end
    return ξ
end

function get_basis(M::AbstractPowerManifold, p, B::AbstractBasis)
    rep_size = representation_size(M.manifold)
    vs = [get_basis(M.manifold, _read(M, rep_size, p, i), B) for i in get_iterator(M)]
    return CachedBasis(B, PowerBasisData(vs))
end
function get_basis(M::AbstractPowerManifold, p, B::DiagonalizingOrthonormalBasis)
    rep_size = representation_size(M.manifold)
    vs = [
        get_basis(
            M.manifold,
            _read(M, rep_size, p, i),
            DiagonalizingOrthonormalBasis(_read(M, rep_size, B.frame_direction, i)),
        ) for i in get_iterator(M)
    ]
    return CachedBasis(B, PowerBasisData(vs))
end
for BT in ManifoldsBase.DISAMBIGUATION_BASIS_TYPES
    if BT == DiagonalizingOrthonormalBasis
        continue
    end
    eval(quote
        @invoke_maker 3 AbstractBasis get_basis(M::AbstractPowerManifold, p, B::$BT)
    end)
end

function get_coordinates(M::AbstractPowerManifold, p, X, B::DefaultOrthonormalBasis)
    rep_size = representation_size(M.manifold)
    vs = [
        get_coordinates(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            B,
        ) for i in get_iterator(M)
    ]
    return reduce(vcat, reshape(vs, length(vs)))
end
function get_coordinates(
    M::AbstractPowerManifold,
    p,
    X,
    B::CachedBasis{<:AbstractBasis,<:PowerBasisData,𝔽},
) where {𝔽}
    rep_size = representation_size(M.manifold)
    vs = [
        get_coordinates(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            _access_nested(B.data.bases, i),
        ) for i in get_iterator(M)
    ]
    return reduce(vcat, reshape(vs, length(vs)))
end
for BT in POWER_BASIS_LIST_CACHED
    eval(quote
        @invoke_maker 4 CachedBasis{<:AbstractBasis,<:PowerBasisData} get_coordinates(
            M::AbstractPowerManifold,
            p,
            X,
            B::$BT,
        )
    end)
end

function get_coordinates!(M::AbstractPowerManifold, Y, p, X, B::DefaultOrthonormalBasis)
    rep_size = representation_size(M.manifold)
    dim = manifold_dimension(M.manifold)
    v_iter = 1
    for i in get_iterator(M)
        # TODO: this view is really suboptimal when `dim` can be statically determined
        get_coordinates!(
            M.manifold,
            view(Y, v_iter:v_iter+dim-1),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            B,
        )
        v_iter += dim
    end
    return Y
end
function get_coordinates!(
    M::AbstractPowerManifold,
    Y,
    p,
    X,
    B::CachedBasis{<:AbstractBasis{ℝ},<:PowerBasisData,ℝ},
)
    TypeTuple = Tuple{
        AbstractPowerManifold,
        Any,
        Any,
        Any,
        CachedBasis{<:AbstractBasis,<:PowerBasisData,ℝ},
    }
    return invoke(get_coordinates!, TypeTuple, M, Y, p, X, B)
end
function get_coordinates!(
    M::AbstractPowerManifold,
    Y,
    p,
    X,
    B::CachedBasis{<:AbstractBasis,<:PowerBasisData,𝔽},
) where {𝔽}
    rep_size = representation_size(M.manifold)
    dim = manifold_dimension(M.manifold)
    v_iter = 1
    for i in get_iterator(M)
        get_coordinates!(
            M.manifold,
            view(Y, v_iter:v_iter+dim-1),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            _access_nested(B.data.bases, i),
        )
        v_iter += dim
    end
    return Y
end

get_iterator(M::PowerManifold{<:Manifold,Tuple{N}}) where {N} = 1:N
@generated function get_iterator(M::PowerManifold{<:Manifold,SizeTuple}) where {SizeTuple}
    size_tuple = size_to_tuple(SizeTuple)
    return Base.product(map(Base.OneTo, size_tuple)...)
end

function get_vector!(
    M::PowerManifold,
    Y,
    p,
    X,
    B::CachedBasis{<:AbstractBasis,<:PowerBasisData},
)
    dim = manifold_dimension(M.manifold)
    rep_size = representation_size(M.manifold)
    v_iter = 1
    for i in get_iterator(M)
        get_vector!(
            M.manifold,
            _write(M, rep_size, Y, i),
            _read(M, rep_size, p, i),
            X[v_iter:v_iter+dim-1],
            _access_nested(B.data.bases, i),
        )
        v_iter += dim
    end
    return Y
end
function get_vector!(M::AbstractPowerManifold, Y, p, X, B::DefaultOrthonormalBasis)
    dim = manifold_dimension(M.manifold)
    rep_size = representation_size(M.manifold)
    v_iter = 1
    for i in get_iterator(M)
        get_vector!(
            M.manifold,
            _write(M, rep_size, Y, i),
            _read(M, rep_size, p, i),
            X[v_iter:v_iter+dim-1],
            B,
        )
        v_iter += dim
    end
    return Y
end

@doc raw"""
    injectivity_radius(M::AbstractPowerManifold[, p])

the injectivity radius on an [`AbstractPowerManifold`](@ref) is for the global case
equal to the one of its base manifold. For a given point `p` it's equal to the
minimum of all radii in the array entries.
"""
function injectivity_radius(M::AbstractPowerManifold, p)
    radius = 0.0
    initialized = false
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        cur_rad = injectivity_radius(M.manifold, _read(M, rep_size, p, i))
        if initialized
            radius = min(cur_rad, radius)
        else
            radius = cur_rad
            initialized = true
        end
    end
    return radius
end
injectivity_radius(M::AbstractPowerManifold) = injectivity_radius(M.manifold)
eval(quote
    @invoke_maker 1 Manifold injectivity_radius(
        M::AbstractPowerManifold,
        rm::AbstractRetractionMethod,
    )
end)
eval(quote
    @invoke_maker 1 Manifold injectivity_radius(
        M::AbstractPowerManifold,
        rm::ExponentialRetraction,
    )
end)

@doc raw"""
    inverse_retract(M::AbstractPowerManifold, p, q, m::InversePowerRetraction)

Compute the inverse retraction from `p` with respect to `q` on an [`AbstractPowerManifold`](@ref) `M`
using an [`InversePowerRetraction`](@ref), which by default encapsulates a inverse retraction
of the base manifold. Then this method is performed elementwise, so the encapsulated inverse
retraction method has to be one that is available on the base [`Manifold`](@ref).
"""
inverse_retract(::AbstractPowerManifold, ::Any...)

function inverse_retract!(M::AbstractPowerManifold, X, p, q, method::InversePowerRetraction)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        inverse_retract!(
            M.manifold,
            _write(M, rep_size, X, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, q, i),
            method.inverse_retraction,
        )
    end
    return X
end

@doc raw"""
    inner(M::AbstractPowerManifold, p, X, Y)

Compute the inner product of `X` and `Y` from the tangent space at `p` on an
[`AbstractPowerManifold`](@ref) `M`, i.e. for each arrays entry the tangent
vector entries from `X` and `Y` are in the tangent space of the corresponding
element from `p`.
The inner product is then the sum of the elementwise inner products.
"""
function inner(M::AbstractPowerManifold, p, X, Y)
    result = zero(number_eltype(X))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result += inner(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            _read(M, rep_size, Y, i),
        )
    end
    return result
end

default_metric_dispatch(::AbstractPowerManifold, ::PowerMetric) = Val(true)

function isapprox(M::AbstractPowerManifold, p, q; kwargs...)
    result = true
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result &= isapprox(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, q, i);
            kwargs...,
        )
    end
    return result
end
function isapprox(M::AbstractPowerManifold, p, X, Y; kwargs...)
    result = true
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result &= isapprox(
            M.manifold,
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            _read(M, rep_size, Y, i);
            kwargs...,
        )
    end
    return result
end

@doc raw"""
    log(M::AbstractPowerManifold, p, q)

Compute the logarithmic map from `p` to `q` on the [`AbstractPowerManifold`](@ref) `M`,
which can be computed using the base manifolds logarithmic map elementwise.
"""
log(::AbstractPowerManifold, ::Any...)

function log!(M::AbstractPowerManifold, X, p, q)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        log!(
            M.manifold,
            _write(M, rep_size, X, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, q, i),
        )
    end
    return X
end

@doc raw"""
    manifold_dimension(M::PowerManifold)

Returns the manifold-dimension of an [`PowerManifold`](@ref) `M`
$=\mathcal N = (\mathcal M)^{n_1,…,n_d}$, i.e. with $n=(n_1,…,n_d)$ the array
size of the power manifold and $d_{\mathcal M}$ the dimension of the base manifold
$\mathcal M$, the manifold is of dimension

````math
\dim(\mathcal N) = \dim(\mathcal M)\prod_{i=1}^d n_i = n_1n_2\cdot…\cdot n_d \dim(\mathcal M).
````
"""
function manifold_dimension(M::PowerManifold{<:Manifold,TSize}) where {TSize}
    return manifold_dimension(M.manifold) * prod(size_to_tuple(TSize))
end

@doc raw"""
    norm(M::AbstractPowerManifold, p, X)

Compute the norm of `X` from the tangent space of `p` on an
[`AbstractPowerManifold`](@ref) `M`, i.e. from the element wise norms the
Frobenius norm is computed.
"""
function norm(M::AbstractPowerManifold, p, X)
    sum_squares = zero(number_eltype(X))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sum_squares +=
            norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))^2
    end
    return sqrt(sum_squares)
end

@doc raw"""
    power_dimensions(M::PowerManifold)

return the power of `M`,
"""
function power_dimensions(M::PowerManifold{<:Manifold,TSize}) where {TSize}
    return size_to_tuple(TSize)
end

function rand(rng::AbstractRNG, d::PowerFVectorDistribution)
    fv = zero_vector(d.type, d.point)
    _rand!(rng, d, fv)
    return fv
end
function rand(rng::AbstractRNG, d::PowerPointDistribution)
    x = allocate_result(d.manifold, rand, d.point)
    _rand!(rng, d, x)
    return x
end

function _rand!(rng::AbstractRNG, d::PowerFVectorDistribution, v::AbstractArray)
    PM = d.type.manifold
    rep_size = representation_size(PM.manifold)
    for i in get_iterator(d.type.manifold)
        copyto!(d.distribution.point, _read(PM, rep_size, d.point, i))
        _rand!(rng, d.distribution, _read(PM, rep_size, v, i))
    end
    return v
end
function _rand!(rng::AbstractRNG, d::PowerPointDistribution, x::AbstractArray)
    M = d.manifold
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        _rand!(rng, d.distribution, _write(M, rep_size, x, i))
    end
    return x
end

@inline function _read(M::AbstractPowerManifold, rep_size::Tuple, x::AbstractArray, i::Int)
    return _read(M, rep_size, x, (i,))
end
@inline function _read(
    ::PowerManifoldMultidimensional,
    rep_size::Tuple,
    x::AbstractArray,
    i::Tuple,
)
    return view(x, rep_size_to_colons(rep_size)..., i...)
end
@inline function _read(
    ::PowerManifoldMultidimensional,
    rep_size::Tuple,
    x::HybridArray,
    i::Tuple,
)
    return x[rep_size_to_colons(rep_size)..., i...]
end
@inline function _read(::PowerManifoldNested, rep_size::Tuple, x::AbstractArray, i::Tuple)
    return view(x[i...], rep_size_to_colons(rep_size)...)
end

@generated function rep_size_to_colons(rep_size::Tuple)
    N = length(rep_size.parameters)
    return ntuple(i -> Colon(), N)
end

function representation_size(M::PowerManifold{<:Manifold,TSize}) where {TSize}
    return (representation_size(M.manifold)..., size_to_tuple(TSize)...)
end

@doc raw"""
    retract(M::AbstractPowerManifold, p, X, method::PowerRetraction)

Compute the retraction from `p` with tangent vector `X` on an [`AbstractPowerManifold`](@ref) `M`
using a [`PowerRetraction`](@ref), which by default encapsulates a retraction of the
base manifold. Then this method is performed elementwise, so the encapsulated retraction
method has to be one that is available on the base [`Manifold`](@ref).
"""
retract(::AbstractPowerManifold, ::Any...)

function retract!(M::AbstractPowerManifold, q, p, X, method::PowerRetraction)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        retract!(
            M.manifold,
            _write(M, rep_size, q, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            method.retraction,
        )
    end
    return q
end

@doc raw"""
    sharp(M::AbstractPowerManifold, p, ξ::FVector{CotangentSpaceType})

Use the musical isomorphism to transform the cotangent vector `ξ` from the tangent space at
`p` on an [`AbstractPowerManifold`](@ref) `M` to a tangent vector.
This can be done elementwise for every entry of `ξ` (and `p`).
"""
sharp(::AbstractPowerManifold, ::Any...)

function sharp!(M::AbstractPowerManifold, X::TFVector, p, ξ::CoTFVector)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sharp!(
            M.manifold,
            FVector(TangentSpace, _write(M, rep_size, X.data, i)),
            _read(M, rep_size, p, i),
            FVector(CotangentSpace, _read(M, rep_size, ξ.data, i)),
        )
    end
    return X
end

function show(
    io::IO,
    mime::MIME"text/plain",
    B::CachedBasis{T,D,𝔽},
) where {T<:AbstractBasis,D<:PowerBasisData,𝔽}
    println(io, "$(T()) for a power manifold")
    for i in Base.product(map(Base.OneTo, size(B.data.bases))...)
        println(io, "Basis for component $i:")
        show(io, mime, _access_nested(B.data.bases, i))
        println(io)
    end
end
function show(io::IO, M::PowerManifold{TM,TSize,ArrayPowerRepresentation}) where {TM,TSize}
    print(io, "PowerManifold($(M.manifold), $(join(TSize.parameters, ", ")))")
end
function show(io::IO, M::PowerManifold{TM,TSize,TPR}) where {TM,TSize,TPR}
    print(io, "PowerManifold($(M.manifold), $(TPR()), $(join(TSize.parameters, ", ")))")
end

support(tvd::PowerFVectorDistribution) = FVectorSupport(tvd.type, tvd.point)
support(d::PowerPointDistribution) = MPointSupport(d.manifold)

@inline function _write(M::AbstractPowerManifold, rep_size::Tuple, x::AbstractArray, i::Int)
    return _write(M, rep_size, x, (i,))
end
@inline function _write(
    M::PowerManifoldMultidimensional,
    rep_size::Tuple,
    x::AbstractArray,
    i::Tuple,
)
    return view(x, rep_size_to_colons(rep_size)..., i...)
end
@inline function _write(M::PowerManifoldNested, rep_size::Tuple, x::AbstractArray, i::Tuple)
    return view(x[i...], rep_size_to_colons(rep_size)...)
end


text:
(s::Symbol)(args...) = WExpr(buildexpr(s), buildexpr.(args))

function buildwith(head, a, b)
    if isa(a, WExpr) && a.head == WSymbol(head)
        WExpr(WSymbol(head), (a.args..., buildexpr(b)))
    else
        WExpr(WSymbol(head), buildexpr.((a, b)))
    end
end

const WTypes = Union{WExpr, WSymbol, MathLink.WInteger, MathLink.WReal, Symbol}

#
# Arithematic operators
#

Base.:+(a::Union{WTypes, T}, b::Union{WTypes, T}) where T <: Number =
    buildwith("Plus", a, b)

Base.:-(a::WTypes) =
    WSymbol("Times")(-1, buildexpr(a))

Base.:-(a::Union{WTypes, T}, b::Union{WTypes, T}) where T <: Number =
    a + (-b)

Base.:*(a::Union{WTypes, T}, b::Union{WTypes, T}) where T <: Number =
    buildwith("Times", a, b)

Base.:/(a::Union{WTypes, T}, b::Union{WTypes, T}) where T <: Number =
    a * WSymbol("Power")(buildexpr(b), -1)

Base.:^(a::Union{WTypes, T}, b::Union{WTypes, T}) where T <: Number =
    WSymbol("Power")(buildexpr(a), buildexpr(b))

#
# Functions
#

Base.:|>(x, f::Symbol) = buildexpr(f)(buildexpr(x))
Base.inv(x::WTypes) = WSymbol("Inverse")(buildexpr(x))


text:
using DataStructures, OPFInstances

function main()

    #############################################################
    # Parameters
    nmax_order1 = 350
    nmax_order2 = 82
    juliapath = joinpath("~", "software", "julia-0.6.3", "bin", "julia")
    juliascript = joinpath(pwd(), "dev", "script_dat_to_hierarchysimple.jl")
    qsubwkdir = pwd()

    date = String(Dates.format(now(), "mm_dd-HHhMM"))
    workdir = joinpath(pwd(), "Mosek_runs", "pararuns", date)

    ispath(workdir) && rm(workdir, recursive=true)
    mkpath(workdir)

    params = OrderedSet()
    tempfolders = Set()

    matpower_path = first(splitdir(getinstancepath("Matpower", "QCQP", "WB2")))
    instances = OrderedSet(sort(readdir(matpower_path), by=x->parse(first(matchall(r"\d+", x)))))

    instances_order1 = filter(x->parse(first(matchall(r"\d+", x))) ≤ nmax_order1, instances)
    instances_order2 = filter(x->parse(first(matchall(r"\d+", x))) ≤ nmax_order2, instances)

    for instance in instances_order1
        tempfolder = randstring(4)
        while tempfolder in tempfolders
            tempfolder = randstring(4)
        end
        push!(tempfolders, tempfolder)

        instancename = first(split(instance, "."))
        push!(params, (instancename, joinpath(matpower_path, instance), 1, joinpath(workdir, tempfolder), 0))
    end

    for instance in instances_order1
        tempfolder = randstring(4)
        while tempfolder in tempfolders
            tempfolder = randstring(4)
        end
        push!(tempfolders, tempfolder)

        instancename = first(split(instance, "."))
        push!(params, (instancename, joinpath(matpower_path, instance), 1, joinpath(workdir, tempfolder), 1))
    end

    for instance in instances_order2
        tempfolder = randstring(4)
        while tempfolder in tempfolders
            tempfolder = randstring(4)
        end
        push!(tempfolders, tempfolder)

        instancename = first(split(instance, "."))
        push!(params, (instancename, joinpath(matpower_path, instance), 2, joinpath(workdir, tempfolder), 0))
    end

    for instance in instances_order2
        tempfolder = randstring(4)
        while tempfolder in tempfolders
            tempfolder = randstring(4)
        end
        push!(tempfolders, tempfolder)

        instancename = first(split(instance, "."))
        push!(params, (instancename, joinpath(matpower_path, instance), 2, joinpath(workdir, tempfolder), 1))
    end

    for param in params
        println(param)
    end

    sleep(3)

    for param in params
        instancename, instancepath, d, logpath, hierarchykind = param
        # println(`qsub -e $logpath/hierarchy_$(instancename).e -o $logpath/hierarchy_$(instancename).o $juliapath dev/script_dat_to_hierarchysimple.jl $instancename $d $logpath`)
        # run(`qsub -e $logpath/hierarchy_$(instancename).e -o $logpath/hierarchy_$(instancename).o -b y julia dev/script_dat_to_hierarchysimple.jl $instancepath $d $logpath`)
        mkpath(logpath)
        open(joinpath(logpath, "runscript.sh"), "w") do f
            println(f, "cd $(pwd())")
            println(f, "qsub -e $logpath/hierarchy_$(instancename).e -o $logpath/hierarchy_$(instancename).o -wd $qsubwkdir -b y $juliapath $juliascript $instancepath $d $logpath $hierarchykind")
        end

        run(`bash $(joinpath(logpath, "runscript.sh"))`)
    end

end

main()


text:
"""
Includes all useful functions for applying GPR to T and wT profiles from Oceananigans.jl simulations.
Uses ProfileData struct to store data and GP struct for performing GPR on the data in ProfileData object.
"""
module GaussianProcess

using LearnConvection.Data

include("kernels.jl")
export  Kernel,
        SquaredExponentialI,
        RationalQuadraticI,
        Matern12I,
        Matern32I,
        Matern52I,
        SpectralMixtureProductI,
        SpectralMixtureProductA,
        kernel_function

include("distances.jl")
export  euclidean_distance,
        derivative_distance,
        antiderivative_distance

include("gp.jl")
export  model_output,
        uncertainty,
        compute_kernel_matrix,
        mean_log_marginal_loss
export  model
        # predict

export get_kernel

"""
```get_kernel(kernel_id::Int64, logγ, logσ, distance; logα=0.0)```

# Description
Returns a Kernel object with the specified parameters.

# Arguments
- `kernel_id`: (Int64). Identification number for the kernel type (see kernel options)
- `logγ`: (Float64). Log(length scale) parameter.
- `logσ`: (Float64). Log(signal variance) parameter.
- `distance`: (Function). Distance function to use in the kernel.

# Keyword Arguments
- `logα`: (Float64). Log(α) parameter if kernel_id = 5.

# kernel options
 1   =>   "Squared exponential kernel:        k(x,x') = σ * exp( - d(x,x')² / 2γ² )",
 2   =>   "Matérn with ʋ=1/2:                 k(x,x') = σ * exp( - d(x,x') / γ )",
 3   =>   "Matérn with ʋ=3/2:                 k(x,x') = σ * (1+c) * exp(-√(3)*d(x,x'))/γ)",
 4   =>   "Matérn with ʋ=5/2:                 k(x,x') = σ * ( 1 + √(5)*d(x,x'))/γ + 5*d(x,x')²/(3*γ^2) ) * exp(-√(5)*d(x,x'))/γ)",
 5   =>   "Rational quadratic kernel:         k(x,x') = σ * (1+(x-x')'*(x-x')/(2*α*(γ²))^(-α)",
"""
# function get_kernel(kernel_id::Int64, logγ, logσ, distance; logα=0.0)
#         # convert from log10 scale
#         γ = 10^logγ
#         σ = 10^logσ
#         α = 10^logα
#         if kernel_id==1; return SquaredExponentialI(γ, σ, distance) end
#         if kernel_id==2; return Matern12I(γ, σ, distance) end
#         if kernel_id==3; return Matern32I(γ, σ, distance) end
#         if kernel_id==4; return Matern52I(γ, σ, distance) end
#         if kernel_id==5; return RationalQuadraticI(γ, σ, α, distance)
#         else; throw(error()) end
# end

function get_kernel(kernel_id::Int64, args...)
        # args = f.(args) # convert from log10 scale

        function g(x)
                if typeof(x) <: Number
                        return 10^x end
                return x
        end

        args = g.(args)

        if kernel_id==1; return SquaredExponentialI(args...) end
        if kernel_id==2; return Matern12I(args...) end
        if kernel_id==3; return Matern32I(args...) end
        if kernel_id==4; return Matern52I(args...) end
        if kernel_id==5; return RationalQuadraticI(args...) end
        if kernel_id==6; return SpectralMixtureProductI(args...)
        else; throw(error()) end
end

end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection