text:
using StatsBase
using Base.Test

X = randn(3, 8)

Z1 = X .- mean(X, 1)
Z2 = X .- mean(X, 2)

w1 = rand(3)
w2 = rand(8)

wv1 = weights(w1)
wv2 = weights(w2)

Z1w = X .- mean(X, wv1, 1)
Z2w = X .- mean(X, wv2, 2)

## reference results

S1 = Z1'Z1
S2 = Z2 * Z2'

Sz1 = X'X
Sz2 = X * X'

S1w = Z1w' * diagm(w1) * Z1w
S2w = Z2w * diagm(w2) * Z2w'

Sz1w = X' * diagm(w1) * X
Sz2w = X * diagm(w2) * X'

## scattermat

@test_approx_eq scattermat(X) S1
@test_approx_eq scattermat(X; vardim=2) S2

@test_approx_eq scattermat(X; mean=0) Sz1
@test_approx_eq scattermat(X; mean=0, vardim=2) Sz2

@test_approx_eq scattermat(X; mean=mean(X,1)) S1
@test_approx_eq scattermat(X; mean=mean(X,2), vardim=2) S2

@test_approx_eq scattermat(X; mean=zeros(1,8)) Sz1
@test_approx_eq scattermat(X; mean=zeros(3), vardim=2) Sz2

## weighted scatter mat

@test_approx_eq scattermat(X, wv1) S1w
@test_approx_eq scattermat(X, wv2; vardim=2) S2w

@test_approx_eq scattermat(X, wv1; mean=0) Sz1w
@test_approx_eq scattermat(X, wv2; mean=0, vardim=2) Sz2w

@test_approx_eq scattermat(X, wv1; mean=mean(X, wv1, 1)) S1w
@test_approx_eq scattermat(X, wv2; mean=mean(X, wv2, 2), vardim=2) S2w

@test_approx_eq scattermat(X, wv1; mean=zeros(1,8)) Sz1w
@test_approx_eq scattermat(X, wv2; mean=zeros(3), vardim=2) Sz2w

# weighted covariance

@test_approx_eq cov(X, wv1) S1w ./ sum(wv1)
@test_approx_eq cov(X, wv2; vardim=2) S2w ./ sum(wv2)

@test_approx_eq cov(X, wv1; mean=0) Sz1w ./ sum(wv1)
@test_approx_eq cov(X, wv2; mean=0, vardim=2) Sz2w ./ sum(wv2)

@test_approx_eq cov(X, wv1; mean=mean(X, wv1, 1)) S1w ./ sum(wv1)
@test_approx_eq cov(X, wv2; mean=mean(X, wv2, 2), vardim=2) S2w ./ sum(wv2)

@test_approx_eq cov(X, wv1; mean=zeros(1,8)) Sz1w ./ sum(wv1)
@test_approx_eq cov(X, wv2; mean=zeros(3), vardim=2) Sz2w ./ sum(wv2)

# mean_and_cov

(m, C) = mean_and_cov(X; vardim=1)
@test m == mean(X, 1)
@test C == cov(X; vardim=1)

(m, C) = mean_and_cov(X; vardim=2)
@test m == mean(X, 2)
@test C == cov(X; vardim=2)

(m, C) = mean_and_cov(X, wv1; vardim=1)
@test m == mean(X, wv1, 1)
@test C == cov(X, wv1; vardim=1)

(m, C) = mean_and_cov(X, wv2; vardim=2)
@test m == mean(X, wv2, 2)
@test C == cov(X, wv2; vardim=2)


text:
# Copyright 2017, Andrey Prokopenko, ORNL

function partition_manual(A, nxyz, mxyz)
    dim = size(nxyz, 2)
    @assert(dim <= 3)

    if dim == 1
        nxyz = [nxyz[1] 0 0]
        mxyz = [mxyz[1] 1 1]
    elseif dim == 2
        nxyz = [nxyz[1] nxyz[2] 0]
        mxyz = [mxyz[1] mxyz[2] 1]
    end
    m = prod(mxyz)

    row_gids = Array{Any}(1,m)
    for K = 1:mxyz[3]
        for J = 1:mxyz[2]
            for I = 1:mxyz[1]
                IJK = (K-1)*mxyz[2]*mxyz[1] + (J-1)*mxyz[1] + I

                (min_x, max_x) = sub1D_bnd(nxyz[1], mxyz[1], I)

                if nxyz[2] != 0
                    (min_y, max_y) = sub1D_bnd(nxyz[2], mxyz[2], J)
                else
                    min_y = 1
                    max_y = 1
                end

                if nxyz[3] != 0
                    (min_z, max_z) = sub1D_bnd(nxyz[3], mxyz[3], K)
                else
                    min_z = 1
                    max_z = 1
                end

                # Calculate inner part othe subdomain
                rows = zeros(Int64, 1, (max_x-min_x+1)*(max_y-min_y+1)*(max_z-min_z+1))
                ind  = 1
                for k = min_z:max_z
                    for j = min_y:max_y
                        for i = min_x:max_x
                            rows(ind) = (k-1)*nxyz[2]*nxyz[1] +
                                        (j-1)*nxyz[1] + i
                            ind = ind+1
                        end
                    end
                end

                row_gids[IJK] = rows
            end
        end
    end

    return row_gids
end

# function [row_gids] = partition_metis(A, m)
# % Requires metismex MATLAB interface to METIS
# % Download here: https://github.com/dgleich/metismex
# map = metismex('PartGraphRecursive', A, m)
#
# row_gids = cell(1,m)
# for I = 1:m
    # row_gids{I} = find(map == (I-1))
# end
# end
#
# function [row_gids] = partition_schwarz(A, row_gids0, schwarz)
#
# % Add overlapped parts
# % The nodes are sorted with respect to the overlap level:
# % first go inner subdomain nodes, then overlap one, then overlap
# % two, and so on
# rows = row_gids0
# for s = 1:schwarz
    # % Find nodes connected to the previous overlap level,
    # % starting from 0
    # [~,cols] = find(A(rows,:))
    # if size(cols,1) == 1
        # cols = cols'
    # end
    # % Update nodes to include those.
    # rows = cat(2, rows, setdiff(cols', rows))
# end
#
# row_gids = rows
# end

