text:
# Aminoacids
# ==========
#
# The amino acid type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"Type representing AminoAcids"
bitstype 8 AminoAcid

# Conversion from/to integers
# ---------------------------

Base.convert(::Type{AminoAcid}, aa::UInt8) = reinterpret(AminoAcid, aa)
Base.convert(::Type{UInt8}, aa::AminoAcid) = reinterpret(UInt8, aa)
Base.convert{T<:Number}(::Type{T}, aa::AminoAcid) = convert(T, UInt8(aa))
Base.convert{T<:Number}(::Type{AminoAcid}, aa::T) = convert(AminoAcid, UInt8(aa))

# Arithmetic and Order
# --------------------

# These methods are necessary when deriving some algorithims
# like iteration, sort, comparison, and so on.
@compat begin
    Base.:-(x::AminoAcid, y::AminoAcid) = Int(x) - Int(y)
    # 0x1c is the size of the amino acid alphabet
    Base.:-(x::AminoAcid, y::Integer) = x + mod(-y, 0x1c)
    Base.:+(x::AminoAcid, y::Integer) = reinterpret(AminoAcid, mod((UInt8(x) + y) % UInt8, 0x1c))
    Base.:+(x::Integer, y::AminoAcid) = y + x
end
Base.isless(x::AminoAcid, y::AminoAcid) = isless(UInt8(x), UInt8(y))

# Amino acid encoding definition
# ------------------------------

"Invalid Amino Acid"
const AA_INVALID = convert(AminoAcid, 0x1c)  # Used during conversion from strings

# lookup table for characters
const char_to_aa = [AA_INVALID for _ in 0x00:0x7f]
const aa_to_char = Vector{Char}(0x1c)

# compatibility bits
const compatbits_aa = Vector{UInt32}(28)

# This set of amino acids is defined by IUPAC-IUB Joint Commission on Biochemical Nomenclature.
# Reference: http://www.insdc.org/documents/feature_table.html#7.4.3

for (aa, doc, code) in [
        ('A', "Alanine",                           0x00),
        ('R', "Arginine",                          0x01),
        ('N', "Asparagine",                        0x02),
        ('D', "Aspartic Acid",                     0x03),
        ('C', "Cysteine",                          0x04),
        ('Q', "Glutamine",                         0x05),
        ('E', "Glutamic Acid",                     0x06),
        ('G', "Glycine",                           0x07),
        ('H', "Histidine",                         0x08),
        ('I', "Isoleucine",                        0x09),
        ('L', "Leucine",                           0x0a),
        ('K', "Lysine",                            0x0b),
        ('M', "Methionine",                        0x0c),
        ('F', "Phenylalanine",                     0x0d),
        ('P', "Proline",                           0x0e),
        ('S', "Serine",                            0x0f),
        ('T', "Threonine",                         0x10),
        ('W', "Tryptophan",                        0x11),
        ('Y', "Tyrosine",                          0x12),
        ('V', "Valine",                            0x13),
        ('O', "Pyrrolysine",                       0x14),  # non-standard
        ('U', "Selenocysteine",                    0x15),  # non-standard
        ('B', "Aspartic Acid or Asparagine",       0x16),  # ambiguous
        ('J', "Leucine or Isoleucine",             0x17),  # ambiguous
        ('Z', "Glutamine or Glutamic Acid",        0x18),  # ambiguous
        ('X', "Unspecified or Unknown Amino Acid", 0x19)]  # ambiguous
    var = Symbol("AA_", aa)
    @eval begin
        @doc $doc const $var = convert(AminoAcid, $code)
        char_to_aa[$(Int(aa)+1)] = char_to_aa[$(Int(lowercase(aa))+1)] = $var
        aa_to_char[$(code)+1] = $aa
    end
    if code ≤ 0x15
        compatbits_aa[code+1] = 1 << code
    elseif code == 0x16
        compatbits_aa[code+1] = 1 << 0x02 | 1 << 0x03
    elseif code == 0x17
        compatbits_aa[code+1] = 1 << 0x09 | 1 << 0x0a
    elseif code == 0x18
        compatbits_aa[code+1] = 1 << 0x05 | 1 << 0x06
    elseif code == 0x19
        compatbits_aa[code+1] = (1 << 0x16) - 1
    end
end

"Terminal"
const AA_Term = convert(AminoAcid, 0x1a)
char_to_aa[Int('*')+1] = AA_Term
aa_to_char[0x1a+1] = '*'
compatbits_aa[0x1a+1] = 1 << 0x1a

"Amino Acid Gap"
const AA_Gap = convert(AminoAcid, 0x1b)
char_to_aa[Int('-') + 1] = AA_Gap
aa_to_char[0x1b+1] = '-'
compatbits_aa[0x1b+1] = 0

Base.colon(start::AminoAcid, stop::AminoAcid) = SymbolRange(start, stop)

Base.isvalid(::Type{AminoAcid}, x::Integer) = 0 ≤ x ≤ 0x1b
Base.isvalid(aa::AminoAcid) = aa ≤ AA_Gap
isambiguous(aa::AminoAcid) = AA_B ≤ aa ≤ AA_X
alphabet(::Type{AminoAcid}) = AA_A:AA_Gap
gap(::Type{AminoAcid}) = AA_Gap

compatbits(aa::AminoAcid) = compatbits_aa[reinterpret(UInt8, aa)+1]

# Conversion from/to Char
# -----------------------

function Base.convert(::Type{AminoAcid}, c::Char)
    @inbounds aa = c <= '\x7f' ? char_to_aa[Int(c)+1] : AA_INVALID
    @assert aa != AA_INVALID error("$(c) is not a valid amino acid")
    return aa
end

Base.convert(::Type{Char}, aa::AminoAcid) = aa_to_char[convert(UInt8, aa) + 1]

# Basic functions
# ---------------

function Base.show(io::IO, aa::AminoAcid)
    if isvalid(aa)
        if aa == AA_Term
            write(io, "AA_Term")
        elseif aa == AA_Gap
            write(io, "AA_Gap")
        else
            write(io, "AA_", Char(aa))
        end
    else
        write(io, "Invalid Amino Acid")
    end
    return
end

function Base.print(io::IO, aa::AminoAcid)
    if !isvalid(aa)
        throw(ArgumentError("invalid amino acid"))
    end
    write(io, Char(aa))
    return
end

# lookup table of 20 standard amino acids
const threeletter_to_aa = Dict(
    "ALA" => AA_A, "ARG" => AA_R, "ASN" => AA_N, "ASP" => AA_D, "CYS" => AA_C,
    "GLN" => AA_Q, "GLU" => AA_E, "GLY" => AA_G, "HIS" => AA_H, "ILE" => AA_I,
    "LEU" => AA_L, "LYS" => AA_K, "MET" => AA_M, "PHE" => AA_F, "PRO" => AA_P,
    "SER" => AA_S, "THR" => AA_T, "TRP" => AA_W, "TYR" => AA_Y, "VAL" => AA_V,
    "ASX" => AA_B, "XLE" => AA_J, "GLX" => AA_Z, "XAA" => AA_X,
    "PYL" => AA_O, "SEC" => AA_U,
)

function Base.parse(::Type{AminoAcid}, s::AbstractString)
    s′ = strip(s)
    if length(s′) == 1
        return convert(AminoAcid, s′[1])
    end
    try
        return threeletter_to_aa[uppercase(s′)]
    catch ex
        if isa(ex, KeyError)
            error("invalid amino acid string: \"$s\" ")
        end
        rethrow()
    end
end

# TODO: tryparse


text:
abstract type TreeIndex{depth} end

Base.length(::TreeIndex{depth}) where {depth} = depth
Base.getindex(index::TreeIndex, i::Int) = (@_propagate_inbounds_meta; index.I[i])
Base.firstindex(index::TreeIndex) = 1
Base.lastindex(index::TreeIndex) = length(index)

TreeIndex(I::Tuple{Vararg{Int}}) = TreeLinearIndex(I)
TreeIndex(I::Tuple{Vararg{CartesianIndex}}) = TreeCartesianIndex(I)
TreeIndex(I...) = TreeIndex(I)

@inline _compute_offsets(::Val{1}, offset_func, x::TreeSize, i::Integer...) = (offset_func(x, i...),)
@inline _compute_offsets(::Val{N}, offset_func, x::TreeSize, i::Integer...) where {N} =
    (offset_func(x, i...), _compute_offsets(Val(N-1), offset_func, Base.tail(x), i...)...)
@inline compute_offsets(offset_func, x::TreeSize, i::Integer...) =
    _compute_offsets(Val(length(x)), offset_func, x, i...)

struct TreeLinearIndex{depth} <: TreeIndex{depth}
    I::NTuple{depth, Int}
end
TreeLinearIndex(I::Int...) = TreeLinearIndex(I)

# offset_linear
@inline function offset_linear(S::TreeSize, i::Vararg{Integer, dim}) where {dim}
    Base._sub2ind(S[1], Tuple(offset_cartesian(S, i...))...)
end
@inline function offset_linear(node::Node, I::Integer...)
    Base._sub2ind(size(node), Tuple(offset_cartesian(node, i...))...)
end

# TreeLinearIndex
@inline TreeLinearIndex(S::TreeSize, i::Integer...) = TreeLinearIndex(compute_offsets(offset_linear, S, i...))
@inline TreeLinearIndex(node::Union{Node, HashNode}, inds::Integer...) = TreeLinearIndex(TreeSize(node), inds...)
function TreeLinearIndex(node::Union{DynamicNode, DynamicHashNode}, inds::Integer...)
    t_child = childtype(node)
    tsize_child = TreeSize(t_child)
    i = Base._sub2ind(size(node), block_index(totalsize(tsize_child), inds...)...)
    I = TreeLinearIndex(tsize_child, inds...).I
    TreeLinearIndex(i, I...)
end

struct TreeCartesianIndex{depth, N} <: TreeIndex{depth}
    I::NTuple{depth, CartesianIndex{N}}
end
TreeCartesianIndex(I::CartesianIndex...) = TreeCartesianIndex(I)

# offset_cartesian
@inline function offset_cartesian(S::TreeSize, I::Integer...)
    dims = totalsize(S)
    dims_child = totalsize(Base.tail(S))
    CartesianIndex(block_index(dims_child, (@. rem(I - 1, dims) + 1)...))
end
@inline offset_cartesian(node::Union{Node, HashNode}, I::Integer...) = offset_cartesian(TreeSize(node), I...)
@inline function offset_cartesian(node::Union{DynamicNode, DynamicHashNode}, I::Integer...)
    @boundscheck checkbounds(CartesianIndices(totalsize(node)), I...)
    dims_child = totalsize(Base.tail(TreeSize(node)))
    CartesianIndex(block_index(dims_child, I...))
end

# TreeCartesianIndex
@inline TreeCartesianIndex(S::TreeSize, i::Integer...) = TreeCartesianIndex(compute_offsets(offset_cartesian, S, i...))
@inline TreeCartesianIndex(node::Union{Node, HashNode}, inds::Integer...) = TreeCartesianIndex(TreeSize(node), inds...)
function TreeCartesianIndex(node::Union{DynamicNode, DynamicHashNode}, inds::Integer...)
    t_child = childtype(node)
    tsize_child = TreeSize(t_child)
    i = CartesianIndex(block_index(totalsize(tsize_child), inds...))
    I = TreeCartesianIndex(tsize_child, inds...).I
    TreeCartesianIndex(i, I...)
end

# totalsize
function totalsize(node::Union{Node, HashNode})
    convert.(Int, totalsize(TreeSize(node)))
end
function totalsize(node::Union{DynamicNode, DynamicHashNode})
    dims = size(node) # don't statically get size to handle the case that rootnode is dynamic node
    dims_child = totalsize(TreeSize(childtype(node)))
    dims .* dims_child
end

ext: jl


text:
# nested package

include("../packages/DeeplyNested/src/DeeplyNested.jl")

@test isdefined(DeeplyNested, :foo)

@test DeeplyNested.foo() == 123

@test isdefined(DeeplyNested, :bar)

@test DeeplyNested.bar() == 404

@test isdefined(DeeplyNested, :baz)

@test DeeplyNested.baz() == 999


