text:
import Base: getindex, size, ndims, convert
export MorletWavelet, MorseWavelet, wavebases, fstd, tstd, ContinuousWaveletTransform, cwt

#
# Mother wavelets, which are convolved with the signal in frequency space
#
abstract MotherWavelet{T<:Real}

# The analytic Morlet wavelet. See:
# Torrence, C., & Compo, G. P. (1998). A practical guide to wavelet
# analysis. Bulletin of the American Meteorological Society, 79(1),
# 61–78.
#
# freq is the frequency in Hz
# k0 is the wave number
immutable MorletWavelet{T} <: MotherWavelet{T}
    freq::Vector{T}
    k0::T
    fourierfactor::T
end
MorletWavelet{T<:Real}(freq::Vector{T}, k0::Real=5.0) =
    MorletWavelet(freq, convert(T, k0), convert(T, (4π)/(k0 + sqrt(2 + k0^2))))

# Generate daughter wavelet (samples x frequencies)
function wavebases{T}(w::MorletWavelet{T}, n::Int, fs::Real=1)
    df = 2π * fs / n
    normconst = sqrt(df) / sqrt(sqrt(π) * n)
    k0 = w.k0

    bases = Array(T, div(n, 2)+1, length(w.freq))
    @inbounds begin
        for k = 1:length(w.freq)
            scale = 1/(w.freq[k] * w.fourierfactor)
            norm = sqrt(scale) * normconst
            bases[1, k] = zero(T)
            for j = 2:size(bases, 1)
                bases[j, k] = norm * exp(-abs2(scale * df * (j-1) - k0)*0.5)
            end
        end
    end
    bases
end

fourierfactor{T}(w::MorletWavelet{T}) = w.fourierfactor
# std in the frequency domain, in units of Hz
fstd{T}(w::MorletWavelet{T}) =
    [(f * w.fourierfactor) / (sqrt(2) * 2π) for f in w.freq]
# std in the time domain, in units of time
tstd{T}(w::MorletWavelet{T}) =
    [1 / (sqrt(2) * f * w.fourierfactor) for f in w.freq]

# Generalized Morse wavelet (first family only). See:
# Olhede, S. C., & Walden, A. T. (2002). Generalized Morse wavelets.
# IEEE Transactions on Signal Processing, 50(11), 2661–2670.
# doi:10.1109/TSP.2002.804066
# Lilly, J. M., & Olhede, S. C. (2009). Higher-order properties of
# analytic wavelets. IEEE Transactions on Signal Processing, 57(1),
# 146–160. doi:10.1109/TSP.2008.2007607
#
# freq is the peak frequency in Hz
# β and γ are wavelet parameters
immutable MorseWavelet{T} <: MotherWavelet{T}
    freq::Vector{T}
    β::T
    γ::T
end
MorseWavelet{T<:Real}(freq::Vector{T}, β::Real, γ::Real) =
    MorseWavelet{T}(freq, convert(T, β), convert(T, γ))

# Generate daughter wavelet (samples x frequencies)
function wavebases{T}(w::MorseWavelet{T}, n::Int, fs::Real=1)
    γ = w.γ
    β = w.β
    # Peak frequency
    ωs = (β/γ)^(1/γ)
    dω0 = fs / n * ωs
    # Olhede & Walden p. 2663, multiplied by sqrt(ωs/(2*pi*f)) for the
    # frequency with an additional correction by 1/n for the
    # unnormalized inverse FFT
    r = (2β+1)/γ
    α0 = 2^(r/2)*sqrt(γ*ωs*fs)/(sqrt(gamma(r))*n)

    bases = Array(T, div(n, 2)+1, length(w.freq))
    for k = 1:length(w.freq)
        f = w.freq[k]
        # Change in angular frequency with each point of the DFT
        dω = dω0/f
        α = α0*f^(-0.5)
        for j = 1:size(bases, 1)
            ω = dω*(j-1)
            bases[j, k] = α * ω^β * exp(-ω^γ)
        end
    end
    bases
end

# Fourier factor based on peak frequency
fourierfactor{T}(w::MorseWavelet{T}) = 2*pi*(w.β/w.γ)^(-1/w.γ)
function fstd{T}(w::MorseWavelet{T})
    γ = w.γ
    β = w.β
    ff = fourierfactor(w)
    σ = sqrt(2^(-2/γ)*(exp(lgamma((2*β+3)/γ)-
                           lgamma((2*β+1)/γ)) -
                       exp(2*(lgamma((2*β+2)/γ)-
                              lgamma((2*β+1)/γ)))))
    [(f * ff * σ) / (2 * pi) for f in w.freq]
end
gammatil(x) = gamma(x)/2^x
function tstd{T}(w::MorseWavelet{T})
    γ = w.γ
    β = w.β
    ff = fourierfactor(w)
    σ = sqrt((β^2*gammatil((2β-1)/γ)+γ^2*
             gammatil((2β+2*γ-1)/γ)-2β*γ*
             gammatil((2β+γ-1)/γ))/gammatil((2β+1)/γ))
    [(σ) / (f * ff) for f in w.freq]
end

#
# Functions for applying wavelets to data
#
immutable ContinuousWaveletTransform{T,S,P1,P2}
    fftin::Vector{T}
    fftout::Vector{S}
    ifftwork::Vector{S}
    bases::Array{T,2}
    coi::Vector{T}
    p1::P1
    p2::P2
end

function ContinuousWaveletTransform{T}(w::MotherWavelet{T}, nfft::Int, fs::Real=1;
                                       coi::Vector=scale!(tstd(w), 2*fs))
    fftin = Array(T, nfft)
    fftout = zeros(Complex{T}, div(nfft, 2)+1)
    ifftwork = zeros(Complex{T}, nfft)
    bases = wavebases(w, nfft, fs)
    length(coi) == size(bases, 2) || isempty(coi) || error("length of coi must match number of frequencies")
    ContinuousWaveletTransform(fftin, fftout, ifftwork, bases, coi, plan_rfft(fftin), plan_bfft!(ifftwork))
end

function evaluate!{T,S<:AbstractFloat}(out::Array{Complex{S}, 2}, t::ContinuousWaveletTransform{T},
                                       signal::AbstractVector{T})
    @inbounds begin
        fftin = t.fftin
        fftout = t.fftout
        ifftwork = t.ifftwork
        bases = t.bases
        coi = t.coi

        nsignal = length(signal)
        nfft = length(fftin)
        nrfft = length(fftout)

        nsignal <= nfft || error("signal exceeds length of transform plan")
        size(out, 1) == length(signal) || error("first dimension of out must match length of signal")
        size(out, 2) == size(bases, 2) || error("second dimension of out must match number of wavelets")

        # Compute signal mean ignoring NaN
        m = 0.0
        for i = 1:nsignal
            v = signal[i]
            m += isnan(v) ? 0.0 : v
        end
        m = convert(T, m/nsignal)
        # Subtract signal mean
        for i = 1:nsignal
            fftin[i] = signal[i] - m
        end
        fftin[nsignal+1:nfft] = zero(T)

        # Interpolate NaN values in input and save their indices to
        # set output to NaN later
        discard_samples = falses(length(signal))
        oldv = zero(eltype(fftin))
        i = 1
        while i <= nsignal
            v = fftin[i]
            if isnan(v)
                nanstart = i

                i += 1
                while i <= nsignal && isnan(fftin[i])
                    i += 1
                end

                nanlen = (i - nanstart)
                st = i > nsignal ? zero(eltype(fftin)) : (fftin[i] - oldv)/(nanlen+1)
                for j = 1:nanlen
                    fftin[nanstart+j-1] = oldv+st*j
                end
                discard_samples[nanstart:i-1] = true
            else
                oldv = v
                i += 1
            end
        end
        # This simplifies replacement with NaN below
        discard_samples[1] = false
        discard_samples[end] = false
        discard_sample_indices = find(discard_samples)

        # Perform FFT of padded signal
        A_mul_B!(fftout, t.p1, fftin)

        for k = 1:size(bases, 2)
            # Multiply by wavelet
            for j = 1:nrfft
                ifftwork[j] = fftout[j] * bases[j, k]
            end

            # We only compute the real FFT, but we may need the imaginary
            # frequencies for some mother wavelets as well
            offset = nrfft + 1 + isodd(nfft)
            for j = nrfft+1:size(bases, 1)
                ifftwork[j] = conj(fftout[offset-j]) * bases[j, k]
            end

            # Zero remaining frequencies
            ifftwork[size(bases, 1)+1:nfft] = zero(Complex{T})

            # Perform FFT
            A_mul_B!(ifftwork, t.p2, ifftwork)

            # Copy to output array
            copy!(out, nsignal*(k-1)+1, ifftwork, 1, nsignal)

            if !isempty(coi)
                # Set NaNs at edges
                coi_length = ceil(Int, coi[k])
                out[1:min(coi_length, nsignal), k] = NaN
                out[max(nsignal-coi_length+1, 1):end, k] = NaN

                # Set NaNs for gaps
                for i in discard_sample_indices
                    out[i, k] = NaN
                    if !discard_samples[i+1]
                        out[i:min(i+coi_length, nsignal), k] = NaN
                    end
                    if !discard_samples[i-1]
                        out[max(i-coi_length, 1):i, k] = NaN
                    end
                end
            end
        end
    end
    out
end

# Friendly interface to ContinuousWaveletTransform
function cwt{T<:Real}(signal::Vector{T}, w::MotherWavelet, fs::Real=1)
    t = ContinuousWaveletTransform(w, nextfastfft(length(signal)), fs)
    evaluate!(Array(Complex{T}, length(signal), size(t.bases, 2)), t, signal)
end


text:
"""
    CFLabelPropagationFinder() <: AbstractCommunityFinder

Represents a label propagation search over the hypergraph `h` that finds
a partition into communities (subsets).
"""
struct CFLabelPropagationFinder <: AbstractCommunityFinder
    max_iter::Int
    seed::Int
end

