text:
cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\PS8-factor")
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables, Optim, HTTP, GLM, ForwardDiff, DataFramesMeta, MultivariateStats

function problemset8()

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS8-factor/nlsy.csv"
df = CSV.read(HTTP.get(url).body)

ols = glm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr), df, Normal())
print(ols)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
asvabMat=convert(Matrix,df[:,r"asvab"])
corr_asvab=cor(asvabMat)
print(corr_asvab)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
ols2 = glm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabAR + asvabCS + asvabMK + asvabNO + asvabPC + asvabWK), df, Normal())
print(ols2)

# The coefficient on the variable schoolt turns negative by inserting asvab variables.
# So, this specification is clearly wrong since increase in 'school time' should increase the wages.
# So there must be a specification bias in this model.

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
M = fit(PCA, asvabMat'; maxoutdim=1)
asvab_principalcomp = MultivariateStats.transform(M, asvabMat')

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
M = fit(FactorAnalysis, asvabMat'; maxoutdim=1)
asvab_factor = MultivariateStats.transform(M, asvabMat')

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 6
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#This is the part where I cannot figure out how to code exactly.
#I think we will start be defining vectors then followed by maximum likelihood/loglikelihood function.

end
problemset8()

max_stars_repo_path: ProblemSets/PS8_Ahmed Chaudhry.jl

max_issues_repo_path: ProblemSets/PS8_Ahmed Chaudhry.jl

max_forks_repo_path: ProblemSets/PS8_Ahmed Chaudhry.jl


text:
"""
    KrylovWrapper(hamiltonian)

Proxy object used in the Krylov iterations, during orbital
improvement. This is useful, since `hamiltonian` may be defined to act
on objects such as vectors living in function spaces (as
e.g. implemented using ContinuumArrays.jl), whereas the SCF iterations
act on the coefficients directly.
"""
mutable struct KrylovWrapper{T,Hamiltonian}
    hamiltonian::Hamiltonian
end

Base.eltype(A::KrylovWrapper{T}) where T = T
Base.size(A::KrylovWrapper{T,<:AbstractMatrix}) where T = size(A.hamiltonian)
"""
    size(::KrylovWrapper)

Returns the dimension of the `KrylovWrapper`. For Hamiltonians which
are not `<:AbstractMatrix`, this needs to be overloaded.
"""
Base.size(A::KrylovWrapper) = size(A.hamiltonian)
Base.size(A::KrylovWrapper, i) = size(A)[i]

"""
    mul!(y, ::KrylovWrapper, x[, α=1, β=0])

Compute the action of the wrapped Hamiltonian on `x` and store it in
`y`
"""
LinearAlgebra.mul!(y, A::KrylovWrapper, x, α::Number=true, β::Number=false) =
    mul!(y, A.hamiltonian, x, α, β)

Base.show(io::IO, kw::KrylovWrapper{T}) where T =
    write(io, "KrylovWrapper{$T} of size $(size(kw))")

"""
    OrthogonalKrylovWrapper(h,A,x′,C,S,tmp)

Proxy object to perform Krylov iterations of the Hamiltonian `h`,
represented by the [`KrylovWrapper`](@ref) `A`, in the subspace
orthogonal to `C`. The component of the test vector orthogonal to `C`
is stored in `x′`, before acting with `A`, and finally orthogonalized
again, i.e. `y = (1-C*C')*A*(1-C*C')*x = (1-C*C')*A*x′`. `S` is the
overlap matrix of the basis functions, used to compute the inner
products.
"""
struct OrthogonalKrylovWrapper{H,K,V,M,SM}
    h::H
    A::K
    x′::V
    C::M
    S::SM
    SC::M
end

function OrthogonalKrylovWrapper(h, A, nC, S)
    m = size(A, 1)
    T = eltype(A)
    x′ = zeros(T, m)
    C = zeros(T, m, nC)
    SC = zeros(T, m, nC)
    OrthogonalKrylovWrapper(h, A, x′, C, S, SC)
end

Base.eltype(A::OrthogonalKrylovWrapper) = eltype(A.A)
Base.size(A::OrthogonalKrylovWrapper) = size(A.A)
Base.size(A::OrthogonalKrylovWrapper, i) = size(A.A, i)

function orthogonalize!(x′, C, SC, x)
    for j = 1:size(C,2)
        c = dot(view(SC, :, j), x)
        BLAS.axpy!(-c, view(C, :, j), x′)
    end
end

function LinearAlgebra.mul!(y, A::OrthogonalKrylovWrapper, x)
    # We cannot simply project `x` into the subspace orthogonal to
    # `C`, since `x` is a test vector from e.g. a Krylov iteration,
    # i.e. it is part of the Krylov space and must not be
    # modified. Hence, we make our own copy, and project out the
    # components parallel to `C` from that one instead.
    copyto!(A.x′, x)
    # This step is necessary for non-orthogonal bases, since the inner
    # product is defined as `a'S*b`, where `S=B'B`, which is
    # non-diagonal for non-orthogonal bases.
    mul!(A.SC, A.S, A.C)

    # x′ <- (1 - C*C')*x
    orthogonalize!(A.x′, A.C, A.SC, x)
    # y <- A*x′
    mul!(reshape(y, :), A.A, A.x′)
    # y <- (1 - C*C')*y
    orthogonalize!(y, A.C, A.SC, y)
end

Base.show(io::IO, okw::OrthogonalKrylovWrapper) =
    write(io, "OrthogonalKrylovWrapper of size $(size(okw)), orthogonal against subspace of $(size(okw.C,2)) vectors")

ext: jl


text:
using Random
include("secrets.jl")
# FLAG is defined in secrets

f(a::Float64) = 4*round(Int64, a,RoundDown)+1

f(a::Int64) = 2*a

function f(a::String)
    s = Vector{UInt8}(a)
    x = Random.randsubseq(1:length(s),0.5)
    y = setdiff(1:length(s),x)

    o = Vector{Any}(undef,length(s))
    o[x] = convert.(Float64,s[x])
    o[y] = convert.(Int64,s[y])

    f.(o)
end

println(f(FLAG))
# [204, 216, 389, 413, 493, 437, 234, 197, 465, 421, 224, 433, 205, 381, 401, 421, 213, 449, 209, 232, 198, 208, 381, 98, 461, 190, 209, 477, 202, 213, 193, 437, 405, 250]
# the above output was observed

max_stars_repo_name: CSEA-IITB/Assignment

max_issues_repo_name: CSEA-IITB/Assignment

max_forks_repo_name: CSEA-IITB/Assignment


text:
##### NEW KEYS #####

struct ShapeKey <: Gadfly.GuideElement
    title::AbstractString
    labels::Vector{String}
    pos::Vector
end
ShapeKey(;title="Shape", labels=[""], pos=Float64[]) = ShapeKey(title, labels, pos)

"""
    Guide.shapekey[(; title="Shape", labels=[""], pos=Float64[])]
    Guide.shapekey(title, labels, pos)

Enable control of the auto-generated shapekey.  Set the key `title` and the item `labels`.
`pos` overrides [Theme(key_position=)](@ref Gadfly) and can be in either
relative (e.g. [0.7w, 0.2h] is the lower right quadrant), absolute (e.g. [0mm,
0mm]), or plot scale (e.g. [0,0]) coordinates.
"""
const shapekey = ShapeKey

function Guide.render(guide::Guide.ShapeKey, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    (theme.key_position == :none) && return Gadfly.Guide.PositionedGuide[]
    gpos = guide.pos
    (theme.key_position == :inside) && (gpos == Float64[]) &&  (gpos = [0.7w, 0.25h])

    # Aesthetics for keys: shape_key_title, shape_label (Function), shape_key_shapes (AbstractDict)
    nshapes = length(unique(aes.shape))
    guide_title = (guide.title!="Shape" || aes.shape_key_title==nothing) ? guide.title : aes.shape_key_title
    shape_key_labels = !(guide.labels==[""]) ? guide.labels : aes.shape_label(1:nshapes)

    colors = [nothing]
    if (aes.shape_key_title !=nothing)  && (aes.color_key_title==aes.shape_key_title)
        colors = collect(keys(aes.color_key_colors))
    end

    title_context, title_width = Guide.render_key_title(guide_title, theme)
    ctxs = render_discrete_key(shape_key_labels, title_context, title_width, theme, shapes=1:nshapes, colors=colors)

    position = right_guide_position
    if gpos != Float64[]
        position = over_guide_position
        ctxs = [compose(context(), (context(gpos[1],gpos[2]), ctxs[1]))]
    elseif theme.key_position == :left
        position = left_guide_position
    elseif theme.key_position == :top
        position = top_guide_position
    elseif theme.key_position == :bottom
        position = bottom_guide_position
    end

    return [Guide.PositionedGuide(ctxs, 0, position)]
end

function render_discrete_key(labels::Vector{String}, title_ctx::Context, title_width::Measure, theme::Gadfly.Theme;
    colors=[nothing], aes_color_label=nothing, shapes=[nothing])

    n = max(length(colors), length(shapes))
    shape1 = shapes[1]
    shapes = (shape1==nothing) ? fill(theme.key_swatch_shape, n) : theme.point_shapes[shapes]
    (colors[1]==nothing) && (colors = fill((theme.key_swatch_color==nothing) ? theme.default_color : theme.key_swatch_color, n))

    # only consider layouts with a reasonable number of columns
    maxcols = theme.key_max_columns < 1 ? 1 : theme.key_max_columns
    maxcols = min(n, maxcols)

    extents = text_extents(theme.key_label_font,
                        theme.key_label_font_size,
                        values(labels)...)

    ypad = 1.0mm
    title_height = title_ctx.box.a[2]
    entry_height = maximum([height for (width, height) in extents]) + ypad
    swatch_size = entry_height / 2

    # return a context with a lyout of numcols columns
    function make_layout(numcols)
        colrows = Array{Int}(undef, numcols)
        m = n
        for i in 1:numcols
            colrows[i] = ceil(Int, m/(1+numcols-i))
            m -= colrows[i]
        end

        xpad = 1mm
        colwidths = Array{Measure}(undef, numcols)
        m = 0
        for (i, nrows) in enumerate(colrows)
            if m == n
                colwidths[i] = 0mm
            else
                colwidth = maximum([width for (width, height) in extents[m+1:m+nrows]])
                colwidth += swatch_size + 2xpad
                colwidths[i] = colwidth
                m += nrows
            end
        end

        ctxwidth = sum(colwidths)
        ctxheight = entry_height * colrows[1] + title_height

        ctxp = ctxpromise() do draw_context
            yoff = 0.5h - ctxheight/2
            outerctx = context()

            compose!(outerctx, (context(xpad, yoff), title_ctx))

            ctx = context(0, yoff + title_height,
                        ctxwidth, ctxheight - title_height,
                        units=UnitBox(0, 0, 1, colrows[1]))

            m = 0
            xpos = 0w
            for (i, nrows) in enumerate(colrows)
                colwidth = colwidths[i]

                x = [0.5cy]
                clrs = colors[m+1:m+nrows]
                shps = shapes[m+1:m+nrows]
                swatches_shapes = [f(x, [y].*cy, [swatch_size/1.5]) for (y,f) in enumerate(shps)]
                sw1 = [(context(), s, fill(c), fillopacity(theme.alphas[1]), stroke(theme.discrete_highlight_color(c)))
                     for (s,c) in zip(swatches_shapes, clrs)]
                swatches = compose!(context(), sw1...)

                swatch_labels = compose!(
                    context(),
                    text([2xpad + swatch_size], [y*cy for y in 1:nrows],
                        collect(values(labels))[m+1:m+nrows], [hleft], [vcenter]),
                    font(theme.key_label_font),
                    fontsize(theme.key_label_font_size),
                    fill(theme.key_label_color))

                col = compose!(context(xpos, yoff), swatches, swatch_labels)
                if aes_color_label != nothing
                    classes = [svg_color_class_from_label(aes_color_label([c])[1]) for c in clrs]
                    #class_jscalls = ["data(\"color_class\", \"$(c)\")" for c in classes]
                    compose!(col,
                        svgclass(classes),
                        jscall(["""
                            data(\"color_class\", \"$(c)\")
                            .click(Gadfly.colorkey_swatch_click)
                            """ for c in classes]))
                end
                compose!(ctx, col)

                m += nrows
                xpos += colwidths[i]
            end

            return compose!(outerctx, ctx,
                            # defeat webkit's asinine default drag behavior
                            jscall("drag(function() {}, function() {}, function() {})"),
                            svgclass("guide colorkey"))
        end

        return compose!(
            context(minwidth=max(title_width, ctxwidth),
                    minheight=ctxheight,
                    units=UnitBox()),
            ctxp)
    end

    return map(make_layout, 1:maxcols)
end


text:
#!/usr/bin/env julia

#=
-------------------------
Log related functionality
-------------------------
=#
struct RheoLogItem
   action      # Nothing, or NamedTuple with fields:
               #     type::Symbol, funct::Symbol, params::NamedTuple, keywords::NamedTuple
   info        # NamedTuple for process output, comments, etc.
end

const RheoLog = Vector{RheoLogItem}

# simply pass text as comment
function RheoLogItem(s::String)
   return(RheoLogItem(Nothing,(comment=s,)))
end

# store info data in log from arbitrary keyword arguments
function RheoLogItem(;kwargs...)
   return(RheoLogItem(Nothing, kwargs.data))
end

#  There is a bug in julia that prevents the proper use of these functions.
#

# function Base.show(io::IO, rli::RheoLogItem)
#     println()
#     print("action = "); println(rli.action)
#     print("info   = "); println(rli.info)
# end
# function Base.display(io::IO, rl::RheoLog)
#     println("Hello!")
#      for rli in rl
#          print(rli)
#      end
# end

#=
-------------------------------
Time data related functionality
-------------------------------
=#
"""
    RheoTimeData(;σ::Vector{T1}, ϵ::Vector{T2}, t::Vector{T3}) where {T1<:Real, T2<:Real, T3<:Real}

`RheoTimeData` struct contains stress, strain and time data.

If preferred, an instance can be generated manually by just providing the three data
vectors in the right order, sampling type will be checked automatically.

# Fields

- `σ`: stress
- `ϵ`: strain
- `t`: time
- `log`: a log of struct's events, e.g. preprocessing
"""
struct RheoTimeData

    σ::Vector{RheoFloat}
    ϵ::Vector{RheoFloat}
    t::Vector{RheoFloat}

    log::Union{RheoLog,Nothing}

end

function RheoTimeData(;ϵ::Vector{T1} = RheoFloat[], σ::Vector{T2} = RheoFloat[], t::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}
    typecheck = check_time_data_consistency(t,ϵ,σ)
    RheoTimeData(convert(Vector{RheoFloat},σ), convert(Vector{RheoFloat},ϵ), convert(Vector{RheoFloat},t),
                    log == nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
end

@enum TimeDataType invalid_time_data=-1 time_only=0 strain_only=1 stress_only=2 strain_and_stress=3

function check_time_data_consistency(t,e,s)
    @assert (length(t)>0)  "Time data empty"

    sdef=(s != RheoFloat[])
    edef=(e != RheoFloat[])
    if (sdef && edef)
        @assert (length(s)==length(t)) && (length(e)==length(t)) "Time data length inconsistent"
        return strain_and_stress
    end

    if (sdef && (!edef))
        @assert (length(s)==length(t)) "Time data length inconsistent"
        return stress_only
    end

    if (edef && (!sdef))
        @assert (length(e)==length(t)) "Time data length inconsistent"
        return strain_only
    end

    if ((!edef) && (!sdef))
        return time_only
    end

    return invalid_time_data

end

function RheoTimeDataType(d::RheoTimeData)
    return check_time_data_consistency(d.t,d.ϵ,d.σ)
end

"""
    hastime(d::RheoTimeData)

returns 'true' if d contains a time array
"""
function hastime(d::RheoTimeData)
    return (d.t != RheoFloat[])
end

"""
    hasstress(d::RheoTimeData)

returns 'true' if d contains stress data
"""
function hasstress(d::RheoTimeData)
    return (d.σ != RheoFloat[])
end

"""
    hasstrain(d::RheoTimeData)

returns 'true' if d contains strain data
"""
function hasstrain(d::RheoTimeData)
    return (d.ϵ != RheoFloat[])
end

@enum LoadingType strain_imposed=1 stress_imposed=2

function Base.show(io::IO, d::RheoTimeData)
    b = length(d.t) > 10
    n = b ? 10 : length(d.t)

    if d.t != RheoFloat[]
        print("t =\t")
        for i = 1:n
            print("$(d.t[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.ϵ != RheoFloat[]
        print("ϵ =\t")
        for i = 1:n
            print("$(d.ϵ[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.σ != RheoFloat[]
        print("σ =\t")
        for i = 1:n
            print("$(d.σ[i])\t")
        end
        println( b ? "..." : "")
    end
end

function +(d1::RheoTimeData, d2::RheoTimeData)

    type1 = RheoTimeDataType(d1)
    type2 = RheoTimeDataType(d2)
    @assert (type1!=invalid_time_data) "Addition error: first parameter invalid"
    @assert (type2!=invalid_time_data) "Addition error: second parameter invalid"
    @assert (type1==type2) "Addition error: parameters inconsistent"
    @assert (type1!=time_only) "Addition error: time only data cannot be added"
    @assert (d1.t == d2.t) "Addition error: timelines inconsistent"

    log =   if (d1.log == nothing) || (d2.log == nothing)
                nothing
            else
                [ RheoLogItem( (type = :source, funct = :+, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
            end

    if (type1==strain_only) && (type2==strain_only)
        return RheoTimeData(RheoFloat[], d1.ϵ+d2.ϵ, d1.t, log)
    end

    if (type1==stress_only) && (type2==stress_only)
        return RheoTimeData(d1.σ+d2.σ, RheoFloat[], d1.t, log)
    end

    if (type1==strain_and_stress) && (type2==strain_and_stress)
        return RheoTimeData(d1.σ+d2.σ, d1.ϵ+d2.ϵ, d1.t, log)
    end

end

function -(d1::RheoTimeData, d2::RheoTimeData)

    type1 = RheoTimeDataType(d1)
    type2 = RheoTimeDataType(d2)
    @assert (type1!=invalid_time_data) "Subtraction error: first parameter invalid"
    @assert (type2!=invalid_time_data) "Subtraction error: second parameter invalid"
    @assert (type1==type2) "Subtraction error: parameters inconsistent"
    @assert (type1!=time_only) "Subtraction error: time only data cannot be added"
    @assert (d1.t == d2.t) "Subtraction error: timelines inconsistent"

    log =   if (d1.log == nothing) || (d2.log == nothing)
                nothing
            else
                [ RheoLogItem( (type = :source, funct = :-, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
            end

    if (type1==strain_only) && (type2==strain_only)
        return RheoTimeData(RheoFloat[], d1.ϵ-d2.ϵ, d1.t, log)
    end

    if (type1==stress_only) && (type2==stress_only)
        return RheoTimeData(d1.σ-d2.σ, RheoFloat[], d1.t, log)
    end

    if (type1==strain_and_stress) && (type2==strain_and_stress)
        return RheoTimeData(d1.σ-d2.σ, d1.ϵ-d2.ϵ, d1.t, log)
    end

end

function -(d::RheoTimeData)

    type = RheoTimeDataType(d)
    @assert (type!=invalid_time_data) "unary - error: parameter invalid"
    @assert (type!=time_only) "unary - error: time only data cannot be manipulated this way"

    log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:-, params=(), keywords=() ), () ) ]

    return RheoTimeData(-d.σ, -d.ϵ, d.t, log)

end

# Required to execute +/- operators for logs

function +(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) + rheologrun(rl2))
end

function -(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) - rheologrun(rl2))
end

function *(operand::Real, d::RheoTimeData)

    type = RheoTimeDataType(d)
    @assert (type!=invalid_time_data) "* error: parameter invalid"
    @assert (type!=time_only) "* error: time only data cannot be manipulated this way"

    log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:*, params=(operand = operand), keywords=() ), () ) ]

    return( RheoTimeData(operand .* d.σ, operand .* d.ϵ, d.t, log) )
end

function *(d::RheoTimeData, operand::Real)
    return(operand * d)
end

function |(d1::RheoTimeData, d2::RheoTimeData)
    # get types and type-check
    type1 = RheoTimeDataType(d1)
    type2 = RheoTimeDataType(d2)

    @assert ((type1==strain_only && type2==stress_only) || (type1==stress_only && type2==strain_only)) "One RheoTimeData must be stress_only and the other must be strain_only"

    # call recursively if not expected way round
    (type1==strain_only && type2==stress_only) && return d2|d1

    # expected case, type1==stress_only, type2==strain_only
    @assert d1.t==d2.t "Data to be combined must have same time-series"

    log = ((d1.log == nothing) || (d2.log == nothing)) ? nothing : [RheoLogItem((type = :source, funct = :|, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ())]

    return RheoTimeData(d1.σ, d2.ϵ, d1.t, log)
end
#=
------------------------------------
Frequency data related functionality
------------------------------------
=#
"""
    RheoFreqData(Gp::Vector{T1}, Gpp::Vector{T2}, ω::Vector{T3}, log::OrderedDict{Any,Any}) where {T1<:Real, T2<:Real, T3<:Real}

`RheoFreqData` contains storage modulus, loss modulus and frequency data.

If preferred, an instance can be generated manually by just providing the three data
vectors in the right order.

# Fields

- `Gp`: storage modulus
- `Gpp`: loss modulus
- `ω`: frequency
- `log`: a log of struct's events, e.g. preprocessing
"""
struct RheoFreqData
    Gp::Vector{RheoFloat}
    Gpp::Vector{RheoFloat}
    ω::Vector{RheoFloat}

    log::Union{RheoLog,Nothing}
end

function RheoFreqData(;Gp::Vector{T1} = RheoFloat[], Gpp::Vector{T2} = RheoFloat[], ω::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}
    typecheck = check_freq_data_consistency(ω,Gp,Gpp)
    RheoFreqData(convert(Vector{RheoFloat},Gp), convert(Vector{RheoFloat},Gpp), convert(Vector{RheoFloat},ω),
    log == nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
end

@enum FreqDataType invalid_freq_data=-1 freq_only=0 with_modulus=1

function check_freq_data_consistency(o,gp,gpp)
    @assert (length(o)>0)  "Freq data empty"

    gpdef=(gp != RheoFloat[])
    gppdef=(gpp != RheoFloat[])
    if (gpdef && gppdef)
        @assert (length(gp)==length(o)) && (length(gpp)==length(o)) "Data length inconsistent"
        return with_modulus
    end

    if ((!gpdef) && (!gppdef))
        return freq_only
    end

    return invalid_freq_data

end

function RheoFreqDataType(d::RheoFreqData)
    return check_freq_data_consistency(d.ω,d.Gp,d.Gpp)
end

function Base.show(io::IO, d::RheoFreqData)
    b = length(d.ω) > 10
    n = b ? 10 : length(d.ω)

    if d.ω != RheoFloat[]
        print("ω =\t")
        for i = 1:n
            print("$(d.ω[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.Gp != RheoFloat[]
        print("Gp =\t")
        for i = 1:n
            print("$(d.Gp[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.Gpp != RheoFloat[]
        print("Gpp =\t")
        for i = 1:n
            print("$(d.Gpp[i])\t")
        end
        println( b ? "..." : "")
    end
end

#=
---------------------------------------
Log-based data generation functionality
---------------------------------------
=#
function rheologrun(rli::RheoLogItem, d=nothing)

      if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
         type=rli.action.type
         if type==:source && d==nothing
            return(eval(rli.action.funct)(rli.action.params...;rli.action.keywords...))
         elseif type==:process && d!=nothing
            return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
         elseif type==:analysis && d!=nothing
            return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
         end
      end
   println(rli.info)
end

"""
    rheologrun(log::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} = nothing)

execute all actions from the log. It applies them to the data `d` provided, or use the log's first action to reload/recreate it otherwise.
"""
function rheologrun(arli::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} = nothing)

  # check first item is a source item
  if d == nothing
        @assert typeof(arli[1].action) <: NamedTuple && :type in keys(arli[1].action) && arli[1].action.type == :source
                "Source missing in RheoLog"
        end

  for rli in arli
      if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
          type=rli.action.type
          if type==:analysis && d!=nothing
              rheologrun(rli, d)
          else  d=rheologrun(rli, d)
          end
      else
          println(rli.info)
      end
  end
  return(d)
end

"""
    showlog(d::Union{RheoTimeData,RheoFreqData})

shows the record of operations on a rheological data.
"""
function showlog(d::Union{RheoTimeData,RheoFreqData})
    if d.log != nothing
        for idx in 1:length(d.log)
            println(idx)
            print("     action = "); println(d.log[idx].action)
            print("     info   = "); println(d.log[idx].info)
        end
    else
        println("No log data available")
    end
end

#=
---------------------------
Model related functionality
---------------------------
=#
"""
    RheoModelClass(name::String, params::Vector{Symbol}, _G::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}, _Ga::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}, _J::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}, _Ja::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}, _Gp::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}, _Gpa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}, _Gpp::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}, _Gppa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}, constraint::FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}}, info::String, expressions::NamedTuple)

`RheoModelClass` contains a model name, it's symbolic parameters and all its moduli (both single-input and array-input versions).

It also contains information about any constraints that must be observed (e.g. the springpot coefficient being inbetween 0 and 1).

Lastly, it also contains additional info about the model which may include a text-art schematic.

Generally, users will want to use the `RheoModelClass` constructor function as shown in the 'Create Your Model' section of the documentation
rather than the default constructor.
"""
struct RheoModelClass

    name::String
    params::Vector{Symbol}

    _G::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
    _Ga::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
    _J::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
    _Ja::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
    _Gp::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
    _Gpa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
    _Gpp::FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
    _Gppa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}

    constraint::FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}}

    info::String
    expressions::NamedTuple

end

function Base.show(io::IO, m::RheoModelClass)
    print(io, "\nModel name: $(m.name)")
    ps=join([string(s) for s in m.params], ", ", " and ")
    print(io, "\n\nFree parameters: $ps\n")
    print(io, m.info)
    return
end

rheoconv(t::Real) = RheoFloat(t)
rheoconv(t::Array{T,1}) where T<:Real = convert(Vector{RheoFloat},t)

invLaplace(f::Function, t::Array{RheoFloat}) = InverseLaplace.talbotarr(f, t)
invLaplace(f::Function, t::RheoFloat) = InverseLaplace.talbot(f, t)

# place holder for undefined moduli/compliance functions
const nanexp = quote NaN end

function RheoModelClass(;name::String,
                         p::Array{Symbol},
                         G::Expr = nanexp,
                         J::Expr = nanexp,
                         Gp::Expr = nanexp,
                         Gpp::Expr = nanexp,
                         constraint::Expr = quote true end,
                         info="", #name * ": model with parameters " * string(join(string.(p), ", "), "."),
                         # flags to avoid bugs related to the FunctionWrappers and MittLeff
                         Ga_safe::Bool = true,
                         Ja_safe::Bool = true
                         )

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))
    expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint, Ga_safe=Ga_safe, Ja_safe=Ja_safe)

    @eval return(RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end)                 |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end)  |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((t,params) -> begin $unpack_expr; $J; end)                 |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end)  |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end)                |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end)               |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        (params -> begin $unpack_expr; $constraint; end)            |> FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}},
        $info, $expressions) )