text:
module UnpackUInt12s

    using SIMD
    using ..UInt24s
    import ..UInt12ArraysBase: UInt12Array, map_idx_to_byte

    include( joinpath("unpack", "unpack_simd_256.jl") )
    include( joinpath("unpack", "unpack_simd_512.jl") )

    # 64 megabtyes when initialized
    const LUT = Ref{Vector{UInt32}}()

    make_lut(x) = make_lut(UInt32(x))
    function make_lut(x::UInt32)
        ( (x & 0xfff000) << 4) | (x & 0x000fff)
    end
    function init_lut()
        if !isdefined(LUT,1) || isempty(LUT[])
            LUT[] = make_lut.( UInt32(0):UInt32(2^24-1) )
        end
    end
    function clear_lut()
        empty!(LUT[])
    end
    function lutConvertToUInt16(A)
        A = reinterpret(UInt24, A)
        init_lut()
        result = LUT[][A .+ 1]
        reinterpret(UInt16, result)
    end

    function Base.convert(::Type{Array{UInt16,N}}, A::UInt12Array{UInt16,B,N}) where {B <: SIMD.FastContiguousArray{UInt8,1}, N}
        len = length(A)
        if len < 64
            if len < 16 || mod(len, 2) == 1
                @debug "Using copy" len
                return copy(A)
            else
                @debug "Using LUT" len
                return lutConvertToUInt16(A.data)
            end
        else
            @debug "Using SIMD" len
            return unpack_uint12_to_uint16(A.data)
        end
    end
    function Base.convert(::Type{Array{UInt16,N}}, A::UInt12Array{UInt16,B,N}) where {B, N}
        len = length(A)
        if len < 16 || mod(len, 2) == 1
            @debug "Using copy" len
            return copy(A)
        else
            @debug "Using LUT" len
            return lutConvertToUInt16(A.data)
        end
    end

    Base.convert(::Type{Array{UInt16}}, A::UInt12Array{UInt16,B,N}) where {B <: SIMD.FastContiguousArray{UInt8,1}, N} =
        Base.convert(Array{UInt16,N}, A)

    function Base.convert(::Type{Array{UInt16,1}}, S::SubArray{UInt16, 1, UInt12Array{UInt16, B, 1}, <: Tuple{UnitRange}}) where {B <: SIMD.FastContiguousArray{UInt8,1}}
        try
            if S |> parentindices |> first |> first |> isodd
                convert(Vector{UInt16}, convert(UInt12Array{UInt16}, S))
            elseif length(S) < 16
                Array{UInt16,1}(S)
            else
                out = Array{UInt16,1}(undef, length(S))
                # Copy first element
                out[1] = S[1]
                # Convert the rest using SIMD
                Arest = convert(UInt12Array{UInt16}, @view S[2:end])
                unpack_uint12_to_uint16(Arest.data, @view out[2:end])
                out
            end
        catch err
            @warn "Unable to convert using SIMD. Defaulting to slower element-wise conversion." err
            Array{UInt16,1}(S)
        end
    end
    Base.convert(::Type{Vector}, S::SubArray{UInt16, 1, UInt12Array{UInt16, B, 1}, <: Tuple{UnitRange}}) where {B <: SIMD.FastContiguousArray{UInt8,1}} =
        Base.convert(Vector{UInt16}, S)
    Base.convert(::Type{Array}, S::SubArray{UInt16, 1, UInt12Array{UInt16, B, 1}, <: Tuple{UnitRange}}) where {B <: SIMD.FastContiguousArray{UInt8,1}} =
        Base.convert(Vector{UInt16}, S)
end


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

mutable struct KeyVaultAndKeyReference <: SwaggerModel
    sourceVault::Any # spec type: Union{ Nothing, SourceVault } # spec name: sourceVault
    keyUrl::Any # spec type: Union{ Nothing, String } # spec name: keyUrl

    function KeyVaultAndKeyReference(;sourceVault=nothing, keyUrl=nothing)
        o = new()
        validate_property(KeyVaultAndKeyReference, Symbol("sourceVault"), sourceVault)
        setfield!(o, Symbol("sourceVault"), sourceVault)
        validate_property(KeyVaultAndKeyReference, Symbol("keyUrl"), keyUrl)
        setfield!(o, Symbol("keyUrl"), keyUrl)
        o
    end
end # type KeyVaultAndKeyReference

const _property_map_KeyVaultAndKeyReference = Dict{Symbol,Symbol}(Symbol("sourceVault")=>Symbol("sourceVault"), Symbol("keyUrl")=>Symbol("keyUrl"))
const _property_types_KeyVaultAndKeyReference = Dict{Symbol,String}(Symbol("sourceVault")=>"SourceVault", Symbol("keyUrl")=>"String")
Base.propertynames(::Type{ KeyVaultAndKeyReference }) = collect(keys(_property_map_KeyVaultAndKeyReference))
Swagger.property_type(::Type{ KeyVaultAndKeyReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_KeyVaultAndKeyReference[name]))}
Swagger.field_name(::Type{ KeyVaultAndKeyReference }, property_name::Symbol) =  _property_map_KeyVaultAndKeyReference[property_name]

function check_required(o::KeyVaultAndKeyReference)
    (getproperty(o, Symbol("sourceVault")) === nothing) && (return false)
    (getproperty(o, Symbol("keyUrl")) === nothing) && (return false)
    true
end

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


text:
#Stochastic collocation towards an exponential bspline and the normal density
using Roots
import AQFED.Math: normcdf, normpdf, norminv
import AQFED.Black: blackScholesFormula, blackScholesVega, impliedVolatility
using LeastSquaresOptim
using BSplines
using SpecialFunctions
struct ExpBSplineCollocation
    g::QuadraticPP #a bspline, can I express price directly in terms of bspline basis?
    forward::Number
end

function solvePositiveStrike(c::ExpBSplineCollocation, expStrike::Number)
    pp = c.g
    n = length(pp.x)
    strike = log(expStrike)
    if strike <= pp.a[1] # is zero
        leftSlope = pp.b[1]
        return (strike - pp.a[1]) / leftSlope + pp.x[1], 0
    elseif strike > pp.a[n]
        rightSlope = pp.b[n] #linear extrapolation a + b x-xn = strike
        return (strike - pp.a[n]) / rightSlope + pp.x[n], n
    end
    i = searchsortedfirst(pp.a, strike)  # x[i-1]<z<=x[i]
    if strike == pp.a[i]
        return pp.x[i], i
    end
    if i > 1
        i -= 1
    end
    x0 = pp.x[i]
    c = pp.c[i]
    b = pp.b[i]
    a = pp.a[i]
    cc = a + x0 * (-b + x0 * c) - strike
    bb = b - 2 * x0 * c
    aa = c
    allck = quadRootsReal(aa, bb, cc)
    for cki in allck
        if cki > pp.x[i] - sqrt(eps(strike)) && cki <= pp.x[i+1] + sqrt(eps(strike))
            return cki, i
        end
    end
    println(allck, " ", pp.x[i], " ", pp.x[i+1], " ", a, " ", b, " ", c, " strike ", strike)
    throw(DomainError(strike, "strike not found"))
end

#QuadRootsReal returns the real roots of a*x^2+ b*x +c =0
function quadRootsReal(a::T, b::T, c::T)::AbstractArray{T} where {T}
    d = b^2 - 4a * c
    if iszero(d)# single root
        return [-b / (2a)]
    elseif d > zero(d)# two real roots
        d = b < zero(d) ? sqrt(d) - b : -sqrt(d) - b
        return [d / (2a), (2c) / d]
    else # two complex roots
        return []
    end
end

function priceEuropean(
    c::ExpBSplineCollocation,
    isCall::Bool,
    strike::Number,
    forward::Number,
    discountDf::Number,
)::Number
    ck, ckIndex = solvePositiveStrike(c, strike)

    useForward = true
    if useForward || isCall
        valuef = hermiteIntegralBounded(c, ck, ckIndex)
        valuek = normcdf(-ck)
        callPrice = valuef - strike * valuek
        putPrice = -(forward - strike) + callPrice
        if isCall
            return callPrice * discountDf
        else
            return putPrice * discountDf
        end
    else #put
        integral = zero(ck)
        pp = c.g
        if ck < pp.x[1]
            integral += firstMomentExtrapolationBounded(p, 1, -300.0, ck)
        else
            integral += firstMomentExtrapolationBounded(p, 1, -300.0, pp.x[1])
        end
        for j = 1:ckIndex-1
            integral += firstMomentBounded(p, j, pp.x[j], pp.x[j+1])
        end
        if ckIndex >= n
            integral += firstMomentExtrapolationBounded(p, n, pp.x[n], ck)
        else
            integral += firstMomentBounded(p, j, pp.x[ckIndex], ck)
        end
        return strike * normcdf(ck) - integral
    end
end

function density(c::ExpBSplineCollocation, strike::Number)::Number
    ck, ckIndex = solvePositiveStrike(c, strike)
    dp = evaluateDerivative(c.g, ck)
    an = normpdf(ck) / (dp * strike)
    return an
end

function adjustForward(lsc::ExpBSplineCollocation)
    theoForward = hermiteIntegral(lsc)
    lnTheo = log(theoForward)
    if isnan(lnTheo) || isinf(lnTheo)
        println("inf forward ", theoForward, " ", lsc.g)
    else
        # println("for ",theoForward)
        lnf = log(lsc.forward)
        lsc.g.a .+= lnf - lnTheo
    end
end

function hermiteIntegral(p::ExpBSplineCollocation)::Number
    return hermiteIntegralBounded(p, -300.0, 1)
end

function firstMomentExtrapolationBounded(p::ExpBSplineCollocation, i::Int, x0t::Number, x1t::Number)::Number
    pp = p.g
    # i = length(pp.x) or i = 1
    x0 = pp.x[i]
    a0 = pp.a[i] + x0 * (-pp.b[i])
    a1 = pp.b[i]
    e = exp(a0 + a1^2 / 2)
    if isinf(e) || e == 0 || isnan(e)
        return zero(a0)
    end
    s = (x0 >= zero(x0)) ? one(x0) : -one(x0)
    ndiff = (normcdf(s * (-x0t + a1)) - normcdf(s * (-x1t + a1)))
    if ndiff == zero(ndiff) #necessary for ForwardDiff
        return ndiff
    else
        return s * e * ndiff
    end
end

function firstMomentBounded(p::ExpBSplineCollocation, i::Int, a::Number, b::Number)::Number
    pp = p.g
    x0 = pp.x[i]
    x1 = pp.x[i+1]
    x0t = max(a, x0)
    x1t = min(b, x1)
    c0 = pp.c[i]
    a0 = pp.a[i] + x0 * (-pp.b[i] + x0 * c0)
    a1 = pp.b[i] - 2 * x0 * c0
    a2 = c0
    onec = 1 - 2 * a2
    lne = a0 + a1^2 / (2 * onec)
    e = exp(lne)
    if iszero(e)
        return e
    elseif isinf(e) || isnan(e)
        return zero(onec)
    elseif abs(onec) < eps(one(onec))
        return e * (x1t - x0t) / sqrt(2 * pi())
    elseif onec > zero(onec)
        sqrtonec = sqrt(onec)
        s = (x0 >= zero(x0)) ? one(x0) : -one(x0)
        ndiff = normcdf(s * (-x0t * sqrtonec + a1 / sqrtonec)) - normcdf(s * (-x1t * sqrtonec + a1 / sqrtonec))
        if ndiff == zero(ndiff) #necessary for ForwardDiff
            return ndiff
        else
            value = s * e * (ndiff) / sqrtonec
            return value
        end
    else
        sqrtonec = sqrt(-onec)
        fip = x1t * sqrtonec + a1 / sqrtonec
        fi = x0t * sqrtonec + a1 / sqrtonec
        s = (fi >= zero(fi)) ? one(x0) : -one(x0)
        useDawson = false
        if !useDawson
            ndiff = (erfi(-s * fi / sqrt(2)) - erfi(-s * fip / sqrt(2)))
            if ndiff == zero(ndiff) #necessary for ForwardDiff
                return ndiff
            else
                value = s * e * ndiff / (2 * sqrtonec)

                if isinf(value) || isnan(value)
                    println(e, " ", fi, " ", fip, " ", sqrtonec)
                    #return zero(x0t)
                end
                return value
            end
        else
            ea = exp(-x0t^2 * onec / 2 + a1 * x0t + a0)
            if iszero(ea)
                return ea
            end
            ebma = exp(-(x1t^2 - x0t^2) * onec / 2 + a1 * (x1t - x0t))
             if isinf(ea) || isinf(ebma) || isnan(ea) || isnan(ebma)
                return zero(ea)
            end
            erfb = dawson(-s * fip / sqrt(2)) * 2 / sqrt(pi)
            erfa = dawson(-s * fi / sqrt(2)) * 2 / sqrt(pi)
            ndiff = iszero(ebma) ? erfa : (erfa - ebma * erfb)
            if ndiff == zero(ndiff) #necessary for ForwardDiff
                return ndiff
            end
            fdif = s * ea * ndiff
            value = fdif / (2 * sqrtonec)
            if isnan(value) || isinf(value)
                println(onec, " ", -x1t^2 * onec / 2 + a1 * x1t + a0, " ", ea, " ", ebma, " dawson ", erfa, " ", erfb)
            end
            return value
        end
    end