function construct_cols(A, row_gids)
    # col_gids is similar to column map
    # It contains a list of columns with nonzero elements in row_gids rows, and
    # is ordered so that the part corresponding to row_gids goes first.
    (_, cols) = findn(A[row_gids,:])
    if size(cols,1) == 1
        cols = cols'
    end
    col_gids = cat(2, row_gids, setdiff(unique(cols'), row_gids))

    return col_gids
end

# function [neigh] = get_neigh_list(sub_index)
# global col_gids0
# global gid2proc
# procs = gid2proc(col_gids0{sub_index})
# neigh = setdiff(unique(procs), sub_index)
# end

function sub1D_bnd(n, m, i)
    @assert(i > 0 && i <= m, "Incorrect i value")

    min_width = trunc(Int64, n/m)
    num_wides = n - m*min_width

    if i <= num_wides
        min = (i-1) *   (min_width+1)
        max = min +      min_width+1
    else
        min = num_wides*(min_width+1) + (i-1 - num_wides)*min_width
        max = min +      min_width
    end
    min = min+1

    return min, max
end

# function [] = plot_stencil(A, coords, color)
# n = size(A,1)
#
# % 1. Separate the matrix into two parts: one with symmetric pattern, and a remainder
# % Get the square part of A (A may be short and wide rectangular matrix)
# S = A(1:n,1:n)
# Z = S
# % Create a pattern matrix
# Z(Z ~= 0) = 1
# % Create the symmetric subpattern
# [r,c] = find((Z+Z') == 2)
# % Create final pattern
# Z = sparse(r,c,ones(1,size(r,1)),n,n)
# S = S .* Z
#
# line_width = 3
# marker = 'o'
#
# % 2. Plot the symmetric part
# [rows,cols] = find(S)
# for k = 1:size(rows)
    # i = rows(k)
    # j = cols(k)
    # if i <= j % need to plot only once
        # plot3([coords(i,1), coords(j,1)], ...
              # [coords(i,2), coords(j,2)], ...
              # [coords(i,3), coords(j,3)], ...
              # 'LineWidth',          line_width, ...
              # 'Marker',             marker, ...
              # 'MarkerFaceColor',    color, ...
              # 'Color',              color)
    # end
# end
#
# % 3. Plot the nonsymmetric part
# [rows,cols] = find(A-S)
# scaling = 0.73
# u = coords(cols,1)-coords(rows,1)
# v = coords(cols,2)-coords(rows,2)
# w = coords(cols,3)-coords(rows,3)
# quiver3(coords(rows,1), coords(rows,2), coords(rows,3), ...
        # scaling*u, scaling*v, scaling*w, ...
        # 'LineWidth',        line_width, ...
        # 'Marker',           marker, ...
        # 'MarkerFaceColor',  color, ...
        # 'AutoScale',        'off', ...
        # 'Color',            color)
#
# axis square
# % axis off
# view(0,90)
#
# set(gcf, 'PaperPositionMode',   'auto')
# set(gcf, 'PaperPosition',       [0 0 100 100])
#
# end
#
# function [alessb, blessa] = compare(ind1, ind2)
# cmp = [ind1 <= ind2]
# alessb = ~sum(find(cmp == 0))
#
# cmp = [ind2 <= ind1]
# blessa = ~sum(find(cmp == 0))
#
# end

ext: jl

max_stars_repo_name: aprokop/ca-preconditioners

max_issues_repo_name: aprokop/ca-preconditioners

max_forks_repo_name: aprokop/ca-preconditioners


text:
type DiscreteMMSE
    rho # snr

    W
    bias

    DiscreteMMSE() = new()
end

function DiscreteMMSE(snr,get_cov;nSamples=1)
    est = DiscreteMMSE()
    est.rho = 10^(0.1*snr)

    est.W = Array{Any}(nSamples)
    est.bias = Array{Any}(nSamples)

    for i in 1:nSamples
        C = get_cov() # get random sample from cov. prior
        nAntennas   = size(C,1)
        est.W[i]    = C/(C + 1/est.rho * eye(nAntennas)) # MMSE filter at sample i
        est.bias[i] = real(logdet(eye(nAntennas) - est.W[i])) # bias term at sample i
    end

    return est
end

"""
    estimate(est::DiscreteMMSE, y)

For each batch `b`, estimate channels `h[:,t,b]` from
observations `y[:,1], ..., y[:,T,b]` by formula

h[:,t,b] = W_est[b] * y[:,t,b]

with

W_est[b] = sum_i{ weights_i[b] W_i }

and

weights_i[b] = softmax( rho*sum_t{ y[:,t,b]' W[i] * y[:,t,b] } + nCoherence*bias[i] )
"""
function estimate(est::DiscreteMMSE, y)
    nAntennas  = size(y,1)
    nCoherence = size(y,2)
    nBatches   = size(y,3);
    nSamples   = length(est.W);

    # calculate weights for all batches
    exps = zeros(nSamples,nBatches)
    for b in 1:nBatches, t in 1:nCoherence, i in 1:nSamples
        exps[i,b] += real(y[:,t,b]'*est.W[i]*y[:,t,b])[1] * est.rho + est.bias[i]
    end
    exps = exps .- maximum(exps,1)

    weights = exp.(exps)
    weights = weights./sum(weights,1)

    # Calculate estimated MMSE filters and channel estimates
    hest = zeros(nAntennas, nCoherence, nBatches) + 1im*zeros(nAntennas, nCoherence, nBatches)
    for b = 1:nBatches
        W_est = zeros(nAntennas, nAntennas)
        for i=1:nSamples
            W_est += est.W[i]*weights[i,b];
        end
        hest[:,:,b] = W_est*y[:,:,b]
    end

    return hest
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
using KernelAbstractions
using Oceananigans.Utils: work_layout
using Oceananigans.Fields
using Oceananigans.Diagnostics: dims_to_result_size
import Oceananigans.Diagnostics: run_diagnostic!

import Oceananigans.Fields: location, total_size, compute!

"""
    Computation{T, R, O, G}

Represents an operation performed over the elements of a field.
"""
struct Computation{T, R, O, G}
      operation :: O
         result :: R
           grid :: G
    return_type :: T
end

"""
    Computation(operation, result; return_type=Array)

Returns a `Computation` representing an `operation` performed over the elements of
`operation.grid` and stored in `result`. `return_type` specifies the output type when the
`Computation` instances is called as a function.
"""
Computation(operation, result; return_type=Array) =
    Computation(operation, result, operation.grid, return_type)

# Utilities for limited field-like behavior
architecture(comp::Computation) = architecture(comp.result)
Base.parent(comp::Computation) = comp # this enables taking a "horizontal average" of a computation
location(comp::Computation) = location(comp.operation)
total_size(comp::Computation) = total_size(comp.result)

"""
    compute!(computation::Computation)

Perform a `computation`. The result is stored in `computation.result`.
"""
function compute!(computation::Computation)
    arch = architecture(computation.result)
    result_data = data(computation.result)
    loc = location(computation.result)

    workgroup, worksize = work_layout(computation.grid, :xyz, include_right_boundaries=true, location=loc)

    compute_kernel! = _compute!(device(arch), workgroup, worksize)

    event = compute_kernel!(result_data, computation.grid, computation.operation; dependencies=Event(device(arch)))

    wait(device(arch), event)

    return nothing
end

"""Compute an `operation` over `grid` and store in `result`."""
@kernel function _compute!(result, grid, operation)
    i, j, k = @index(Global, NTuple)
    @inbounds result[i, j, k] = operation[i, j, k]
end

"""
    (computation::Computation)(args...)

Performs the `compute(computation)` and returns the result if `isnothing(return_type)`,
or the result after being converted to `return_type`.
"""
function (computation::Computation)(args...)
    compute!(computation)
    return computation.return_type(interior(computation.result))
end

function (computation::Computation{<:Nothing})(args...)
    compute!(computation)
    return computation.result
end

#####
##### Functionality for using computations with Average
#####

"""
    Average(op::AbstractOperation, result; dims, kwargs...)

Returns the representation of an `Average` over the operation `op`, using `result` as
a temporary array to store the result of `operation` computed on `op.grid`.
"""
function Average(op::AbstractOperation, result; dims, kwargs...)
    computation = Computation(op, result)
    return Average(computation; dims=dims, kwargs...)
end

"""
    Average(op::AbstractOperation, model; dims, kwargs...)

Returns the representation of an `Average` over the operation `op`, using
`model.pressures.pHY′` as a temporary array to store the result of `operation` computed on
`op.grid`.
"""
Average(op::AbstractOperation, model::AbstractModel; dims, kwargs...) =
    Average(op, model.pressures.pHY′; dims=dims, kwargs...)

"""
    Average(field::Field; dims, iteration_interval=nothing, time_interval=nothing, return_type=Array)

Construct an `Average` of `field` along the dimensions specified by the tuple `dims`.

After the average is computed it will be stored in the `result` property.

The `Average` can be used as a callable object that computes and returns the average.

An `iteration_interval` or `time_interval` (or both) can be passed to indicate how often to
run this diagnostic if it is part of `simulation.diagnostics`. `iteration_interval` is a
number of iterations while `time_interval` is a time interval in units of `model.clock.time`.

A `return_type` can be used to specify the type returned when the `Average` is
used as a callable object. The default `return_type=Array` is useful when running a GPU
model and you want to save the output to disk by passing it to an output writer.
"""
function Average(field::Computation; dims, iteration_interval=nothing, time_interval=nothing, return_type=Array,
                 with_halos=true)

    dims isa Union{Int, Tuple} || error("Average dims must be an integer or tuple!")
    dims isa Int && (dims = tuple(dims))

    length(dims) == 0 && error("dims is empty! Must average over at least one dimension.")
    length(dims) > 3  && error("Models are 3-dimensional. Cannot average over 4+ dimensions.")
    all(1 <= d <= 3 for d in dims) || error("Dimensions must be one of 1, 2, 3.")

    arch = architecture(field)
    result_size = dims_to_result_size(field, dims, field.grid)
    result = zeros(arch, field.grid, result_size...)

    return Average(field, dims, result, iteration_interval, time_interval, 0.0,
                   return_type, with_halos)
end

"""Compute the average of a computation."""
function run_diagnostic!(avg::Average{<:Computation}, model)
    compute!(avg.field)
    zero_halo_regions!(parent(avg.field.result), (Cell, Cell, Cell), model.grid)
    sum!(avg.result, parent(avg.field.result))
    normalize_sum!(avg)
    return nothing
end


text:
s = """
throw(MethodError("blah"))
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
MethodError("blah")
"""
msgs = lintstr(s)
@test msgs[1].code == :W448
@test occursin("MethodError is an Exception but it is not enclosed in a throw", msgs[1].message)


text:
## gzip file io ##
module GZip
using Libdl
using Base.Libc
import Base: show, fd, close, flush, truncate, seek,
             seekend, skip, position, eof, read,
             readline, write, unsafe_write, peek

export
  GZipStream,
  show,

# io functions
# open,  ## not exported; use as GZip.open(...)
  gzopen,
  gzdopen,
  fd,
  close,
  flush,
  truncate,
  seek,
  skip,
  position,
  eof,
  read,
  readline,
  write,
  unsafe_write,
  peek,

# lower-level io functions
  gzgetc,
  gzungetc,
  gzgets,
  gzputc,
  gzwrite,
  gzread,
  gzbuffer,

# File offset
  ZFileOffset,

# GZError, ZError, related constants (zlib_h.jl)
  GZError,
  ZError,
  Z_OK,
  Z_STREAM_END,
  Z_NEED_DICT,
  Z_ERRNO,
  Z_STREAM_ERROR,
  Z_DATA_ERROR,
  Z_MEM_ERROR,
  Z_BUF_ERROR,
  Z_VERSION_ERROR,

# Compression constants (zlib_h.jl)
  Z_NO_COMPRESSION,
  Z_BEST_SPEED,
  Z_BEST_COMPRESSION,
  Z_DEFAULT_COMPRESSION,

# Compression strategy (zlib_h.jl)
  Z_FILTERED,
  Z_HUFFMAN_ONLY,
  Z_RLE,
  Z_FIXED,
  Z_DEFAULT_STRATEGY,

# Default buffer sizes
  Z_DEFAULT_BUFSIZE,
  Z_BIG_BUFSIZE

include("zlib_h.jl")

const GZLIB_VERSION = unsafe_string(ccall((:zlibVersion, GZip._zlib), Ptr{UInt8}, ()))

# Expected line length for strings
const GZ_LINE_BUFSIZE = 256

# Wrapper around gzFile
mutable struct GZipStream <: IO
    name::AbstractString
    gz_file::Ptr{Cvoid}
    buf_size::Int

    _closed::Bool

    function GZipStream(name::AbstractString, gz_file::Ptr{Cvoid}, buf_size::Int)
        x = new(name, gz_file, buf_size, false)
        finalizer(close, x)
        x
    end
end
GZipStream(name::AbstractString, gz_file::Ptr{Cvoid}) = GZipStream(name, gz_file, Z_DEFAULT_BUFSIZE)

# gzerror
function gzerror(err::Integer, s::GZipStream)
    e = Int32[err]
    if !s._closed
        msg_p = ccall((:gzerror, _zlib), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Int32}),
                      s.gz_file, e)
        msg = (msg_p == C_NULL ? "" : unsafe_string(msg_p))
    else
        msg = "(GZipStream closed)"
    end
    (e[1], msg)
end
gzerror(s::GZipStream) = gzerror(0, s)

mutable struct GZError <: Exception
    err::Int32
    err_str::AbstractString

    GZError(e::Integer, str::AbstractString) = new(Int32(e), str)
    GZError(e::Integer, s::GZipStream) = (a = gzerror(e, s); new(a[1], a[2]))
    GZError(s::GZipStream) = (a = gzerror(s); new(a[1], a[2]))
end

# show
show(io::IO, s::GZipStream) = print(io, "GZipStream(", s.name, ")")

macro test_eof_gzerr(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val))
            if eof($(esc(s)))  throw(EOFError())  else  throw(GZError($(esc(s))))  end
        end
        ret
    end
end

macro test_eof_gzerr2(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val)) && !eof($(esc(s))) throw(GZError($(esc(s)))) end
        ret
    end
end

macro test_gzerror(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val)) throw(GZError(ret, $(esc(s)))) end
        ret
    end
end

macro test_gzerror0(s, cc)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret <= 0 throw(GZError(ret, $(esc(s)))) end
        ret
    end
end

macro test_z_ok(cc)
    quote
        ret = $(esc(cc))
        if (ret != Z_OK) throw(ZError(ret)) end
        ret
    end
end

# Easy access to gz reading/writing functions (Internal)
gzgetc(s::GZipStream) =
    @test_eof_gzerr(s, ccall((:gzgetc, _zlib), Int32, (Ptr{Cvoid},), s.gz_file), -1)

gzgetc_raw(s::GZipStream) = ccall((:gzgetc, _zlib), Int32, (Ptr{Cvoid},), s.gz_file)

gzungetc(c::Integer, s::GZipStream) =
    @test_eof_gzerr(s, ccall((:gzungetc, _zlib), Int32, (Int32, Ptr{Cvoid}), c, s.gz_file), -1)

gzgets(s::GZipStream, a::Array{UInt8}) =
    @test_eof_gzerr2(s,
                     ccall((:gzgets, _zlib), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}, Int32),
                           s.gz_file, a, Int32(length(a))),
                     C_NULL)

gzgets(s::GZipStream, p::Ptr{UInt8}, len::Integer) =
    @test_eof_gzerr2(s,
                     ccall((:gzgets, _zlib), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}, Int32),
                           s.gz_file, p, Int32(len)),
                     C_NULL)