end

# Cool replacement function inspired from
# https://stackoverflow.com/questions/29778698/julia-lang-metaprogramming-turn-expression-into-function-with-expression-depend
# This probably could go in its own module as this is very useful.
expr_replace!(ex, s, v) = ex
expr_replace!(ex::Symbol, s, v) = s == ex ? v : ex

function expr_replace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = expr_replace!(ex.args[i], s, v)
    end
    return ex
end

expr_replace(ex, s, v) = expr_replace!(copy(ex), s, v)

function expr_replace(ex, nt)
    e = copy(ex)
    k=keys(nt)
    v=values(nt)
    for i in 1:length(nt)
        expr_replace!(e, k[i], v[i])
    end
    return e
end

function modulusexists(modulus_single_input)
   # 5.0 is not meaningful, just an input to check that the function does not return NaN
    # at a non-zero value, which is the default function output for moduli that have
    # not been defined.
    return !isnan(modulus_single_input(5.0))
end

#
#  Function used to check and convert any named tupple of parameters into an array ordered as expected by the moduli functions of the RheoModelClass
#

function check_and_reorder_parameters(nt::NamedTuple, params_order::Vector{Symbol}; err_string::String="calling function")
    # check that every parameter in m exists in the named tuple nt
    @assert all(i->i in keys(nt),params_order) "Missing parameter(s) in " * err_string
    # check that no extra parameters have been provided
    @assert length(params_order) == length(nt) "Mismatch number of model parameters and parameters provided in " * err_string

    p = map(i->RheoFloat(nt[i]), params_order)