end

function hermiteIntegralBounded(p::ExpBSplineCollocation, ck::Number, ckIndex::Int)::Number
    pp = p.g
    n = length(pp.x)
    i = ckIndex
    if i > length(pp.x)
        i -= 1
    end
    integral = zero(ck)
    if ck < pp.x[1]
        integral += firstMomentExtrapolationBounded(p, 1, ck, pp.x[1])
        i = 0
    elseif ck < pp.x[n]
        #include logck to x[i] with i-1 coeffs
        integral += firstMomentBounded(p, i, ck, pp.x[i+1])
    end
    for j = i+1:n-1
        integral += firstMomentBounded(p, j, pp.x[j], pp.x[j+1])
    end
    integral += firstMomentExtrapolationBounded(p, n, max(pp.x[n], ck), 300.0)
    if isnan(integral) || isinf(integral)
        println(ck, " ", ckIndex, " inf integral ", p.g)
        throw(DomainError("infinite integral"))
    end
    return integral
end

function makeExpBSplineCollocationGuess(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    size = 0,
    minSlope = eps(Float64)
) where {T}
    b = fitExpBSplineBachelier(strikes, callPrices, weights, τ, forward, discountDf, size = size, minSlope = minSlope)
end

function makeExpBSplineCollocation(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    minSlope = 1e-4,
    penalty = 0.0,
    size = 0,
    rawFit = false,
)::Tuple{ExpBSplineCollocation,Number} where {T} #return collocation and error measure
    strikesf, pricesf, weightsf = filterConvexPrices(
        strikes,
        callPrices ./ discountDf,
        weights,
        forward,
        tol = minSlope + sqrt(eps(one(minSlope))),
    )
    isoc = makeExpBSplineCollocationGuess(strikesf, pricesf, weightsf, τ, forward, 1.0, size = size, minSlope = minSlope)
    isoc, m =
        rawFit ? fit(isoc, strikes, callPrices, weights, forward, discountDf, minSlope = minSlope, penalty = penalty) :
        fit(isoc, strikesf, pricesf, weightsf, forward, discountDf, minSlope = minSlope, penalty = penalty)
    return isoc, m
end

using FastGaussQuadrature
function fitExpBSplineBachelier(strikes, prices, weights, τ, forward, discountDf; size::Int = 0, minSlope::Float64 = eps(Float64))
    m = length(strikes)
    i = findfirst(x -> x > forward, strikes)
    if i == nothing
        i = m
    elseif i == 1
        i = 2
    end
    price = (prices[i] * (forward - strikes[i-1]) + prices[i-1] * (strikes[i] - forward)) / (strikes[i] - strikes[i-1])
    strike = forward
    #bvol = Bachelier.bachelierImpliedVolatility(price, true, strike, τ, forward, discountDf)
    # isoc = IsotonicCollocation(Polynomials.Polynomial([sqrt(bvol * sqrt(τ))]), Polynomials.Polynomial([0.0]), forward)
    vol = impliedVolatility(true, price, forward, strike, τ, discountDf)
    σ = vol * sqrt(τ)
    # need to use black because ys need to be > 0. xs < 0 ok.
    strikesf, pif, xf = makeXFromUndiscountedPrices(strikes, prices, slopeTolerance = minSlope)
    mindx = minimum(xf[2:end] - xf[1:end-1])
    if mindx < zero(xf[1])
        throw(DomainError(mindx, "dx negative, x is decreasing"))
    end

    local x::Vector{Float64}
    local n::Int
    if size == 0
        n = length(xf)
        x = copy(xf)
        #@. x[1:n-1] = xf[1:n-1]+xf[2:n])/2
        # x[1:n-1] = xf[1:n-1]
        x[n] = max(xf[end] * 1.2, 2.0)
        #x[1] = xf[2] - (xf[2] - xf[1]) / 2
    else
        n = max(size, 3) + 1
        # x = collect(range(min(xf[1] * 1.2, -3.0), stop = max(xf[end] * 1.2, 3.0), length = n))
        x = gausshermite(n)[1]
        x = x .* (max(-xf[1] * 1.2, xf[end] * 1.2, 1.5) / x[end])
    end
    a = @. (-σ^2 / 2 + log(forward) + σ * x)
    b = @. (σ * ones(Float64, n))
    c = zeros(Float64, n)
    pp = QuadraticPP(a, b, c, x)
    isoc = ExpBSplineCollocation(pp, forward)
    # adjustForward(isoc)
    return isoc
end

Base.length(p::ExpBSplineCollocation) = Base.length(p.g.x)
Base.broadcastable(p::ExpBSplineCollocation) = Ref(p)

abstract type Bijection end
struct IdentityTransformation{T} <: Bijection end
(f::IdentityTransformation{T})(x::T) where {T} = x
inv(f::IdentityTransformation{T}, x::T) where {T} = x

struct ClosedTransformation{T} <: Bijection
    minValue::T
    maxValue::T
end

function (f::ClosedTransformation{T})(x) where {T}
    return (f.minValue + f.maxValue) / 2 + (f.minValue - f.maxValue) / 2 * cos(x)
end

function inv(f::ClosedTransformation{T}, y) where {T}
    x = (y - f.minValue + (f.minValue - f.maxValue) / 2) / (f.minValue - f.maxValue) * 2
    return acos(x)
end

struct SquareMinTransformation{T} <: Bijection
    minValue::T
end

function (f::SquareMinTransformation{T})(x) where {T}
    return x^2 + f.minValue
end

function inv(f::SquareMinTransformation{T}, y) where {T}
    return sqrt(y - f.minValue)
end

struct ExpMinTransformation{T} <: Bijection
    minValue::T
end

function (f::ExpMinTransformation{T})(x) where {T}
    return exp(x) + f.minValue
end

function inv(f::ExpMinTransformation{T}, y) where {T}
    return log(y - f.minValue)
end

using ForwardDiff

function fit(isoc::ExpBSplineCollocation, strikes, prices, weights, forward, discountDf; minSlope = 1e-6, penalty = 0.0)
    iter = 0
    basis = BSplineBasis(3, isoc.g.x)
    t = BSplines.knots(basis)
    c = zeros(Float64, length(basis)) # = length x + 1
    spl = convert(BSplines.Spline, isoc.g)
    ct = zeros(Float64, length(spl.coeffs) - 1)
    minValue = max(1e-8 * (spl.coeffs[end] - spl.coeffs[1]), minSlope)
    maxValue = 4 * (spl.coeffs[end] - spl.coeffs[1])
    #transform = ExpMinTransformation(minValue)
    transform = ClosedTransformation(minValue, maxValue)
    for i = 1:length(ct)
        ct[i] = inv(transform, min(max(spl.coeffs[i+1] - spl.coeffs[i], minValue), maxValue))
        # println(i, " ", spl.coeffs[i+1] - spl.coeffs[i]," ct ",ct[i])
    end
    function obj(ct0::AbstractArray{T}) where {T}

        ct = @. transform(ct0)
        α = zeros(T, length(basis)) # = length x + 1
        α[1] = -sum(ct) #balance out theoretical forward such that it does not explode
        for i = 2:length(α)
            α[i] = ct[i-1] + α[i-1]
        end
        # println("iteration ",iter," ", ForwardDiff.value.(ct0))
        spl = BSplines.Spline(basis, α)
        pp = convert(QuadraticPP, spl)
        lsc = ExpBSplineCollocation(pp, forward)
        # ck = solvePositiveStrike(lsc, strikes[1])[1]
        # verr = firstMomentBounded(lsc, 1, pp.x[1], pp.x[2])
        # println(verr)
        adjustForward(lsc)
        iter += 1
        verr = @. weights * (priceEuropean(lsc, true, strikes, forward, discountDf) - prices)
        if penalty > 0
            vpen = @. ((1 / lsc.g.b[2:end] - 1 / lsc.g.b[1:end-1]) * penalty) #more appropriate if transform is unbounded > 0
            # vpen = @. (lsc.g.c[1:end] * penalty) #bounded transform?
            return vcat(verr, vpen)
        else
            return vcat(verr)
        end
    end
        function obj!(fvec, x)
        fvec[:] = obj(x)
        fvec
    end
#     cfg = ForwardDiff.JacobianConfig(obj, ct)
# function jac!(fvec, x)
#         fvec[:] = ForwardDiff.jacobian(obj, x, cfg)
#         #println("jac ",ForwardDiff.value.(fvec))
#         if isnan(ForwardDiff.value.(fvec[1]))
#             println("NaN in fvec for x ", x)
#         end
#         fvec
#     end
    # xerr = [0.39995341844904403, 0.5701199797586093, 0.5719276151194864, 0.5749687387719313, 0.5792870114145262, 0.5849459867261593, 0.5920318162549387, 0.6006571022625484, 0.6109662554578916, 0.6231428961116537, 0.637420114896599, 0.6540948511536501, 0.6735483718820394, 0.6962760710263824, 0.7229319988442942, 0.7543975847042282, 0.7918919014142847, 0.8373318965074066, 0.8941117058099329, 0.968128639742918, 1.0708664149833151, 1.2290358686740526, 1.528853236357424, 3.1467136183063698, 0.6495297531849857, 0.2032851901865146, 0.1649241311004164, 0.16488058052073906, 0.1945626540660544, 0.4238521533062266, 0.3181125312715509, 0.20542436305187858, 0.1820776017058426, 0.16554249062832802, 0.15359025795483686, 0.14744503512937399, 0.14343801501723907, 0.13812943948290135, 0.13187658168464514, 0.12548380500682077, 0.11896539357695363, 0.1132340096687979, 0.10887520969143545, 0.10616838922548384, 0.10466025306165908, 0.10356814434805951, 0.10255067548487518, 0.10182025744343283, 0.10145458270349705, 0.10093016161063079, 0.09983396460674955, 0.09814047452830171, 0.09601978776914528, 0.09376538827320721, 0.09175270079823829, 0.0903275648714271, 0.08961473802430946, 0.08937007223655123, 0.08926136226468968, 0.08902365495660793, 0.08855590378343575, 0.08796143487378127, 0.0874257083869379, 0.08698206005804007, 0.08658184258148705, 0.08616569498269007, 0.08573478298521706, 0.08531502458771066, 0.08496648156850824, 0.08477127394096373, 0.08483486149431202, 0.08511029161991776, 0.08549169345727906, 0.08587987096959282, 0.08622352366061321, 0.08651601223406952, 0.08677746695661273, 0.08701303156828352, 0.08722439396408083, 0.08741327232368615, 0.08758138688337708, 0.08773043542525312, 0.08786207262280761, 0.08797789320757941, 0.08807941878729295, 0.0881680880488771, 0.08824525001618294, 0.08831215999677243, 0.08836997783949553, 0.08841976812961458, 0.08846250196607737, 0.08849905999232086, 0.08853023638406257, 0.08855674353241154, 0.08857921719602688, 0.08859822193060779, 0.0886142566364226, 0.08862776009442899, 0.08863911638821935, 0.08864866013245495, 0.06267943175334381]
    # fvec = ForwardDiff.jacobian(obj, xerr, cfg)
    # if isnan(fvec[1])
    #     throw(DomainError(fvec[1]))
    # end
    # throw(DomainError(0))
    outlen = length(strikes)
    if penalty > 0
        outlen += length(isoc.g.b) - 1
    end
    #fit = optimize(obj, ct, LevenbergMarquardt(); autodiff = :forward, show_trace = false, iterations = 1024)
    fit = optimize!(
        LeastSquaresProblem(x = ct, f! = obj!, autodiff = :forward,
        # g! = jac!, #useful to debug issue with ForwardDiff NaNs
        output_length = outlen),
        LevenbergMarquardt();
        iterations = 1024,
    )
    # fit = fsolve(obj!, jac!, ct, outlen; show_trace=false, method=:lm, tol=1e-8) #fit.x
    println(iter, " fit ", fit, obj(fit.minimizer)) #obj(fit.x))  #fit.f
    ct0 = fit.minimizer

    ct = @. transform(ct0)
    c[1] = -sum(ct)
    for i = 2:length(c)
        c[i] = ct[i-1] + c[i-1]
    end
    spl = BSplines.Spline(basis, c)
    pp = convert(QuadraticPP, spl)
    measure = fit.ssr #sqrt(sum(x -> x^2, fit.f)/length(fit.f)) #fit.ssr
    lsc = ExpBSplineCollocation(pp, forward)
    adjustForward(lsc)
    return lsc, measure