"""
    findcommunities(h::Hypergraph, method::CFLabelPropagationFinder)

Implements the label propagation algorithm over a hypergraph `h`.

NOTE
The algorithm works on a single connected component.
An AssertionError is thrown otherwise.

This algorithm generalizes the one proposed for graphs by Raghavan et al.
(Raghavan, U. N., Albert, R., and Kumara, S. Near linear time algorithm to detect
community structures in large-scale networks. Physical review. E, Statistical,
nonlinear, and soft matter physics 76 (2007).)

The proposed algorithm modifies the propagation rule, splitting it into two phases:
hyperedge labeling and vertex labeling.

For more information see `Section 4` in the paper
Alessia Antelmi, Gennaro Cordasco, Bogumił Kamiński, Paweł Prałat,
Vittorio Scarano, Carmine Spagnuolo, Przemyslaw Szufel
*Analyzing, Exploring, and Visualizing Complex Networks via Hypergraphs Using SimpleHypergraphs.jl.*
Journal Internet Mathematics (2020). https://doi.org/10.24166/im.01.2020
"""
function findcommunities(h::Hypergraph, method::CFLabelPropagationFinder)
    @assert length(get_connected_components(h)) == 1

    rng = MersenneTwister(method.seed)
    vlabels = Dict{Int, Int}()
    helabels = Dict{Int, Int}()

    for v in 1:nhv(h)
        push!(vlabels, v=>v)
    end

    vertices = collect(1:nhv(h))
    edges = collect(1:nhe(h))

    stop = false
    iter = 0

    while !stop && iter < method.max_iter
        stop = true
        shuffle!(rng, edges)

        for e in edges
            length(getvertices(h, e)) == 0 && continue
            l = SimpleHypergraphs.compute_edge_label(h, e, vlabels, helabels, rng)
            push!(helabels, e => l)
        end

        shuffle!(rng, vertices)
        for v in vertices
            l = SimpleHypergraphs.compute_vertex_label(h, v, vlabels, helabels, rng)
            if l != vlabels[v]
                stop = false
                push!(vlabels, v => l)
            end
        end

        iter+=1
    end

    np_vertices = unique(values(vlabels))
    np_edges = unique(values(helabels))

    comms_vertices = Dict{Int, Set{Int}}()
    comms_hyperedges = Dict{Int, Set{Int}}()

    for pv in vlabels
        push!(
            get!(comms_vertices, pv[2], Set{Int}()),
            pv[1]
        )
    end

    for pe in helabels
        push!(
            get!(comms_hyperedges, pe[2], Set{Int}()),
            pe[1]
        )
    end

    labels = Array{Int}(undef, nhv(h))
    for i in 1:nhv(h)
        labels[i] = vlabels[i]
    end

    hlabels = Array{Int}(undef, nhe(h))
    for i in 1:nhe(h)
        if length(getvertices(h, i)) == 0
            hlabels[i] = -1
        else
            hlabels[i] = helabels[i]
        end
    end

    (np=collect(values(comms_vertices)), hep=collect(values(comms_hyperedges)), vlabels=labels, helabels=hlabels, iter=iter)
end

"""
    compute_vertex_label(h::Hypergraph, v::Int, vlabels::Dict{Int,Int}, helabels::Dict{Int,Int}, rng::MersenneTwister)

Vertices labeling phase. Computes the label of each vertex according to the most
frequent label among the hyperedges it belongs to.
"""
function compute_vertex_label(h::Hypergraph, v::Int, vlabels::Dict{Int,Int}, helabels::Dict{Int,Int}, rng::MersenneTwister)
    hesᵥ = gethyperedges(h, v)
    vL = Dict{Int,Int}()

    max = 0
    maxL = Set{Int}()

    for e in shuffle!(rng, collect(keys(hesᵥ)))
        l = helabels[e]

        if !haskey(vL, l)
            push!(vL, l=>0)
        end

        push!(
            vL,
            l => vL[l] + (length(getvertices(h, e)))
        )

        if vL[l] == max
            push!(maxL, l)
        elseif vL[l] > max
            max = vL[l]
            maxL = Set{Int}()
            push!(maxL, l)
        end
    end

    if in(vlabels[v], maxL)
        return vlabels[v]
    end

    return first(maxL)
end

"""
    compute_edge_label(h::Hypergraph, e::Int, vlabels::Dict{Int,Int}, helabels::Dict{Int,Int}, rng::MersenneTwister)

Hyperedges labeling phase. Computes the labels of the hyperedges according  to
the  most frequent label among the vertices contained in that hyperedge.
"""
function compute_edge_label(h::Hypergraph, e::Int, vlabels::Dict{Int,Int}, helabels::Dict{Int,Int}, rng::MersenneTwister)
    vₑ = getvertices(h,e)
    eL = Dict{Int,Int}()

    max = 0
    maxL = Set{Int}()

    for v in shuffle!(rng, collect(keys(vₑ)))
        l = vlabels[v]

        if !haskey(eL, l)
            push!(eL, l=>0)
        end

        push!(
            eL,
            l => eL[l]+1
        )

        if eL[l] == max
            push!(maxL, l)
        elseif eL[l] > max
            max = eL[l]
            maxL = Set{Int}()
            push!(maxL, l)
        end
    end

    if haskey(helabels, e) && in(helabels[e], maxL)
        return helabels[e]
    end

    return first(maxL)
end

ext: jl


text:
"""
structs helping managing and storing data
"""
module DataStructs
using ..BasicStructs
using Parameters,  Dictionaries
export WindowControlStruct,AnnotationStruct,DataToScrollDims,valueForMasToSetStruct,SimpleLineTextStruct, CalcDimsStruct,RawDataToDisp,TwoDimRawDat, ThreeDimRawDat, DataToDisp,FullScrollableDat,SingleSliceDat,SimpleLineTextStruct

"""
hold raw Data that can be send to be displayed
"""
abstract type RawDataToDisp end

"""
2 dimensional ata for displaying single slice
struct is mutable becouse in case of the masks data can be changed multiple times and rapidly
"""
@with_kw mutable struct TwoDimRawDat{T} <: RawDataToDisp
   type::Type{T}= UInt8# easy access to type
   name::String=""#associated name
   dat::AbstractArray{T, 2}=ones(type,2,2)# raw pixel data
end#2DimRawDat

"""
3 dimensional data for displaying single slice
struct is mutable becouse in case of the masks data can be changed multiple times and rapidly
"""
@with_kw mutable struct ThreeDimRawDat{T} <: RawDataToDisp
   type::Type{T}= UInt8# easy access to type
   name::String=""#associated name
   dat::AbstractArray{T, 3}=ones(type,2,2,2)# raw voxel data
end#2DimRawDat

"""
given Vector of tuples where first is string and second is RawDataToDisp
it creates dictionary where keys are those strings - names and values are indicies where they are found
"""
function getLocationDict(listt)::Dictionary{String, Int64}
   return Dictionary(map(it->it.name,listt),collect(eachindex(listt)))

end#getLocationDict

"""
hold Data that can be send to be displayed with required metadata
"""
abstract type DataToDisp end

"""
Struct holding line of text with some text metadata
"""
@with_kw struct SimpleLineTextStruct
    text::String = "" #text to be displayed
    fontSize::Int = 110 # size of letters
    extraLineSpace::Int = 1 # if left to 1 we will get standard line spacing if more distance to line below will increase
end#simpleTextStruct

"""
stores additional data about full dimensions of scrollable dat - this is necessery for switching slicing plane orientation efficiently
"""
@with_kw mutable  struct  DataToScrollDims
    imageSize::Tuple{Int64, Int64, Int64} = (1,1,1)#amount of voxels in each dimensions
    voxelSize::Tuple{Float64, Float64, Float64}  = (1.0,1.0,1.0)#physical size of each pixel
    dimensionToScroll::Int= 3 # by which dimension we should scroll so for example if set to 3 one and we have slice number x we will get data A by A[:,:,x] if dimensionToScroll = 2 ->A[:,x,:]...

end#DataToScrollDims

"""
Data that can be displayed and scrolled (so we have multiple slices)
struct is mutable becouse in case of the masks data can be changed multiple times and rapidly
"""
@with_kw mutable struct FullScrollableDat<: DataToDisp
    dataToScrollDims::DataToScrollDims=DataToScrollDims()#stores additional data about full dimensions of scrollable dat - this is necessery for switching slicing plane orientation efficiently
    dimensionToScroll::Int= 3 # by which dimension we should scroll so for example if set to 3 one and we have slice number x we will get data A by A[:,:,x] if dimensionToScroll = 2 ->A[:,x,:]...
    dataToScroll::Vector{ThreeDimRawDat}=[ThreeDimRawDat()] # tuples where first entry is name of image that we given in configuration, and second entry is data that we want to pass
# data to display in form of a list Of tuples where first entry will be used as headtitle for the data that is an value ;second entry -  value is a vector where each entry will be displayed in separate line
    mainTextToDisp::Vector{SimpleLineTextStruct}=[] # data about text that is to be displayed - text that will be displayd for all data (scrolling will not affect it)
    sliceTextToDisp::Vector{Vector{SimpleLineTextStruct}}=[] # text that will be associated with given slice - length of this needs to be the same as  number of slices we want to scroll through
    #all metrics that were not measured and are possible in ResultMetrics struct will have associated value = -1
    segmMetr::ResultMetrics=ResultMetrics() #results of metrics for whole 3d image
    segmMetrs::Vector{ResultMetrics}=[] #results of metrics for each slice - array needs to be of the same size as number of slices in passed data
    nameIndexes::Dictionary{String, Int64}= getLocationDict(dataToScroll)  #gives a way of efficient querying by supplying dictionary where key is a name we are intrested in and a key is index where it is located in our array
    slicesNumber::Int32=1# number of available slices
end #fullScrollableDat

"""
Data for displaying single slice
struct is mutable becouse in case of the masks data can be changed multiple times and rapidly
"""
@with_kw mutable struct SingleSliceDat<: DataToDisp
    listOfDataAndImageNames::Vector{TwoDimRawDat}=[TwoDimRawDat()]   # tuples where first entry is name of image that we given in configuration, and second entry is data that we want to pass
    textToDisp::Vector{SimpleLineTextStruct}=[] # data about text that is to be displayed
    segmMetr::ResultMetrics=ResultMetrics() #results metrics associated with this slice
    nameIndexes::Dictionary{String, Int64}= getLocationDict(listOfDataAndImageNames)  #gives a way of efficient querying by supplying dictionary where key is a name we are intrested in and a key is index where it is located in our array
    sliceNumber::Int=1 # if we want it to be tamporarly  associated with some slice in scrollable data
end #SingleSliceDat