gzputc(s::GZipStream, c::Integer) =
    @test_gzerror(s,
                  ccall((:gzputc, _zlib), Int32, (Ptr{Cvoid}, Int32),
                        s.gz_file, Int32(c)),
                  -1)

gzwrite(s::GZipStream, p::Ptr, len::Integer) =
    len == 0 ? Int32(0) :
               @test_gzerror0(s, ccall((:gzwrite, _zlib), Int32, (Ptr{Cvoid}, Ptr{Cvoid}, UInt32),
                                       s.gz_file, p, len))

gzread(s::GZipStream, p::Ptr, len::Integer) =
    @test_gzerror(s,
                  ccall((:gzread, _zlib), Int32, (Ptr{Cvoid}, Ptr{Cvoid}, UInt32),
                        s.gz_file, p, len),
                  -1)

let _zlib_h = Libdl.dlopen(_zlib)
    global gzbuffer, _gzopen, _gzseek, _gztell, _gzrewind, _gzdirect, _gzoffset

    # Doesn't exist in zlib 1.2.3 or earlier
    if Libdl.dlsym_e(_zlib_h, :gzbuffer) != C_NULL
        gzbuffer(gz_file::Ptr, gz_buf_size::Integer) =
           ccall((:gzbuffer, _zlib), Int32, (Ptr{Cvoid}, UInt32), gz_file, gz_buf_size)
    else
        gzbuffer(gz_file::Ptr, gz_buf_size::Integer) = Int32(-1)
    end

    #####

    # Use 64-bit functions if available

    if Libdl.dlsym_e(_zlib_h, :gzopen64) != C_NULL && (z_off_t_sz == 8 || !Sys.iswindows())
        const _gzopen = :gzopen64
        const _gzseek = :gzseek64
        const _gztell = :gztell64
        const _gzoffset = :gzoffset64
    else
        const _gzopen = :gzopen
        const _gzseek = :gzseek
        const _gztell = :gztell
        const _gzoffset = :gzoffset
    end
    const _gzrewind = :gzrewind
    const _gzdirect = :gzdirect
end

function gzopen(fname::AbstractString, gzmode::AbstractString, gz_buf_size::Integer)
    # gzmode can contain extra characters specifying
    # * compression level (0-9)
    # * strategy ('f' => filtered data, 'h' -> Huffman-only compression,
    #             'R' -> run-length encoding, 'F' -> fixed code compression)
    #
    # '+' is also not allowed

    # For windows, force binary mode; doesn't hurt on unix
    if !('b' in gzmode)
        gzmode *= "b"
    end

    gz_file = ccall((_gzopen, _zlib), Ptr{Cvoid}, (Ptr{UInt8}, Ptr{UInt8}), fname, gzmode)
    if gz_file == C_NULL
        throw(GZError(-1, "gzopen failed"))
    end
    if gz_buf_size != Z_DEFAULT_BUFSIZE
        if gzbuffer(gz_file, gz_buf_size) == -1
            # Generally a non-fatal error, although it shouldn't happen here
            gz_buf_size = Z_DEFAULT_BUFSIZE
        end
    end
    s = GZipStream(fname, gz_file, gz_buf_size)
    peek(s) # Set EOF-bit for empty files
    return s
end
gzopen(fname::AbstractString, gzmode::AbstractString) = gzopen(fname, gzmode, Z_DEFAULT_BUFSIZE)
gzopen(fname::AbstractString) = gzopen(fname, "rb", Z_DEFAULT_BUFSIZE)
open(args...) = gzopen(args...)

function gzopen(f::Function, args...)
    io = gzopen(args...)
    try f(io)
    finally close(io)
    end
end

function gzdopen(name::AbstractString, fd::Integer, gzmode::AbstractString, gz_buf_size::Integer)
    if !('b' in gzmode)
        gzmode *= "b"
    end

    # Duplicate the file descriptor, since we have no way to tell gzclose()
    # not to close the original fd
    dup_fd = Libc.dup(Libc.RawFD(fd))

    gz_file = ccall((:gzdopen, _zlib), Ptr{Cvoid}, (Int32, Ptr{UInt8}), dup_fd, gzmode)
    if gz_file == C_NULL
        throw(GZError(-1, "gzdopen failed"))
    end
    if gz_buf_size != Z_DEFAULT_BUFSIZE
        if gzbuffer(gz_file, gz_buf_size) == -1
            # Generally a non-fatal error, although it shouldn't happen here
            gz_buf_size = Z_DEFAULT_BUFSIZE
        end
    end
    s = GZipStream(name, gz_file, gz_buf_size)
    peek(s) # Set EOF-bit for empty files
    return s
end
gzdopen(fd::Integer, gzmode::AbstractString, gz_buf_size::Integer) = gzdopen(string("<fd ",fd,">"), fd, gzmode, gz_buf_size)
gzdopen(fd::Integer, gz_buf_size::Integer) = gzdopen(fd, "rb", gz_buf_size)
gzdopen(fd::Integer, gzmode::AbstractString) = gzdopen(fd, gzmode, Z_DEFAULT_BUFSIZE)
gzdopen(fd::Integer) = gzdopen(fd, "rb", Z_DEFAULT_BUFSIZE)
gzdopen(s::IOStream, args...) = gzdopen(fd(s), args...)

fd(s::GZipStream) = error("fd is not supported for GZipStreams")

function close(s::GZipStream)
    if s._closed
        return Z_STREAM_ERROR
    end
    s._closed = true

    s.name *= " (closed)"

    ret = (@test_z_ok ccall((:gzclose, _zlib), Int32, (Ptr{Cvoid},), s.gz_file))

    return ret
end

flush(s::GZipStream, fl::Integer) =
    @test_z_ok ccall((:gzflush, _zlib), Int32, (Ptr{Cvoid}, Int32), s.gz_file, Int32(fl))
flush(s::GZipStream) = flush(s, Z_SYNC_FLUSH)

truncate(s::GZipStream, n::Integer) = throw(MethodError(truncate, (GZipStream, Integer)))

# Note: seeks to byte position within uncompressed data stream
function seek(s::GZipStream, n::Integer)
    # Note: band-aid to avoid a bug occurring on uncompressed files under Windows
    @static if Sys.iswindows()
        if (ccall((_gzdirect, _zlib), Cint, (Ptr{Cvoid},), s.gz_file)) == 1
            ccall((_gzrewind, _zlib), Cint, (Ptr{Cvoid},), s.gz_file)!=-1 ||
                error("seek (gzseek) failed")
        end
    end
    ccall((_gzseek, _zlib), ZFileOffset, (Ptr{Cvoid}, ZFileOffset, Int32),
           s.gz_file, n, SEEK_SET)!=-1 || # Mimic behavior of seek(s::IOStream, n)
        error("seek (gzseek) failed")
end

# Note: skips bytes within uncompressed data stream
skip(s::GZipStream, n::Integer) =
    (ccall((_gzseek, _zlib), ZFileOffset, (Ptr{Cvoid}, ZFileOffset, Int32),
           s.gz_file, n, SEEK_CUR)!=-1 ||
     error("skip (gzseek) failed")) # Mimic behavior of skip(s::IOStream, n)

if GZLIB_VERSION > "1.2.3.9"
position(s::GZipStream, raw::Bool=false) = raw ?
    ccall((_gzoffset, _zlib), ZFileOffset, (Ptr{Cvoid},), s.gz_file) :
      ccall((_gztell, _zlib), ZFileOffset, (Ptr{Cvoid},), s.gz_file)
else
position(s::GZipStream, raw::Bool=false) =
      ccall((_gztell, _zlib), ZFileOffset, (Ptr{Cvoid},), s.gz_file)
end

eof(s::GZipStream) = Bool(ccall((:gzeof, _zlib), Int32, (Ptr{Cvoid},), s.gz_file))

function peek(s::GZipStream)
    c = gzgetc_raw(s)
    if c != -1
        gzungetc(c, s)
    end
    c
end

# Mimics read(s::IOStream, a::Array{T})
function read(s::GZipStream, a::Array{T}) where {T}
    if isbits(T)
        nb = length(a)*sizeof(T)
        # Note: this will overflow and succeed without warning if nb > 4GB
        ret = ccall((:gzread, _zlib), Int32,
                    (Ptr{Cvoid}, Ptr{Cvoid}, UInt32), s.gz_file, a, nb)
        if ret == -1
            throw(GZError(s))
        end
        if ret < nb
            throw(EOFError())  # TODO: Do we have/need a way to read without throwing an error near the end of the file?
        end
        peek(s) # force eof to be set
        a
    else
        invoke(read!, Tuple{IO,Array}, s, a)
    end
end

function read(s::GZipStream, ::Type{UInt8})
    ret = gzgetc(s)
    if ret == -1
        throw(GZError(s))
    end
    peek(s) # force eof to be set
    UInt8(ret)
end

# For this function, it's really unfortunate that zlib is
# not integrated with ios
function read(s::GZipStream, ::Type{String}; bufsize::Int = Z_BIG_BUFSIZE)
    buf = Array{UInt8}(undef, bufsize)
    len = 0
    while true
        ret = gzread(s, pointer(buf)+len, bufsize)
        if ret == 0
            # check error status to make sure stream was not truncated
            # (we won't normally get an error until the close, because it's
            # possible that the file is still being written to.)

            ## *** Disabled, to allow the function to return the buffer ***
            ## *** Truncation error will be generated on gzclose... ***

            #(err, msg) = gzerror(s)
            #if err != Z_OK
            #    throw(GZError(err, msg))
            #end

            # Resize buffer to exact length
            if length(buf) > len
                resize!(buf, len)
            end
            return String(copy(buf))
        end
        len += ret
        # Grow the buffer so that bufsize bytes will fit
        resize!(buf, bufsize+len)
    end
end

function readline(s::GZipStream)
    buf = Array{UInt8}(undef, GZ_LINE_BUFSIZE)
    pos = 1

    if gzgets(s, buf) == C_NULL      # Throws an exception on error
        return ""
    end

    while(true)
        # since gzgets didn't return C_NULL, there must be a \0 in the buffer
        # eos = search(buf, '\0', pos)
        eos = findnext(x->x==UInt8('\0'), buf, pos)::Int
        # @assert eos ≢ nothing
        if eos == 1 || buf[eos-1] == UInt8('\n')
            return String(copy(resize!(buf, eos-1)))
        end

        # If we're at the end of the file, return the string
        if eof(s)
            return String(copy(resize!(buf, eos-1)))
        end

        # Otherwise, append to the end of the previous buffer

        # Grow the buffer so that there's room for GZ_LINE_BUFSIZE chars
        add_len = GZ_LINE_BUFSIZE - (length(buf)-eos+1)
        resize!(buf, add_len+length(buf))
        pos = eos

        # Read in the next chunk
        if gzgets(s, pointer(buf)+pos-1, GZ_LINE_BUFSIZE) == C_NULL
            # eof(s); remove extra buffer space
            return String(copy(resize!(buf, length(buf)-add_len)))
        end
    end