end


text:
###
### LongSequence specific specializations of src/biosequence/indexing.jl
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# assumes `i` is positive and `bitsof(A)` is a power of 2

function Base.getindex(seq::LongSequence, part::UnitRange{<:Integer})
	@boundscheck checkbounds(seq, part)
	newseq = typeof(seq)(length(part))
	return copyto!(newseq, 1, seq, first(part), length(part))
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::AbstractVector{<:Integer}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    return unsafe_setindex!(seq, other, locs)
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::AbstractVector{Bool}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    return unsafe_setindex!(seq, other, locs)
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::UnitRange{<:Integer}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    return copyto!(seq, locs.start, other, 1, length(locs))
end

function Base.setindex!(seq::LongSequence{A},
			            other::LongSequence{A}, ::Colon) where {A}
    return setindex!(seq, other, 1:lastindex(seq))
end

function Base.setindex!(seq::LongSequence{A}, x, ::Colon) where {A}
    return setindex!(seq, x, 1:lastindex(seq))
end

@inline function encoded_setindex!(s::LongSequence, v::Unsigned, i::BitIndex)
    vi, off = i
    data = encoded_data(s)
    bits = @inbounds data[vi]
	v_ = v % encoded_data_eltype(s)
    @inbounds data[vi] = (v_ << off) | (bits & ~(bindata_mask(s) << off))
    return s
end

@inline function unsafe_setindex!(seq::LongSequence, x, locs::AbstractVector{<:Integer})
	bin = encode(Alphabet(seq), convert(eltype(seq), x))
    for i in locs
        encoded_setindex!(seq, bin, bitindex(seq, i))
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence, x, locs::AbstractVector{Bool})
    bin = encode(Alphabet(seq), convert(eltype(seq), x))
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence{A}, other::LongSequence{A},
	                      locs::AbstractVector{<:Integer}) where {A <: Alphabet}
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence{A},
                          other::LongSequence{A},
                          locs::AbstractVector{Bool}) where {A <: Alphabet}
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        unsafe_setindex!(seq, other[j += 1], i)
    end
    return seq
end

@inline function encoded_setindex!(seq::LongSequence{A},
				   bin::Unsigned, i::Integer) where {A <: Alphabet}
    j, r = bitindex(seq, i)
	bin_ = bin % encoded_data_eltype(seq)
    data = encoded_data(seq)
    @inbounds data[j] = (bin_ << r) | (data[j] & ~(bindata_mask(seq) << r))
    return seq
end


text:
## algebra
#include("fit_ss.jl")

function pow0(z::ComplexF64,k::Int=0)::ComplexF64
  if k==0
    return 1
  elseif k>0
    return 0
  else
    return complex(NaN)
  end
end

pow0(z::Symbol)=""

function pow1(z::ComplexF64,k::Int=0)::ComplexF64
  if k==0
    return z
  elseif k==1
    return 1
  elseif k>1
    return 0
  else
    return complex(NaN)#force NaN if a negative derivative is requested
  end
end

pow1(z::Symbol)="$z"

function pow2(z::ComplexF64,k::Int=0)::ComplexF64
  if k==0
    return z^2
  elseif k==1
    return 2*z
  elseif k==2
    return 2
  elseif k>2
    return 0
  else
    return complex(NaN) #force NaN if a negative derivative is requested
  end
end
pow2(z::Symbol)="$z^2"

#TODO turn this into a macro
function pow(z::ComplexF64,k::Int64,a::Int64)::ComplexF64
  if k>a>0
    f= 0
  elseif k>=0
    f=1
    i=a
    for j=0:k-1
      f*=i
      i-=1
    end
    f*=z^(a-k)
  else
    f=complex(NaN)
  end
  return f
end

function pow(z::ComplexF64,k::Int64,a::Number)::ComplexF64
  if k>=0
    f=1
    i=a
    for j=0:k-1
      f*=i
      i-=1
    end
    f*=z^(a-k)
  else
    f=complex(NaN)
  end
  return f
end

function pow_a(a::Number)

  function f(z::ComplexF64,k::Int64=0)::ComplexF64
    let a=a
      return pow(z,k,a)
    end#let
  end
  #This code fails incremental compilation.
  #TODO: check why
  # if a==0
  #   f(z::Symbol)=""
  # elseif a==1
  #   f(z::Symbol)="$z"
  # else
  #   f(z::Symbol)="$z^$a"
  # end
  function f(z::Symbol)
    let a=a
      if a==0
        return ""
      elseif a==1
        return "$z"
      else
        return
        "$z^$a"
      end
    end
  end
  return f
end

function generate_exp_az(a)
  function f(z::ComplexF64,k::Int)::ComplexF64
    let a=a
      if k>=0
        return a^k*exp(a*z)
      else
        return NaN+NaN*1im
      end
    end
  end

  function f(z::Symbol)::String
    let a=a
      return "exp(a$z)"
    end
  end
  return f
end

function exp_az(z::ComplexF64,a::ComplexF64,k::Int)::ComplexF64
  if k>=0
    return a^k*exp(a*z)
  else
    return nothing
  end
end

function exp_delay(ω::ComplexF64,τ::ComplexF64,m::Int,n::Int)::ComplexF64
  a=-1.0im
  u(z,l)=pow(z,l,m)
  f=0
  for i = 0:n
    f+=binomial(n,i)*u(τ,i)*(a*ω)^(n-i)
  end
  f*=a^m*exp(a*ω*τ)
  return f
end
function exp_delay(ω::Symbol,τ::Symbol)::String
  return "exp(-i$ω$τ)"
end

###

#TODO: improve closures by using let blocks or fastclosure package

function generate_stsp_z(A,B,C,D)
  function stsp_z(z,n)
    f = (-im)^n*factorial(n)*C*(im*z*LinearAlgebra.I-A)^(-n-1)*B
    if n==0
      f = f+D
    end
    return f[1]
  end
  return stsp_z
end

function generate_z_g_z(g)
  function z_g_z(z,n)
    if n == 0
      f = z*g(z,0)
    else
      f = z*g(z,n)+n*g(z,n-1)
    end
    return f
  end
  return z_g_z
end

tau_delay=exp_delay

#TODO create a macro for generating functions

# function powa(z::ComplexF64,k::Int,a::Int)
#   if
# end
function z_exp_iaz(z::ComplexF64, a::ComplexF64,m=0,n=0)::ComplexF64
  if m==n==0
    return z*exp(1im*a*z)
  elseif m==1
    (1im*a*z +1)*exp(1im*a*z)
  elseif n==1
    return 1im*z^2*exp(1im*a*z)
  else
    print("Argh!!!!!")
  end
end

function z_exp__iaz(z::ComplexF64, a::ComplexF64,m=0,n=0)::ComplexF64
  if m==n==0
    return z*exp(-1im*a*z)
  elseif m==1
    (-1im*a*z +1)*exp(-1im*a*z)
  elseif n==1
    return- 1im*z^2*exp(-1im*a*z)
  else
    print("Argh!!!!!")
  end
end

function exp_pm(s)
a=s*1.0im
  function exp_delay(ω,τ,m,n)
    u(z,l)=pow(z,l,m)
    f=0
    for i = 0:n
      f+=binomial(n,i)*u(τ,i)*(a*ω)^(n-i)
    end
    f*=a^m*exp(a*ω*τ)
    return f
  end
return exp_delay
end

function exp_ax2(z,a,n)
  if a==0.0+0im
    if n==0
      return 1.0+0.0im
    else
      return 0.0+0.0im
    end
  end

  f=0.0+0.0im
  A=a^n
  Z=z^n
  cnst=2^n*factorial(n)
  for k=0:n÷2
    coeff=0.0+0.0im
    coeff+=cnst*4.0^(-k)/factorial(k)/factorial(n-2*k)
    coeff*=A
    coeff*=Z
    f+=coeff
    A/=a
    Z/=z^2
  end
  f*=exp(a*z^2)
  return f
end

function exp_az2mzit(z,τ,a,m,n,k)
  f=pow_a(n+2*k)
  g(z,l)=exp_ax2(z,a,l)
  h(z,l)=exp_delay(z,τ,l,0)
  coeff=0.0+0.0im
  multinomialcoeff=factorial(m)
  for ii=0:m
    multi_ii=factorial(ii)
    coeff_ii=h(z,ii)
    for jj=0:m-ii
      multi_jj=factorial(jj)
      coeff_jj=g(z,jj)
      kk=m-jj-ii
      multi=multinomialcoeff/multi_ii/multi_jj/factorial(kk)
      coeff+=multi*f(z,kk)*coeff_jj*coeff_ii
    end
  end
    coeff*=(-1.0im)^n
    return coeff
end

function generate_Σy_exp_ikx(y)
  N=length(y)
  function Σy_exp_ikx(z::ComplexF64,n::Int)
    f=0.0+0.0im
    for (k,y) in enumerate(y)
      k-=1 #zero-based counting of wave-number
      f+=k^n*y*exp(2*pi*1.0im*k/N*z)
    end
    f*=(2*pi*1.0im/N)^n
    return f
  end
  return Σy_exp_ikx
end

function generate_gz_hz(g,h)
  function func(z::ComplexF64,k::Int)
    f=0.0+0.0im
    for i = 0:k
      f+=binomial(k,i)*h(z,k-i)*g(z,i)
    end
    return f
  end
  return func
end

function generate_1_gz(g)
  function func(z::ComplexF64,k::Int)
    if k==0
      return 1-g(z,k)
    else
      return -g(z,k)
    end
  end
  return func
end

function Σnexp_az2mzit(args...)
  J=(length(args)-2)÷6
  f=0.0+0.0im
  z=args[1] #argument is started by eigenfrequency...
  m=args[J*3+1+1]
  i=2:4
  for j=0:J-1
    nn,τ,a=args[i.+j*3] #... and continues with repeated sets of n, τ and a
    l,n,k=args[i.+3*j.+3*J.+1]
    f+=pow1(nn,l)*exp_az2mzit(z,τ,a,m,n,k)
  end
  return f
end


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