"""
struct holding  data needed for calculating proper mouse position ,
getting proper size for the texture depending on image dimensions
getting into account  proportions of diffrent parts of display
usefull stats for proper text display
"""
@with_kw  struct CalcDimsStruct
    #imageDims = texture dimensions of main image texture
    imageTextureWidth::Int32=Int32(1)
    imageTextureHeight::Int32=Int32(1)
    #exture dimensions of texture for displaying text
    textTexturewidthh::Int32=Int32(1)
    textTextureheightt::Int32=Int32(1)
    textTextureZeros::Matrix{UInt8}= zeros(UInt8,textTexturewidthh, textTextureheightt ) # used in order to refresh the text texture
    #windowDims
    windowWidth::Int64=1
    windowHeight::Int64=1
    #required ratios - set ratios we want to achieve
    fractionOfMainIm::Float32= Float32(1.0) # needed for controlling  how much space we want for text
    heightToWithRatio::Float32= Float32(1.0) # needed to take into account proper height to width ratio, so main texture would not get srtetched unnatrurally
    #some statistics of a window
    avWindWidtForMain::Int32= Int32(round(windowWidth*fractionOfMainIm ))# how many units we have for the main image in window in width - in case of height  main image will have all height it requires
    avWindHeightForMain::Int32 = windowHeight # how much of height is available for main image
    avMainImRatio::Float32= Float32(windowHeight/avWindWidtForMain )# without any change it will give ratio between available height and width of place for main image - later we will compare it to actual required ratio of image
    #basic data for corrections required for display and calculating mouse positions
    correCtedWindowQuadHeight::Int32  = avMainImRatio# corrected main quad height in window coordinate system
    correCtedWindowQuadWidth::Int32  = avWindWidtForMain# corrected main quad width in window coordinate system
    # now we still need ratio of the resulting quad window size after corrections relative to  total window size
    quadToTotalHeightRatio::Float32= correCtedWindowQuadHeight/avWindHeightForMain
    quadToTotalWidthRatio::Float32= correCtedWindowQuadWidth/avWindWidtForMain
    #opengl corrections - used from top bottom left and right - remember that in open gl max height and width is 2!
    widthCorr::Float32=0.0
    heightCorr::Float32=0.0
    #and analogical in window coordinates
    windowWidthCorr::Int32=Int32(round( (widthCorr/2)*windowWidth))
    windowHeightCorr::Int32= Int32(round((heightCorr/2)*windowHeight))

    #calculated vertex positions for main quads
    mainImageQuadVert::Vector{Float32}= [] # vertex positions of main quad
    wordsImageQuadVert::Vector{Float32} =[]# vertex positions of text display quad

    #calculated sizes for glBufferData
    mainQuadVertSize::Int64 = sizeof(mainImageQuadVert)
    wordsQuadVertSize::Int64 = sizeof(wordsImageQuadVert)

end#CalcDimsStruct

"""
simple struct that when passed is giving information about what should be current value we are setting to the mask
"""
@with_kw struct valueForMasToSetStruct
    value::Int64=1# value that will be used to set  pixels where we would interact with mouse
    text::SimpleLineTextStruct = SimpleLineTextStruct(text= "value of mask to set is  $(value)")
end#valueForMasToSetStruct

"""
struct holding data usefull to controll display window
"""
@with_kw struct WindowControlStruct
    letterCode::String="F2" # shorthand to activate predefined window
    #values represented in shader and controlling window of display
    min_shown_white::Int32 =400
    max_shown_black::Int32 =-200
    toIncrease::Bool= false#true if we intend to increse a treshold value
    toDecrease::Bool= false#true if we intend to decrese a treshold value
    upper::Bool= false#true if we intend to modify upper treshold
    lower::Bool= false#true if we intend to modify lower treshold
    maskContributionToChange::Bool= false#true if we intend to modify mask contribution
end#windowControlStruct

"""
struct holding data usefull to controll mouse interaction with image
    - so stroke width etc.
"""
@with_kw struct AnnotationStruct
    strokeWidthChange::Int32=0 #controls how widely we want to increase visible set of points around point clicked
end#windowControlStruct

end# DataStructs


text:
################
# Comment Type #
################