end

write(s::GZipStream, b::UInt8) = gzputc(s, b)
write(s::GZipStream, a::Array{UInt8}) = gzwrite(s, pointer(a), sizeof(a))
unsafe_write(s::GZipStream, p::Ptr{UInt8}, nb::UInt) = gzwrite(s, p, nb)

function write(s::GZipStream, a::SubArray{T,N,Array}) where {T,N}
    if !isbits(T) || stride(a,1)!=1
        return invoke(write, Tuple{Any,AbstractArray}, s, a)
    end
    colsz = size(a,1)*sizeof(T)
    if N==1
        write(s, pointer(a, 1), colsz)
    else
        cartesian_map((idxs...)->write(s, pointer(a, idxs), colsz),
                      tuple(1, size(a)[2:end]...))
    end
end

end # module GZip


text:
XFloat16(x::Float16)  = reinterpret(XFloat16, Float32(x))
XFloat16(x::Float64)  = reinterpret(XFloat16, Float32(x))
XFloat16(x::BigFloat) = reinterpret(XFloat16, Float32(x))

XFloat32(x::Float16)  = reinterpret(XFloat32, Float64(x))
XFloat32(x::Float32)  = reinterpret(XFloat32, Float64(x))
XFloat32(x::BigFloat) = reinterpret(XFloat32, Float64(x))

Float16(x::XFloat16)  = Float16(reinterpret(Float32, x))
Float64(x::XFloat16)  = Float64(reinterpret(Float32, x))
BigFloat(x::XFloat16) = BigFloat(reinterpret(Float32, x))

Float16(x::XFloat32)  = Float16(reinterpret(Float64, x))
Float32(x::XFloat32)  = Float32(reinterpret(Float64, x))
BigFloat(x::XFloat32) = BigFloat(reinterpret(Float64, x))

for T in (:Int8, :Int16, :Int32, :Int64, :Int128, :BigInt,
          :UInt8, :UInt16, :UInt32, :UInt64, :UInt128)
  @eval begin
    XFloat16(x::$T) = reinterpret(XFloat16, Float32(x))
    XFloat32(x::$T) = reinterpret(XFloat32, Float64(x))
    $T(x::XFloat16) = $T(reinterpret(Float32,x))
    $T(x::XFloat32) = $T(reinterpret(Float64,x))
  end
end

XFloat16(x::Bool) = reinterpret(XFloat16, Float32(x))
XFloat32(x::Bool) = reinterpret(XFloat32, Float64(x))
Bool(x::XFloat16) = Bool(reinterpret(Float32, x))
Bool(x::XFloat32) = Bool(reinterpret(Float64, x))

XFloat16(x::Irrational{S}) where {S} = XFloat16(Float32(x))
XFloat32(x::Irrational{S}) where {S} = XFloat32(Float64(x))

function Base.frexp(x::XFloat16)
  fr,xp = frexp(Float32(x))
  return XFloat16(fr), xp
end

function Base.frexp(x::XFloat32)
  fr,xp = frexp(Float64(x))
  return XFloat32(fr), xp
end

function Base.ldexp(x::XFloat16, xp::I) where {I<:Integer}
   fr = Float32(x)
   ldxp = ldexp(fr, xp)
   return XFloat16(ldxp)
end

function Base.ldexp(x::XFloat32, xp::I) where {I<:Integer}
   fr = Float64(x)
   ldxp = ldexp(fr, xp)
   return XFloat32(ldxp)
end


text:
function make_fqparser( filename; forcegzip=false )
   fopen = open( filename, "r" )
   if isgzipped( filename ) || forcegzip
      to_open = ZlibInflateInputStream( fopen, reset_on_end=true )
   else
      to_open = BufferedInputStream( fopen )
   end
   FASTQ.Reader( to_open, fill_ambiguous=DNA_A )
end

# modified for Bio v0.2 with tryread_bool!
@inline function read_chunk!( chunk, parser )
   i = 1
   while i <= length(chunk) && !eof(parser)
      read!( parser, chunk[i] )
      i += 1
   end
   while i <= length(chunk)
      pop!(chunk) # clean up if we are at the end
   end
   parser
end

function allocate_chunk( parser; size=10000 )
  chunk = Vector{eltype(parser)}( size )
  for i in 1:length(chunk)
     chunk[i] = eltype(parser)()
  end
  chunk
end

function allocate_fastq_records( size::Int=10000 )
   chunk = Vector{FASTQRecord}( undef, size )
   for i in 1:length(chunk)
      chunk[i] = FASTQRecord()
   end
   chunk
end

function process_reads!( parser, param::AlignParam, lib::GraphLib, quant::GraphLibQuant,
                         multi::MultiMapping{SGAlignSingle}, mod::B;
                         bufsize=150, sam=false, qualoffset=33 ) where B <: BiasModel

   reads  = allocate_fastq_records( bufsize )
   mean_readlen = 0.0
   total        = 0
   mapped       = 0
   if sam
      stdbuf = BufferedOutputStream( stdout )
      write_sam_header( stdbuf, lib )
   end

   fqfile = open("BTLA-A1.fastq", "w")
   fqwriter = FASTQ.Writer(fqfile)

   while length(reads) > 0
      read_chunk!( reads, parser )
      total += length(reads)
#      println(stderr, "$(length(reads)) reads/$total............")
      @inbounds for i in 1:length(reads)
#= Lots of these 'reads' objects have an 'EMPTY SEQUENCE' and the ones that do not
   appear to have a sequence that does not match the .sequence data.  Here's one.

Whippet.FASTQRecord(< EMPTY SEQUENCE >, UInt8[], FASTX.FASTQ.Record:
   identifier: HISEQ:587:C6RMKANXX:1:1307:3191:12232/1
  description: <missing>
     sequence: CTCAGAGTGAAAGGCTGGAAAACAAATTTCCAAGCAAATGGTCTGAAGAAACAAGCTGGAGTAGCCATTCTAATATCGAATAAAATCGACTTCCAACCCAAAGTTATCAAAAAAGACAAGGAGGGA
      quality: UInt8[...])
=#
         fill!( reads[i], qualoffset )
#= The 'fill' this refers to appears to be that it fills in the EMPTY SEQUENCE with
   the sequence in the '.sequence' field.
=#
         align = ungapped_align( param, lib, reads[i] )
         if !isnull( align )
            ## Builds us a little description string to contain scores of the matches
            ## and mismatches for this gene.
            descString = ""
	    for k in 1:length(align.value)
               for m in 1:length(align.value[k].path)
                  descString *= lib.names[align.value[k].path[m].gene] * ":";
                  descString *= string(align.value[k].path[m].node) * "(";
                  descString *= string(align.value[k].path[m].score.matches) * ",";
                  descString *= string(align.value[k].path[m].score.mismatches) * ",";
                  descString *= string(align.value[k].path[m].score.mistolerance) * ")/";
               end
            end
            ## Writes a FASTQ file for any reads that match the given gene.
            if lib.names[first(first(align.value).path).gene] == "ENSMUSG00000052013"
               write(fqwriter,
                     FASTQ.Record(identifier(reads[i].raw),
                                  descString,
                                  sequence(reads[i].raw),
                                  quality(reads[i].raw); offset=33));
            end
            biasval = count!( mod, reads[i].sequence )
            if length( align.value ) > 1
               push!( multi, align.value, biasval, quant, lib )
               sam && write_sam( stdbuf, reads[i], align.value, lib, qualoffset=qualoffset )
            else
               count!( quant, align.value[1], biasval )
               sam && write_sam( stdbuf, reads[i], align.value[1], lib, qualoffset=qualoffset )
            end
            mapped += 1
            @fastmath mean_readlen += (length(reads[i].sequence) - mean_readlen) / mapped
         end
      end
      if total % 100000 == 0
         GC.gc()
      end
   end # end while

   close(fqfile)

   if sam
      close(stdbuf)
   end
   mapped,total,mean_readlen
end

function process_paired_reads!( fwd_parser, rev_parser, param::AlignParam,
                                lib::GraphLib, quant::GraphLibQuant,
                                multi::MultiMapping{SGAlignPaired}, mod::B;
                                bufsize=50, sam=false, qualoffset=33 ) where B <: BiasModel

   fwd_reads  = allocate_fastq_records( bufsize )
   rev_reads  = allocate_fastq_records( bufsize )
   mean_readlen = 0.0
   total        = 0
   mapped       = 0
   if sam
      stdbuf = BufferedOutputStream( stdout )
      write_sam_header( stdbuf, lib )
   end
   while length(fwd_reads) > 0 && length(rev_reads) > 0
      read_chunk!( fwd_reads, fwd_parser )
      read_chunk!( rev_reads, rev_parser )
      total += length(fwd_reads)
      @inbounds for i in 1:length(fwd_reads)
         fill!( fwd_reads[i], qualoffset )
         fill!( rev_reads[i], qualoffset )
         fwd_aln,rev_aln = ungapped_align( param, lib, fwd_reads[i], rev_reads[i] )
         if !isnull( fwd_aln ) && !isnull( rev_aln )
            biasval = count!( mod, fwd_reads[i].sequence, rev_reads[i].sequence )
            if length( fwd_aln.value ) > 1
               push!( multi, fwd_aln.value, rev_aln.value, biasval, quant, lib )
               sam && write_sam( stdbuf, fwd_reads[i], rev_reads[i], fwd_aln.value, rev_aln.value, lib,
                                 paired=true, is_pair_rc=param.is_pair_rc, qualoffset=qualoffset )
            else
               count!( quant, fwd_aln.value[1], rev_aln.value[1], biasval )
               sam && write_sam( stdbuf, fwd_reads[i], fwd_aln.value[1], lib,
                                 paired=true, fwd_mate=true, is_pair_rc=param.is_pair_rc,
                                 qualoffset=qualoffset )
               sam && write_sam( stdbuf, rev_reads[i], rev_aln.value[1], lib,
                                 paired=true, fwd_mate=false, is_pair_rc=param.is_pair_rc,
                                 qualoffset=qualoffset )
            end
            mapped += 1
            @fastmath mean_readlen += (length(fwd_reads[i].sequence) - mean_readlen) / mapped
         end
      end
   end # end while
   if sam
      close(stdbuf)
   end
   mapped,total,mean_readlen
end


text:
using Test
using SmartShockFinder

EXAMPLES_DIR = joinpath(pathof(SmartShockFinder) |> dirname |> dirname, "examples")

@test_nowarn rm(joinpath("datasets"), recursive=true, force=true)

@time @testset "SmartShockFinder" begin
  @test_nowarn mkpath(joinpath("datasets", "1d"))

  @testset "examples/1d/traindata_NNPPh.jl" begin
    @test_nowarn include(joinpath(EXAMPLES_DIR, "1d", "traindata_NNPPh.jl"))

    @test isfile(joinpath("datasets", "1d", "traindata1dNNPPh.h5"))
  end
end

@test_nowarn rm(joinpath("datasets"), recursive=true, force=true)


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