mutable struct NetworkConfigurationDiagnosticProfile <: SwaggerModel
    direction::Any # spec type: Union{ Nothing, Direction } # spec name: direction
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    source::Any # spec type: Union{ Nothing, String } # spec name: source
    destination::Any # spec type: Union{ Nothing, String } # spec name: destination
    destinationPort::Any # spec type: Union{ Nothing, String } # spec name: destinationPort

    function NetworkConfigurationDiagnosticProfile(;direction=nothing, protocol=nothing, source=nothing, destination=nothing, destinationPort=nothing)
        o = new()
        validate_property(NetworkConfigurationDiagnosticProfile, Symbol("direction"), direction)
        setfield!(o, Symbol("direction"), direction)
        validate_property(NetworkConfigurationDiagnosticProfile, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(NetworkConfigurationDiagnosticProfile, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        validate_property(NetworkConfigurationDiagnosticProfile, Symbol("destination"), destination)
        setfield!(o, Symbol("destination"), destination)
        validate_property(NetworkConfigurationDiagnosticProfile, Symbol("destinationPort"), destinationPort)
        setfield!(o, Symbol("destinationPort"), destinationPort)
        o
    end
end # type NetworkConfigurationDiagnosticProfile

const _property_map_NetworkConfigurationDiagnosticProfile = Dict{Symbol,Symbol}(Symbol("direction")=>Symbol("direction"), Symbol("protocol")=>Symbol("protocol"), Symbol("source")=>Symbol("source"), Symbol("destination")=>Symbol("destination"), Symbol("destinationPort")=>Symbol("destinationPort"))
const _property_types_NetworkConfigurationDiagnosticProfile = Dict{Symbol,String}(Symbol("direction")=>"Direction", Symbol("protocol")=>"String", Symbol("source")=>"String", Symbol("destination")=>"String", Symbol("destinationPort")=>"String")
Base.propertynames(::Type{ NetworkConfigurationDiagnosticProfile }) = collect(keys(_property_map_NetworkConfigurationDiagnosticProfile))
Swagger.property_type(::Type{ NetworkConfigurationDiagnosticProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkConfigurationDiagnosticProfile[name]))}
Swagger.field_name(::Type{ NetworkConfigurationDiagnosticProfile }, property_name::Symbol) =  _property_map_NetworkConfigurationDiagnosticProfile[property_name]

function check_required(o::NetworkConfigurationDiagnosticProfile)
    (getproperty(o, Symbol("direction")) === nothing) && (return false)
    (getproperty(o, Symbol("protocol")) === nothing) && (return false)
    (getproperty(o, Symbol("source")) === nothing) && (return false)
    (getproperty(o, Symbol("destination")) === nothing) && (return false)
    (getproperty(o, Symbol("destinationPort")) === nothing) && (return false)
    true
end

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


text:
"""
    NLNS(n_steps, x₀)

A nonlinear, nonstationary process.

## Equations
The systems is given by the following stochastic model:

```math
x(t+1) = x(t) + a(t) a(t-1)
```

with ``a(t) \textasciitilde N(0,1)`` being stationary. The nonlinearity lies in
the ``a(t) a(t-1)`` term.
This is, in essence, a nonlinear random walk. `x₀` sets the initial condition.

This example is from [Lucio et al. (2012)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202).

## References

Lucio et al., Phys. Rev. E *85*, 056202 (2012), after J. Timmer,
Phys. Rev. E *58*, 5153 (1998). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)

"""
function NLNS(n_steps, x₀)
    a = rand(Uniform(0, 1), n_steps)

    x = Vector{Float64}(n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = x[i-1] + a[i]*a[i-1]
    end

    x
end

"""
    AR1(n_steps, x₀, k)

Simple AR(1) model with no static transformation from [Lucio et al. (2012)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202).

## Equations
The system is given by the following map [1]:

```math
x(t+1) = k x(t) + a(t),
```

where ``a(t)`` is a draw from a normal distribution with zero mean and unit
variance. `x₀` sets the initial condition and `k` is the tunable parameter in
the map.

# References

Lucio et al., Phys. Rev. E *85*, 056202 (2012). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function AR1(n_steps, x₀, k)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = k*x[i-1] + a[i]
    end
    x
end

"""
    NSAR2(n_steps, x₀, x₁)

Cyclostationary AR(2) process from [Lucio et al. (2012)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202).

## References

Lucio et al., Phys. Rev. E *85*, 056202 (2012), after J. Timmer,
Phys. Rev. E *58*, 5153 (1998). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)

"""
function NSAR2(n_steps, x₀, x₁)
    T₀ = 50.0
    τ = 10.0
    M = 5.5
    Tmod = 250
    σ₀ = 1.0
    a_1_0 = 2*cos(2*π/T₀)*exp(-1/τ)
    a₂ = -exp(-2 / τ)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(n_steps)
    x[1], x[2] = x₀, x₁

    for i = 3:n_steps
        T_t = T₀ + M*sin(2*π*i / Tmod)
        a₁_t = 2*cos(2*π/T_t) * exp(-1 / τ)

        σ = σ₀^2/(1 - a_1_0^2 - a₂^2 - (2*(a_1_0^2) * a₂)/(1 - a₂)) *
            (1 - a₁_t - a₂^2 - (2*(a₁_t^2) * a₂)/(1 - a₂))
        x[i] = a₁_t * x[i-1] + a₂*x[i-2] + rand(Normal(0, sqrt(σ)))
    end
    x
end

"""
    randomwalk(n_steps, x₀)

Linear random walk (AR(1) process with a unit root) from [Lucio et al. (2012)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202).
This is an example of a nonstationary linear process.

# References

Lucio et al., Phys. Rev. E *85*, 056202 (2012), after J. Timmer,
Phys. Rev. E *58*, 5153 (1998). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function randomwalk(n_steps, x₀)
    a = rand(Normal(), n_steps)

    x = Vector{Float64}(n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = x[i-1] + a[i]
    end
    x
end

"""
    SNLST(n_steps, x₀, k)

Dynamically linear process transformed by a strongly nonlinear static
transformation (SNLST), from
[Lucio et al., 2012](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202).

## Equations
The system is by the following map:

```math
x(t) = k x(t-1) + a(t)
```

with the transformation ``s(t) = x(t)^3``.

# References

Lucio et al., Phys. Rev. E *85*, 056202 (2012), after J. Timmer,
Phys. Rev. E *58*, 5153 (1998).
[https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)

"""
function SNLST(n_steps, x₀, k)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = k*x[i-1] + a[i]
    end
    x.^3
end

# Keyword versions of the functions
AR1(;n_steps = 500, x₀ = rand(), k = rand())    = AR1(n_steps, x₀, k)
NLNS(;n_steps = 500, x₀ = rand())               = NLNS(n_steps, x₀)
SNLST(;n_steps = 500, x₀ = rand(), k = rand())  = SNLST(n_steps, x₀, k)
randomwalk(;n_steps = 500, x₀ = rand())         = randomwalk(n_steps, x₀)
NSAR2(;n_steps = 500, x₀ = rand(), x₁ = rand()) = NSAR2(n_steps, x₀, x₁)

export
    AR1,
    NLNS,
    SNLST,
    randomwalk,
    NSAR2


text:
#=На прямоугольном поле
произвольных размеров расставить
маркеры в виде "шахматных" клеток,
начиная с юго-западного угла поля,
когда каждая отдельная "шахматная" клетка имеет
размер n x n клеток поля (n - это параметр функции).
Начальное положение Робота - произвольное, конечное - совпадает с начальным.
Клетки на севере и востоке могут получаться "обрезанными" - зависит от соотношения
размеров поля и "шахматных" клеток.
(Подсказка: здесь могут быть полезными
две глобальных переменных, в которых будут
содержаться текущие декартовы координаты Робота
относительно начала координат в левом нижнем
углу поля, например)=#
X=0
Y=0
SIZE=0
function special_chesser(r::Robot, n::Int) #n-параметр функции
    global SIZE, Y, X
    SIZE=n
    vertical=moves!(r,Sud)
    horizontal=moves!(r,West)
    marker_now(r)
    go_back(r,Sud,Y-vertical)
    go_back(r,West,X-horizontal)
end

function to_point(r::Robot, side::HorizonSide) #Смещение по координатам
    global X, Y
    if side==Nord
        Y+=1
    elseif side==Sud
        Y-=1
    elseif side==Ost
        X+=1
    elseif side==West
        X-=1
    end
    move!(r,side)
end

function marker_now(r::Robot) #Маркировка
    mark_line(r,Ost)
    side=West
    while isborder(r,Nord)==false
        to_point(r,Nord)
        mark_line(r,side)
        side=inverse(side)
    end
end

function mark_line(r::Robot,side::HorizonSide) #Маркировка линии
    term_mark(r,X,Y,SIZE)
    while isborder(r,side)==false
        to_point(r,side)
        term_mark(r,X,Y,SIZE)
    end
end

function term_mark(r::Robot, x::Int, y::Int, size::Int) #Условие маркировки
    if ((mod(x,size+size) in 0:size-1)&&(mod(y,size+size) in 0:size-1))==true
        putmarker!(r)
    end
    if ((mod(x,size+size) in size:size+size-1)&&(mod(y,size+size) in size:size+size-1))==true
        putmarker!(r)
    end
end

function moves!(r::Robot,side::HorizonSide) #Описание в предыдущих задачах
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function go_back(r::Robot,side::HorizonSide,num_steps::Int) #Возвращение в исходную клетку
    if (num_steps<0)==true
        side=inverse(side)
        num_steps*=-1
    end
    for a in 1:num_steps
        move!(r,side)
    end
end

function inverse(side::HorizonSide) #Описание в предыдущих задачах и в библиотечном файле
    HorizonSide(mod(Int(side)+2,4))
end
#Данная задача основана на Задаче 7 - добавлены улучшения и новые функции

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.cluster;
#
# import typedb.client.api.TypeDBOptions;
# import typedb.client.api.TypeDBSession;
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.database.Database;
# import typedb.client.core.CoreClient;
# import typedb.client.core.CoreSession;
# import org.slf4j.Logger;
# import org.slf4j.LoggerFactory;
#
# public class ClusterSession implements TypeDBSession {
#     private static final Logger LOG = LoggerFactory.getLogger(TypeDBSession.class);
#     private final ClusterClient clusterClient;
#     private final TypeDBOptions.Cluster options;
#     private CoreClient coreClient;
#     private CoreSession coreSession;
#
#     public ClusterSession(ClusterClient clusterClient, String serverAddress, String database, Type type, TypeDBOptions.Cluster options) {
#         this.clusterClient = clusterClient;
#         this.coreClient = clusterClient.coreClient(serverAddress);
#         LOG.debug("Opening a session to '{}'", serverAddress);
#         this.coreSession = coreClient.session(database, type, options);
#         this.options = options;
#     }
#
#     @Override
#     public TypeDBTransaction transaction(TypeDBTransaction.Type type) {
#         return transaction(type, TypeDBOptions.cluster());
#     }
#
#     @Override
#     public TypeDBTransaction transaction(TypeDBTransaction.Type type, TypeDBOptions options) {
#         TypeDBOptions.Cluster clusterOpt = options.asCluster();
#         if (clusterOpt.readAnyReplica().isPresent() && clusterOpt.readAnyReplica().get()) {
#             return transactionAnyReplica(type, clusterOpt);
#         } else {
#             return transactionPrimaryReplica(type, options);
#         }
#     }
#
#     private TypeDBTransaction transactionPrimaryReplica(TypeDBTransaction.Type type, TypeDBOptions options) {
#         return transactionFailsafeTask(type, options).runPrimaryReplica();
#     }
#
#     private TypeDBTransaction transactionAnyReplica(TypeDBTransaction.Type type, TypeDBOptions.Cluster options) {
#         return transactionFailsafeTask(type, options).runAnyReplica();
#     }
#
#     private FailsafeTask<TypeDBTransaction> transactionFailsafeTask(TypeDBTransaction.Type type, TypeDBOptions options) {
#         return new FailsafeTask<TypeDBTransaction>(clusterClient, database().name()) {
#
#             @Override
#             TypeDBTransaction run(ClusterDatabase.Replica replica) {
#                 return coreSession.transaction(type, options);
#             }
#
#             @Override
#             TypeDBTransaction rerun(ClusterDatabase.Replica replica) {
#                 if (coreSession != null) coreSession.close();
#                 coreClient = clusterClient.coreClient(replica.address());
#                 coreSession = coreClient.session(database().name(), ClusterSession.this.type(), ClusterSession.this.options());
#                 return coreSession.transaction(type, options);
#             }
#         };
#     }
#
#     @Override
#     public TypeDBSession.Type type() {
#         return coreSession.type();
#     }
#
#     @Override
#     public TypeDBOptions.Cluster options() {
#         return options;
#     }
#
#     @Override
#     public boolean isOpen() {
#         return coreSession.isOpen();
#     }
#
#     @Override
#     public void close() {
#         coreSession.close();
#     }
#
#     @Override
#     public Database database() {
#         return coreSession.database();
#     }
# }


text:
#=
interp1.jl
emulate matlab's interp1()
=#

export interp1

using Interpolations

"""
    yi = interp1(x, y, xi ; how=Gridded(Linear()), extrap=0)

1D interpolation of `y = f(x)` at points `xi`

In:
- `x::AbstractVector{<:Real}`
- `y::AbstractVector{<:Number}`

Option:
- `how::Interpolations.InterpolationType` default `Gridded(Linear())`
- `extrap::Any` how to extrapolate, e.g., `Flat()`; default `0`
other options from Interpolations.jl are `Line()` `Periodic()` `Reflect()` `Throw()`

Output is same size as input `xi`
"""
function interp1(x::AbstractVector{<:Real}, y::AbstractVector{<:Number}, xi ;
		how::Interpolations.InterpolationType = Gridded(Linear()),
		extrap::Any = 0)
#		extrap::Union{<:Interpolations.BoundaryCondition, <:Number} = 0 # fails?

	fun = interpolate((x,), y, how)
	if !isnothing(extrap)
		fun = extrapolate(fun, extrap)
	end
	fun.(xi)
end


text:
macro unless(condition, body)
  :(if !$condition
    $(esc(body))
  end)
end

x, y = 0, 1
@unless (x == 5) begin
  y = 3                  # should execute
end
@unless (x == 0) begin
  y = 5                  # should not execute
end

@assert y == 3

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using SEML, Test

@time begin include("test1.jl") end

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
windoptions() = Dict(
    :gisregion => "Europe8",            # "Europe8", "Eurasia38", "Scand3"
    :filenamesuffix => "",              # e.g. "_landx2" to save high land availability data as "GISdata_solar2018_Europe8_landx2.mat"

    :onshore_density => 5,              # about 30% of existing farms have at least 5 W/m2, will become more common
    :offshore_density => 8,             # varies a lot in existing parks (4-18 W/m2)
                                        # For reference: 10D x 5D spacing of 3 MW turbines (with 1D = 100m) is approximately 6 MW/km2 = 6 W/m2
    :area_onshore => .08,               # area available for onshore wind power after the masks have been applied
    :area_offshore => .33,              # area available for offshore wind power after the masks have been applied

    :distance_elec_access => 150,       # max distance to grid [km] (for wind classes of category B and offshore)
    :persons_per_km2 => 150,            # not too crowded, max X persons/km2
                                        # US census bureau requires 1000 ppl/mile^2 = 386 ppl/km2 for "urban" (half in Australia)
                                        # roughly half the people of the world live at density > 300 ppl/km2
    :max_depth => 40,                   # max depth for offshore wind [m]
    :min_shore_distance => 5,           # minimum distance to shore for offshore wind [km]
    :exclude_landtypes => [0,11,13],    # exclude water, wetlands and urban areas. See codes in table below.
    :protected_codes => [1,2,3,4,5,8],  # IUCN codes to be excluded as protected areas. See codes in table below.

    :scenarioyear => "ssp2_2050",       # default scenario and year for population and grid access datasets
    :era_year => 2018,                  # which year of the ERA5 time series to use
    :rescale_to_wind_atlas => true,     # rescale the ERA5 time series to fit annual wind speed averages from the Global Wind Atlas

    :res => 0.01,                       # resolution of auxiliary datasets [degrees per pixel]
    :erares => 0.28125,                 # resolution of ERA5 datasets [degrees per pixel]

    :onshoreclasses_min => [2,5,6,7,8],     # lower bound on annual onshore wind speeds for class X    [0:0.25:12.25;]
    :onshoreclasses_max => [5,6,7,8,99],    # upper bound on annual onshore wind speeds for class X    [0.25:0.25:12.5;]
    :offshoreclasses_min => [3,6,7,8,9],    # lower bound on annual offshore wind speeds for class X
    :offshoreclasses_max => [6,7,8,9,99],   # upper bound on annual offshore wind speeds for class X

    :downsample_masks => 1,     # set to 2 or higher to scale down mask sizes to avoid GPU errors in Makie plots for large regions
    :classB_threshold => 0.001  # minimum share of pixels within distance_elec_access km that must have grid access
                                # for a pixel to be considered for wind class B.
)
    # Land types
    #     0      'Water'
    #     1      'Evergreen Needleleaf Forests'
    #     2      'Evergreen Broadleaf Forests'
    #     3      'Deciduous Needleleaf Forests'
    #     4      'Deciduous Broadleaf Forests'
    #     5      'Mixed Forests'
    #     6      'Closed Shrublands'
    #     7      'Open Shrublands'
    #     8      'Woody Savannas'
    #     9      'Savannas'
    #    10      'Grasslands'
    #    11      'Permanent Wetlands'
    #    12      'Croplands'
    #    13      'Urban'
    #    14      'Cropland/Natural'
    #    15      'Snow/Ice'
    #    16      'Barren'

    # Protected areas (IUCN codes from the WDPA)
    #    1      'Ia'                'Strict Nature Reserve'
    #    2      'Ib'                'Wilderness Area'
    #    3      'II'                'National Park'
    #    4      'III'               'Natural Monument'
    #    5      'IV'                'Habitat/Species Management'
    #    6      'V'                 'Protected Landscape/Seascape'
    #    7      'VI'                'Managed Resource Protected Area'
    #    8      'Not Reported'      'Not Reported'
    #    9      'Not Applicable'    'Not Applicable'
    #    10     'Not Assigned'      'Not Assigned'

mutable struct WindOptions
    gisregion               ::String
    filenamesuffix          ::String
    onshore_density         ::Float64           # W/m2
    offshore_density        ::Float64           # W/m2
    area_onshore            ::Float64           # share [0-1]
    area_offshore           ::Float64           # share [0-1]
    distance_elec_access    ::Float64           # km
    persons_per_km2         ::Float64           # persons/km2
    max_depth               ::Float64           # m
    min_shore_distance      ::Float64           # km
    exclude_landtypes       ::Vector{Int}
    protected_codes         ::Vector{Int}
    scenarioyear            ::String
    era_year                ::Int
    rescale_to_wind_atlas   ::Bool
    res                     ::Float64           # degrees/pixel
    erares                  ::Float64           # degrees/pixel
    onshoreclasses_min      ::Vector{Float64}
    onshoreclasses_max      ::Vector{Float64}
    offshoreclasses_min     ::Vector{Float64}
    offshoreclasses_max     ::Vector{Float64}
    downsample_masks        ::Int
    classB_threshold        ::Float64
end

WindOptions() = WindOptions("","",0,0,0,0,0,0,0,0,[],[],"",0,false,0,0,[],[],[],[],0,0.0)

function WindOptions(d::Dict{Symbol,Any})
    options = WindOptions()
    for (key,val) in d
        setproperty!(options, key, val)
    end
    return options
end

function GISwind(; savetodisk=true, plotmasks=false, optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, era_year, filenamesuffix, downsample_masks = options

    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land, protected, lonrange, latrange =
                read_datasets(options)

    mask_onshoreA, mask_onshoreB, mask_offshore =
        create_wind_masks(options, regions, offshoreregions, gridaccess, popdens, topo, land, protected, lonrange, latrange,
                            plotmasks=plotmasks, downsample=downsample_masks)

    plotmasks == :onlymasks && return nothing

    windatlas, meanwind, windspeed = read_wind_datasets(options, lonrange, latrange)

    windCF_onshoreA, windCF_onshoreB, windCF_offshore, capacity_onshoreA, capacity_onshoreB, capacity_offshore =
        calc_wind_vars(options, windatlas, meanwind, windspeed, regions, offshoreregions, regionlist,
                mask_onshoreA, mask_onshoreB, mask_offshore, lonrange, latrange)

    if savetodisk
        mkpath(in_datafolder("output"))
        matopen(in_datafolder("output", "GISdata_wind$(era_year)_$gisregion$filenamesuffix.mat"), "w", compress=true) do file
            write(file, "CFtime_windonshoreA", windCF_onshoreA)
            write(file, "CFtime_windonshoreB", windCF_onshoreB)
            write(file, "CFtime_windoffshore", windCF_offshore)
            write(file, "capacity_onshoreA", capacity_onshoreA)
            write(file, "capacity_onshoreB", capacity_onshoreB)
            write(file, "capacity_offshore", capacity_offshore)
        end
    end

    nothing
    # return windCF_onshoreA, windCF_onshoreB, windCF_offshore, capacity_onshoreA, capacity_onshoreB, capacity_offshore
end

function read_datasets(options)
    @unpack res, gisregion, scenarioyear = options

    println("\nReading auxiliary datasets...")
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    gridaccess = JLD.load(in_datafolder("gridaccess_$scenarioyear.jld"), "gridaccess")[lonrange,latrange]
    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")[lonrange,latrange]
    topo = JLD.load(in_datafolder("topography.jld"), "topography")[lonrange,latrange]
    land = JLD.load(in_datafolder("landcover.jld"), "landcover")[lonrange,latrange]
    protected = JLD.load(in_datafolder("protected.jld"), "protected")[lonrange,latrange]

    popdens = pop ./ cellarea'

    # drawmap(log.(pop))

    return regions, offshoreregions, regionlist, gridaccess, popdens, topo, land, protected, lonrange, latrange
end

function read_wind_datasets(options, lonrange, latrange)
    @unpack res, erares, era_year = options
    windatlas = getwindatlas()[lonrange,latrange]

    println("Reading ERA5 wind speeds and calculating capacity factors...")
    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)

    @time meanwind, windspeed = h5open(in_datafolder("era5wind$era_year.h5"), "r") do file
        if length(eralonranges) == 1
            file["meanwind"][eralonranges[1], eralatrange],
                file["wind"][:, eralonranges[1], eralatrange]
        else
            [file["meanwind"][eralonranges[1], eralatrange]; file["meanwind"][eralonranges[2], eralatrange]],
                [file["wind"][:, eralonranges[1], eralatrange] file["wind"][:, eralonranges[2], eralatrange]]
        end
    end
    return windatlas, meanwind, windspeed
end

function create_wind_masks(options, regions, offshoreregions, gridaccess, popdens, topo, land, protected, lonrange, latrange; plotmasks=false, downsample=1)
    @unpack res, gisregion, exclude_landtypes, protected_codes, distance_elec_access, persons_per_km2,
                min_shore_distance, max_depth, classB_threshold, filenamesuffix = options

    println("Creating masks...")

    goodland = (regions .> 0)
    for i in exclude_landtypes
        goodland[land .== i] .= false
    end
    protected_area = zeros(Bool, size(protected))
    for i in protected_codes
        protected_area[protected .== i] .= true
    end

    # Pixels with electricity access for onshore wind A
    gridA = (gridaccess .> 0)

    # Pixels with electricity access for onshore wind B and offshore wind
    km_per_degree = π*2*6371/360
    disk = diskfilterkernel(distance_elec_access/km_per_degree/res)
    gridB = (imfilter(gridaccess, disk) .> max(1e-9, classB_threshold)) # avoid artifacts if classB_threshold == 0

    # println("MAKE SURE MASKS DON'T OVERLAP! (regions & offshoreregions, mask_*)")

    # all mask conditions
    mask_onshoreA = gridA .& (popdens .< persons_per_km2) .& goodland .& .!protected_area
    mask_onshoreB = (gridB .& .!gridA) .& (popdens .< persons_per_km2) .& goodland .& .!protected_area

    # shoreline mask for offshore wind
    disk = diskfilterkernel(min_shore_distance/km_per_degree/res)
    shore = (imfilter(regions .> 0, disk) .> 1e-6)

    # all mask conditions
    mask_offshore = gridB .& .!shore .& (topo .> -max_depth) .& (offshoreregions .> 0) .& .!protected_area

    if plotmasks != false   # can == :onlymasks as well
        # drawmap(land)
        isregion = (regions .> 0) .& (regions .!= NOREGION)

        # mask values refer to colors in ColorBrewer Set2_7:
        # https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#colorbrewer-1
        masks = zeros(Int16, size(regions))
        masks[(masks .== 0) .& (popdens .> persons_per_km2)] .= 2
        masks[(masks .== 0) .& protected_area] .= 3
        masks[(masks .== 0) .& .!gridA .& .!gridB] .= 4
        masks[(masks .== 0) .& .!goodland] .= 1
        masks[(masks .== 0) .& .!gridA .& gridB] .= 8
        masks[(masks .== 0) .& isregion] .= 7
        masks[regions .== 0] .= 0
        masks[regions .== NOREGION] .= NOREGION
        legendtext = ["bad land type", "high population", "protected area", "no grid", "", "", "wind plant A", "wind plant B"]
        maskmap("$(gisregion)_masks_wind$filenamesuffix", masks, legendtext, lonrange, latrange; legend=true, downsample=downsample)

        # drawmap(.!shore .& (offshoreregions .> 0))
        # drawmap((topo .> -max_depth) .& (offshoreregions .> 0))
        # drawmap(.!protected_area .& (offshoreregions .> 0))
        # drawmap(mask_offshore)
    end

    return mask_onshoreA, mask_onshoreB, mask_offshore
end

# 0 - 29 m/s
const windparkcurve = [
    0.0, 0.0014, 0.0071, 0.0229, 0.0545, 0.1067, 0.1831, 0.2850, 0.4085, 0.5434,
    0.6744, 0.7847, 0.8614, 0.9048, 0.9266, 0.9353, 0.9373, 0.9375, 0.9375, 0.9375,
    0.9375, 0.9375, 0.9375, 0.9311, 0.8683, 0.6416, 0.2948, 0.0688, 0.0063, 0.0
]

function speed2capacityfactor(windspeed)
    if windspeed >= 29 || windspeed < 0
        return 0.0
    end
    fw = floor(Int, windspeed)
    frac = windspeed - fw
    return (1-frac).*windparkcurve[fw+1] + frac.*windparkcurve[ceil(Int, windspeed)+1]
end

function increment_windCF!(cf::AbstractVector{<:AbstractFloat}, speed_or_cf::AbstractVector{<:AbstractFloat}, factor, rescale::Bool)
    if rescale
        @inbounds for i = 1:length(cf)
            cf[i] += speed2capacityfactor(speed_or_cf[i]*factor)
        end
    else
        @inbounds for i = 1:length(cf)
            cf[i] += speed_or_cf[i]
        end
    end
end

function getclasses0(windatlas, classes_min, classes_max)
    class = zeros(Int16, size(windatlas))
    for c = 1:length(classes_min)
        class[(windatlas .>= classes_min[c]) .& (windatlas .< classes_max[c])] .= c
    end
    return class
end

function getclasses(windatlas, classes_min, classes_max)
    class = similar(windatlas, Int16)
    nclasses = length(classes_min)
    @inbounds for i = 1:length(windatlas)
        wa = windatlas[i]
        if wa < classes_min[1] || wa > classes_max[end]
            class[i] = 0
            continue
        end
        for c = 1:nclasses
            if wa < classes_max[c]
                class[i] = c
                break
            end
        end
    end
    return class
end

function makewindclasses(options, windatlas)
    println("Allocating pixels to classes using the Global Wind Atlas...")
    # println("CHANGE TO CAPACITY FACTOR LATER!")

    @unpack onshoreclasses_min, onshoreclasses_max, offshoreclasses_min, offshoreclasses_max = options

    onshoreclass = getclasses(windatlas, onshoreclasses_min, onshoreclasses_max)
    offshoreclass = getclasses(windatlas, offshoreclasses_min, offshoreclasses_max)

    return onshoreclass, offshoreclass
end

function calc_wind_vars(options, windatlas, meanwind, windspeed, regions, offshoreregions, regionlist,
                mask_onshoreA, mask_onshoreB, mask_offshore, lonrange, latrange)

    println("Calculating GW potential and hourly capacity factors for each region and wind class...")
    println("Interpolate ERA5 wind speeds later (maybe 4x runtime).")

    onshoreclass, offshoreclass = makewindclasses(options, windatlas)
    eralons, eralats, lonmap, latmap, cellarea = eralonlat(options, lonrange, latrange)

    _, _, meanwind_allyears, _ = annualwindindex(options)
    @assert size(meanwind_allyears) == size(meanwind)

    @unpack onshoreclasses_min, offshoreclasses_min, rescale_to_wind_atlas, res, erares,
                onshore_density, area_onshore, offshore_density, area_offshore = options

    numreg = length(regionlist)
    nonshoreclasses, noffshoreclasses = length(onshoreclasses_min), length(offshoreclasses_min)
    yearlength, nlons, nlats = size(windspeed)

    capacity_onshoreA = zeros(numreg,nonshoreclasses)
    capacity_onshoreB = zeros(numreg,nonshoreclasses)
    capacity_offshore = zeros(numreg,noffshoreclasses)
    windCF_onshoreA = zeros(yearlength,numreg,nonshoreclasses)
    windCF_onshoreB = zeros(yearlength,numreg,nonshoreclasses)
    windCF_offshore = zeros(yearlength,numreg,noffshoreclasses)
    count_onshoreA = zeros(Int,numreg,nonshoreclasses)
    count_onshoreB = zeros(Int,numreg,nonshoreclasses)
    count_offshore = zeros(Int,numreg,noffshoreclasses)

    if rescale_to_wind_atlas
        println("\nRescaling ERA5 wind speeds to match annual wind speeds from the Global Wind Atlas.")
        # println("This will increase run times by an order of magnitude (since GWA has very high spatial resolution).")
    end

    # Run times vary wildly depending on geographical area (because of far offshore regions with mostly zero wind speeds).
    # To improve the estimated time of completing the progress bar, iterate over latitudes in random order.
    Random.seed!(1)
    updateprogress = Progress(nlats, 1)
    for j in randperm(nlats)
        eralat = eralats[j]
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for i = 1:nlons
            meanwind[i,j] == 0 && continue
            wind = rescale_to_wind_atlas ? windspeed[:, i, j] : speed2capacityfactor.(windspeed[:, i, j])
            eralon = eralons[i]
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]

            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                offreg = offshoreregions[r,c]
                area = cellarea[c]

                # @views is needed to make sure increment_windCF!() works with matrix
                # slice. It's also faster since it avoids making copies.
                # We divide by meanwind_allyears instead of meanwind in order to
                # account for more and less windy years, while still rescaling
                # hourly ERA5 wind speeds to match annual averages from the Global
                # Wind Atlas.
                if reg > 0 && reg != NOREGION
                    class = onshoreclass[r,c]
                    @views if reg > 0 && class > 0 && mask_onshoreA[r,c] > 0
                        capacity_onshoreA[reg,class] += 1/1000 * onshore_density * area_onshore * area
                        increment_windCF!(windCF_onshoreA[:,reg,class], wind, windatlas[r,c] / meanwind_allyears[i,j], rescale_to_wind_atlas)
                        count_onshoreA[reg,class] += 1
                    elseif reg > 0 && class > 0 && mask_onshoreB[r,c] > 0
                        capacity_onshoreB[reg,class] += 1/1000 * onshore_density * 2 * area_onshore * area
                        increment_windCF!(windCF_onshoreB[:,reg,class], wind, windatlas[r,c] / meanwind_allyears[i,j], rescale_to_wind_atlas)
                        count_onshoreB[reg,class] += 1
                    end
                end
                if offreg > 0 && offreg != NOREGION
                    offclass = offshoreclass[r,c]
                    @views if offreg > 0 && offclass > 0 && mask_offshore[r,c] > 0
                        capacity_offshore[offreg,offclass] += 1/1000 * offshore_density * area_offshore * area
                        increment_windCF!(windCF_offshore[:,offreg,offclass], wind, windatlas[r,c] / meanwind_allyears[i,j], rescale_to_wind_atlas)
                        count_offshore[offreg,offclass] += 1
                    end
                end
            end
        end
        next!(updateprogress)
    end

    for y = 1:yearlength
        windCF_onshoreA[y,:,:] ./= count_onshoreA
        windCF_onshoreB[y,:,:] ./= count_onshoreB
        windCF_offshore[y,:,:] ./= count_offshore
    end

    return windCF_onshoreA, windCF_onshoreB, windCF_offshore, capacity_onshoreA, capacity_onshoreB, capacity_offshore
end

# returns a wind speed index for 1979-2019, i.e. normalized average wind speed over
# each pixel in each region for each year
# also return a matrix of annual wind speeds in each ERA5 pixel over the region
function annualwindindex(options; resource=:wind, sites_quantile=1.0, aggregateregions=[])
    @unpack gisregion, era_year, res, erares, filenamesuffix = options
    regions, _, regionlist, lonrange, latrange = loadregions(gisregion)
    smallregions = resize_categorical(regions, regionlist, lonrange, latrange;
                skipNOREGION=true)
    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)
    varname = (resource == :wind) ? "annualwind" : "annualssrd"
    annualwind = getmonthlywind(:annual, resource, eralonranges, eralatrange, filenamesuffix)
    meanwind = meandrop(annualwind, dims=1)
    nyears, nreg = size(annualwind, 1), length(regionlist)
    aggregateregions = isempty(aggregateregions) ? [[r] for r = 1:nreg] : aggregateregions
    regionmasks = [in.(smallregions, Ref(aggreg)) for aggreg in aggregateregions]
    cutoffs = [quantile(vec(meanwind[(meanwind.>0) .& regmask]), 1-sites_quantile)
            for regmask in regionmasks]
    masks = [(meanwind .> cutoffs[i]) .& regionmasks[i] for i = 1:length(regionmasks)]
    index = [sum(annualwind[y,:,:] .* mask) ./ sum(mask) for y = 1:nyears, mask in masks]
    return index, index ./ mean(index, dims=1), meanwind, annualwind[era_year-1979+1,:,:]