@ghdef mutable struct Comment
    body::Union{String, Nothing}
    path::Union{String, Nothing}
    diff_hunk::Union{String, Nothing}
    original_commit_id::Union{String, Nothing}
    commit_id::Union{String, Nothing}
    id::Union{Int, Nothing}
    original_position::Union{Int, Nothing}
    position::Union{Int, Nothing}
    line::Union{Int, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    issue_url::Union{HTTP.URI, Nothing}
    pull_request_url::Union{HTTP.URI, Nothing}
    user::Union{Owner, Nothing}
end

Comment(id::Real) = Comment(Dict("id" => id))

namefield(comment::Comment) = comment.id

kind_err_str(kind) = ("Error building comment request: :$kind is not a valid kind of comment.\n"*
                      "The only valid comment kinds are: :issue, :review, :commit")
###############
# API Methods #
###############

@api_default function comment(api::GitHubAPI, repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_get_json(api, path; options...))
end

@api_default function comments(api::GitHubAPI, repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/$(name(item))/comments"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/$(name(item))/comments"
    elseif kind == :commit
        path = "/repos/$(name(repo))/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end
    results, page_data = gh_get_paged_json(api, path; options...)
    return map(Comment, results), page_data
end

@api_default function create_comment(api::GitHubAPI, repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/$(name(item))/comments"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/$(name(item))/comments"
    elseif kind == :commit
        path = "/repos/$(name(repo))/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_post_json(api, path; options...))
end

@api_default function edit_comment(api::GitHubAPI, repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_patch_json(api, path; options...))
end

@api_default function delete_comment(api::GitHubAPI, repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return gh_delete(api, path; options...)
end


text:
using HPAT

@acc hpat function cumsum_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    df1 = aggregate(df, :id, :sx=sum(:x), :sy=sum(:y))
    return sum(df1[:sx])
end

cumsum_test(ENV["HOME"]*"/tmp/udf.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
using Mimi
using MimiFUND

m = MimiFUND.get_model()
run(m)

explore(m)


text:
module NCRegressionTests

using Test
import OrderedCollections
import NCDatasets
const NC = NCDatasets
import PrettyTables

function get_nc_data(ds, var::String)
    if haskey(ds, var)
        return ds[var]
    else
        for key in keys(ds.group)
            if haskey(ds.group[key], var)
                return ds.group[key][var]
            end
        end
    end
    error("No key $var for mse computation.")
    return nothing
end

"""
    compute_mse(;
        job_name::String,
        reference_mse::OrderedCollections.OrderedDict,
        ds_filename_computed::String,
        ds_filename_reference::String,
        varname::Function = s -> s,
        debug_print = false,
        compare_mse = default_compare_mse,
    )

Returns a `Dict` of mean-squared errors between
`NCDataset`s `ds_filename_computed` and
`ds_filename_reference` for all keys in `reference_mse`.
Keys in `reference_mse` may directly map to keys in
the `NCDataset`s, or they may be mapped to the keys
via `varname`.
"""
function compute_mse(;
    job_name::String,
    reference_mse::OrderedCollections.OrderedDict,
    ds_filename_computed::String,
    ds_filename_reference::String,
    varname::Function = s -> s,
    debug_print = false,
    compare_mse = default_compare_mse,
)
    @assert isfile(ds_filename_computed)
    @assert isfile(ds_filename_reference)

    computed_mse = NC.Dataset(ds_filename_computed, "r") do ds_computed
        NC.Dataset(ds_filename_reference, "r") do ds_reference
            _compute_mse(job_name, reference_mse, ds_computed, ds_reference, varname, debug_print, compare_mse)
        end
    end
end

function _compute_mse(job_name, reference_mse, ds_computed, ds_reference, varname, debug_print, compare_mse)

    mse = OrderedCollections.OrderedDict()
    # Ensure z_tcc and fields are consistent lengths:
    best_keys = keys(reference_mse)
    n_keys = length(best_keys)
    variables = Vector{String}(undef, n_keys)
    variables .= varname.(best_keys)
    computed_mse = zeros(n_keys)
    table_reference_mse = zeros(n_keys)
    mse_reductions = zeros(n_keys)
    data_scales_now = zeros(n_keys)
    data_scales_ref = zeros(n_keys)

    for (i, key) in enumerate(best_keys)
        str_key = varname(key)
        if debug_print
            @info "Computing mse for `$key`, `$str_key`"
        end
        data_computed_arr = vec(get_nc_data(ds_computed, str_key))
        data_reference_arr = vec(get_nc_data(ds_reference, str_key))

        # Interpolate data
        # Compute data scale
        data_scale_now = sum(abs.(data_computed_arr)) / length(data_computed_arr)
        data_scale_ref = sum(abs.(data_reference_arr)) / length(data_reference_arr)
        data_scales_now[i] = data_scale_now
        data_scales_ref[i] = data_scale_ref

        # Compute mean squared error (mse)
        mse_single_var = sum((data_computed_arr .- data_reference_arr) .^ 2)
        scaled_mse = mse_single_var / data_scale_ref^2 # Normalize by data scale

        mse[key] = scaled_mse
        mse_reductions[i] = (reference_mse[key] - mse[key]) / reference_mse[key] * 100
        table_reference_mse[i] = reference_mse[key]
        computed_mse[i] = mse[key]
    end

    # Tabulate output
    header = (
        ["Variable", "Data scale", "Data scale", "MSE", "MSE", "MSE"],
        ["", "Computed", "Reference", "Computed", "Reference", "Reduction (%)"],
    )

    table_data = hcat(variables, data_scales_now, data_scales_ref, computed_mse, table_reference_mse, mse_reductions)

    hl_worsened_mse = PrettyTables.Highlighter(
        (data, i, j) -> !compare_mse(data[i, 4], data[i, 5]) && j == 4,
        PrettyTables.crayon"red bold",
    )
    hl_worsened_mse_reduction = PrettyTables.Highlighter(
        (data, i, j) -> !compare_mse(data[i, 4], data[i, 5]) && j == 6,
        PrettyTables.crayon"red bold",
    )
    hl_improved_mse = PrettyTables.Highlighter(
        (data, i, j) -> compare_mse(data[i, 4], data[i, 5]) && j == 6,
        PrettyTables.crayon"green bold",
    )
    @info "Regression tables for `$job_name`"
    PrettyTables.pretty_table(
        table_data;
        header,
        formatters = PrettyTables.ft_printf("%.16e", 4:5),
        header_crayon = PrettyTables.crayon"yellow bold",
        subheader_crayon = PrettyTables.crayon"green bold",
        highlighters = (hl_worsened_mse, hl_improved_mse, hl_worsened_mse_reduction),
        crop = :none,
    )

    return mse
end

default_compare_mse(computed_mse::Real, reference_mse::Real) =
    reference_mse - sqrt(eps()) ≤ computed_mse ≤ reference_mse + sqrt(eps())

# Sometimes we don't have variables to compare against, so
# we have default behavior for handling these cases.
# Users may want to overload these methods for their own use-case.
default_compare_mse(computed_mse::String, reference_mse::String) = true
default_compare_mse(computed_mse::Real, reference_mse::String) = true
default_compare_mse(computed_mse::String, reference_mse::Real) = true

function test_mse(computed_mse, reference_mse, compare_mse = default_compare_mse)
    @testset "Regression tests" begin
        for key in keys(reference_mse)
            mse_not_regressed = compare_mse(computed_mse[key], reference_mse[key])
            if !mse_not_regressed
                @info "Regression failed for `$key`. (computed | best) mse: `$(computed_mse[key])` | `$(reference_mse[key])`"
            end
            @test mse_not_regressed
        end
    end
end

end # module


text:
#
# Define a global problem and ist constructors
#
# ---
import Random: randperm
export getGradient, getCost, getProximalMap
export Problem, HessianProblem

"""
    Problem
Specify properties (values) and related functions for computing
a certain optimization problem.
"""
abstract type Problem end
#
# 1) Function defaults / Fallbacks
#
"""
    getCost(p,x)

evaluate the cost function `F` stored within a [`Problem`](@ref) at the [`MPoint`](@ref) `x`.
"""
function getCost(p::P,x::MP) where {P <: Problem, MP <: MPoint}
  return p.costFunction(x)
end
getGradient(p::Pr,x::P) where {Pr <: Problem, P <: MPoint} =
    throw(ErrorException("no gradient found in $(typeof(p)) to evaluate for a $(typeof(x))."))
getProximalMap(p::Pr,λ,x::P,i) where {Pr <: Problem, P <: MPoint} =
    throw(ErrorException("No proximal map No. $(i) found in $(typeof(p)) to evaluate for $(typeof(x)) with $(typeof(λ))."))
getSubGradient(p::Pr,x::P) where {Pr <: Problem, P <: MPoint} =
        throw(ErrorException("no subgradient found in $(typeof(p)) to evaluate for a $(typeof(x))."))

"""
    HessianProblem <: Problem
For now this is just a dummy problem to carry information about a Problem also providing a Hessian
"""
mutable struct HessianProblem{mT <: Manifold} <: Problem
    M::mT
    costFunction::Function
    Heassian::Function
end


text:
# # Example gallery
#
# Semi-curated collection of examples.
#
# ## Lines and markers
#
# ### A simple scatter plot

using SplitApplyPlot, CairoMakie

df = (x=rand(100), y=rand(100))
data(df) * mapping(:x, :y) |> plot
AbstractPlotting.save("simplescatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](simplescatter.svg)
#
# ### A simple lines plot

x = range(-π, π, length=100)
y = sin.(x)
df = (; x, y)
data(df) * mapping(:x, :y) * visual(Lines) |> plot
AbstractPlotting.save("simplelines.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](simplelines.svg)
#
# ### Lines and scatter combined plot

x = range(-π, π, length=100)
y = sin.(x)
df = (; x, y)
data(df) * mapping(:x, :y) * (visual(Scatter) + visual(Lines)) |> plot
AbstractPlotting.save("simplescatterlines1.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](simplescatterlines1.svg)
#
x = range(-π, π, length=100)
y = sin.(x)
df1 = (; x, y)
df2 = (x=rand(10), y=rand(10))
m = mapping(:x, :y)
geoms = data(df1) * visual(Lines) + data(df2) * visual(Scatter)
plot(m * geoms)
AbstractPlotting.save("simplescatterlines2.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](simplescatterlines2.svg)
#
# ### Linear regression on a scatter plot

df = (x=rand(100), y=rand(100), z=rand(100))
m = data(df) * mapping(:x, :y)
geoms = linear() + visual(Scatter) * mapping(color=:z)
plot(m * geoms)
AbstractPlotting.save("linefit.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](linefit.svg)
#
# ### Overload default geometry

df = (x=rand(100), y=rand(100), z=rand(100))
m = data(df) * mapping(:x, :y)
geoms = visual(Scatter) * linear() + visual(Scatter) * mapping(color=:z)
plot(m * geoms)
AbstractPlotting.save("linefitscatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](linefitscatter.svg)
#
# ## Faceting
#
# The "facet style" is only applied with an explicit call to `facet!`.
#
# ### Facet grid

df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
data(df) * mapping(:x, :y, col=:i, row=:j) |> plot |> facet!
AbstractPlotting.save("facetscatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](facetscatter.svg)
#
# ### Facet wrap

df = (x=rand(100), y=rand(100), l=rand(["a", "b", "c", "d", "e"], 100))
data(df) * mapping(:x, :y, layout=:l) |> plot |> facet!
AbstractPlotting.save("facetwrapscatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](facetwrapscatter.svg)
#

# ### Embedding facets
#
# All SplitApplyPlot plots can be inserted in any figure position, where the rest
# of the figure is managed by vanilla Makie.
# For example

df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
resolution = (1200, 600)
fig = Figure(; resolution)
ax = Axis(fig[1, 1], title="Some plot")
layer = data(df) * mapping(:x, :y, col=:i, row=:j)
subfig = fig[1, 2:3]
ag = plot!(subfig, layer)
facet!(subfig, ag)
for ae in ag
    Axis(ae).xticklabelrotation[] = π/2
end
display(fig)
AbstractPlotting.save("nestedfacet.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](nestedfacet.svg)

# ### Adding traces to only some subplots

df1 = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
df2 = (x=[0, 1], y=[0.5, 0.5], i=fill("a", 2), j=fill("e", 2))
m = mapping(:x, :y, col=:i, row=:j)
geoms = data(df1) * visual(Scatter) + data(df2) * visual(Lines)
m * geoms |> plot |> facet!
AbstractPlotting.save("facetscatterlines.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](facetscatterlines.svg)
#
# ## Statistical analyses
#
# ### Density plot

df = (x=randn(1000), c=rand(["a", "b"], 1000))
data(df) * mapping(:x, color=:c) * SplitApplyPlot.density(bandwidth=0.5) |> plot
AbstractPlotting.save("density.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](density.svg)
#
# Using the recipe from AbstractPlotting also works (let us try to figure out whether we need an analysis or not).

df = (x=randn(1000), c=rand(["a", "b"], 1000))
data(df) * mapping(:x, col=:c) * visual(AbstractPlotting.Density) |> plot |> facet!
AbstractPlotting.save("densityvisual.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](densityvisual.svg)
#

df = (x=randn(1000), c=rand(["a", "b"], 1000))
layer = data(df) * mapping(:x, color=:c) * SplitApplyPlot.density(bandwidth=0.5) *
    visual(orientation=:vertical)
"Not yet supported" # hide

# ## Discrete scales
#
# By default categorical ticks, as well as names from legend entries, are taken from the
# value of the variable converted to a string. Scales can be equipped with labels to
# overwrite that

df = (x=rand(["a", "b", "c"], 100), y=rand(100))
data(df) * mapping(:x, :y) * visual(BoxPlot) |> plot
AbstractPlotting.save("boxplot.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](boxplot.svg)
#

df = (x=rand(["a", "b", "c"], 100), y=rand(100))
layer = data(df) *
    mapping(
        :x => renamer("a" => "label1", "b" => "label2", "c" => "label3"),
        :y
    ) * visual(BoxPlot)
plot(layer)
AbstractPlotting.save("relabel.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](relabel.svg)
#
# The order can also be changed by tweaking the scale
layer = data(df) *
    mapping(
        :x => renamer("b" => "label b", "a" => "label a", "c" => "label c"),
        :y
    ) * visual(BoxPlot)
plot(layer)
AbstractPlotting.save("reorder.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](reorder.svg)
#
# ## Continuous scales

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
layer = data(df) *
    mapping(
        :x,
        :y => log => "√x + 20x + 100 (log scale)",
    ) * visual(Lines)
plot(layer)
AbstractPlotting.save("logscale1.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](logscale1.svg)
#

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
layer = data(df) *
    mapping(
        :x,
        :y => "√x + 20x + 100",
    ) * visual(Lines)
plot(layer, axis=(yscale=log,))
AbstractPlotting.save("logscale2.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](logscale2.svg)
#
# ## Custom scales
#
# Sometimes, there is no default palettes for a specific attribute. In that
# case, the user can pass their own.

using Colors
x=repeat(1:20, inner=20)
y=repeat(1:20, outer=20)
u=cos.(x)
v=sin.(y)
c=rand(Bool, length(x))
d=rand(Bool, length(x))
df = (; x, y, u, v, c, d)
colors = [colorant"#E24A33", colorant"#348ABD"]
heads = ['▲', '●']
layer = data(df) *
    mapping(:x, :y, :u, :v) *
    mapping(arrowhead=:c => nonnumeric) *
    mapping(arrowcolor=:d => nonnumeric) *
    visual(Arrows, arrowsize=10, lengthscale=0.3)
plot(layer; palettes=(arrowcolor=colors, arrowhead=heads))
AbstractPlotting.save("arrows.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](arrows.svg)
#
# ## Axis and figure keywords
#
# ### Axis tweaking
#
# To tweak one or more axes, simply use the `axis` keyword when plotting. For example
#
df = (x=rand(100), y=rand(100), z=rand(100))
m = data(df) * mapping(:x, :y)
geoms = linear() + visual(Scatter) * mapping(color=:z)
plot(m * geoms, axis=(aspect=1,))
AbstractPlotting.save("axis.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](axis.svg)
#
# ### Figure tweaking
#
df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
m = data(df) * mapping(:x, :y, layout=:c)
geoms = linear() + visual(Scatter) * mapping(color=:z)
fg = plot(m * geoms, axis=(aspect=1,), figure=(resolution=(1200, 600),))
facet!(fg)
AbstractPlotting.save("figure.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](figure.svg)
#
# ## Multiple selection
#
# Selecting multiple columns at once can have two possible applications. One is
# "wide data", the other is on-the-fly creating of novel columns.
#
# ### Wide data
#

df = (a=randn(100), b=randn(100), c=randn(100))
m = data(df) * mapping((:a, :b, :c) .=> "some label") * mapping(color=dims(1))
plot(m * SplitApplyPlot.density())
AbstractPlotting.save("widedensity.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](widedensity.svg)
#

df = (a=rand(100), b=rand(100), c=rand(100), d=rand(100))
m = data(df) * mapping(1, 2:4, color=dims(1))
geoms = linear() + visual(Scatter)
fg = plot(m * geoms)
AbstractPlotting.save("widescatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](widescatter.svg)
#
# The wide format is combined with broadcast semantics.

df = (sepal_length=rand(100), sepal_width=rand(100), petal_length=rand(100), petal_width=rand(100))
xvars = ["sepal_length", "sepal_width"]
yvars = ["petal_length" "petal_width"]
m = data(df) * mapping(
    xvars .=> "sepal",
    yvars .=> "petal",
    row=dims(1) => c -> split(xvars[c], '_')[2],
    col=dims(2) => c -> split(yvars[c], '_')[2],
)
geoms = linear() + visual(Scatter)
facet!(plot(m * geoms))
AbstractPlotting.save("widefacet.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](widefacet.svg)
#
# ### New columns on the fly
#
df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
m = data(df) * mapping(:x, (:x, :y, :z) => (+) => "x + y + z", layout=:c)
geoms = linear() + visual(Scatter) * mapping(color=:z)
fg = plot(m * geoms)
facet!(fg)
AbstractPlotting.save("columnonthefly.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](columnonthefly.svg)
#


text:
using PyPlot

"""
    y = ℓ(c, j, x)

Computes `j`th Lagrange polynomial with respect to the centres `c[k]`.
Thus, if `1 ≤ k ≤ n` then `ℓ` has degree `n-1` and `1 ≤ j ≤ n`.

"""
function ℓ(c, j, x)
    n = length(c)
    @assert 1 ≤ j ≤ n
    y = 1.0
    for k = 1:j-1
        y *= ( x - c[k] ) / ( c[j] - c[k] )
    end
    for k = j+1:n
        y *= ( x - c[k] ) / ( c[j] - c[k] )
    end
    return y
end

function ℓ(c, j, x::AbstractArray)
    n = length(x)
    y = fill(1.0, n)
    for k = 1:n
        y[k] = ℓ(c, j, x[k])
    end
    return y
end

n = 6
x = range(-1; stop=1, length=201)
c = range(-1; stop=1, length=n)

figure(1)
for j = 1:n
    plot(x, ℓ(c, j, x))
end
for j = 1:n
    plot(c[j:j], [1.0], "ok", c[j:j], [0.0], "ok")
end
grid(true)
title("Lagrange polynomials of degree $(n-1)")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
# functions related to mesh warping/coordinate updating

"""
  This function updates the coordinates of the vertices of a given element.
  In parallel, it is the responsibility of the caller to update the coordinates
  consistenly for shared elements.  Note that the user must not change
  the topology of the mesh (ie. the adjacency of mesh entities), only
  their positions.

  **Inputs**

   * elnum: the element number to update
   * coords_new: the `mesh.dim` x `coord_numNodesPerElement` array of new
            coordinates.
   * snap: if true, snap the coordinates to the geometry (if the geometric
           model supports snapping).  Setting this to false is strongly
           discouraged and can cause inaccurate coordinate and derivative
           evaluations in other functions.  If true, coords_new will be
           overwritten with the snapped coordinates

  For `coords_new`, each column contains the coordinate of a coordinate node,
  starting
  with the vertices, followed by the edges (if quadratic).
  The entities must be ordered the same as the downward adjacencies of
  the elements.  The mesh.topo field combined with mesh.element_vertnums
  can be used to figure this out.

  Note that even if `coords_new` is complex valued, this function only uses
  the real part.  If you need to propigate complex number through the
  metric calculation, see [`recalcCoordinatesAndMetrics`](@ref).

  After completeing all calls to update_coords, users *must* call
  commit_coords().

  This function is supported for backwards compatability, see the other
  methods of `update_coords` which work on the entire mesh at once, which may
  be more efficient.
"""
function update_coords(mesh::PumiMesh, elnum::Integer,  coords_new::AbstractMatrix, snap::Bool=true)

#  @assert size(coords_new, 1) == mesh.dim
#  @assert (size(coords_new, 2) == mesh.dim + 1 || size(coords_new, 2) == mesh.dim+1 + mesh.numTypePerElement[2]) # number of verts or number of verts + edges

  verts = Array{Ptr{Void}}(12)
  coords_j = zeros(Float64, 3)
#  newx_j = zeros(Float64, 3)
#  xi_j = zeros(Float64, 2)
#  g = apf.getModel(mesh.m_ptr)
#  can_eval = gmi.can_eval(g)
#  snap = snap && can_eval

  el_i = mesh.elements[elnum]
  apf.getDownward(mesh.m_ptr, el_i, 0, verts)

  if mesh.coord_order == 1
    @assert (size(coords_new, 2) >= mesh.numTypePerElement[1])
  elseif mesh.coord_order == 2
    @assert (size(coords_new, 2) == (mesh.numTypePerElement[1] + mesh.numTypePerElement[2]))
  else
    throw(ErrorException("Unsupported coordinate order"))
  end

  #TODO: make this a proper master loop construct
  for j=1:(mesh.dim + 1)  # loop over verts
    for k=1:mesh.dim
      coords_j[k] = real(coords_new[k, j])
    end
    setCoords(mesh, verts[j], 0, coords_j, snap)
    for k=1:mesh.dim
      coords_new[k, j] = coords_j[k]
    end
  end

  if apf.hasNodesIn(mesh.coordshape_ptr, 1)
    offset = mesh.dim + 1
    nedges = apf.getDownward(mesh.m_ptr, el_i, 1, verts)
    for j=1:nedges
      for k=1:mesh.dim
        coords_j[k] = real(coords_new[k, j + offset])
      end

      setCoords(mesh, verts[j], 0, coords_j, snap)
      for k=1:mesh.dim
        coords_new[k, j + offset] = coords_j[k]
      end
    end
  end

  return nothing
end

#------------------------------------------------------------------------------
# Internal API for updating xyz and xi coordinates of individual nodes

"""
  Set the parametric coordinates for a node.  Users should *not* call
  apf.setParam directly, use this function instead.  This function
  updates both the parametric and xyz coordinates.

  Note that this function stores the parametric coordinates as defined by
  the CAD engine, not the hybrid xyz-xi format used by `interface_geo.jl`.
  In particular, for 2D meshes, the xi coordinates for a mesh entity classified
  on a model region should be the CAD parametric coordinates, not the xyz.

  **Inputs**

   * mesh: mesh object
   * entity: apf::MeshEntity*
   * node: local node number (0 based)
   * xi: vector of length 3 containing the parametric coordinates.

  **Inputs/Outputs**

   * coords: if provided, will be overwritten with the xyz coordinates
             of the point
"""
function setCoordsXi(mesh::PumiMesh, entity::Ptr{Void}, node::Integer,
                     xi::AbstractVector{Float64},
                     coords::AbstractArray{Float64}=NullArray{Float64}(0))

  g = apf.getModel(mesh.m_ptr)
  me = apf.toModel(mesh.m_ptr, entity)
  me_dim = apf.getModelType(mesh.m_ptr, me)

  @assert me_dim < 3  # there are no parametric coordinates for regions
  @assert me_dim > 0  # vertex parametric coordinates cannot be modified
  @assert mesh.geoNums.can_eval

  # check xi range
  xi_range = Array{Float64}(2)
  for i=1:me_dim
    gmi.range(g, me, i-1, xi_range)
    @assert xi[i] >= xi_range[1]
    @assert xi[i] <= xi_range[2]
  end

  e_dim = apf.getDimension(mesh.m_ptr, entity)
  if e_dim == 0
    apf.setParam(mesh.m_ptr, entity, xi)
  else  # higher dimension xi coordinates are stored in a field
    apf.setComponents(mesh.geoNums.xi_fptr, entity, node, xi)
  end

  coords2 = Array{Float64}(3)  #TODO: cache this
  gmi.geval(g, me, xi, coords2)
  apf.setPoint(mesh.m_ptr, entity, node, coords2)

  if length(coords) > 0
    for i=1:3
      coords[i] = coords2[i]
    end
  end

  return nothing
end

"""
  Retrieve the parametric coordinates for a mesh node, and optionally the
  xyz coordinates as well.  This works for *any* MeshEntity, not only
  vertices like `apf.getParam`.  This function should be preferred in
  all cases over `apf.getParam`.

  **Inputs**

   * mesh: mesh object
   * entity: apf::MeshEntity*
   * node: local node number (0 based)

  **Inputs/Outputs**

   * xi: vector of length 3 to be overwritten with the parametric coordinates.
   * coords: if provided, will be overwritten with the xyz coordinates
             of the point
"""
function getCoordsXi(mesh::PumiMesh, entity::Ptr{Void}, node::Integer,
                     xi::AbstractVector{Float64},
                     coords::AbstractArray{Float64}=NullArray{Float64}(0))

  g = apf.getModel(mesh.m_ptr)
  me = apf.toModel(mesh.m_ptr, entity)
  me_dim = apf.getModelType(mesh.m_ptr, me)

  @assert me_dim < 3  # there are no parametric coordinates for regions
  #@assert me_dim > 0  # vertex parametric coordinates cannot be modified
  @assert mesh.geoNums.can_eval

  e_dim = apf.getDimension(mesh.m_ptr, entity)
  if e_dim == 0
    apf.getParam(mesh.m_ptr, entity, xi)
  else  # higher dimension xi coordinates are stored in a field
    apf.getComponents(mesh.geoNums.xi_fptr, entity, node, xi)
  end

  if length(coords) > 0
    apf.getPoint(mesh.m_ptr, entity, node, coords)
  end

  return nothing
end

"""
  This function updates the xyz coordinates for a given node, and also
  the xi coordinates, which are computed from the xyz coordinates.

  **Inputs**

   * mesh
   * node: local node number (0 based)
   * coords: vector of length 3 containing the xyz coordinates
   * snap: Bool, if true, the coordinates will be snapped onto the geometry
           and the `coords` overwritten with the new coordinates, default
           true.  Setting this to false can result in inaccurate derivative
           calculations.  If the geometric model does not support evaluation
           points, this option has no effect.
"""
function setCoords(mesh::PumiMesh, entity::Ptr{Void}, node::Integer,
                   coords::AbstractArray{Float64}, snap::Bool=true)

  geonums = mesh.geoNums
  newx = Array{Float64}(3)
  xi = Array{Float64}(3)

  me = apf.toModel(mesh.m_ptr, entity)
  me_dim = apf.getModelType(mesh.m_ptr, me)

  if snap && geonums.can_eval
    # update xyz and xi coordinates consistently
    has_xi = getSnappedCoords(mesh, entity, true, coords, newx, xi)
    if has_xi
      if me_dim > 0
        setCoordsXi(mesh, entity, node, xi, newx)
      end
      for i=1:3
        coords[i] = newx[i]
      end
    else
      apf.setPoint(mesh.m_ptr, entity, node, coords)
    end

  elseif !snap && geonums.can_eval
    # update the xyz and xi coordinates, but don't make them consistent
    has_xi = getSnappedCoords(mesh, entity, false, coords, newx, xi)
    if has_xi && me_dim > 0
      setCoordsXi(mesh, entity, node, xi)
    end
    apf.setPoint(mesh.m_ptr, entity, node, coords)  # overwrite coordinates from
                                               # setCoordsXi
  else  # can't eval
    apf.setPoint(mesh.m_ptr, entity, node, coords)
  end

  return nothing
end

"""
  Retrieve the xyz coordinates of a node

  **Inputs**

   * mesh
   * entity: the MeshEntity
   * node: the local node number (0 based)

  **Inputs/Outputs**

   * coords: vector of length 3 to be overwritten with coordinates
"""
function getCoords(mesh::PumiMesh, entity::Ptr{Void}, node::Integer,
                   coords::AbstractArray{Float64})

  apf.getPoint(mesh.m_ptr, entity, node, coords)

  return nothing
end

"""
  Internal function to snap coordinates to the geometry.  Only call this
  function if the geometric model supports evaluating coordinates.

  **Inputs**

   * mesh
   * entity: the MeshEntity
   * snap: if the xyz coordinates should be snapped
   * coords: vector of length 3 containing xyz coordinates

  **Inputs/Outputs**

   * newx: new xyz coordinates, snapped to the geometry if `snap` is true,
           same as `coords` if false
   * xi: the CAD parametric coordinates

  **Outputs**

   * has_xi: Bool, if true, the model entity that `entity` is classified on
             has parametric coordinates.  If false, the value of `xi` on exit
             is undefined.
"""
function getSnappedCoords(mesh::PumiMesh, entity::Ptr{Void}, snap::Bool,
                          coords::AbstractVector, newx::AbstractVector,
                          xi::AbstractVector)

  g = apf.getModel(mesh.m_ptr)
  me = apf.toModel(mesh.m_ptr, entity)
  me_dim = apf.getModelType(mesh.m_ptr, me)
  has_xi = me_dim != 3  # only snap if the parametric coordinates exist

  if has_xi
    if me_dim == 0  # Simmetrix doesn't have closest_point for vertices
      fill!(xi, 0)
      gmi.geval(g, me, xi, newx)
    else
      gmi.closest_point(g, me, coords, newx, xi)
    end
  end

  if !snap || !has_xi
    for k=1:mesh.dim
      newx[k] = coords[k]
    end
  end

  return has_xi
end

"""
  This must be called after all calls to update_coords are complete. It
  update all fields of the mesh that are derived from the coordinates
  (dxidx, jac, node coordinate etc.), given that update_coords does not
  change the mesh topology.

  The verify keyword argument determines whether or not to run the Pumi
  verifier on the updated mesh.  It can be expensive, so it is recommended
  to only run it on small meshes.

  The `write_vis` keyword controls whether a vtk file, named `pre_commit` is
  written just before the metrics are recalculated (if an element is turned
  inside out, the metric recalculation is where the error gets thrown).
"""
function commit_coords(mesh::PumiMesh, sbp, opts; verify=true, write_vis=false)
# users must call this function when they have finished updating the coordinates
# the verify kwarg determines if the Pumi verifier on the new mesh
# that should check for negative volumes

  apf.acceptChanges(mesh.m_ptr)
  if verify
    apf.Verify(mesh.m_ptr)
    # TODO; call our verification code too
  end

  if write_vis
    writeVisFiles(mesh, "pre_commit")
  end

  getAllCoordinatesAndMetrics(mesh, sbp, opts, verify=verify)

end


text:
using NamedTupleTools: namedtuple

ntkeys(::Type{NamedTuple{K,V}}) where {K, V} = K
ntvaltype(::Type{NamedTuple{K,V}}) where {K, V} = V

export schema

"""
    schema(::Type)

`schema` turns a type into a value that's easier to work with.

Example:

    julia> nt = (a=(b=[1,2],c=(d=[3,4],e=[5,6])),f=[7,8]);

    julia> NT = typeof(nt)
    NamedTuple{(:a, :f),Tuple{NamedTuple{(:b, :c),Tuple{Array{Int64,1},NamedTuple{(:d, :e),Tuple{Array{Int64,1},Array{Int64,1}}}}},Array{Int64,1}}}

    julia> schema(NT)
    (a = (b = Array{Int64,1}, c = (d = Array{Int64,1}, e = Array{Int64,1})), f = Array{Int64,1})
"""
function schema end

schema(::NamedTuple{(), Tuple{}}) = NamedTuple()
schema(::Type{NamedTuple{(), Tuple{}}}) = NamedTuple()

function schema(NT::Type{NamedTuple{names, T}}) where {names, T}
    return namedtuple(ntkeys(NT), schema(ntvaltype(NT)))
end

function schema(TT::Type{T}) where {T <: Tuple}
    return schema.(Tuple(TT.types))
end

schema(t::T) where {T <: Tuple} = schema(T)

schema(t::T) where {T <: NamedTuple} = schema(T)

schema(T) = T


text:
export MultiMPIFile

mutable struct MultiMPIFile <: MPIFile
  files::Vector{MPIFile}

  function MultiMPIFile(filenames::Vector{String})
    return new([MPIFile(f) for f in filenames])
  end

  function MultiMPIFile(f::Vector{T}) where T<:MPIFile
    return new(f)
  end

end

getindex(f::MultiMPIFile, index::Integer) = f.files[index]

length(f::MultiMPIFile) = length(f.files)

start_(f::MultiMPIFile) = 1
next_(f::MultiMPIFile,state) = (f[state],state+1)
done_(f::MultiMPIFile,state) = state > length(f.files)
iterate(f::MultiMPIFile, s=start_(f)) = done_(f, s) ? nothing : next_(f, s)

function show(io::IO, f::MultiMPIFile)
  print(io, "Multi MPI File:\n")
  for (i,file) in enumerate(f.files)
    print(io,"$i. ", file)
  end
end

acqNumPeriodsPerFrame(f::MultiMPIFile) = length(f.files)*acqNumFrames(f.files[1])
acqNumFrames(f::MultiMPIFile) = 1

for op in [:filepath, :version, :uuid, :time, :studyName, :studyNumber, :studyTime, :studyUuid, :studyDescription,
            :experimentName, :experimentNumber, :experimentUuid, :experimentDescription,
            :experimentSubject, :experimentHasMeasurement,
            :experimentIsSimulation, :experimentIsCalibration, :experimentHasProcessing,
            :tracerName, :tracerBatch, :tracerVendor, :tracerVolume, :tracerConcentration,
            :tracerSolute, :tracerInjectionTime,
            :scannerFacility, :scannerOperator, :scannerManufacturer, :scannerName,
            :scannerTopology, :acqNumBGFrames,
            :acqStartTime,
            :dfNumChannels, :dfBaseFrequency, :dfDivider,
            :dfCycle, :dfWaveform, :rxNumChannels, :acqNumAverages, :rxBandwidth,
            :rxNumSamplingPoints, :rxTransferFunction, :rxTransferFunctionFileName,
            :rxInductionFactor, :rxUnit, :rxDataConversionFactor]
  @eval $op(f::MultiMPIFile) = $op(f.files[1])
end

rxHasTransferFunction(f::MultiMPIFile) = !(any(rxHasTransferFunction.(f) .== 0))

for op in [ :dfStrength, :dfPhase ]
  @eval begin function $op(f::MultiMPIFile)
       tmp = $op(f.files[1])
       newVal = similar(tmp, size(tmp,1), size(tmp,2),
                        acqNumFrames(f.files[1]),length(f.files))
       for c=1:length(f.files)
         tmp = $op(f.files[c])
         for y=1:acqNumFrames(f.files[1])
           for a=1:size(tmp,1)
             for b=1:size(tmp,2)
               newVal[a,b,y,c] = tmp[a,b]
             end
           end
         end
       end
      return reshape(newVal,size(newVal,1),size(newVal,2),:)
    end
  end
end

function acqOffsetField(f::MultiMPIFile)
   tmp = acqOffsetField(f.files[1])
   newVal = similar(tmp, 3, acqNumFrames(f.files[1]),length(f.files))
   for c=1:length(f.files)
     tmp = acqOffsetField(f.files[c])
     for b=1:acqNumFrames(f.files[1])
       for a=1:3
           newVal[a,b,c] = tmp[a,1,1,1]
       end
     end
   end
  return reshape(newVal,3,1,:)
end

function acqGradient(f::MultiMPIFile)
   tmp = acqGradient(f.files[1])
   newVal = similar(tmp, 3, 3, acqNumFrames(f.files[1]),length(f.files))
   for c=1:length(f.files)
     tmp = acqGradient(f.files[c])
     for b=1:acqNumFrames(f.files[1])
       for a=1:3
         for d=1:3
             newVal[a,d,b,c] = tmp[a,d,1,1]
         end
       end
     end
   end
  return reshape(newVal,3,3,1,:)
end

for op in [:measIsFourierTransformed, :measIsTFCorrected,
           :measIsBGCorrected,
           :measIsFastFrameAxis, :measIsFramePermutation, :measIsFrequencySelection,
           :measIsSpectralLeakageCorrected,
           :measFramePermutation, :measIsBGFrame]
  @eval $op(f::MultiMPIFile) = $op(f.files[1])
end

experimentHasReconstruction(f::MultiMPIFile) = false

## Achtung hack in der Schleife acqNumFrames(fi) statt acqNumFrames(f)
# notwendig, da hier Sprung zwischen MultiMPIFile und MPIFile
function measData(f::MultiMPIFile, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),
                  receivers=1:rxNumChannels(f))
  data = zeros(Float32, rxNumSamplingPoints(f), length(receivers),
                        length(frames),length(periods))
  #for (i,p) in enumerate(periods)
  #  data[:,:,:,i,:] = measData(f.files[p], frames, 1, receivers)
  #end
  for (i,fi) in enumerate(f.files)
    fr_fi=acqNumFrames(fi)
    data[:,:,:,fr_fi*(i-1)+1:fr_fi*i] = measData(fi, 1:fr_fi, 1, receivers)
  end
  return reshape(data,size(data,1),size(data,2),:,1)
end

function measDataTDPeriods(f::MultiMPIFile, periods=1:acqNumPeriods(f),
              receivers=1:rxNumChannels(f))

  data = zeros(Float32, rxNumSamplingPoints(f), length(receivers), length(periods))
  for (i,p) in enumerate(periods)
    l = divrem(p-1, acqNumPeriods(f.files[1]) )
    data[:,:,i] = measDataTDPeriods(f.files[l[1]+1], l[2]+1, receivers)
  end

  return data
end

# TODO: define functions for multi calibration data
#  if experimentIsCalibration(f)
#    for op in [:calibSNR, :calibFov, :calibFovCenter,
#               :calibSize, :calibOrder, :calibPositions, :calibOffsetField,
#               :calibDeltaSampleSize, :calibMethod]
#      setparam!(params, string(op), eval(op)(f))
#    end
#  end


text:
module Grid

using Base.Iterators

"""
    cell_vertices(d)

Return the local multi-indices of vertices of the `d`-dimensional brick element.

This function returns an object `𝓛` of type `CartesianIndices`. An
element `v ∈ 𝓛` represents the vertex with coordinates `(x[1], …, x[d])`

```
x[i] = (-1)^v[i] * h[i] / 2,    i = 1, …, d
```

See [Geometry of the reference brick element](@ref).

"""
cell_vertices(d::Int) = CartesianIndices(Tuple(fill(1:2, d)))

"""
    cell_vertices(p, N)

Return the global multi-indices of the vertices of a grid-cell.

The cell is specified through its `CartesianIndex`, `p`; the size of the grid is
defined by `N`.

This function returns a `d`-dimensional array `𝒢` of `CartesianIndex`, such that
`𝒢[l]` is the global index of the vertex with local index
`l ∈ CartesianIndices(1:2, …, 1:2)`.

Note that periodic boundary conditions are used.

"""
function cell_vertices(p::CartesianIndex{d}, N::NTuple{d,Int}) where {d}
    p_ = collect(Tuple(p))
    𝓛 = map(collect, product(fill(0:1, d)...))
    return [CartesianIndex(((p_ .+ q .- 1) .% N .+ 1)...) for q in 𝓛]
end

export cell_vertices

end


text:
# Functions for chemostatting reactants,
# and producing simplified reaction nets.

"""
    chemostatted!(reactions::Vector{Reaction}, spec::Vector{Int}, z)

Make a chemostatted version of the reaction system,
setting species spec to concentration z
"""
function chemostatted!(reactions::Vector{Reaction}, spec::Vector{Int}, z)
    for i = 1:length(reactions)
        r = reactions[i]

        nr = length(r.reactants)
        j=1
        while (j <= nr)
            idx = findfirst(x->x==r.reactants[j], spec)
            if idx != nothing
                nr -= 1
                r.kf *= z[idx]^r.stoichr[j]

                deleteat!(r.reactants, j)
                deleteat!(r.stoichr,   j)
            else
                j += 1
            end
        end

        np = length(r.products)
        j=1
        while (j <= np)
            idx = findfirst(x->x==r.products[j], spec)
            if idx != nothing
                np -= 1
                r.kr *= z[idx]^r.stoichp[j]

                deleteat!(r.products, j)
                deleteat!(r.stoichp,  j)
            else
                j += 1
            end
        end
    end
end

"""
    relabel!(reactions::Vector{Reaction}, prm)

re-label species according to permutation prm,
i.e. [1,3,2] maps [1,2,3]->[1,3,2]
"""
function relabel!(reactions::Vector{Reaction}, prm)
    for i = 1:length(reactions)
        r = reactions[i]

        nr = length(r.reactants)
        for j=1:nr
            r.reactants[j] = prm[r.reactants[j]]
        end

        np = length(r.products)
        for j=1:np
            r.products[j] = prm[r.products[j]]
        end
    end
end


text:
module ThinFilmsTools

include("Utils.jl")
export Utils

include("CommonStructures.jl")
using .CommonStructures
export PlaneWave,
       LayerTMMO,
       ModelFit,
       BoundariesFit

include("RefractiveIndicesDB.jl")
export RIdb

include("RefractiveIndicesModels.jl")
export RI

include("SpectraDBExamples.jl")
export SpectraDB

include("TransferMatrixMethod.jl")
using .TransferMatrixMethod
export tmm_optics, TMMOptics

include("FTFSpectrum.jl")
using .FTFSpectrum
export FitTMMOptics,
       Reflectance,
       Transmittance,
       Ellipsometry,
       NoAlpha,
       UseAlpha,
       MeanAbs,
       SumAbs,
       SumMeanAbs,
       space_solution_ema,
       theoretical_spectrum,
       normalize_reflectance,
       fit_tmm_optics

include("ThreeOmegaMethod.jl")
using .ThreeOmegaMethod
export ThreeOmega,
       LayerTOM,
       HeaterGeometry,
       Source,
       three_omega

include("PlottingTools.jl")
using .PlottingTools
export TOMPlot,
       Spectrum1D,
       SpectrumAngle1D,
       Spectrum2D,
       EMF2D,
       EMFAngle2D,
       FitSpectrum,
       FitSpectrumEllip,
       SpaceSolution,
       RIprofile,
       PBGDispersion1D,
       PBGDispersion1Dalt,
       PBGDispersion1Dimre,
       PBGDispersion2D,
       PBGDispersion2Dalt

include("FitBellCurves.jl")
using .FitBellCurves
export FitCurveModel, fit_curve_model

end # module


text:
"""
    sliceview(A, code)
    slicecopy(A, code)

Slice array `A` according to `code`, a tuple of length `ndims(A)`,
in which `:` indicates a dimension of the slices, and `*` a dimension separating them.
For example if `code = (:,*,:)` then slices are either `view(A, :,i,:)`
or `A[:,i,:]` with `i=1:size(A,2)`.
"""
function sliceview(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    if code == (:,*)
        collect(eachcol(A))
    elseif code == (*,:)
        collect(eachrow(A))
    elseif count(isequal(*), code) == 1
        collect(eachslice(A, dims = findfirst(isequal(*), code)))
    else
        iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
        [ view(A, i...) for i in iter ]
    end
end

@doc @doc(sliceview)
function slicecopy(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
    [ A[i...] for i in iter ]
end

sliceview(A::AbstractArray{T,N}) where {T,N} = sliceview(A, ntuple(i-> i==N ? (*) : (:), N))
slicecopy(A::AbstractArray{T,N}) where {T,N} = slicecopy(A, ntuple(i-> i==N ? (*) : (:), N))

"""
    glue!(B, A, code)
    copy_glue(A, code) = glue!(Array{T}(...), A, code)

Copy the contents of an array of arrays into one larger array,
un-doing `sliceview` / `slicecopy` with the same `code`.
Also called `stack` or `align` elsewhere.

    cat_glue(A, code)
    red_glue(A, code)

The same result, but calling either things like `hcat(A...)`
or things like `reduce(hcat, A)`.
The code must be sorted like `(:,:,:,*,*)`, except that `(*,:)` is allowed.

    glue(A)
    glue(A, code)

If `code` is omitted, the default is like `(:,:,*,*)`
with `ndims(first(A))` colons first, then `ndims(A)` stars.
If the inner arrays are `StaticArray`s (and the code is sorted) then it calls `static_glue`.
Otherwise it will call `red_glue`, unless code is unsuitable for that, in which case `copy_glue`.
"""
function glue(A::AbstractArray{<:AbstractArray{T,IN},ON}, code::Tuple=defaultcode(IN,ON)) where {T,IN,ON}
    if iscodesorted(code) || code == (*,:)
        red_glue(A, code)
    else
        copy_glue(A, code)
    end
end

glue(A::AbstractArray{<:Number,ON}, code::Tuple=defaultcode(0,ON)) where {ON} = A

defaultcode(IN::Int, ON::Int) = ntuple(d-> d<=IN ? (:) : (*), IN+ON)

@inline function red_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    if code == (:,*)
        reduce(hcat, A) # this is fast, specially optimised
    elseif code == (*,:)
        # mapreduce(transpose, vcat, A) # this is slow
        reduce(vcat, PermuteDims.(A)) # now avoiding transpose
    # elseif count(isequal(*), code) == 1 && code[end] == (*)
    #     reduce((x,y) -> cat(x,y; dims = length(code)), A) # this is slow
    elseif iscodesorted(code)
        flat = reduce(hcat, vec(vec.(A)))
        reshape(flat, gluedsize(A, code))
    else
        throw(ArgumentError("can't glue code = $code with reduce(cat...)"))
    end
end

@inline function cat_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    if code == (:,*)
        hcat(A...)
    elseif code == (*,:)
        vcat(PermuteDims.(A)...)
    elseif count(isequal(*), code) == 1 && code[end] == (*)
        cat(A...; dims = length(code))
    elseif iscodesorted(code)
        flat = cat(vec(A)...; dims = length(code)-ndims(A)+1)
        reshape(flat, gluedsize(A, code))
    else
        throw(ArgumentError("can't glue code = $code with cat(A...)"))
    end
end

@inline function copy_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    B = Array{eltype(first(A))}(undef, gluedsize(A, code))
    glue!(B, A, code)
end

function glue!(B::AbstractArray{T,N}, A::AbstractArray{IT,ON}, code::Tuple) where {T,N,IT,ON}
    gluecodecheck(A, code)
    N == ndims(A) + ndims(first(A))  || throw(DimensionMismatch("wrong size target"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(B,d) : Ref(:), Val(N))...)
    for i in iter
        copyto!(view(B, i...), A[decolonise(i)...] )
    end
    B
end

function gluecodecheck(A::AbstractArray, code::Tuple)
    colons = countcolons(code)
    inner = ndims(first(A))
    outer = ndims(A)
    colons == inner || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(first(A)) == $inner not cannot be glued with code = $(pretty(code))"))
    length(code) - colons == outer || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(A) == $outer cannot be glued with code = $(pretty(code))"))
end

@generated function decolonise(i::Tuple)
    ind = Int[]
    for k in 1:length(i.parameters)
        if i.parameters[k] != Colon
            push!(ind, k)
        end
    end
    Expr(:tuple, [Expr(:ref, :i, k) for k in ind]...)
end

@generated function countcolons(code::Tuple)
    n = 0
    for s in code.parameters
        if s == Colon
            n += 1
        end
    end
    n
end

@generated function iscodesorted(code::Tuple) # true if all colons before all stars
    colons = true
    sorted = true
    for s in code.parameters
        if s != Colon && colons # then we're at transition, change flag
            colons = false
        elseif s == Colon && !colons # then a : is following a *
            sorted = false
        end
    end
    sorted
end

@generated function gluedsize(A::AbstractArray{IT, N}, code::Tuple) where {IT, N}
    list = Any[]
    dout = 1
    din = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(first(A),$din) ))
            din += 1
        else
            push!(list, :( size(A,$dout) ))
            dout += 1
        end
    end
    :( ($(list...),) )
end

"""
    B = orient(A, code)

Reshapes `A` such that its nontrivial axes lie in the directions where `code` contains a `:`,
by inserting axes on which `size(B, d) == 1` as needed.
Throws an error if `ndims(A) != length(code)`.
"""
@generated function orient(A::AbstractArray, code::Tuple)
    list = Any[]
    pretty = Any[] # just for error
    d = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(A,$d) ))
            push!(pretty, ":")
            d += 1
        else
            push!(list, 1)
            push!(pretty, "*")
        end
    end
    str = join(pretty, ", ")
    d-1 == ndims(A) || throw(ArgumentError(
        "orient(A, ($str)) expeceted $(d-1) dimensions, got ndims(A) = $(ndims(A))"))
    :(reshape(A, ($(list...),)))
end

# performance: transpose faster than reshape? https://github.com/JuliaArrays/LazyArrays.jl/issues/16
orient(A::AbstractVector{T}, ::Tuple{typeof(*),Colon}) where {T <: Number} = PermuteDims(A)
#=
V = rand(500);
@time @reduce W[i] := sum(j,k) V[i]*V[j]*V[k] lazy;
# was 0.140318 seconds seconds, now 0.030957 seconds, factor 4
=#

# performance: avoid reshaping lazy transposes, as this is very slow in broadcasting
orient(A::Union{PermutedDimsArray, LinearAlgebra.Transpose, LinearAlgebra.Adjoint}, code::Tuple) =
    orient(collect(A), code)
#=
A = rand(10,100); B = rand(100,100); C = rand(100,10);
@time @reduce D[a,d] := sum(b,c) A[a,b] * B[b,c] * C[c,d] lazy;
# was 0.142307 seconds 76.296 MiB, now 0.001855 seconds 10.516 KiB, factor 75
=#

# performance: avoid that if you can just extract A.parent. TODO make @generated perhaps?
const LazyRowVec = Union{
    LinearAlgebra.Transpose{<:Any,<:AbstractVector},
    LinearAlgebra.Adjoint{<:Real,<:AbstractVector} }
orient(A::LazyRowVec, ::Tuple{typeof(*),Colon,Colon}) = orient(A.parent, (*,*,:))
orient(A::LazyRowVec, ::Tuple{Colon,typeof(*),Colon}) = orient(A.parent, (*,*,:))
orient(A::LazyRowVec, ::Tuple{typeof(*),typeof(*),Colon,Colon}) = orient(A.parent, (*,*,*,:))
orient(A::LazyRowVec, ::Tuple{typeof(*),Colon,typeof(*),Colon}) = orient(A.parent, (*,*,*,:))
orient(A::LazyRowVec, ::Tuple{Colon,typeof(*),typeof(*),Colon}) = orient(A.parent, (*,*,*,:))

# tidier but not strictly necc: for @cast Z[i] = A[3] + B[i]
orient(A::AbstractArray{<:Any, 0}, ::Tuple{typeof(*)}) = first(A)

"""
	rview(A, :,1,:) ≈ (@assert size(A,2)==1; view(A, :,1,:))

This simply reshapes `A` so as to remove a trivial dimension where indicated.
Throws an error if size of `A` is not 1 in the indicated dimension.

Will fail silently if given a `(:,3,:)` or `(:,\$c,:)`,
for which `needview!(code) = true`, so the macro should catch such cases.
"""
function rview(A::AbstractArray, code...)
	# This nice error has to happen at runtime, and thus is slow
	# all(decolonise(code) .== 1) ||
	#     throw(ArgumentError("can't use rview(A, code...) with code = $(pretty(code)), it accepts only 1 and :"))
	rview(A, code)
end

@generated function rview(A::AbstractArray, code::Tuple)
    list = Any[]
    pretty = Any[] # just for error
    ex = quote end
    d = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(A,$d) ))
            push!(pretty, ":")
            d += 1
        elseif s == Int
            push!(pretty, "1")
            str = "rview(A) expected size(A,$d) == 1"
            push!(ex.args, :(size(A,$d)==1|| throw(DimensionMismatch($str)) ) )
            d += 1
        end
    end
    str = join(pretty, ", ")
    d-1 == ndims(A) || throw(ArgumentError(
        "rview(A, ($str)) expeceted $(d-1) dimensions, got ndims(A) = $(ndims(A))"))
    push!(ex.args, :(reshape(A, ($(list...),))) )
    ex
end

rview(A::LazyRowVec, ::Tuple{Int,Colon}) = A.parent

rview(A::Union{PermutedDimsArray, LinearAlgebra.Transpose, LinearAlgebra.Adjoint}, code::Tuple) =
    view(A, code)
# cannot use rview(collect(A), code) as this may be on LHS, e.g. sum!(rview(),...)

"""
    rvec(x) = [x]
    rvec(A) = vec(A)

Really-vec... extends `LinearAlgebra.vec` to work on scalars too.
"""
rvec(x::Number) = [x]
rvec(A) = LinearAlgebra.vec(A)

"""
    PermuteDims(A::Matrix)
    PermuteDims(A::Vector)

Lazy like `transpose`, but not recursive:
calls `PermutedDimsArray` unless `eltype(A) <: Number`.
"""
PermuteDims(A::AbstractMatrix) = PermutedDimsArray(A, (2,1))
PermuteDims(A::AbstractMatrix{T}) where {T<:Number} = transpose(A)

PermuteDims(A::AbstractVector) = reshape(A,1,:)
PermuteDims(A::AbstractVector{T}) where {T<:Number} = transpose(A)

"""
    diagview(M)

Like `diag(M)` but makes a view.
"""
function diagview(A::AbstractMatrix)
    r,c = size(A)
    view(A, 1:(r+1):r*c)
end

diagview(A::LinearAlgebra.Diagonal) = A.diag


text:
#-------------------------------------------------------------------------------
# SimplestLogger  - from source code

"""
    SimplestLogger(stream=stderr, min_level=Info)

Simplistic logger for logging all messages with level greater than or equal to
`min_level` to `stream`.
"""

struct SimplestLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel
    message_limits::Dict{Any,Int}
end

SimplestLogger(stream::IO=stderr, level=Info) = SimplestLogger(stream, level, Dict{Any,Int}())

Logging.shouldlog(logger::SimplestLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

Logging.min_enabled_level(logger::SimplestLogger) = logger.min_level

Logging.catch_exceptions(logger::SimplestLogger) = false

function Logging.handle_message(logger::SimplestLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end
    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)

    println(iob, message)

    write(logger.stream, take!(buf))
    nothing
end


text:
function randWeight(a)

    n = a.n
    (ai,aj) = findnz(a)
    m = length(ai)

    # potentials or edge-based

    if (rand(1)[1] < .3)
        w = rand(m)

    else
        v = randn(a.n)

        # mult by matrix ?
        if (rand(1)[1] < .5)

            invdeg = spdiagm(1./(a*ones(size(a)[1])))
            if (rand(1)[1] < .5)
                for i in 1:10
                    v = a * (invdeg * v)
                    v = v - mean(v)
                end
            else
                for i in 1:10
                    v = v - a * (invdeg * v)
                    v = v - mean(v)
                end
            end
        end

        w = abs(v[ai]-v[aj])

    end

    # reciprocate or not?

    if (rand(1)[1] < .5)
        w = 1./w
    end

    w = w / mean(w)

    ar = sparse(ai,aj,w,n,n)
    ar = ar + ar';
    return ar
end


text:
#=
S. Boldo and J.M. Muller,
Exact and Approximated Error of the FMA,
IEEE Transactions on Computers, 60(2011), 157–164.

function fma_with_error(a::Float64, b::Float64, c::Float64)
    hi  = fma(a, b, c)
    u1, u2 = two_prod(a, b)
    v1, v2 = two_sum(c, u2)
    w1, w2 = two_sum(u1, v1)
    t = w1 - hi + w2
    md, lo = two_hilo_sum(t, v2)
    return hi, md, lo
end

  (a * b) + c == hi + md + lo
   fma(a,b,c) == hi
   |md + lo| = unit_roundoff_error(|hi|)/2
   |lo| = unit_roundoff_error(|md|)/2

function two_fma(a::Float64, b::Float64, c::Float64)
    hi  = fma(a, b, c)
    u1, u2 = two_prod(a, b)
    v1, v2 = two_sum(c, u2)
    w1, w2 = two_sum(u1, v1)
    md = (w1 - hi + w2) + v2
    return hi, md
end

=#

function Base.fma(xhi::T, xlo::T, yhi::T, ylo::T, zhi::T, zlo::T) where {T<:Float64}
   chi, c1 = two_prod(xhi, yhi)
   t0 = xlo * ylo
   t1 = fma(xhi, ylo, t0)
   c2 = fma(xlo, yhi, t1)
   c3 = c1 + c2
   dhi, dlo = two_hilo_sum(chi, c3)

   shi, slo = two_sum(zhi, dhi)
   thi, tlo = two_sum(zlo, dlo)
   c = slo + thi
   vhi, vlo = two_hilo_sum(shi, c)
   w = tlo + vlo
   hi, lo = two_hilo_sum(vhi, w)
   return Double64((hi, lo))
end

@inline function Base.fma(x::T, y::T, z::T) where {T<:Double64}
   return fma(Hi(x), Lo(x), Hi(y), Lo(y), Hi(z), Lo(z))
end

for T in (:Double64, :ComplexD64)
   @eval @inline function Base.muladd(x::$T, y::$T, z::$T)
      return x*y + z
   end
end

#=
   Quadruple-precision BLAS using Bailey's arithmetic with FMA instruction
   by S. Yamada, t. Ina, N. Sasa, Y.Idomura, M. Machida, T. Imamura
   2017 IEEE International Parallel and Distributed Processing Symposium Workshops
=#

function Base.:(*)(ahi::T, alo::T, bhi::T, blo::T) where {T<:Float64}
   p1 = ahi * bhi
   p2 = fma(ahi, bhi, -p1)
   p2 += alo*bhi + ahi*blo
   chi = p1 + p2
   clo = p2 - (chi - p1)
   return Double64((chi, clo))
end

@inline function Base.:(*)(a::T, b::T) where {T<:Double64}
   return (*)(Hi(a), Lo(a), Hi(b), Lo(b))
end