#
# 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

@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 +(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) + rheologrun(rl2))
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 -(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) - rheologrun(rl2))
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

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

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

#=
---------------------------
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)

function RheoModelClass(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Expr = quote NaN end,
                         J::Expr = quote NaN end,
                         Gp::Expr = quote NaN end,
                         Gpp::Expr = quote NaN end,
                         constraint::Expr = quote true end,
                         info=name)

    # 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)

    @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 model_parameters(nt::NamedTuple, params::Vector{Symbol}, err_string::String)
    # check that every parameter in m exists in the named tuple nt
    @assert all(i->i in keys(nt),params) "Missing parameter(s) in " * err_string
    # check that no extra parameters have been provided
    @assert length(params) == length(nt) "Mismatch number of model parameters and parameters provided in " * err_string

    p = map(i->RheoFloat(nt[i]), params)
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)}( ( G, J, Gp, Gpp ) )

    @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 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

"""
    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 it's 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 = model_parameters(nt0, m.params,"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)}( ( G, J, Gp, Gpp ) )

    @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


text:
# Script to compute value added losses due to brexit

# so I dont have to open the scripts and run them, should be replaced by a function in the end
include("src/leontief_matrix.jl")
include("src/mfn_tariffs.jl")
include("src/trade_elasticities.jl")

# ---------- Transform dataframes into matrices for computation ------------------------------------------------------------------------------------------------------

soft_brexit = all_tariffs.NTB_soft # N*S×1
hard_brexit = all_tariffs.avg_tariff .+ all_tariffs.NTB_hard # N*S×1
elasticities = all_elasticities.sigma # N*S×1

begin
    EU27 = ["AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA", "LVA",
    "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN", "ESP", "SWE"]
end

ctry = sort([EU27; "GBR"]) # all vectors/matrices are sorted alphabetically
position_GBR = findfirst(x -> x == "GBR", ctry) # find position of GBR

exports_to_GBR = exports[:, position_GBR] # N*S×1

# take out GBR rows
N = 27 # EU27 (no GBR)
S = 56

soft_brexit = soft_brexit[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1
hard_brexit = hard_brexit[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1
elasticities = elasticities[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1
exports_to_GBR = exports_to_GBR[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1
#v_add_share = v_add_share[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1
v_add_share = v_add_share[Not(position_GBR*S-S+1:position_GBR*S)] # N*S×1

L = L[Not(position_GBR*S-S+1:position_GBR*S), Not(position_GBR*S-S+1:position_GBR*S)] # N*S×N*S

# ---------- Computation -------------------------------------------------------------------------------------------------------------------------------------------

M_soft_brexit = similar(L) # N*S×N*S
M_hard_brexit = similar(L) # N*S×N*S

for i in 1:size(L, 1)
    for j in 1:size(L, 2)
        M_soft_brexit[i, j] = (elasticities[i] - 1) * soft_brexit[i]/100 * L[i, j] * exports_to_GBR[i] # N*S×N*S
        M_hard_brexit[i, j] = (elasticities[i] - 1) * hard_brexit[i]/100 * L[i, j] * exports_to_GBR[i] # N*S×N*S
    end
end

replace!(v_add_share, NaN => 0.0) # remove all NaN, mostly sectors U, T
# does value added need to stay within (0,1)?
v_add_share = [x < 1.0 && x > 0.0 ? x : 0.0 for x in v_add_share] # remove negative values and values above 1

repeat(1:S:size(M_soft_brexit, 2), inner=S)

# sum over all columns, i.e. per reporting country-sector pair (first S sectors are direct effects) but we need to select the right partner country (domestic sectors)!
index = hcat(1:size(M_soft_brexit, 1), repeat(1:S:size(M_soft_brexit, 2), inner=S)) # index for iterating over

direct_soft = [sum(M_soft_brexit[index[i, 1], index[i, 2]:index[i, 2]+S-1]) for i in 1:size(index, 1)] # N*S×1
direct_hard = [sum(M_hard_brexit[index[i, 1], index[i, 2]:index[i, 2]+S-1]) for i in 1:size(index, 1)] # N*S×1

indirect_soft = [sum(M_soft_brexit[index[i, 1], Not(index[i, 2]:index[i, 2]+S-1)]) for i in 1:size(index, 1)] # N*S×1
indirect_hard = [sum(M_hard_brexit[index[i, 1], Not(index[i, 2]:index[i, 2]+S-1)]) for i in 1:size(index, 1)] # N*S×1

# value added loss per country-sector pair
direct_loss_soft = -v_add_share .* direct_soft # N*S×1
direct_loss_hard = -v_add_share .* direct_hard # N*S×1

indirect_loss_soft = -v_add_share .* indirect_soft # N*S×1
indirect_loss_hard = -v_add_share .* indirect_hard # N*S×1

total_loss_soft = direct_loss_soft .+ indirect_loss_soft # N*S×1
total_loss_hard = direct_loss_hard .+ indirect_loss_hard # N*S×1

# value added loss per country
direct_ctry_loss_soft = [sum(direct_loss_soft[i:i+S-1]) for i in 1:S:S*N] # N×1
direct_ctry_loss_hard = [sum(direct_loss_hard[i:i+S-1]) for i in 1:S:S*N] # N×1

indirect_ctry_loss_soft = [sum(indirect_loss_soft[i:i+S-1]) for i in 1:S:S*N] # N×1
indirect_ctry_loss_hard = [sum(indirect_loss_hard[i:i+S-1]) for i in 1:S:S*N] # N×1

total_ctry_loss_soft = [sum(total_loss_soft[i:i+S-1]) for i in 1:S:S*N] # N×1
total_ctry_loss_hard = [sum(total_loss_hard[i:i+S-1]) for i in 1:S:S*N] # N×1

# some plots and tables

ctry_loss_wide = hcat(sort(EU27), direct_ctry_loss_soft, direct_ctry_loss_hard, indirect_ctry_loss_soft, indirect_ctry_loss_hard, total_ctry_loss_soft, total_ctry_loss_hard)
ctry_loss_wide = DataFrame(ctry_loss_wide, :auto)
rename!(ctry_loss_wide, [:iso3, :direct_soft, :direct_hard, :indirect_soft, :indirect_hard, :total_soft, :total_hard])

XLSX.writetable("clean/all_losses.xlsx", ctry_loss_wide, overwrite=true)

ctry_loss_long = stack(ctry_loss_wide, Not(:iso3))
ctry_loss_long = hcat(ctry_loss_long, DataFrame(reduce(vcat, permutedims.(split.(ctry_loss_long.variable, "_"))), [:channel, :scenario]))
select(ctry_loss_long, Not(:variable))

using StatsPlots

function grouped_bar_plot(scenario, channel, type, name)
    df = filter(row -> row.channel in channel && row.scenario == scenario, ctry_loss_long)
    plot = groupedbar(df[:, :iso3], -df.value/1000, group=df.channel, bar_width=0.8, bar_position=type, xrotation=90, color=[:red :blue])
    xlabel!("EU27")
    ylabel!("billion euro")
    title!("Loss in value added in a $scenario Brexit scenario")
    savefig(plot, "images/$name.png")
    #savefig(plot, "clean/$name.html")
    return plot
end

# direct vis-a-vis indirect loss per country
grouped_bar_plot("soft", ["direct", "indirect"], :dodge, "soft_split")
grouped_bar_plot("hard", ["direct", "indirect"], :dodge, "hard_split")

# total loss per country
grouped_bar_plot("soft", ["direct", "indirect"], :dodge, "soft_total")
grouped_bar_plot("hard", ["direct", "indirect"], :dodge, "hard_total")


text:
using DFTK
using Test
include("testcases.jl")

@testset "Forces on semiconductor (using total energy)" begin
    function energy(pos)
        Ecut = 5
        Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
        atoms = [Si => pos]
        model = model_DFT(silicon.lattice, atoms, :lda_xc_teter93)
        basis = PlaneWaveBasis(model, Ecut, kgrid=[2, 2, 2])

        is_converged = DFTK.ScfConvergenceDensity(1e-10)
        scfres = self_consistent_field(basis; is_converged=is_converged)
        sum(values(scfres.energies)), forces(scfres)
    end

    # symmetrical positions, forces should be 0
    pos0 = [(ones(3)) / 8, -ones(3) / 8]
    _, F1 = energy(pos0)
    @test norm(F1) < 1e-4

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8]  # displace a bit from equilibrium
    disp = rand(3)
    ε = 1e-8
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    E1, F1 = energy(pos1)
    E2, _ = energy(pos2)

    diff_findiff = -(E2 - E1) / ε
    diff_forces = dot(F1[1][1], disp)
    @test abs(diff_findiff - diff_forces) < 1e-4
end

@testset "Forces on metal (using free energy)" begin
    function silicon_energy_forces(pos; smearing=Smearing.FermiDirac())
        Ecut = 4
        Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
        atoms = [Si => pos]
        model = model_DFT(silicon.lattice, atoms, :lda_xc_teter93;
                          temperature=0.03, smearing=smearing)
        # TODO Put kshift=[1/2, 0, 0] here later
        basis = PlaneWaveBasis(model, Ecut, kgrid=[2, 1, 2])

        n_bands = 10
        is_converged = DFTK.ScfConvergenceDensity(5e-10)
        scfres = self_consistent_field(basis, n_bands=n_bands,
                                       is_converged=is_converged,
                                      )
        sum(values(scfres.energies)), forces(scfres)
    end

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8] # displace a bit from equilibrium
    disp = rand(3)
    ε = 1e-6
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    for (tol, smearing) in [(0.001, Smearing.FermiDirac), (5e-5, Smearing.Gaussian)]
        E1, F1 = silicon_energy_forces(pos1; smearing=smearing())
        E2, _ = silicon_energy_forces(pos2; smearing=smearing())

        diff_findiff = -(E2 - E1) / ε
        diff_forces = dot(F1[1][1], disp)

        @test abs(diff_findiff - diff_forces) < tol
    end
end


text:
# Constants used in the device API.

# Constant to represent the transmit direction
const SOAPY_SDR_TX = 0

# Constant to represent the receive direction
const SOAPY_SDR_RX = 1

# Indicate end of burst for transmit or receive.
# For write, end of burst if set by the caller.
# For read, end of burst is set by the driver.
const SOAPY_SDR_END_BURST  = (1 << 1)

# Indicates that the time stamp is valid.
# For write, the caller must set has time when timeNs is provided.
# For read, the driver sets has time when timeNs is provided.
const SOAPY_SDR_HAS_TIME = (1 << 2)

# Indicates that stream terminated prematurely.
# This is the flag version of an overflow error
# that indicates an overflow with the end samples.
const SOAPY_SDR_END_ABRUPT = (1 << 3)

# Indicates transmit or receive only a single packet.
# Applicable when the driver fragments samples into packets.
# For write, the user sets this flag to only send a single packet.
# For read, the user sets this flag to only receive a single packet.
const SOAPY_SDR_ONE_PACKET = (1 << 4)

# Indicate that this read call and the next results in a fragment.
# Used when the implementation has an underlying packet interface.
# The caller can use this indicator and the SOAPY_SDR_ONE_PACKET flag
# on subsequent read stream calls to re-align with packet boundaries.
const SOAPY_SDR_MORE_FRAGMENTS = (1 << 5)

# Indicate that the stream should wait for an external trigger event.
# This flag might be used with the flags argument in any of the
# stream API calls. The trigger implementation is hardware-specific.
const SOAPY_SDR_WAIT_TRIGGER = (1 << 6)


text:
# Helper functions for training Flux models
using Flux: train!, @epochs, params, SGD, mse

function train_and_eval!(model, tr_data, val_data, n_epochs)
    loss(x, y) = mse(model(x), y)
    opt = SGD(params(model), η=0.001)
    @epochs n_epochs train!(loss, tr_data, opt)

end

max_stars_repo_name: dillondaudert/neural-alu

max_issues_repo_name: dillondaudert/neural-alu

max_forks_repo_name: dillondaudert/neural-alu


text:
import Base.Libdl: dlopen
export stinger_core_lib, stinger_alg_lib

if "STINGER_LIB_PATH" in keys(ENV)
    const stinger_core_lib = dlopen(joinpath(ENV["STINGER_LIB_PATH"],"libstinger_core"))
    const stinger_alg_lib = dlopen(joinpath(ENV["STINGER_LIB_PATH"],"libstinger_alg"))
else
    basedir = joinpath(dirname(@__DIR__), "deps", "usr", "lib")
    const stinger_core_lib = dlopen(joinpath(basedir, "libstinger_core"))
    const stinger_core_alg = dlopen(joinpath(basedir, "libstinger_alg"))
end


text:
# Load X and y variable
using JLD
data = load("multiData.jld")
(X,y,Xtest,ytest) = (data["X"],data["y"],data["Xtest"],data["ytest"])

# Data is already roughly standardized, but let's add bias
n = size(X,1)
X = [ones(n,1) X]

# Do the same transformation to the test data
t = size(Xtest,1)
Xtest = [ones(t,1) Xtest]

# Fit one-vs-all logistic regression model
#include("logReg.jl")
#model = logRegOnevsAll(X,y)
include("softmaxClassifier.jl")
model = softmaxClassifier(X,y)

# Compute training and validation error
yhat = model.predict(X)
trainError = mean(yhat .!= y)
@show(trainError)
yhat = model.predict(Xtest)
validError = mean(yhat .!= ytest)
@show(validError)

# Plot results
k = maximum(y)
include("plot2Dclassifier.jl")
plot2Dclassifier(X,y,model,Xtest=Xtest,ytest=ytest,biasIncluded=true,k=5)
show()


text:
function basis(elem::Quad, N, r, s)
    Np = convert(Int, (N + 1) * (N + 1))
    sk = 1
    V, Vr, Vs = ntuple(x->zeros(length(r), Np), 3)
    for j = 0:N
        P_j = jacobiP(s, 0, 0, j)
        for i = 0:N
            P_i = jacobiP(r, 0, 0, i)
            V[:, sk]  = P_i .* P_j
            Vr[:, sk] = grad_jacobiP(r, 0, 0, i) .* P_j
            Vs[:, sk] = P_i .* grad_jacobiP(s, 0, 0, j)
            sk += 1
        end
    end
    return V,Vr,Vs
end

function nodes(elem::Quad, N)
    r1D, w1D = gauss_lobatto_quad(0, 0, N)
    s, r = meshgrid(r1D)
    return r[:], s[:]
end

function equi_nodes(elem::Quad, N)
    r1D = LinRange(-1, 1, N + 1)
    s, r = meshgrid(r1D)
    return r[:], s[:]
end

function quad_nodes(elem::Quad,N)
    r1D, w1D = gauss_quad(0, 0, N)
    s, r = meshgrid(r1D)
    ws, wr = meshgrid(w1D)
    w = @. wr * ws
    return r[:], s[:], w[:]
end


text:
module UnderwaterAcoustics

using DSP
using LinearAlgebra
using SignalAnalysis

include "basic.jl"
include "sources.jl"
include "pm_core.jl"
include "pm_3ray.jl"

end # module


text:
abstract type AbstractGrid{dim, T} end

const Vec = Ferrite.Vec

"""
```
struct RectilinearGrid{dim, T, N, M, TG<:Ferrite.Grid{dim, <:Ferrite.Cell{dim,N,M}, T}} <: AbstractGrid{dim, T}
    grid::TG
    nels::NTuple{dim, Int}
    sizes::NTuple{dim, T}
    corners::NTuple{2, Vec{dim, T}}
    white_cells::BitVector
    black_cells::BitVector
    constant_cells::BitVector
end
```