end

function annualwindindex(; resource=:wind, sites_quantile=1.0, aggregateregions=[], optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    return annualwindindex(options, resource=resource, sites_quantile=sites_quantile,
                    aggregateregions=aggregateregions)
end

function seasonalwindprofile(; resource=:wind, sites_quantile=1.0, aggregateregions=[], optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, era_year, res, erares, filenamesuffix = options
    regions, _, regionlist, lonrange, latrange = loadregions(gisregion)
    smallregions = resize_categorical(regions, regionlist, lonrange, latrange;
                skipNOREGION=true)
    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)
    annualwind = getmonthlywind(:annual, resource, eralonranges, eralatrange, filenamesuffix)
    monthlywind = getmonthlywind(:monthly, resource, eralonranges, eralatrange, filenamesuffix)
    meanwind = meandrop(annualwind, dims=1)
    cutoff = quantile(vec(meanwind[meanwind.>0]), 1-sites_quantile)
    nyears, nreg = size(monthlywind, 1), length(regionlist)
    aggregateregions = isempty(aggregateregions) ? [[r] for r = 1:nreg] : aggregateregions
    profile = [sum(meandrop(monthlywind[m:12:end,:,:], dims=1) .*
        (mask = (meanwind .> cutoff) .& in.(smallregions, Ref(aggreg)))) ./ sum(mask)
                    for m = 1:12, aggreg in aggregateregions]
end

function getmonthlywind(time, resource, eralonranges, eralatrange, filenamesuffix)
    h5open(in_datafolder("era5monthly$resource$filenamesuffix.h5"), "r") do file
        varname = (resource == :wind) ? "$(time)wind" : "$(time)ssrd"
        if length(eralonranges) == 1
            file[varname][:, eralonranges[1], eralatrange]
        else
            [file[varname][:, eralonranges[1], eralatrange] file[varname][:, eralonranges[2], eralatrange]]
        end
    end
end

# Quick and ugly copy/paste hack to create resource maps for wind classes combined with masks.
function GISwindmap(; optionlist...)
    options = WindOptions(merge(windoptions(), optionlist))
    @unpack gisregion, era_year, filenamesuffix = options

    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land, protected, lonrange, latrange =
                read_datasets(options)
    windatlas, meanwind, windspeed = read_wind_datasets(options, lonrange, latrange)

    mask_onshoreA, mask_onshoreB, mask_offshore =
        create_wind_masks(options, regions, offshoreregions, gridaccess, popdens, topo, land, protected, lonrange, latrange, plotmasks=true)

    onshoremap, offshoremap = calc_wind_map(options, windatlas, meanwind, windspeed, regions, offshoreregions, regionlist,
                mask_onshoreA, mask_onshoreB, mask_offshore, lonrange, latrange)

    return onshoremap, offshoremap
end

function calc_wind_map(options, windatlas, meanwind, windspeed, regions, offshoreregions, regionlist,
                mask_onshoreA, mask_onshoreB, mask_offshore, lonrange, latrange)

    println("Calculating GW potential and hourly capacity factors for each region and wind class...")
    # println("Interpolate ERA5 wind speeds later (maybe 4x runtime).")

    onshoreclass, offshoreclass = makewindclasses(options, windatlas)
    eralons, eralats, lonmap, latmap, cellarea = eralonlat(options, lonrange, latrange)

    @unpack onshoreclasses_min, offshoreclasses_min, rescale_to_wind_atlas, res, erares,
                onshore_density, area_onshore, offshore_density, area_offshore = options

    numreg = length(regionlist)
    nonshoreclasses, noffshoreclasses = length(onshoreclasses_min), length(offshoreclasses_min)
    yearlength, nlons, nlats = size(windspeed)

    onshoremap = zeros(Int16, size(regions))
    offshoremap = zeros(Int16, size(regions))

    if rescale_to_wind_atlas
        println("\nRescaling ERA5 wind speeds to match annual wind speeds from the Global Wind Atlas.")
        println("This will increase run times by an order of magnitude (since GWA has very high spatial resolution).")
    end

    # Run times vary wildly depending on geographical area (because of far offshore regions with mostly zero wind speeds).
    # To improve the estimated time of completing the progress bar, iterate over latitudes in random order.
    Random.seed!(1)
    updateprogress = Progress(nlats, 1)
    @inbounds for j in randperm(nlats)
        eralat = eralats[j]
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for i = 1:nlons
            meanwind[i,j] == 0 && continue
            wind = rescale_to_wind_atlas ? windspeed[:, i, j] : speed2capacityfactor.(windspeed[:, i, j])
            eralon = eralons[i]
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]

            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                area = cellarea[c]
                class = onshoreclass[r,c]
                offreg = offshoreregions[r,c]
                offclass = offshoreclass[r,c]

                # can't use elseif here, probably some overlap in the masks
                # @views is needed to make sure increment_windCF!() works with matrix slices
                # also faster since it avoids making copies
                @views if reg > 0 && class > 0 && mask_onshoreA[r,c] > 0
                    onshoremap[r,c] = class
                elseif reg > 0 && class > 0 && mask_onshoreB[r,c] > 0
                    # onshoremap[r,c] = class
                elseif offreg > 0 && offclass > 0 && mask_offshore[r,c] > 0
                    offshoremap[r,c] = class
                end
            end
        end
        next!(updateprogress)
    end

    return onshoremap, offshoremap