end

"""
    freeze_params(m::RheoModelClass, nt0::NamedTuple)

Return a new `RheoModelClass` with some of the parameters frozen to specific values

# Arguments

- `m`: original `RheoModelClass`
- `nt0`: named tuple with values for each parameter to freeze

# Example
```@example
julia> SLS2_mod = freeze_params( SLS2, (G₀=2,η₂=3.5))
[...]

julia> SLS2.G(1,[2,1,2,3,3.5])
3.8796492

julia> SLS2_mod.G(1,[1,2,3])
3.8796492
```

"""
function freeze_params(m::RheoModelClass, nt0::NamedTuple)

    # check that every parameter in m exists in the named tuple nt
    @assert all( i-> i in m.params,keys(nt0)) "A parameter to freeze is not involved in the model"
    # convert values format for consistency
    nt=NamedTuple{keys(nt0)}([RheoFloat(i) for i in nt0])

    # create array of remaining variables
    p = filter(s -> !(s in keys(nt)),m.params)

    name="$(m.name) with set parameters: $nt"
    info = m.info

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))

    # This section creates moduli expressions with material parameters
    # replaced by specific values.

    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)
    constraint = expr_replace(m.expressions.constraint, nt)

    expressions=NamedTuple{(:G, :J, :Gp, :Gpp, :constraint, :Ga_safe, :Ja_safe)}( ( G, J, Gp, Gpp, constraint, m.expressions.Ga_safe, m.expressions.Ja_safe ) )

    @eval return( RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((t,params) -> begin $unpack_expr; $J; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        (params -> begin $unpack_expr; $constraint; end) |> FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}},
        $info, $expressions)   )