A type that represents a rectilinear grid with corner points `corners`.

- `dim`: dimension of the problem
- `T`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `M`: number of faces in a cell of the grid
- `grid`: a Ferrite.Grid struct
- `nels`: number of elements in every dimension
- `sizes`: dimensions of each rectilinear cell
- `corners`: 2 corner points of the rectilinear grid
- `white_cells`: cells fixed to be void during optimization
- `black_cells`: cells fixed to have material during optimization
- `constant_cells`: cells fixed to be either void or have material during optimization
"""
@params struct RectilinearGrid{dim, T, N, M, TG<:Ferrite.Grid{dim, <:Ferrite.Cell{dim,N,M}, T}} <: AbstractGrid{dim, T}
    grid::TG
    nels::NTuple{dim, Int}
    sizes::NTuple{dim, T}
    corners::NTuple{2, Vec{dim,T}}
    white_cells::BitVector
    black_cells::BitVector
    constant_cells::BitVector
end

"""
    RectilinearGrid(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim,T}) where {dim, T, CellType}

Constructs an instance of [`RectilinearGrid`](@ref).

- `dim`: dimension of the problem
- `T`: number type for coordinates
- `nels`: number of elements in every dimension
- `sizes`: dimensions of each rectilinear cell

Example:

```
rectgrid = RectilinearGrid((60,20), (1.0,1.0))
```
"""
function RectilinearGrid(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim,T}) where {dim, T, CellType}
    if dim === 2
        if CellType === :Linear
            geoshape = Quadrilateral
        else
            geoshape = QuadraticQuadrilateral
        end
    else
        geoshape = Hexahedron
    end
    corner1 = Vec{dim}(fill(T(0), dim))
    corner2 = Vec{dim}((nels .* sizes))
    grid = generate_grid(geoshape, nels, corner1, corner2);

    N = nnodes(geoshape)
    M = Ferrite.nfaces(geoshape)
    ncells = prod(nels)
    return RectilinearGrid(grid, nels, sizes, (corner1, corner2), falses(ncells), falses(ncells), falses(ncells))
end

nnodespercell(::RectilinearGrid{dim,T,N,M}) where {dim, T, N, M} = N
nfacespercell(::RectilinearGrid{dim,T,N,M}) where {dim, T, N, M} = M

left(rectgrid::RectilinearGrid, x) = x[1] ≈ rectgrid.corners[1][1]
right(rectgrid::RectilinearGrid, x) = x[1] ≈ rectgrid.corners[2][1]
bottom(rectgrid::RectilinearGrid, x) = x[2] ≈ rectgrid.corners[1][2]
top(rectgrid::RectilinearGrid, x) = x[2] ≈ rectgrid.corners[2][2]
back(rectgrid::RectilinearGrid, x) = x[3] ≈ rectgrid.corners[1][3]
front(rectgrid::RectilinearGrid, x) = x[3] ≈ rectgrid.corners[2][3]
middlex(rectgrid::RectilinearGrid, x) = x[1] ≈ (rectgrid.corners[1][1] + rectgrid.corners[2][1]) / 2
middley(rectgrid::RectilinearGrid, x) = x[2] ≈ (rectgrid.corners[1][2] + rectgrid.corners[2][2]) / 2
middlez(rectgrid::RectilinearGrid, x) = x[3] ≈ (rectgrid.corners[1][3] + rectgrid.corners[2][3]) / 2

nnodes(cell::Type{Ferrite.Cell{dim,N,M}}) where {dim, N, M} = N
nnodes(cell::Ferrite.Cell) = nnodes(typeof(cell))

"""
    LGrid(::Type{Val{CellType}}, ::Type{T}; length = 100, height = 100, upperslab = 50, lowerslab = 50) where {T, CellType}
    LGrid(::Type{Val{CellType}}, nel1::NTuple{2,Int}, nel2::NTuple{2,Int}, LL::Vec{2,T}, UR::Vec{2,T}, MR::Vec{2,T}) where {CellType, T}

Constructs a `Ferrite.Grid` that represents the following L-shaped grid.

```
        upperslab   UR
       ............
       .          .
       .          .
       .          .
height .          .                     MR
       .          ......................
       .                               .
       .                               . lowerslab
       .                               .
       .................................
     LL             length

```

Examples:

```
LGrid(upperslab = 30, lowerslab = 70)
LGrid(Val{:Linear}, (2, 4), (2, 2), Vec{2,Float64}((0.0,0.0)), Vec{2,Float64}((2.0, 4.0)), Vec{2,Float64}((4.0, 2.0)))
```
"""
function LGrid(::Type{Val{CellType}}, ::Type{T}; length = 100, height = 100, upperslab = 50, lowerslab = 50) where {T, CellType}
    @assert length > upperslab
    @assert height > lowerslab
    LGrid(Val{CellType}, (upperslab, height), (length-upperslab, lowerslab),
        Vec{2,T}((0.0,0.0)), Vec{2,T}((T(upperslab), T(height))),
        Vec{2,T}((T(length), T(lowerslab))))
end
function LGrid(::Type{Val{CellType}}, nel1::NTuple{2,Int}, nel2::NTuple{2,Int},
    LL::Vec{2,T}, UR::Vec{2,T}, MR::Vec{2,T}) where {CellType, T}

    if CellType === :Linear
        return _LinearLGrid(nel1, nel2, LL, UR, MR)
    else
        return _QuadraticLGrid(nel1, nel2, LL, UR, MR)
    end
end

function _LinearLGrid(nel1::NTuple{2,Int}, nel2::NTuple{2,Int},
    LL::Vec{2,T}, UR::Vec{2,T}, MR::Vec{2,T}) where {T}

    @assert nel1[2] > nel2[2]

    midpointindy = round(Int, nel2[2]/2) + 1
    nodes = Node{2,T}[]
    cells = Quadrilateral[]
    boundary = Tuple{Int,Int}[]
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["right"]  = Set{Tuple{Int,Int}}()
    facesets["top"]    = Set{Tuple{Int,Int}}()
    nodesets = Dict{String, Set{Int}}()
    nodesets["load"] = Set{Int}()

    # Lower left rectangle
    nel_x1 = nel1[1]; nel_y1 = nel2[2];
    n_nodes_x1 = nel_x1 + 1; n_nodes_y1 = nel_y1 + 1
    n_nodes1 = n_nodes_x1 * n_nodes_y1

    _LR = Vec{2,T}((UR[1], LL[2]))
    _UL = Vec{2,T}((LL[1], MR[2]))
    _UR = Vec{2,T}((UR[1], MR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x1, n_nodes_y1, LL, _LR, _UR, _UL)

    node_array1 = reshape(collect(1:n_nodes1), (n_nodes_x1, n_nodes_y1))
    for j in 1:nel_y1, i in 1:nel_x1
        push!(cells, Quadrilateral((node_array1[i,j], node_array1[i+1,j], node_array1[i+1,j+1], node_array1[i,j+1])))
        if i == 1
            push!(boundary, (length(cells), 4))
        end
        if j == 1
            push!(boundary, (length(cells), 1))
        end
    end

    # Lower right rectangle
    offsetstep = (MR[1] - _LR[1])/nel2[1]
    indexoffset = length(nodes)

    nel_x2 = nel2[1] - 1; nel_y2 = nel2[2]
    n_nodes_x2 = nel_x2 + 1; n_nodes_y2 = nel_y2 + 1
    n_nodes2 = n_nodes_x2 * n_nodes_y2

    _LL = Vec{2,T}((_LR[1] + offsetstep, _LR[2]))
    _LR = Vec{2,T}((MR[1], LL[2]))
    _UL = Vec{2,T}((_UR[1] + offsetstep, MR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x2, n_nodes_y2, _LL, _LR, MR, _UL)

    node_array2 = reshape(collect(indexoffset+1:indexoffset+n_nodes2), (n_nodes_x2, n_nodes_y2))
    for j in 1:nel_y2
        push!(cells, Quadrilateral((node_array1[end,j], node_array2[1,j], node_array2[1,j+1], node_array1[end,j+1])))
        j == 1 && push!(boundary, (length(cells), 1))
        j == nel_y2 && push!(boundary, (length(cells), 3))
        if nel_x2 == 1
            push!(boundary, (length(cells), 2))
            push!(facesets["right"], (length(cells), 2))
        end
        for i in 1:nel_x2
            push!(cells, Quadrilateral((node_array2[i,j], node_array2[i+1,j], node_array2[i+1,j+1], node_array2[i,j+1])))
            if i == nel_x2
                push!(boundary, (length(cells), 2))
                push!(facesets["right"], (length(cells), 2))
            end
            j == 1 && push!(boundary, (length(cells), 1))
            j == nel_y2 && push!(boundary, (length(cells), 3))
        end
    end

    push!(nodesets["load"], node_array2[end, midpointindy])

    # Upper left rectangle
    offsetstep = (UR[2] - MR[2])/(nel1[2] - nel2[2])
    indexoffset = length(nodes)

    nel_x3 = nel1[1]; nel_y3 = nel1[2] - nel2[2] - 1
    n_nodes_x3 = nel_x3 + 1; n_nodes_y3 = nel_y3 + 1
    n_nodes3 = n_nodes_x3 * n_nodes_y3

    _LL = Vec{2,T}((LL[1], MR[2] + offsetstep))
    _LR = Vec{2,T}((UR[1], MR[2] + offsetstep))
    _UL = Vec{2,T}((LL[1], UR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x3, n_nodes_y3, _LL, _LR, UR, _UL)

    # Generate cells
    node_array3 = reshape(collect(indexoffset+1:indexoffset+n_nodes3), (n_nodes_x3, n_nodes_y3))

    for i in 1:nel_x3
        push!(cells, Quadrilateral((node_array1[i,end], node_array1[i+1,end], node_array3[i+1,1], node_array3[i,1])))
        i == 1 && push!(boundary, (length(cells), 4))
        i == nel_x3 && push!(boundary, (length(cells), 2))
    end
    for j in 1:nel_y3, i in 1:nel_x3
        push!(cells, Quadrilateral((node_array3[i,j], node_array3[i+1,j], node_array3[i+1,j+1], node_array3[i,j+1])))
        i == 1 && push!(boundary, (length(cells), 4))
        i == nel_x3 && push!(boundary, (length(cells), 2))
        if j == nel_y3
            push!(boundary, (length(cells), 3))
            push!(facesets["top"], (length(cells), 3))
        end
    end

    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)

    return Grid(cells, nodes, facesets=facesets, nodesets=nodesets,
        boundary_matrix=boundary_matrix)
end

function _QuadraticLGrid(nel1::NTuple{2,Int}, nel2::NTuple{2,Int},
    LL::Vec{2,T}, UR::Vec{2,T}, MR::Vec{2,T}) where {T}

    @assert nel1[2] > nel2[2]

    midpointindy = round(Int, nel2[2]/2) + 1
    nodes = Node{2,T}[]
    cells = QuadraticQuadrilateral[]
    boundary = Tuple{Int,Int}[]
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["right"]  = Set{Tuple{Int,Int}}()
    facesets["top"]    = Set{Tuple{Int,Int}}()
    nodesets = Dict{String, Set{Int}}()
    nodesets["load"] = Set{Int}()

    # Lower left rectangle
    nel_x1 = nel1[1]; nel_y1 = nel2[2];
    n_nodes_x1 = 2*nel_x1 + 1; n_nodes_y1 = 2*nel_y1 + 1
    n_nodes1 = n_nodes_x1 * n_nodes_y1

    _LR = Vec{2,T}((UR[1], LL[2]))
    _UL = Vec{2,T}((LL[1], MR[2]))
    _UR = Vec{2,T}((UR[1], MR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x1, n_nodes_y1, LL, _LR, _UR, _UL)

    node_array1 = reshape(collect(1:n_nodes1), (n_nodes_x1, n_nodes_y1))
    for j in 1:nel_y1, i in 1:nel_x1
        push!(cells, QuadraticQuadrilateral((node_array1[2*i-1,2*j-1], node_array1[2*i+1,2*j-1],
                                             node_array1[2*i+1,2*j+1],node_array1[2*i-1,2*j+1],
                                             node_array1[2*i,2*j-1],node_array1[2*i+1,2*j],
                                             node_array1[2*i,2*j+1],node_array1[2*i-1,2*j],
                                             node_array1[2*i,2*j])))
        if i == 1
            push!(boundary, (length(cells), 4))
        end
        if j == 1
            push!(boundary, (length(cells), 1))
        end
    end

    # Lower right rectangle
    offsetstep = (MR[1] - _LR[1])/nel2[1]/2
    indexoffset = length(nodes)

    nel_x2 = nel2[1] - 1; nel_y2 = nel2[2]
    n_nodes_x2 = 2*nel_x2 + 2; n_nodes_y2 = 2*nel_y2 + 1
    n_nodes2 = n_nodes_x2 * n_nodes_y2

    _LL = Vec{2,T}((_LR[1] + offsetstep, _LR[2]))
    _LR = Vec{2,T}((MR[1], LL[2]))
    _UL = Vec{2,T}((_UR[1] + offsetstep, MR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x2, n_nodes_y2, _LL, _LR, MR, _UL)

    node_array2 = reshape(collect(indexoffset+1:indexoffset+n_nodes2), (n_nodes_x2, n_nodes_y2))
    for j in 1:nel_y2
        push!(cells, QuadraticQuadrilateral((node_array1[end,2*j-1], node_array2[2,2*j-1],
                                             node_array2[2,2*j+1], node_array1[end,2*j+1],
                                             node_array2[1,2*j-1], node_array2[2,2*j],
                                             node_array2[1,2*j+1], node_array1[end,2*j],
                                             node_array2[1,2*j])))
        j == 1 && push!(boundary, (length(cells), 1))
        j == nel_y2 && push!(boundary, (length(cells), 3))
        if nel_x2 == 1
            push!(boundary, (length(cells), 2))
            push!(facesets["right"], (length(cells), 2))
        end
        for i in 1:nel_x2
            push!(cells, QuadraticQuadrilateral((node_array2[2*i,2*j-1],
                                                 node_array2[2*i+2,2*j-1],
                                                 node_array2[2*i+2,2*j+1],
                                                 node_array2[2*i,2*j+1],
                                                 node_array2[2*i+1,2*j-1],
                                                 node_array2[2*i+2,2*j],
                                                 node_array2[2*i+1,2*j+1],
                                                 node_array2[2*i,2*j],
                                                 node_array2[2*i+1,2*j])))
            if i == nel_x2
                push!(boundary, (length(cells), 2))
                push!(facesets["right"], (length(cells), 2))
            end
            j == 1 && push!(boundary, (length(cells), 1))
            j == nel_y2 && push!(boundary, (length(cells), 3))
        end
    end

    push!(nodesets["load"], node_array2[end, midpointindy])

    # Upper left rectangle
    offsetstep = (UR[2] - MR[2])/(nel1[2] - nel2[2])/2
    indexoffset = length(nodes)

    nel_x3 = nel1[1]; nel_y3 = nel1[2] - nel2[2] - 1
    n_nodes_x3 = 2*nel_x3 + 1; n_nodes_y3 = 2*nel_y3 + 2
    n_nodes3 = n_nodes_x3 * n_nodes_y3

    _LL = Vec{2,T}((LL[1], MR[2] + offsetstep))
    _LR = Vec{2,T}((UR[1], MR[2] + offsetstep))
    _UL = Vec{2,T}((LL[1], UR[2]))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x3, n_nodes_y3, _LL, _LR, UR, _UL)

    # Generate cells
    node_array3 = reshape(collect(indexoffset+1:indexoffset+n_nodes3), (n_nodes_x3, n_nodes_y3))

    for i in 1:nel_x3
        push!(cells, QuadraticQuadrilateral((node_array1[2i-1,end], node_array1[2i+1,end],
                                             node_array3[2i+1,2], node_array3[2i-1,2],
                                             node_array1[2i,end], node_array3[2i+1,1],
                                             node_array3[2i,2], node_array3[2i-1,1],
                                             node_array3[2i,1])))
        i == 1 && push!(boundary, (length(cells), 4))
        i == nel_x3 && push!(boundary, (length(cells), 2))
    end
    for j in 1:nel_y3, i in 1:nel_x3
        push!(cells, QuadraticQuadrilateral((node_array3[2i-1,2j],
                                             node_array3[2i+1,2j],
                                             node_array3[2i+1,2j+2],
                                             node_array3[2i-1,2j+2],
                                             node_array3[2i,2j],
                                             node_array3[2i+1,2j+1],
                                             node_array3[2i,2j+2],
                                             node_array3[2i-1,2j+1],
                                             node_array3[2i,2j+1])))

        i == 1 && push!(boundary, (length(cells), 4))
        i == nel_x3 && push!(boundary, (length(cells), 2))
        if j == nel_y3
            push!(boundary, (length(cells), 3))
            push!(facesets["top"], (length(cells), 3))
        end
    end

    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)

    return Grid(cells, nodes, facesets=facesets, nodesets=nodesets,
        boundary_matrix=boundary_matrix)
end

function TieBeamGrid(::Type{Val{CellType}}, ::Type{T}=Float64, refine=1) where {T, CellType}
    if CellType === :Linear
        return _LinearTieBeamGrid(T, refine)
    else
        return _QuadraticTieBeamGrid(T, refine)
    end
end

function _LinearTieBeamGrid(::Type{T}=Float64, refine=1) where {T}
    nodes = Node{2,T}[]
    cells = Quadrilateral[]
    boundary = Tuple{Int,Int}[]
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["leftfixed"]  = Set{Tuple{Int,Int}}()
    facesets["toproller"]    = Set{Tuple{Int,Int}}()
    facesets["rightload"]  = Set{Tuple{Int,Int}}()
    facesets["bottomload"] = Set{Tuple{Int,Int}}()

    # Lower left rectangle
    nel_x1 = 32 * refine; nel_y1 = 3 * refine;
    n_nodes_x1 = nel_x1 + 1; n_nodes_y1 = nel_y1 + 1
    n_nodes1 = n_nodes_x1 * n_nodes_y1

    LL = Vec{2,T}((0, 0))
    LR = Vec{2,T}((T(nel_x1 / refine), T(0)))
    UR = Vec{2,T}((T(nel_x1 / refine), T(nel_y1 / refine)))
    UL = Vec{2,T}((T(0), T(nel_y1 / refine)))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x1, n_nodes_y1, LL, LR, UR, UL)

    node_array1 = reshape(collect(1:n_nodes1), (n_nodes_x1, n_nodes_y1))
    for j in 1:nel_y1, i in 1:nel_x1
        push!(cells, Quadrilateral((node_array1[i,j], node_array1[i+1,j],
                                    node_array1[i+1,j+1],node_array1[i,j+1])))
        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
            push!(facesets["leftfixed"], (cidx, 4))
        end
        if i == nel_x1
            cidx = length(cells)
            push!(boundary, (cidx, 2))
            push!(facesets["rightload"], (cidx, 2))
        end
        if j == 1
            cidx = length(cells)
            push!(boundary, (cidx, 1))
            if i == 31
                push!(facesets["bottomload"], (cidx, 1))
            end
        end
        if j == nel_y1 && i != 31
            cidx = length(cells)
            push!(boundary, (cidx, 3))
        end
    end

    nel_x2 = 1 * refine; nel_y2 = 3 * refine + refine - 1
    n_nodes_x2 = nel_x2 + 1; n_nodes_y2 = nel_y2 + 1
    n_nodes2 = n_nodes_x2 * n_nodes_y2
    indexoffset = length(nodes)
    LL = Vec{2,T}((T(30), (nel_y1 + T(1)) / refine))
    LR = Vec{2,T}((T(31), (nel_y1 + T(1)) / refine))
    UR = Vec{2,T}((T(31), nel_y1 / refine + T(4)))
    UL = Vec{2,T}((T(30), nel_y1 / refine + T(4)))

    Ferrite._generate_2d_nodes!(nodes, n_nodes_x2, n_nodes_y2, LL, LR, UR, UL)
    node_array2 = reshape(collect(indexoffset+1:indexoffset+n_nodes2), (n_nodes_x2, n_nodes_y2))

    t = 30
    for i in 1:refine
        push!(cells, Quadrilateral((node_array1[t*refine+i, nel_y1 + 1],
                                    node_array1[t*refine+i+1, nel_y1 + 1],
                                    node_array2[i+1,1],
                                    node_array2[i,1])))
        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
        end
        if i == refine
            cidx = length(cells)
            push!(boundary, (cidx, 2))
        end
    end

    for j in 1:nel_y2, i in 1:nel_x2
        push!(cells, Quadrilateral((node_array2[i,j], node_array2[i+1,j],
                                    node_array2[i+1,j+1], node_array2[i,j+1])))
        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
        end
        if i == nel_x2
            cidx = length(cells)
            push!(boundary, (cidx, 2))
        end
        if j == nel_y2
            cidx = length(cells)
            push!(boundary, (cidx, 3))
            push!(facesets["toproller"], (cidx, 3))
        end
    end

    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)
    return Grid(cells, nodes, facesets=facesets,
        boundary_matrix=boundary_matrix)
end

function _QuadraticTieBeamGrid(::Type{T}=Float64, refine = 1) where {T}
    nodes = Node{2,T}[]
    cells = QuadraticQuadrilateral[]
    boundary = Tuple{Int,Int}[]
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["leftfixed"]  = Set{Tuple{Int,Int}}()
    facesets["toproller"]    = Set{Tuple{Int,Int}}()
    facesets["rightload"]  = Set{Tuple{Int,Int}}()
    facesets["bottomload"] = Set{Tuple{Int,Int}}()

    # Lower left rectangle
    nel_x1 = 32*refine; nel_y1 = 3*refine;
    n_nodes_x1 = 2*nel_x1 + 1; n_nodes_y1 = 2*nel_y1 + 1
    n_nodes1 = n_nodes_x1 * n_nodes_y1

    LL = Vec{2,T}((0, 0))
    LR = Vec{2,T}((T(nel_x1/refine), T(0)))
    UR = Vec{2,T}((T(nel_x1/refine), T(nel_y1/refine)))
    UL = Vec{2,T}((T(0), T(nel_y1/refine)))
    Ferrite._generate_2d_nodes!(nodes, n_nodes_x1, n_nodes_y1, LL, LR, UR, UL)

    node_array1 = reshape(collect(1:n_nodes1), (n_nodes_x1, n_nodes_y1))
    for j in 1:nel_y1, i in 1:nel_x1
        push!(cells, QuadraticQuadrilateral((node_array1[2*i-1,2*j-1], node_array1[2*i+1,2*j-1],
                                             node_array1[2*i+1,2*j+1],node_array1[2*i-1,2*j+1],
                                             node_array1[2*i,2*j-1],node_array1[2*i+1,2*j],
                                             node_array1[2*i,2*j+1],node_array1[2*i-1,2*j],
                                             node_array1[2*i,2*j])))
        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
            push!(facesets["leftfixed"], (cidx, 4))
        end
        if i == nel_x1
            cidx = length(cells)
            push!(boundary, (cidx, 2))
            push!(facesets["rightload"], (cidx, 2))
        end
        if j == 1
            cidx = length(cells)
            push!(boundary, (cidx, 1))
            if i == 31
                push!(facesets["bottomload"], (cidx, 1))
            end
        end
        if j == nel_y1 && i != 31
            cidx = length(cells)
            push!(boundary, (cidx, 3))
        end
    end

    nel_x2 = 1*refine; nel_y2 = 3*refine + refine - 1
    n_nodes_x2 = 2*nel_x2 + 1; n_nodes_y2 = 2*nel_y2 + 2
    n_nodes2 = n_nodes_x2 * n_nodes_y2
    indexoffset = length(nodes)
    LL = Vec{2,T}((T(30), (nel_y1 + T(0.5)) / refine))
    LR = Vec{2,T}((T(31), (nel_y1 + T(0.5)) / refine))
    UR = Vec{2,T}((T(31), nel_y1/refine + T(4)))
    UL = Vec{2,T}((T(30), nel_y1/refine + T(4)))

    Ferrite._generate_2d_nodes!(nodes, n_nodes_x2, n_nodes_y2, LL, LR, UR, UL)
    node_array2 = reshape(collect(indexoffset+1:indexoffset+n_nodes2), (n_nodes_x2, n_nodes_y2))

    t = 30
    for i in 1:refine
        push!(cells, QuadraticQuadrilateral((node_array1[2*(refine*t+i-1)+1, 2*nel_y1+1],
                                            node_array1[2*(refine*t+i-1)+3, 2*nel_y1+1],
                                            node_array2[1+2i, 2],
                                            node_array2[2i-1, 2],
                                            node_array1[2*(refine*t+i-1)+2, 2*nel_y1+1],
                                            node_array2[1+2i, 1],
                                            node_array2[2i, 2],
                                            node_array2[2i-1, 1],
                                            node_array2[2i, 1])))

        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
        end
        if i == refine
            cidx = length(cells)
            push!(boundary, (cidx, 2))
        end
    end

    for j in 1:nel_y2, i in 1:nel_x2
        push!(cells, QuadraticQuadrilateral((node_array2[2*i-1,2*j], node_array2[2*i+1,2*j],
                                             node_array2[2*i+1,2*j+2],node_array2[2*i-1,2*j+2],
                                             node_array2[2*i,2*j],node_array2[2*i+1,2*j+1],
                                             node_array2[2*i,2*j+2],node_array2[2*i-1,2*j+1],
                                             node_array2[2*i,2*j+1])))
        if i == 1
            cidx = length(cells)
            push!(boundary, (cidx, 4))
        end
        if i == nel_x2
            cidx = length(cells)
            push!(boundary, (cidx, 2))
        end
        if j == nel_y2
            cidx = length(cells)
            push!(boundary, (cidx, 3))
            push!(facesets["toproller"], (cidx, 3))
        end
    end

    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)
    return Grid(cells, nodes, facesets=facesets,
        boundary_matrix=boundary_matrix)
end


text:
using Test, Pkg
using MPI

if haskey(Pkg.installed(), "CuArrays")
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 32

send_mesg = ArrayType{Float64}(undef,N)
recv_mesg = ArrayType{Float64}(undef,N)
recv_mesg_expected = ArrayType{Float64}(undef,N)

fill!(send_mesg, Float64(rank))
fill!(recv_mesg_expected, Float64(src))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

stats = MPI.Waitall!([sreq, rreq])
@test rreq isa MPI.Request
@test sreq isa MPI.Request
@test MPI.Get_source(stats[2]) == src
@test MPI.Get_tag(stats[2]) == src+32
@test recv_mesg == recv_mesg_expected

(done, stats) = MPI.Testall!([sreq, rreq])
@test done

if rank == 0
    MPI.send(send_mesg, dst, rank+32, comm)
    recv_mesg = recv_mesg_expected
elseif rank == size-1
    (recv_mesg, _) = MPI.recv(src, src+32, comm)
else
    (recv_mesg, _) = MPI.recv(src, src+32, comm)
    MPI.send(send_mesg, dst, rank+32, comm)
end
@test recv_mesg == recv_mesg_expected

if rank == 0
    MPI.Send(Float64(rank), dst, rank+32, comm)
    recv_val = Float64(src)
elseif rank == size-1
    (recv_val, _) = MPI.Recv(Float64, src, src+32, comm)
else
    (recv_val, _) = MPI.Recv(Float64, src, src+32, comm)
    MPI.Send(Float64(rank), dst, rank+32, comm)
end
@test recv_val == Float64(src)

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

(inds, stats) = MPI.Waitsome!([sreq, rreq])
req_arr = [sreq,rreq]
for i in inds
    (done, status) = MPI.Test!( req_arr[i] )
    @test done
end

sreq = MPI.Isend(send_mesg, dst, rank+32, comm)
MPI.Cancel!(sreq)
@test sreq.buffer == nothing

GC.gc()

MPI.Finalize()
@test MPI.Finalized()