end

max_stars_repo_name: niklasjakobsson/GlobalEnergyGIS

max_issues_repo_name: niklasjakobsson/GlobalEnergyGIS

max_forks_repo_name: niklasjakobsson/GlobalEnergyGIS


text:
"""
`SolutionPODE`: Solution of a partitioned ordinary differential equation

Contains all fields necessary to store the solution of an PODE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nt`: number of time steps to store
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ni]` with `q[:,0,:]` the initial conditions
* `p`:  solution `p[nd, nt+1, ni]` with `p[:,0,:]` the initial conditions
* `ntime`: number of time steps to compute
* `nsave`: store every nsave'th time step (default: 1)
* `nwrite`: save data to disk after every nwrite'th time step (default: ntime)
* `offset`: counter for file offset
* `counter`: counter for copied solution entries

### Constructors

```julia
SolutionPODE(equation, Δt, ntimesteps; nsave=DEFAULT_NSAVE, nwrite=DEFAULT_NWRITE, filename=nothing)
SolutionPODE(t::TimeSeries, q::DataSeries, p::DataSeries, ntimesteps)
SolutionPODE(h5io::SolutionHDF5)
SolutionPODE(file::String)
```

The constructors `SSolutionPODE` create a `SolutionPODE` with internal data structures
for serial simulations (i.e., standard arrays), while the constructors `PSolutionPODE`
create a `SolutionPODE` with internal data structures for parallel simulations (i.e.,
shared arrays).

The usual way to initialise a `Solution` is by passing an equation, which for
`SolutionPODE` has to be an [`PODE`](@ref), [`HODE`](@ref), [`IODE`](@ref) or
[`LODE`](@ref), a time step `Δt` and the number of time steps `ntimesteps`.
The optional parameters `nsave` and `nwrite` determine the intervals for
storing the solution and writing to file, i.e., if `nsave > 1` only every
`nsave`'th solution is actually stored, and every `nwrite`'th time step the
solution is stored to disk.

The other constructors, either passing a `TimeSeries` and two `DataSeries` or a
filename are used to read data from previous simulations.

"""
mutable struct SolutionPODE{dType,tType,N} <: DeterministicSolution{dType,tType,N}
    nd::Int
    nt::Int
    ni::Int
    t::TimeSeries{tType}
    q::DataSeries{dType,N}
    p::DataSeries{dType,N}
    ntime::Int
    nsave::Int
    nwrite::Int
    offset::Int
    counter::Vector{Int}
    periodicity::dType

    function SolutionPODE{dType,tType,N}(nd, nt, ni, t, q, p, ntime, nsave, nwrite, periodicity = zero(q[begin])) where {dType<:Union{Number,AbstractArray},tType<:Real,N}
        new(nd, nt, ni, t, q, p, ntime, nsave, nwrite, 0, zeros(Int, ni), periodicity)
    end