end

function freeze_params(m::RheoModelClass; kwargs...)
    return(freeze_params(m,kwargs.data))
end

"""
    RheoModel(_G::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}, _Ga::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}, _J::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}, _Ja::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}, _Gp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}, _Gpa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}, _Gpp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}, _Gppa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}, expressions::NamedTuple)

`RheoModel` contains all known moduli of a particular model, as for a `RheoModelClass` model name. However, a `RheoModel`
has all its parameters fixed to known values.

Generally, users will begin with a defined `RheoModelClass` (e.g. SLS) and then specialise the parameters,
rather than calling the default constructor explicitly.
"""
struct RheoModel

    _G::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Ga::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _J::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Ja::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _Gp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Gpa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _Gpp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Gppa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}

    expressions::NamedTuple
    params::NamedTuple
    info::String    # do we need this?
    # log::OrderedDict{Any,Any}

end

function RheoModel(m::RheoModelClass; kwargs...)
    return(RheoModel(m,kwargs.data))
end

function RheoModel(m::RheoModelClass, nt0::NamedTuple)

    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(nt0, m.params, err_string = "model definition")
    nt=NamedTuple{Tuple(m.params)}(p)

    # string printed
    info = string("\nModel: $(m.name)\n\nParameter values: $nt \n",m.info)

    # This attaches expressions to the models with parameters substituted by values.
    # Future development: this could be disabled with a key word if processing time is an issue
    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)

    expressions=NamedTuple{(:G, :J, :Gp, :Gpp, :Ga_safe, :Ja_safe)}( ( G, J, Gp, Gpp, m.expressions.Ga_safe, m.expressions.Ja_safe ) )

    @eval return( RheoModel(
    (t -> begin $G; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ta -> begin [$G for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (t -> begin $J; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ta -> begin [$J for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (ω -> begin $Gp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ωa -> begin [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (ω -> begin $Gpp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ωa -> begin [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    $expressions, $nt, $info) )
end

function Base.show(io::IO, m::RheoModel)
    print(io,m.info)
end

#
# These are internal functions to access the array versions of the relaxation and creep moduli
# They are used in processing.jl as well as in the user level moduli functions
#

function _Ga(m::RheoModelClass)
    if m.expressions.Ga_safe
        m._Ga
    else
        (ta,p) -> [m._G(t,p) for t in ta]
    end
end

function _Ja(m::RheoModelClass)
    if m.expressions.Ja_safe
        m._Ja
    else
        (ta,p) -> [m._J(t,p) for t in ta]
    end
end

function _Ga(m::RheoModel)
    if m.expressions.Ga_safe
        m._Ga
    else
        ta -> [m._G(t) for t in ta]
    end
end

function _Ja(m::RheoModel)
    if m.expressions.Ja_safe
        m._Ja
    else
        ta -> [m._J(t) for t in ta]
    end
end

#
# These are functions enabling users to access the different moduli functions for each model
#

"""
    relaxmod(m[, t, params])

provide access to the relaxation modulus (G) of a given model m. It can be used with a broad range of inputs.

* When a time value is provided (or an array of time values), the function returns the corresponding value(s) of the relaxation modulus.

    relaxmod(m::RheoModel, time (single value or array))

    relaxmod(m::RheoModelClass, time (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    relaxmod(Maxwell, 1, k=1., η=1)

    relaxmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, relaxmod returns the relaxation function itself.

    relaxmod(m::RheoModel)

    relaxmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    G = relaxmod(m)
    G([0,1,2])
"""
function relaxmod(m::RheoModel, t::Number)
    m._G(t)
end

function relaxmod(m::RheoModel, ta::Vector{T}) where T <: Number
    _Ga(m)(ta)
end

relaxmod(m) = x -> relaxmod(m,x)

function relaxmod(m::RheoModelClass, t::Number, params::Vector{T}) where T <: Number
    m._G(t, params)
end

function relaxmod(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    _Ga(m)(ta, params)
end

function relaxmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.params, err_string = "relaxmod")
    relaxmod(m, x, p)
end

function relaxmod(m::RheoModelClass, x; kwargs...)
    relaxmod(m, x, kwargs.data)
end

relaxmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> relaxmod(m, x, params)
relaxmod(m::RheoModelClass, params::NamedTuple) =  x -> relaxmod(m, x, params)
relaxmod(m::RheoModelClass; kwargs...) =  x -> relaxmod(m, x, kwargs.data)

"""
    creepcomp(m[, t, params])

provide access to the creep compliance (J) of a given model m. It can be used with a broad range of inputs.

* When a time value is provided (or an array of time values), the function returns the corresponding value(s) of the creep compliance.

    creepcomp(m::RheoModel, time (single value or array))

    creepcomp(m::RheoModelClass, time (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    creepcomp(Maxwell, 1, k=1., η=1)

    creepcomp(Maxwell, [1,2,3], [1,1])

* When no time value is provided, creepcomp returns the creep compliance function itself.

    creepcomp(m::RheoModel)

    creepcomp(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    J = creepcomp(m)
    J([0,1,2])
"""
function creepcomp(m::RheoModel, t::Number)
    m._J(t)
end

function creepcomp(m::RheoModel, ta::Vector{T}) where T <: Number
    _Ja(m)(ta)
end

creepcomp(m) = x -> creepcomp(m,x)

function creepcomp(m::RheoModelClass, t::Number, params::Vector{T}) where T <: Number
    m._J(t, params)
end

function creepcomp(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    _Ja(m)(ta, params)
end

function creepcomp(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.params, err_string = "relaxmod")
    creepcomp(m, x, p)
end

function creepcomp(m::RheoModelClass, x; kwargs...)
    creepcomp(m, x, kwargs.data)
end

creepcomp(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> creepcomp(m, x, params)
creepcomp(m::RheoModelClass, params::NamedTuple) =  x -> creepcomp(m, x, params)
creepcomp(m::RheoModelClass; kwargs...) =  x -> creepcomp(m, x, kwargs.data)

"""
    storagemod(m[, ω, params])

provide access to the storage modulus (G') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the storage modulus.

    storagemod(m::RheoModel, frequency (single value or array))

    storagemod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    storagemod(Maxwell, 1, k=1., η=1)

    storagemod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, storagemod returns the storage modulus function itself.

    storagemod(m::RheoModel)

    storagemod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gp = storagemod(m)
    Gp([0,1,2])
"""
function storagemod(m::RheoModel, ω::Number)
    m._Gp(ω)
end

function storagemod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gpa(ωa)
end

storagemod(m) = x -> storagemod(m,x)

function storagemod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gp(ω, params)
end

function storagemod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
      m._Gpa(ωa, params)
end

function storagemod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.params, err_string = "relaxmod")
    storagemod(m, x, p)
end

function storagemod(m::RheoModelClass, x; kwargs...)
    storagemod(m, x, kwargs.data)
end

storagemod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> storagemod(m, x, params)
storagemod(m::RheoModelClass, params::NamedTuple) =  x -> storagemod(m, x, params)
storagemod(m::RheoModelClass; kwargs...) =  x -> storagemod(m, x, kwargs.data)

"""
    lossmod(m[, ω, params])

provide access to the loss modulus (G'') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the loss modulus.

    lossmod(m::RheoModel, frequency (single value or array))

    lossmod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    lossmod(Maxwell, 1, k=1., η=1)

    lossmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, lossmod returns the loss modulus function itself.

    lossmod(m::RheoModel)

    lossmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gpp = lossmod(m)
    Gpp([0,1,2])
"""
function lossmod(m::RheoModel, ω::Number)
    m._Gpp(ω)
end

function lossmod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gppa(ωa)
end

lossmod(m) = x -> lossmod(m,x)

function lossmod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gpp(ω, params)
end

function lossmod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Gppa(ωa, params)
end

function lossmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.params, err_string = "relaxmod")
    lossmod(m, x, p)
end

function lossmod(m::RheoModelClass, x; kwargs...)
    lossmod(m, x, kwargs.data)
end

lossmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> lossmod(m, x, params)
lossmod(m::RheoModelClass, params::NamedTuple) =  x -> lossmod(m, x, params)
lossmod(m::RheoModelClass; kwargs...) =  x -> lossmod(m, x, kwargs.data)

"""
    dynamicmod(m[, ω, params])

provide access to the complex dynamic modulus (G' + i G'') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the complex dynamic modulus.

    dynamicmod(m::RheoModel, frequency (single value or array))

    dynamicmod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    dynamicmod(Maxwell, 1, k=1., η=1)

    dynamicmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, lossmod returns the complex dynamic modulus function itself.

    dynamicmod(m::RheoModel)

    dynamicmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gpp = dynamicmod(m)
    Gpp([0,1,2])

Note: use abs() and angle() to get the magnitude and phase of the complex modulus.
"""
function dynamicmod(m::RheoModel, ω::Number)
    m._Gp(ω) + m._Gpp(ω) * im
end

function dynamicmod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gpa(ωa) + m._Gppa(ωa) * im
end

dynamicmod(m) = x -> dynamicmod(m,x)

function dynamicmod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gp(ω) + m._Gpp(ω) * im
end

function dynamicmod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Gpa(ωa, params) + m._Gppa(ωa, params) * im
end

function dynamicmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.params, err_string = "relaxmod")
    dynamicmod(m, x, p)
end

function dynamicmod(m::RheoModelClass, x; kwargs...)
    dynamicmod(m, x, kwargs.data)
end

dynamicmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> dynamicmod(m, x, params)
dynamicmod(m::RheoModelClass, params::NamedTuple) =  x -> dynamicmod(m, x, params)
dynamicmod(m::RheoModelClass; kwargs...) =  x -> dynamicmod(m, x, kwargs.data)


text:
module Hydra

using Random, LinearAlgebra, Distances, NamedArrays, StatsModels,Statistics,TexTables,LoopVectorization
include("hydra2.jl")
include("output.jl")

export hydra2,
    permute,
    Euclidean,
    SqEuclidean,
    PeriodicEuclidean,
    Cityblock,
    TotalVariation,
    Chebyshev,
    Minkowski,
    Jaccard,
    BrayCurtis,
    RogersTanimoto,
    @formula

end


text:
using DiffEqUncertainty, DiffEqBase, OrdinaryDiffEq, DiffEqProblemLibrary
using Test

function g(du,u,p,t)
  σ,ρ,β = p
  x,y,z = u
  du[1] = σ*(y-x)
  du[2] = x*(ρ-z) - y
  du[3] = x*y - β*z
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,10.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(g,u0,tspan,p)

cb = ProbIntsUncertainty(1e4,5)
solve(prob,Tsit5())
monte_prob = EnsembleProblem(prob)
sim = solve(monte_prob,Tsit5(),trajectories=10,callback=cb,adaptive=false,dt=1/10)

#using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

function fitz(du,u,p,t)
  V,R = u
  du[1] = 3.0*(V - V^3/3 + R)
  du[2] = -(1/3.0)*(V -  0.2 - 0.2*R)
end
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
prob = ODEProblem(fitz,u0,tspan)

cb = ProbIntsUncertainty(0.1,1)
sol = solve(prob,Euler(),dt=1/10)
monte_prob = EnsembleProblem(prob)
sim = solve(monte_prob,Euler(),trajectories=100,callback=cb,adaptive=false,dt=1/10)

#using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

cb = AdaptiveProbIntsUncertainty(5)
sol = solve(prob,Tsit5())
monte_prob = EnsembleProblem(prob)
sim = solve(monte_prob,Tsit5(),trajectories=100,callback=cb,abstol=1e-3,reltol=1e-1)

#using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)


text:
struct AsymptoticCMETest{K<:Kernel,V,M,S} <: HypothesisTests.HypothesisTest
    """Kernel."""
    kernel::K
    """Number of observations."""
    nsamples::Int
    """Number of test locations."""
    ntestsamples::Int
    """UCME estimate."""
    estimate::V
    """Mean deviation for each test location."""
    mean_deviations::M
    """Test statistic."""
    statistic::S
end

function AsymptoticCMETest(
    estimator::UCME, predictions::AbstractVector, targets::AbstractVector
)
    # determine number of observations and test locations
    nsamples = length(predictions)
    testpredictions = estimator.testpredictions
    ntestsamples = length(testpredictions)

    # create matrix of deviations (observations × test locations)
    kernel = estimator.kernel
    deviations =
        CalibrationErrors.unsafe_ucme_eval.(
            (kernel,),
            predictions,
            targets,
            permutedims(testpredictions),
            permutedims(estimator.testtargets),
        )

    # compute UCME estimate
    mean_deviations = Statistics.mean(deviations; dims=1)
    estimate = sum(abs2, mean_deviations) / ntestsamples

    # compute test statistic
    C = LinearAlgebra.Symmetric(Statistics.covm(deviations, mean_deviations))
    x = vec(mean_deviations)
    statistic = nsamples * LinearAlgebra.dot(x, C \ x)

    return AsymptoticCMETest(kernel, nsamples, ntestsamples, estimate, x, statistic)
end

# HypothesisTests interface

HypothesisTests.default_tail(::AsymptoticCMETest) = :right

## have to specify and check keyword arguments in `pvalue` and `confint` to
## force `tail = :right` due to the default implementation in HypothesisTests

function HypothesisTests.pvalue(test::AsymptoticCMETest; tail=:right)
    if tail === :right
        StatsFuns.chisqccdf(test.ntestsamples, test.statistic)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

HypothesisTests.testname(::AsymptoticCMETest) = "Asymptotic CME test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::AsymptoticCMETest)
    return "Mean vector", zero(test.mean_deviations), test.mean_deviations
end

function HypothesisTests.show_params(io::IO, test::AsymptoticCMETest, ident="")
    println(io, ident, "number of observations: ", test.nsamples)
    println(io, ident, "number of test locations: ", test.ntestsamples)
    println(io, ident, "UCME estimate: ", test.estimate)
    return println(io, ident, "test statistic: ", test.statistic)
end


text:
# Test module for functions
module Tst
    include("../src/interpol.jl")
    import .Interpol
    #using .Interpol

    # Interpol.cartesian_interpolation(inMat1,inMat2,steps)
    # say_hello()

    # your other test code here
end

max_stars_repo_name: dasha-shchep/LIIC

max_issues_repo_name: dasha-shchep/LIIC

max_forks_repo_name: dasha-shchep/LIIC


text:
using Invest
using Test

@testset "American" begin
    american = AmericanCall(1.0, 2.0, 1.1, 0.9, 0.01, 5 / 12)
    usa = AmericanPut(1.0, 2.0, 1.1, 0.9, 0.01, 5 / 12, 1)
    @test isa(american, AmericanCall)
    @test isa(usa, AmericanPut)
end

@testset "European" begin
    european = EuropeanCall(1.0, 2.0, 0.1, 0.3, 5 / 12)
    euro = EuropeanPut(1.0, 2.0, 0.1, 0.3, 5 / 12, 1)
    @test isa(european, EuropeanCall)
    @test isa(euro, EuropeanPut)
end

# TODO: Add tests for American options with calculated σ

@testset "Pricing" begin
    e = EuropeanCall(90.0, 85.0, 0.1, 0.3, 6 / 12)
    @test isa(e, EuropeanCall)
    C, P = price(e)
    @test round(C; digits=2) == 12.69
    @test round(P; digits=2) == 3.54
    @test round(price(EuropeanCall(437.5, 460.0, 0.0125, 0.42, 12 / 12))[2]; digits=2) ==
          82.85
    @test round(price(EuropeanCall(20.0, 20.0, 0.06, 0.2, 3 / 12))[1]; digits=5) == 0.94938
    @test round(price(EuropeanCall(20.0, 20.0, 0.06, 0.2, 3 / 12))[2]; digits=5) == 0.65162
    @test round(price(EuropeanCall(50.0, 50.0, 0.10, 0.3, 3 / 12))[2]; digits=6) == 2.375941
end

# TODO: Add tests for underlying derivative price of put option

@testset "Greeks" begin
    e = EuropeanPut(40.0, 45.0, 0.08, 0.3, 4 / 12)
    @test isa(e, EuropeanPut)
    @test !isat(e)
    @test isat(EuropeanPut(2.0, 2.0, 0.01, 0.2, 1 / 12))
    @test round(Δ(e)[1]; digits=2) == 0.33
    @test round(Γ(e); digits=3) == 0.052
    @test round(Θ(e); digits=4) == -0.0129
    C, _ = dprice(e, 3, 1 / 52)
    @test round((C); digits=2) == 2.6
    @test round(Δ(EuropeanPut(1.0, 1.0, 0.10, 0.25, 0.5))[1]; digits=3) == 0.645
end


text:
using Turing, ContinuousBenchmarks.TuringTools, ContinuousBenchmarks.Reporter

data = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

log_report("Dummy Benchmark started!")

@model constrained_test(obs) = begin
    p ~ Beta(2,2)
    for i = 1:length(obs)
        obs[i] ~ Bernoulli(p)
    end
    p
end

log_report("Dummy model constructed!")

bench_res = @tbenchmark(HMC(1000, 1.5, 3), constrained_test, data)

log_report("Dummy benchmark finished!")

# bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
LOG_DATA = build_log_data("Dummy-Benchmark", bench_res...)
print_log(LOG_DATA)

log_report("Dummy benchmark reported!")


text:
"""
inspectQNITensor is a developer-level debugging tool
to look at internals or properties of QNITensors
"""
function inspectQNITensor(T::ITensor,is::QNIndexSet)
  #@show T.store.blockoffsets
  #@show T.store.data
  println("Block fluxes:")
  for b in nzblocks(T)
    @show flux(T,b)
  end
end
inspectQNITensor(T::ITensor,is::IndexSet) = nothing
inspectQNITensor(T::ITensor) = inspectQNITensor(T,inds(T))

"""
    pause()

Pauses execution until a key (other than 'q') is pressed.
Entering 'q' exits the program. The `pause()` function
is useful for inspecting output of programs at certain
points while giving the option to continue.
"""
function pause()
  print(stdout, "(Paused) ")
  c = read(stdin, 1)
  c == UInt8[0x71] && exit(0)
  return
end


text:
struct TimeGridPoint
  idx::Int64
  val::BranchPoint
end

struct TimeGrid <: AbstractVector{TimeGridPoint}
  contour::Contour
  points::Vector{TimeGridPoint}
  step::Vector{ComplexF64}
  branch_bounds::Vector{Tuple{TimeGridPoint, TimeGridPoint}}

  function TimeGrid(c::Contour; npts_real = 0, npts_imag = 0)
    points = TimeGridPoint[]
    step = ComplexF64[]
    branch_bounds = Tuple{TimeGridPoint, TimeGridPoint}[]

    for b in c.branches
      npts = (b.domain == imaginary_branch) ? npts_imag : npts_real
      @assert npts > 1 "every branch must have at least two points (b = $b, npts = $npts)"

      indices = (1:npts) .+ length(points)
      branch_points = b.(range(0, 1, length=npts))
      time_grid_points = TimeGridPoint.(indices, branch_points)

      append!(points, time_grid_points)
      append!(step, (b.max_val - b.min_val) / (npts - 1))
      push!(branch_bounds, (time_grid_points[1], time_grid_points[end]))
    end
    return new(c, points, step, branch_bounds)
  end
end

### AbstractArray Interface ###
IndexStyle(::Type{<:TimeGrid}) = IndexLinear()
size(grid::TimeGrid) = size(grid.points)
getindex(grid::TimeGrid, i::Int) = grid.points[i]
setindex!(grid::TimeGrid, v::TimeGridPoint, i::Int) = grid.points[i] = v

"get ith point on branch b of grid"
function getindex(grid::TimeGrid, b::BranchEnum, i::Int)
  bounds = branch_bounds(grid, b)
  return grid[bounds[1].idx + i - 1]
end

"get all points on branch b of grid"
function getindex(grid::TimeGrid, b::BranchEnum)
  bounds = branch_bounds(grid, b)
  return grid[bounds[1].idx:bounds[2].idx]
end

function branch_bounds(grid::TimeGrid, b::BranchEnum)
  return grid.branch_bounds[grid.contour.branch_indices[Int(b)]]
end

function length(grid::TimeGrid, b::BranchEnum)
  bounds = branch_bounds(grid, b)
  return bounds[2].idx - bounds[1].idx + 1
end

function step(grid::TimeGrid, b::BranchEnum)
  return grid.step[grid.contour.branch_indices[Int(b)]]
end

function integrate(f, grid::TimeGrid, t1::TimeGridPoint, t2::TimeGridPoint)
  # wrap around
  (t1.idx < t2.idx) && return integrate(f, grid, t1, grid[1]) + integrate(f, grid, grid[end], t2)

  integral = zero(fieldtype(BranchPoint, :val)(0.0) * f(t2)) # extra evaluation just to get the type

  first_branch_idx = grid.contour.branch_indices[Int(t2.val.domain)]
  last_branch_idx =  grid.contour.branch_indices[Int(t1.val.domain)]

  first = t2
  for b in first_branch_idx:last_branch_idx
    Δt = grid.step[b]
    last = (b == last_branch_idx ? t1 : grid.branch_bounds[b][2])
    @assert first.val.domain == last.val.domain
    if (last.idx != first.idx) #trapezoid rule
      branch_integral = 0.5 * (f(grid[first.idx]) + f(grid[last.idx]))
      for i in (first.idx+1):(last.idx-1)
        t = @inbounds grid[i]
        branch_integral += f(t)
      end
      integral += (Δt * branch_integral)
    end
    b == last_branch_idx || (first = grid[last.idx + 1])
  end
  return integral
end

integrate(f, grid::TimeGrid) = integrate(f, grid, grid[end], grid[1])

realtimes(grid::TimeGrid) = map(t -> real(t.val.val), grid[forward_branch])
imagtimes(grid::TimeGrid) = map(t -> real(1.0im * t.val.val), grid[imaginary_branch]) # τ = it


text:
# An abstract type for fixed-point acceleration methods
# Fixed point problem x = g(x) with residual f(x) = x - g(x)
export AbstractAccelerator
"""
    AbstractAccelerator

Abstract supertype for acceleration objects that can be used to speed up a fixed-point iterations g = g(x) of a nonexpansive operator `g`.
They must implement the following methods to communicate with the fixed-point algorithm:
  - update!(aa::AbstractAccelerator, g, x, num_iter) #stores the fixed-point iterates
  - accelerate!(g::AbstractVector, x, aa::AbstractAccelerator, num_iter) #recombines past iterates to determine an accelerated point and overwrites `g`
  - restart!(aa::AbstractAccelerator, args...; kwargs...) # algorithm tells the accelerator to restart

The algorithm has to be able to query the following information:
  - was_successful(aa::AbstractAccelerator) --> Bool #indicate whether accelerate! was succesful at the last iteration

The following method is optional:
  - log!(aa::AbstractAccelerator, args...; kwargs...) # algorithm tells accelerator to log certain information for debugging
"""
abstract type AbstractAccelerator end

# update!(aa::AbstractAccelerator, g, x, iter) = nothing
# accelerate!(g, x, aa::AbstractAccelerator, iter) = nothing
# was_successful(aa::AbstractAccelerator) = false
# restart!(aa::AbstractAccelerator) = nothing

log!(aa::AbstractAccelerator, args...; kwargs...) = false

max_issues_repo_name: oxfordcontrol/COSMOAccelerators

max_forks_repo_name: oxfordcontrol/COSMOAccelerators


text:
using Test, MLToolkit

@testset "Neural SBC" begin
    i_dim = 20
    k_init = 10
    batch_size = 50

    @testset "$sbc" for sbc in [MeanFieldSBC(i_dim, k_init), StructuredSBC(i_dim, k_init)]
        is_all_true = true
        @test begin
            for k_out = 1:k_init
                x = AT(rand(FT, i_dim, batch_size))
                dist_kuma, dist_gumbel = sbc(x, k_out)
                is_all_true = is_all_true && dist_kuma isa BatchKumaraswamy
                is_all_true = is_all_true && dist_gumbel isa BatchGumbelBernoulliLogit
            end
            is_all_true
        end
    end

    @warn "`neural/sbc.jl` is not tested for online initialisation yet."
end


text:
####### RAYLEIGH FLOW ########

HeatFlux(h01::StagnationEnthalpy,h02::StagnationEnthalpy) = HeatFlux(h02-h01)

function HeatFlux(T01::StagnationTemperature,T02::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas)
  h01 = StagnationEnthalpy(T01;gas=gas)
  h02 = StagnationEnthalpy(T02;gas=gas)
  return HeatFlux(h01,h02)
end

function T0OverT0Star(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return TemperatureRatio(((γ+1)*M^2*(2+(γ-1)*M^2))/(1+γ*M^2)^2)
end

function TOverTStar(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return TemperatureRatio(((γ+1)^2*M^2)/(1+γ*M^2)^2)
end

function MachNumber(T0_over_T0star::TemperatureRatio,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
   value(T0_over_T0star) <= 1.0 || error("T0/T0* must be 1 or smaller")

    Msub = find_zero(x -> T0OverT0Star(MachNumber(x),RayleighFlow,gas=gas)-T0_over_T0star,(0.001,1),order=16)
    max_T0ratio = T0OverT0Star(MachNumber(1e10),RayleighFlow,gas=gas)

    if value(T0_over_T0star) < value(max_T0ratio)
        return MachNumber(Msub)
    else
        Msup = find_zero(x -> T0OverT0Star(MachNumber(x),RayleighFlow,gas=gas)-T0_over_T0star,(1,1e10),order=16)
        return MachNumber(Msub), MachNumber(Msup)
    end
end

function POverPStar(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return PressureRatio((γ+1)/(1+γ*M^2))
end

function VOverVStar(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return VelocityRatio(((γ+1)*M^2)/(1+γ*M^2))
end

function ρOverρStar(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  return DensityRatio(1/VOverVStar(M,RayleighFlow,gas=gas))
end

function P0OverP0Star(M::MachNumber,::Type{RayleighFlow};gas::PerfectGas=DefaultPerfectGas)
  γ = SpecificHeatRatio(gas)
  return PressureRatio((γ+1)/(1+γ*M^2)*((2+(γ-1)*M^2)/(γ+1))^(γ/(γ-1)))
end


text:
module gRPC

import ProtoBuf: call_method

using BitFlags
using CodecZlib
using Nghttp2
using ProtoBuf
using Sockets

export gRPCChannel, gRPCController, gRPCServer
export SerializeStream, DeserializeStream
export handle_request, call_method

"""
    gRPC status error codes.

    https://grpc.github.io/grpc/core/md_doc_statuscodes.html
"""
@enum(StatusCode::Int32,
    # Not an error; returned on success.
    STATUS_OK = 0,
    # The operation was cancelled, typically by the caller.
    STATUS_CANCELLED = 1,
    # Unknown error. For example, this error may be returned when a Status value received from
    # another address space belongs to an error space that is not known in this address space.
    STATUS_UNKNOWN = 2,
    # The client specified an invalid argument.
    STATUS_INVALID_ARGUMENT = 3,
    # The deadline expired before the operation could complete.
    STATUS_DEADLINE_EXCEEDED = 4,
    # Some requested entity (e.g., file or directory) was not found.
    STATUS_NOT_FOUND = 5,
    # The entity that a client attempted to create (e.g., file or directory) already exists.
    STATUS_ALREADY_EXISTS = 6,
    # The caller does not have permission to execute the specified operation.
    # PERMISSION_DENIED must not be used for rejections caused by exhausting some resource.
    STATUS_PERMISSION_DENIED = 7,
    # Some resource has been exhausted, perhaps a per-user quota, or perhaps the entire file system is out of space.
    STATUS_RESOURCE_EXHAUSTED = 8,
    # The operation was rejected because the system is not in a state required for the operation's execution.
    STATUS_FAILED_PRECONDITION = 9,
    # The operation was aborted, typically due to a concurrency issue such as a sequencer check failure or transaction abort.
    STATUS_ABORTED = 10,
    # The operation was attempted past the valid range.
    STATUS_OUT_OF_RANGE = 11,
    # The operation is not implemented or is not supported/enabled in this service.
    STATUS_UNIMPLEMENTED = 12,
    # Internal errors. This means that some invariants expected by the underlying system have been broken.
    # This error code is reserved for serious errors.
    STATUS_INTERNAL = 13,
    # The service is currently unavailable. This is most likely a transient condition,
    # which can be corrected by retrying with a backoff.
    # Note that it is not always safe to retry non-idempotent operations.
    STATUS_UNAVAILABLE = 14,
    # Unrecoverable data loss or corruption.
    STATUS_DATA_LOSS = 15,
    # The request does not have valid authentication credentials for the operation.
    STATUS_UNAUTHENTICATED = 16)

const STATUS_CODES = Dict(string(Int(i)) => i for i in instances(StatusCode))

"""
    gRPC error.
"""
struct gRPCError <: Exception
    status_code::StatusCode
    msg::AbstractString

    gRPCError(status_code::StatusCode, msg::AbstractString) = new(status_code, msg)
end

"""
    gRPC channel.
"""
mutable struct gRPCChannel <: ProtoRpcChannel
    session::Http2ClientSession

    gRPCChannel(session::Http2ClientSession) = new(session)
end

"""
    gRPC Server implementation.
"""
mutable struct gRPCServer
    proto_services::Dict{String,ProtoService}
    is_running::Bool

    gRPCServer(routes::Dict{String,ProtoService}) = new(routes, true)
end

"""
    gRPC Controller.
"""
struct gRPCController <: ProtoRpcController end

# https://github.com/grpc/grpc/blob/master/doc/PROTOCOL-HTTP2.md

"""
    gRPC Http2 responces.
"""
const DEFAULT_STATUS_200 = [
    ":status" => "200",
    "content-type" => "application/grpc",
    "grpc-encoding" => "gzip"]

const DEFAULT_GZIP_ENCRYPTION_STATUS_200 = [
    ":status" => "200",
    "content-type" => "application/grpc",
    "grpc-encoding" => "gzip"]

const DEFAULT_TRAILER = ["grpc-status" => "0"]

"""
    Serialize a collection of proto objects into a IO stream.
"""
mutable struct SerializeStream <: IO
    itr::Base.Iterators.Enumerate
    next::Any
    buffer::IOBuffer
    eof::Bool

    SerializeStream(itr) = new(itr, iterate(itr), IOBuffer(), false)
end

function internal_read(serialize_stream::SerializeStream)::Bool
    if serialize_stream.eof
        # No more elements available.
        return false
    end

    if !isnothing(serialize_stream.next)
        (i, state) = serialize_stream.next
        (_, element) = i

        # Write to the buffer.
        iob = serialize_object(element)

        current_position = position(serialize_stream.buffer)
        seekend(serialize_stream.buffer)
        write(serialize_stream.buffer, iob)
        seek(serialize_stream.buffer, current_position)

        # Get the next element from the iterator.
        serialize_stream.next = iterate(serialize_stream.itr, state)
    else
        serialize_stream.eof = true
    end

    return !serialize_stream.eof
end

function ensure_in_buffer(serialize_stream::SerializeStream, nb::Integer)
    should_read = true

    # Process the enumeration until there is no more elements available in the collection.
    while should_read
        if bytesavailable(serialize_stream.buffer) >= nb || serialize_stream.eof
            should_read = false
        end

        if should_read && internal_read(serialize_stream)
            continue
        end

        # The read stopped.
        break
    end
end

Base.eof(serialize_stream::SerializeStream)::Bool = serialize_stream.eof && eof(serialize_stream.buffer)

function Base.read(serialize_stream::SerializeStream, nb::Integer)::Vector{UInt8}
    ensure_in_buffer(serialize_stream, nb)

    return read(serialize_stream.buffer, nb)
end

function Base.read(serialize_stream::SerializeStream, ::Type{UInt8})
    ensure_in_buffer(serialize_stream, 1)

    return read(serialize_stream.buffer, UInt8)
end

"""
    Deserialize a collection of proto objects from the IO stream.
"""
struct DeserializeStream{T} <: AbstractChannel{T}
    io::IO

    function DeserializeStream{T}() where {T<:ProtoType}
        return new(devnull)
    end

    function DeserializeStream{T}(io::IO) where {T<:ProtoType}
        return new(io)
    end
end

function Base.Iterators.Enumerate{T}() where {T<:ProtoType}
    return DeserializeStream{T}()
end

function Base.iterate(stream::DeserializeStream{T}, s=nothing) where {T<:ProtoType}
    io = stream.io

    if eof(io)
        return nothing
    end

    local compressed::Bool
    try
        compressed = read(io, UInt8)
    catch e
        if isa(e, EOFError)
            return nothing
        else
            rethrow()
        end
    end

    data_length = ntoh(read(io, UInt32))

    instance::T = T()

    if data_length != 0
        io = IOBuffer(read(io, data_length))

        if compressed == 1
            io = GzipDecompressorStream(io)
        end

        readproto(io, instance)

        if compressed == 1
            finalize(io)
        end
    end

    return (instance, 1)
end

"""
    Deserialize a proto object instance from the io stream.
"""
function deserialize_object(io::IO, instance_type::Type{T}) where {T<:ProtoType}
    instance = instance_type()
    is_compressed = read(io, UInt8)
    data_length = ntoh(read(io, UInt32))

    if data_length != 0
        io = IOBuffer(read(io, data_length))

        if is_compressed == 1
            io = GzipDecompressorStream(io)
        end

        readproto(io, instance)

        if is_compressed == 1
            finalize(io)
        end
    end

    return instance
end

"""
    Deserialize a stream of proto objects from the io stream.
"""
function deserialize_object(io::IO, ::Type{AbstractChannel{T}}) where {T<:ProtoType}
    deserialize_stream = DeserializeStream{T}(io)
    return deserialize_stream
end

"""
    Serialize the instance of the proto object into the io buffer.
"""
function serialize_object(instance::ProtoType)::IO
    iob = IOBuffer()

    # Compression.
    #    write(iob, UInt8(1))

    #    iob_proto = IOBuffer()
    #    data_len = writeproto(iob_proto, instance)
    #    seek(iob_proto, 0)

    #    compressed = transcode(GzipCompressor, read(iob_proto))
    #    @show compressed

    #    @show compressed

    #    write(iob, hton(UInt32(length(compressed))))
    #    write(iob, compressed)

    #    @show iob

    # No compression
    write(iob, UInt8(0))
    # Placeholder for the serialized object length.
    write(iob, hton(UInt32(0)))
    data_length = writeproto(iob, instance)
    seek(iob, 1)
    write(iob, hton(UInt32(data_length)))

    #
    seek(iob, 0)
    return iob
end

function serialize_object(instances)::IO
    serialize_stream = SerializeStream(enumerate(instances))
    return serialize_stream
end

"""
    Process server request.
"""
function handle_request(http2_server_session::Http2ServerSession, controller::gRPCController, server::gRPCServer)
    request_stream::Http2Stream = recv(http2_server_session)

    headers = request_stream.headers

    method = headers[":method"]
    path = headers[":path"]
    path_components = split(path, "/"; keepempty=false)

    if length(path_components) != 2
        # Missing or invalid path in request's header.
        return nothing
    end

    service_name, method_name = path_components

    proto_service::ProtoService = server.proto_services[service_name]

    method = find_method(proto_service, method_name)

    request_type = get_request_type(proto_service, method)

    request_argument = deserialize_object(request_stream, request_type)

    response = call_method(proto_service, method, controller, request_argument)

    io = serialize_object(response)

    return submit_response(request_stream, io, DEFAULT_STATUS_200, DEFAULT_TRAILER)
end

"""
    Client request.
"""
function call_method(
    channel::ProtoRpcChannel,
    service::ServiceDescriptor,
    method::MethodDescriptor,
    controller::ProtoRpcController,
    request)
    path = "/" * service.name * "/" * method.name
    headers = [
        ":method" => "POST",
        ":path" => path,
        ":authority" => "localhost:5000",
        ":scheme" => "http",
        "user-agent" => "grpc-julia",
        "accept-encoding" => "identity,gzip",
        "content-type" => "application/grpc",
        "grpc-accept-encoding" => "identity,deflate,gzip",
        "te" => "trailers"]

    iob = gRPC.serialize_object(request)

    response_stream = submit_request(channel.session, iob, headers)

    response_status_code::StatusCode = STATUS_OK

    # Check grps-status if it is available.
    if haskey(response_stream.headers, "grpc-status")
        grpc_status = response_stream.headers["grpc-status"]

        # Get the response code.
        if haskey(STATUS_CODES, grpc_status)
            response_status_code = STATUS_CODES[grpc_status]
        else
            # Server returned an invalid responose code.
            response_status_code = STATUS_UNKNOWN
        end
    end

    if response_status_code != STATUS_OK
        local response_message::String

        if haskey(response_stream.headers, "grpc-message")
            response_message = response_stream.headers["grpc-message"]
        else
            # Server did not include a response in the message.
            response_message = "unknown"
        end

        throw(gRPCError(response_status_code, response_message))
    end

    response_type = get_response_type(method)

    instance = deserialize_object(response_stream, response_type)

    return instance
end

end # module gRPC


text:
segy_file_1   = string(path, "/SampleFiles/SEGY/03.334.12.09.00.0362.1")
segy_file_2   = string(path, "/SampleFiles/SEGY/test_rev1.segy")
segy_file_3   = string(path, "/SampleFiles/SEGY/04.322.18.06.46.92C3.1.coords1.segy")
segy_file_4   = string(path, "/SampleFiles/SEGY/04.322.18.06.46.92C3.1.coords2.segy")
segy_be_file  = string(path, "/SampleFiles/SEGY/02.104.00.00.08.7107.2")
segy_fpat     = string(path, "/SampleFiles/SEGY/02.104.00.00.08.7107.2")
segy_nodal    = string(path, "/SampleFiles/SEGY/FORGE_78-32_iDASv3-P11_UTC190428135038.sgy")

printstyled("  SEG Y\n", color=:light_green)
printstyled("    Helper functions\n", color=:light_green)
printstyled("      IBM-float\n", color=:light_green)

#=
  These bit representations are from Wikipedia's page on IBM-float,
  https://en.wikipedia.org/wiki/IBM_hexadecimal_floating_point
=#
y = [-118.625, +7.2370051e75, +5.397605e-79, 0.0, 0.1]
for (i, str) in enumerate( ["11000010011101101010000000000000",
                            "01111111111111111111111111111111",
                            "00000000000100000000000000000000",
                            "00000000000000000000000000000000",
                            "01000000000110011001100110011010"] )
  u = hton(parse(UInt32, str, base=2))
  @test isapprox(ibmfloat(u), y[i], rtol=10*eps(Float32))
end

printstyled("      auto_coords\n", color=:light_green)
sc = Int16[1,1]
lat0 = -39.4209
lon0 = -71.94061
ele = 2856.3
xy = [round(Int32, lon0*1.0e4), round(Int32, lat0*1.0e4)]
y,x = auto_coords(xy, sc)
@test latlon2xy(y,x) == xy

xy2 = latlon2xy(lat0, lon0)
auto_coords(xy2, sc)
@test isapprox([auto_coords(xy2, sc)...], [lat0, lon0], rtol=eps(Float32))

S1 = read_data(segy_file_3)
S2 = read_data(segy_file_4)
@test abs(S1.loc[1].lat-S2.loc[1].lat) < 1.0e-4
@test abs(S1.loc[1].lon-S2.loc[1].lon) < 1.0e-4
@test abs(S1.loc[1].el-S2.loc[1].el) < 1.0e-4
@test S1.x[1] == S2.x[1]

printstyled("      trid\n", color=:light_green)
@test trid(Int16(2), S1.fs[1], 1.0) == "EHZ"
@test trid(Int16(2), S1.fs[1], 1/30) == "HHZ"

printstyled("    SEG Y rev 1\n", color=:light_green)
if has_restricted
  SEG = verified_read_data("segy", segy_file_2)
  redirect_stdout(out) do
    segyhdr(segy_file_2)
  end

  printstyled("      data integrity\n", color=:light_green)
  x = collect(1:1:8)
  for i in [1,2,3,5,8]
    fname = string(path, "/SampleFiles/SEGY/test_rev1_", i, ".segy")
    S = read_data("segy", fname)
    @test isapprox(S.x[1], x)
  end
  @test_throws ErrorException read_data("segy", string(path, "/SampleFiles/SEGY/test_rev1_4.segy"))

  printstyled("      header accuracy\n", color=:light_green)
  S = read_data("segy", segy_file_2)
  @test S.n == 36
  @test length(unique([i[1:6] for i in S.id])) == 12
  for i in 1:S.n
    @test S.fs[i] ≈ 2000.0
    @test S.gain[i] ≈ 320.0
    @test S.units[i] == "m/s"
    @test length(S.x[i]) == 32000
    @test size(S.t[i], 1) == 2
  end
  @test length(findall([endswith(i, "GHZ") for i in S.id])) == 12
  @test length(findall([endswith(i, "GH1") for i in S.id])) == 12
  @test length(findall([endswith(i, "GH2") for i in S.id])) == 12
else
  printstyled("    Skipped SEG Y rev 1\n", color=:red)
end

printstyled("    PASSCAL/NMT SEG Y\n", color=:light_green)
SEG = verified_read_data("passcal", segy_file_1, full=true)

printstyled("      header integrity\n", color=:light_green)

SEG = verified_read_data("passcal", segy_file_1, full=true)
@test SEG.misc[1]["gain_const"] == 32
@test SEG.gain[1] ≈ SEG.misc[1]["scale_fac"]
@test isapprox(1.0/SEG.gain[1], 4.47021e-07/SEG.misc[1]["gain_const"], atol=eps(Float32))
@test SEG.fs[1] == 100.0 == 1.0e6 / SEG.misc[1]["delta"]
@test lastindex(SEG.x[1]) == 247698
@test SEG.misc[1]["trace_seq_line"] == 3
@test SEG.misc[1]["trace_seq_file"] == 3
@test SEG.misc[1]["rec_no"] == 1
@test SEG.misc[1]["channel_no"] == 2
@test SEG.misc[1]["trace_id_code"] == 3
@test SEG.misc[1]["h_units_code"] == 2
@test SEG.misc[1]["nx"] == 32767
@test SEG.misc[1]["samp_rate"] == 10000
@test SEG.misc[1]["gain_type"] == 1
@test SEG.misc[1]["year"] == 2003
@test SEG.misc[1]["day"] ==  334
@test SEG.misc[1]["hour"] == 12
@test SEG.misc[1]["minute"] == 9
@test SEG.misc[1]["second"] == 0
@test SEG.misc[1]["ms"] == 5
@test SEG.misc[1]["time_code"] == 2
@test SEG.misc[1]["trigyear"] == 2003
@test SEG.misc[1]["trigday"] == 334
@test SEG.misc[1]["trighour"] == 12
@test SEG.misc[1]["trigminute"] == 9
@test SEG.misc[1]["trigsecond"] == 0
@test SEG.misc[1]["trigms"] == 5
@test SEG.misc[1]["data_form"] == 1
@test SEG.misc[1]["inst_no"] == 0x016a # 0362
@test strip(SEG.misc[1]["sensor_serial"]) == "UNKNOWN"
@test strip(SEG.misc[1]["station_name"]) == "362"

# Location
printstyled("      sensor position\n", color=:light_green)
h_sc = Float64(get(SEG.misc[1], "h_sc", 1.0))
h_sc = abs(h_sc)^(h_sc < 0.0 ? -1 : 1)
z_sc = Float64(get(SEG.misc[1], "z_sc", 1.0))
z_sc = abs(z_sc)^(z_sc < 0.0 ? -1 : 1)
x = get(SEG.misc[1], "rec_x", 0.0)
y = get(SEG.misc[1], "rec_y", 0.0)
z = get(SEG.misc[1], "rec_ele", 0.0)
# @test SEG.loc[1].lat == y*h_sc == 45.2896      # 45.2896 in wash.sta
# @test SEG.loc[1].lon == x*h_sc == -121.7915    # 121.791496 in wash.sta
# @test SEG.loc[1].el == z*z_sc == 1541.0       # 1541.0 in wash.sta

printstyled("      data integrity\n", color=:light_green)
@test Float64(SEG.misc[1]["max"]) == maximum(SEG.x[1]) == 396817
@test Float64(SEG.misc[1]["min"]) == minimum(SEG.x[1]) == -416512
# @test ≈(SEG.x[1][1:10], [47.0, 46.0, 45.0, 44.0, 51.0, 52.0, 57.0, 59.0, 40.0, 34.0])
@test length(SEG.x[1]) == SEG.misc[1]["num_samps"] == 247698

redirect_stdout(out) do
  segyhdr(segy_file_1, passcal=true)
end

printstyled("      big-endian support\n", color=:light_green)

SEG = verified_read_data("passcal", segy_be_file, full=true, swap=true)
@test SEG.n == 1
@test SEG.id ==  ["...spn"] # lol, WHY BOB
@test isapprox(1.0/SEG.gain[1], 5.92346875e-8, atol=eps(Float32))
@test SEG.misc[1]["trigyear"] == SEG.misc[1]["year"] == 2002

printstyled("    wildcard support\n", color=:light_green)

SEG = verified_read_data("passcal", segy_fpat, full=true, swap=true)
@test SEG.n == 1
@test Float64(SEG.misc[1]["max"]) == maximum(SEG.x[1]) == 49295.0
@test Float64(SEG.misc[1]["min"]) == minimum(SEG.x[1]) == -54454.0
@test ≈(SEG.x[1][1:5], [-615.0, -3994.0, -4647.0, -895.0, 190.0])
@test length(SEG.x[1]) == SEG.misc[1]["num_samps"] == 180027

printstyled("    keyword ll\n", color=:light_green)
for ll in (0x01, 0x04, 0x07)
  S = read_data("segy", segy_nodal, ll=ll)
  @test S.n == 33
  @test first(S.id) == ".0.01.YYY"
  @test last(S.id) == ".0.0X.YYY"
end

for ll in (0x02, 0x05, 0x06)
  S = read_data("segy", segy_nodal, ll=ll)
  @test S.n == 1
  @test S.id[1] == ".0.00.YYY"
end

S = read_data("segy", segy_nodal)
@test S.n == 1
@test S.id[1] == ".0..YYY"
@test_throws InexactError read_data("segy", segy_nodal, ll=0x03)