end

function SolutionPODE(equation::Union{PODE{DT,TT1,AT},HODE{DT,TT1,AT},IODE{DT,TT1,AT},LODE{DT,TT1,AT}}, Δt::TT2, ntimesteps::Int;
    nsave::Int = DEFAULT_NSAVE, nwrite::Int = DEFAULT_NWRITE, filename = nothing) where {DT,TT1,TT2,AT}
    @assert nsave > 0
    @assert ntimesteps == 0 || ntimesteps ≥ nsave
    @assert nwrite == 0 || nwrite ≥ nsave
    @assert mod(ntimesteps, nsave) == 0

    if nwrite > 0
        @assert mod(nwrite, nsave) == 0
        @assert mod(ntimesteps, nwrite) == 0
    end

    TT = promote_type(TT1, TT2)
    N = (nsamples(equation) > 1 ? 2 : 1)
    nd = ndims(equation)
    ni = nsamples(equation)
    nt = div(ntimesteps, nsave)
    nt = (nwrite == 0 ? nt : div(nwrite, nsave))
    nw = (nwrite == 0 ? ntimesteps : nwrite)

    @assert nd > 0
    @assert ni > 0
    @assert nt ≥ 0
    @assert nw ≥ 0

    t = TimeSeries{TT}(nt, Δt, nsave)
    q = DataSeries(equation.q₀, nt, ni)
    p = DataSeries(equation.p₀, nt, ni)
    s = SolutionPODE{AT,TT,N}(nd, nt, ni, t, q, p, ntimesteps, nsave, nw, periodicity(equation))
    set_initial_conditions!(s, equation)

    if !isnothing(filename)
        isfile(filename) ? @warn("Overwriting existing HDF5 file.") : nothing
        create_hdf5!(s, filename)
    end

    return s
end

function SolutionPODE(t::TimeSeries{TT}, q::DataSeries{AT,N}, p::DataSeries{AT,N}, ntimesteps::Int) where {AT,TT,N}
    @assert ndims(q) == ndims(p)
    @assert ntime(q) == ntime(p)
    @assert nsamples(q) == nsamples(p)

    # extract parameters
    nd = ndims(q)
    ni = nsamples(q)
    nt = t.n
    ns = div(ntimesteps, nt)

    @assert mod(ntimesteps, nt) == 0

    # create solution
    SolutionPODE{AT,TT,N}(nd, nt, ni, t, q, p, ntimesteps, ns, 0)
end

function SolutionPODE(h5io::SolutionHDF5)
    # get hdf5 file descriptor
    h5 = hdf5(h5io)

    # read attributes
    ntimesteps = read(attributes(h5)["ntime"])
    nsave = read(attributes(h5)["nsave"])
    nsamples = read(attributes(h5)["nsamples"])

    # reading data arrays
    t = TimeSeries(read(h5["t"]), nsave)
    q = fromarray(DataSeries, read(h5["q"]), nsamples)
    p = fromarray(DataSeries, read(h5["p"]), nsamples)

    # create solution
    SolutionPODE(t, q, p, ntimesteps)
end

function SolutionPODE(file::String)
    load(file) do h5io
        SolutionPODE(h5io::SolutionHDF5)
    end
end

Base.:(==)(sol1::SolutionPODE, sol2::SolutionPODE) = (
    sol1.nd == sol2.nd
    && sol1.nt == sol2.nt
    && sol1.ni == sol2.ni
    && sol1.t == sol2.t
    && sol1.q == sol2.q
    && sol1.p == sol2.p
    && sol1.ntime == sol2.ntime
    && sol1.nsave == sol2.nsave
    && sol1.nwrite == sol2.nwrite
    && sol1.offset == sol2.offset
    && sol1.counter == sol2.counter
    && sol1.periodicity == sol2.periodicity)

@inline GeometricBase.counter(sol::SolutionPODE) = sol.counter
@inline GeometricBase.offset(sol::SolutionPODE) = sol.offset
@inline GeometricBase.lastentry(sol::SolutionPODE) = sol.ni == 1 ? sol.counter[1] - 1 : sol.counter .- 1
@inline GeometricBase.nsamples(sol::SolutionPODE) = sol.ni
@inline GeometricBase.nsave(sol::SolutionPODE) = sol.nsave
@inline GeometricBase.ntime(sol::SolutionPODE) = sol.ntime
@inline GeometricBase.timesteps(sol::SolutionPODE) = sol.t
@inline GeometricBase.eachtimestep(sol::SolutionPODE) = 1:sol.nt*sol.nsave
@inline GeometricBase.periodicity(sol::SolutionPODE) = sol.periodicity

function set_initial_conditions!(sol::SolutionPODE, equ::Union{PODE,HODE,IODE,LODE})
    set_initial_conditions!(sol, equ.t₀, equ.q₀, equ.p₀)
end

function set_initial_conditions!(sol::SolutionPODE{AT}, t₀::Real, q₀::AT, p₀::AT) where {AT}
    set_data!(sol.q, q₀, 0)
    set_data!(sol.p, p₀, 0)
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function set_initial_conditions!(sol::SolutionPODE{AT}, t₀::Real, q₀::AbstractVector{AT}, p₀::AbstractVector{AT}) where {AT}
    for i in eachindex(q₀, p₀)
        set_data!(sol.q, q₀[i], 0, i)
        set_data!(sol.p, p₀[i], 0, i)
    end
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function get_initial_conditions!(sol::SolutionPODE{AT,TT}, asol::AtomicSolutionPODE{DT,TT,AT}, k, n = 1) where {DT,TT,AT<:AbstractArray{DT}}
    get_solution!(sol, asol.q, asol.p, n - 1, k)
    asol.t = sol.t[n-1]
    asol.q̃ .= 0
    asol.p̃ .= 0
end

function get_initial_conditions!(sol::SolutionPODE{AT}, q::AT, p::AT, k, n = 1) where {AT}
    get_solution!(sol, q, p, n - 1, k)
end

function get_initial_conditions(sol::SolutionPODE, k, n = 1)
    get_solution(sol, n - 1, k)
end

function get_solution!(sol::SolutionPODE{AT}, q::AT, p::AT, n, k = 1) where {AT}
    get_data!(sol.q, q, n, k)
    get_data!(sol.p, p, n, k)
end

function get_solution(sol::SolutionPODE{AT,TT,1}, n, k = 1) where {AT,TT}
    @assert k == 1
    (sol.t[n], sol.q[n], sol.p[n])
end

function get_solution(sol::SolutionPODE{AT,TT,2}, n, k = 1) where {AT,TT}
    (sol.t[n], sol.q[n, k], sol.p[n, k])
end

function set_solution!(sol::SolutionPODE, t, q, p, n, k)
    set_solution!(sol, q, p, n, k)
end

function set_solution!(sol::SolutionPODE{AT,TT}, asol::AtomicSolutionPODE{DT,TT,AT}, n, k = 1) where {DT,TT,AT<:AbstractArray{DT}}
    set_solution!(sol, asol.t, asol.q, asol.p, n, k)
end

function set_solution!(sol::SolutionPODE{AT,TT}, q::AT, p::AT, n, k = 1) where {AT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        set_data!(sol.p, p, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function GeometricBase.reset!(sol::SolutionPODE)
    reset!(sol.q)
    reset!(sol.p)
    compute_timeseries!(sol.t, sol.t[end])
    sol.counter .= 1
    sol.offset += sol.nt
end

function Base.close(solution::SolutionPODE)
    close(solution.h5)
end


text:
struct DTWHelper{D,T <: PrimitiveModel,S} <: AbstractDistance
    q::Vector{T}
    y::Vector{T}
    d::D
    storage::S
end

DTWHelper(d::AbstractDistance, q::TimeVaryingAR, y::TimeVaryingAR) = DTWHelper(
    map(m -> PrimitiveModel(d, m), q.models),
    map(m -> PrimitiveModel(d, m), y.models),
    d,
    nothing,
)

function DTWHelper(d::OptimalTransportRootDistance, q::TimeVaryingAR, y::TimeVaryingAR)
    Q = map(m -> PrimitiveModel(d, m), q.models)
    Y = map(m -> PrimitiveModel(d, m), y.models)
    n = length(Q[1].weights)
    T = eltype(Q[1].weights)
    C = zeros(T, n, n)
    DTWHelper(Q, Y, d, (SinkhornLogWorkspace(T, n, n), C))
end

function evaluate(d::DTWHelper{<:EuclideanRootDistance}, m1::PrimitiveModel, m2::PrimitiveModel; kwargs...)
    evaluate(d.d, m1.roots, m2.roots, m1.weights, m2.weights)
end

function evaluate(d::DTWHelper{<:OptimalTransportRootDistance}, m1::PrimitiveModel, m2::PrimitiveModel; kwargs...)
    d.d.divergence === nothing ||
        throw(ArgumentError("Divergences are not supported for efficient DTW"))
    w,C = d.storage
    distmat_euclidean!(C, m1.roots, m2.roots, d.d.p)
    Γ = sinkhorn_log!(w, C, m1.weights, m2.weights; β = d.d.β, kwargs...)[1]
    dot(Γ, C)
end

function DynamicAxisWarping.dtw_cost(d::DTWHelper, q, y, args...; kwargs...)
    DynamicAxisWarping.dtw_cost(q, y, d.d, args...; kwargs...)
end

"""
    DynamicAxisWarping.dtwnn(q::TimeVaryingAR, y::TimeVaryingAR, dist, rad::Int, args...; kwargs...)

Wrapper for `dtwnn`. To save allocations between multiple calls to `dtwnn`, you may manually create a helper object and call

```julia
using DynamicAxisWarping
h = DTWHelper(dist, q, y)
dtwnn(h, rad::Int, args...; kwargs...)
```
"""
function DynamicAxisWarping.dtwnn(q::TimeVaryingAR, y::TimeVaryingAR, dist, rad::Int, args...; kwargs...)
    h = DTWHelper(dist, q, y)
    DynamicAxisWarping.dtwnn(h, rad, args...; kwargs...)
end

function DynamicAxisWarping.dtwnn(d::DTWHelper, rad::Int, args...; kwargs...)
    w = DynamicAxisWarping.DTWWorkspace(d.q, d, rad)
    DynamicAxisWarping.dtwnn(w, d.y, args...; kwargs...)
end