text:
using p7zip_jll
using Colors
using Mmap

"""
    readptx(
        fn::AbstractString,
        scale::EnergyScale, # Energy scale for hyperspectrum
        props::Dict{Symbol,Any},
        nch::Int; # Number of channels in hyperspectrum
        blocksize = 1,  # Size of blocks to sum to create averaged hyperspectrum
        dets=[true,true,true,true], # Which detectors to include
        frames=1:typemax(Int)) # which frames to include in hyperspectrum
    readptx(
        fn::AbstractString,
        scale::EnergyScale, # Energy scale for hyperspectrum
        nch::Int; # Number of channels in hyperspectrum
        blocksize = 1,  # Size of blocks to sum to create averaged hyperspectrum
        dets=[true,true,true,true], # Which detectors to include
        frames=1:typemax(Int)) # which frames to include in hyperspectrum

Read a SEMantics .ptx file into a HyperSpectrum.
"""
function readptx(
    fn::AbstractString,
    scale::EnergyScale,
    props::Dict{Symbol,<:Any},
    nch::Int;
    blocksize = 1,
    frames = 1:typemax(Int),
    dets = ( true, true, true, true ),
    datatype = UInt16
)
    function readmsg(io::IO)
        hdr = Array{UInt8}(undef, 32)
        if readbytes!(io, hdr) == 32
            name = String(hdr[1:findfirst(c -> c == 0x0, hdr)-1])
            len = reinterpret(UInt32, hdr[17:20])[1]
            msg = Array{UInt8}(undef, len)
            if readbytes!(io, msg) == len
                return (name, msg)
            end
        end
        return ("End", nothing)
    end
    mktempdir() do tmp_dir
        p7zip() do p7z_exe
            # Use pkzip to unzip the contents to a temporary directory
            run(`$p7z_exe -bsp0 -bso0 e "$fn" -o"$tmp_dir"`)
        end
        # Read and interpret the header
        xml = read(joinpath(tmp_dir,"Header"))
        hdr = readxml(IOBuffer(xml))
        nprops = Dict{Symbol,Any}(props)
        nprops[:Name] = splitpath(fn)[end]
        dwell = 1.0e-6 * parse(Float64, findfirst("/OpticState/mDwell", hdr).content)
        fov = 0.1 * parse(Float64, findfirst("/OpticState/mFieldOfView", hdr).content)
        width = parse(Int, findfirst("/OpticState/mWidth", hdr).content)
        height = parse(Int, findfirst("/OpticState/mHeight", hdr).content)
        rx = parse(Int, findfirst("/OpticState/mRaster/x", hdr).content)
        ry = parse(Int, findfirst("/OpticState/mRaster/y", hdr).content)
        rw = parse(Int, findfirst("/OpticState/mRaster/width", hdr).content)
        rh = parse(Int, findfirst("/OpticState/mRaster/height", hdr).content)
        hss = HyperSpectrum(
            scale, # EnergyScale
            nprops, # Properties
            ((rh - ry)÷blocksize, (rw - rx)÷blocksize), # Pixel dimensions
            nch, # Channel depth
            datatype, # Default UInt16,
            axisnames = ( "X", "Y" ),
            fov=(fov*(rw/width),fov*(rh/height)) # Image scale
        )
        open(joinpath(tmp_dir, "Data"),"r", lock=true) do data
            realtime, deadtime = zero(UInt64), zero(UInt64)
            maxchannel, maxframeid, minframeid = zero(Int32), zero(Int32), typemax(Int32)
            # elapsetime, prevelapse = zero(UInt64), typemax(UInt32)
            seen = Set{String}()
            images = Dict{Tuple{Int32,Int32},Array}()
            cmd, msg = readmsg(data)
            bpp=8
            while cmd != "End"
                # println("$cmd[$(isnothing(msg) ? "nothing" : length(msg))]")
                if cmd == "GetMapData"
                    # The output array consists of list of datagrams, each has got five double-word (32-bit) values (DW0 to
                    # DW4). Little endian is storage is applied. Most significant byte (bits 31-24) of W0 tells the event type:
                    # err = reinterpret(Int32, msg[1:4])[1]
                    (reinterpret(UInt32, msg[5:8])[1]!=0) && warn("A buffer overflow occured during the acquisition of the map data.")
                    # usage = reinterpret(Float32, msg[9:12])[1]
                    # @assert reinterpret(UInt32, msg[13:16])[1]==0
                    ws = reinterpret(UInt32, msg[13:end])
                    frame = zero(eltype(ws))
                    for i = 1:length(ws)÷5
                        w0, w1, w2, w3, w4 = ws[5*(i-1)+1:5*i]
                        # println("$(ws[5*(i-1)+1:5*i])")
                        # msg[2] # timestamp in 10 ns units
                        msgid = (w0 >>> 24) & 0xFF
                        if msgid == 0x80 # dead time
                            # W0 = 0x80xxxxxx – dead-time measurement
                            # W0 – bits 0-7: filter index (1= filter 1, 2 = filter 2, …) bits 8-11: detector index (0=1st detector, 1=2nd detector, …)
                            # W1 – real-time timer value in 10 ns
                            # W2 – dead-time timer value in 10 ns
                            # W3 – event counter
                            x_det = 1 + ((w0 >>> 8) & 0x0F)
                            if (x_det>=1) && (x_det<=length(dets)) && dets[x_det] && frame in frames
                                realtime += w1
                                deadtime += w2
                            end
                        elseif msgid == 0x00 # measured energy
                            # W0 = 0x00xxxxxx – measured energy
                            # W0 – bits 0-7: filter index (1= filter 1, 2 = filter 2, …)  bits 8-11: detector index (0=1st detector, 1=2nd detector, ...)
                            # W1 – energy in eV (24.8 signed fixed point number)
                            # W2 – timestamp in 10ns
                            # W3 – pixel index
                            # W4 – frame index
                            x_det, frame = 1 + ((w0 >>> 8) & 0x0F), w4
                            if (x_det>=1) && (x_det<=length(dets)) && dets[x_det] && frame in frames
                                ch = channel(Float64(w1 >>> 8), scale)
                                r, c = Int(w3 ÷ rw - ry + 1)÷blocksize, Int(w3 % rw - rx + 1)÷blocksize
                                #elapsetime = ( prevelapse == typemax(UInt64) ? 0 : elapsetime + (w2 >= prevelapse ? UInt64(w2 - prevelapse) : zero(typeof(elapsetime))) # w2 + (prevelapse - typemax(UInt32)))
                                #prevelapse = w2
                                if checkbounds(Bool, hss.counts, ch, r, c)
                                    @inbounds hss.counts[ch, r, c] += 1
                                end
                            end
                        end
                        # W0 = 0x10xxxxxx – line sync
                        # W0 – bit 0: “Line” signal state bits 8-11: detector index (0=1st detector, 1=2nd detector, ...)
                        # W2 – timestamp in 10ns
                        # W4 – frame index
                        # W0 = 0x11xxxxxx – pileup
                        # W0 – bits 8-11: detector index (0=1st detector, 1=2nd detector, …)
                        # W2 – timestamp in 10ns
                        # W4 – frame index
                        # W0 = 0x13xxxxxx – preamplifier reset
                        # W0 – bits 8-11: detector index (0=1st detector, 1=2nd detector, …)
                        # W2 – timestamp in 10ns
                        # W4 – frame index
                    end
                elseif cmd == "ScData"
                    # See SharkSEM API ScScanLine documentation
                    # ScData(in unsigned int frameid, in int channel, in unsigned int index, in int bpp, in char[] data);
                    frameid = reinterpret(Int32, msg[1:4])[1]+1
                    maxframeid = max(frameid, maxframeid)
                    minframeid = min(frameid, minframeid)
                    ch = reinterpret(Int32, msg[5:8])[1]+1
                    maxchannel = max(ch, maxchannel)
                    index  = reinterpret(Int32, msg[9:12])[1]
                    bpp = reinterpret(Int32, msg[13:16])[1]
                    @assert bpp==8 || bpp==16
                    if !haskey(images,(frameid, ch))
                        images[(frameid, ch)] = (bpp == 8 ? zeros(UInt8, (rh - ry, rw - rx)) : zeros(UInt16, (rh - ry, rw - rx)))
                    end
                    img=images[(frameid, ch)]
                    imgdata = (bpp==8 ? reinterpret(UInt8, msg[21:end]) : reinterpret(UInt16, msg[21:end]))
                    img[index-ry+1:index-rx+length(imgdata)] .= imgdata
                else
                    if !(cmd in seen)
                        println("Skipping: $cmd")
                        push!(seen, cmd)
                    end
                end
                cmd, msg = readmsg(data)
            end
            if (maxchannel>0) && (maxframeid>0)
                for frameid in minframeid:maxframeid
                    image = zeros(eltype(images[(minframeid,1)]),  (rh - ry, rw - rx, 1:maxchannel))
                    for ch in 1:maxchannel
                        if haskey(images, (frameid, ch))
                            image[:,:,ch] .= images[(frameid,ch)]
                        end
                    end
                    hss[Symbol("Image$(maxframeid-minframeid+1)")] = Gray.(reinterpret(bpp==8 ? N0f8 : N0f16, permutedims(image, (2, 1, 3))))
                end
            end
            hss[:Dwell] = dwell * blocksize^2
            # Elapse time is broken!!!
            # hss[:Elapse] = 1.0e-8 * elapsetime # second
            hss[:RealTime] = (1.0e-8 * blocksize^2 / ((rh - ry) * (rw - rx))) * realtime  # seconds
            hss.livetime .= (1.0e-8 * blocksize^2 / ((rh - ry) * (rw - rx))) * (realtime - deadtime)  # seconds
            hss[:DeadFraction] = deadtime / realtime
            return hss
        end # close the file
    end # delete the temporary directory
end
readptx(
    fn::AbstractString,
    scale::EnergyScale,
    nch::Int;
    blocksize = 1,
    frames = 1:typemax(Int),
    dets = ( true, true, true, true ),
    datatype = UInt16
) = readptx(fn, scale, Dict{Symbol,Any}(), nch, blocksize=blocksize, frames=frames, dets=dets, datatype=datatype)

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
function horizontal_projection(v::StaticVector{3})
    @inbounds return SVector(v[1], v[2])
end

function horizontal_projection(r::RotMatrix{3})
    @inbounds return RotMatrix(r[SVector(1, 2), SVector(1, 2)])
end

critically_damped_gains(k::Number) = PDGains(k, 2 * sqrt(k))

"""
``A x \\le b``
"""
struct HRep{T, TA<:AbstractMatrix{T}, Tb<:AbstractVector{T}}
    A::TA
    b::Tb
end

function HRep(hull::ConvexHull)
    HRep(hrep(hull)...)
end

function Base.:*(M::StaticMatrix, hrep::HRep)
    HRep(M * hrep.A, M * hrep.b)
end

function Base.:+(hrep::HRep, c::StaticVector)
    HRep(hrep.A, hrep.b + hrep.A * c)
end

function Base.in(x::AbstractVector, hrep::HRep)
    all(hrep.A * x .<= hrep.b)
end

function RigidBodyDynamics.transform(fhull::Framed{C}, tf::Transform3D) where C<:ConvexHull
    @framecheck fhull.frame tf.from
    R = horizontal_projection(rotation(tf))
    p = horizontal_projection(translation(tf))
    return in_frame(tf.to, C(map(x -> R * x + p, vertices(unwrap(fhull)))))
end

const SHRep{N, M, T, L} = HRep{T, SMatrix{N, M, T, L}, SVector{N, T}}
const MHRep{N, M, T, L} = HRep{T, MMatrix{N, M, T, L}, MVector{N, T}}

Base.zero(::Type{SHRep{N, M, T, L}}) where {N, M, T, L} = HRep(zero(SMatrix{N, M, T, L}), zero(SVector{N, T}))
Base.zero(::Type{MHRep{N, M, T, L}}) where {N, M, T, L} = HRep(zero(MMatrix{N, M, T, L}), zero(MVector{N, T}))

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Plots
x = rand(10^6)
b = 10.0 .^ (-5:0)
histogram(x; alpha=0.3, label="", bin=b, norm=true)

# %%
using Plots, StatsBase
x = rand(10^6)
b = 10.0 .^ (-5:0)
h = fit(Histogram{Float64}, x, b)
h.weights ./= diff(h.edges[1])
plot(h; alpha=0.3, label="")

# %%
using Plots
n = 10^6
f(x) = n * x * log(10)
x = rand(n)
b = 10.0 .^ (-5:0)
histogram(x; alpha=0.3, label="", bin=b, xscale=:log10, xlim=extrema(b))
x = 10.0 .^ (-5:0.01:0)
plot!(x, f; label="n × pdf on log-scaled xaxis", legend=:topleft)

# %%
using Plots
n = 10^6
f(x) = n * x * log(10)
x = rand(n)
b = 10.0 .^ (-5:0.2:0)
histogram(x; alpha=0.3, label="", bin=b, xscale=:log10, xlim=extrema(b))
x = 10.0 .^ (-5:0.01:0)
plot!(x, f; label="n × pdf on log-scaled xaxis", legend=:topleft)

# %%
using Plots
n = 10^6
f(x) = n * x * log(10)
x = rand(n)
b = 10.0 .^ (-5:0.2:0)
histogram(x; alpha=0.3, label="", bin=b, xscale=:log10, xlim=extrema(b))
x = 10.0 .^ (-5:0.01:0)
plot!(x, f; label="n × pdf on log-scaled xaxis", legend=:topleft)

# %%
using Plots, StatsBase
n = 10^6
f(x) = n * x * log(10)
x = rand(n)
b = 10.0 .^ (-5:0.2:0)
h = fit(Histogram{Float64}, x, b)
h.weights ./= diff(log10.(h.edges[1]))
plot(h; alpha=0.3, label="", xscale=:log10, xlim=extrema(b))
x = 10.0 .^ (-5:0.01:0)
plot!(x, f; label="n × pdf on log-scaled xaxis", legend=:topleft)

# %%
using Plots, StatsBase, Random
X = 3randexp(10^6)
bin = [0; [2^(k/2) - 2^(-1/2) for k in 0:8]]
@show round.(bin; digits=2)
h = fit(Histogram{Float64}, X, bin)
h.weights ./= diff(h.edges[1])
plot(h; alpha=0.3, label="")

# %%
using Plots, StatsBase, Random
n = 10^6
f(x) = n * exp(-x/3)/3 * x * log(10)
X = 3randexp(10^6)
bin = [1e-3; [2^(k/2) - 2^(-1/2) for k in 0:8]]
@show round.(bin; digits=2)
h = fit(Histogram{Float64}, X, bin)
h.weights ./= diff(log10.(h.edges[1]))
plot(h; alpha=0.3, label="", xscale=:log10, xlim=extrema(bin))
x = exp.(range(log.(extrema(bin))...; length=500))
plot!(x, f; label="n × pdf on log-scaled xaxis", legend=:topleft)

# %%

max_stars_repo_path: 0014/manual plotting of a histogram.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/manual plotting of a histogram.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/manual plotting of a histogram.jl

max_forks_repo_name: genkuroki/public


text:
using StatsBase
using Catch22
import StatsBase.std
using LinearAlgebra

function reStandardise(F::AbstractArray)
    idxs = vec(nanrows(F) .| constantrows(F))
    𝛔 = StatsBase.std(F, dims=2)
    𝛍 = mean(F, dims=2)
    𝛔[idxs] .= Inf
    return X -> normalise(X, 𝛍, 𝛔, standardise, 2)
end
export reStandardise

function reZero(F::AbstractArray, α::Float64=10.0)
    # We want to adjust a baseline data set so that features have a variance of 0 if the constant baseline is similar to the high dimensional baseline  (otherwise, the baseline does not change variance).
    # For this, a hyperbolic function of constant and high dimensional baseline variances is used to scale the output
    𝛔 = vec(StatsBase.std(F, dims=2))
    f(X) =  begin
                #C = 1 .- min.((1,), 𝛔./(vec(StatsBase.std(X, dims=2))))
                C = exp.(-α.*𝛔./(vec(StatsBase.std(X, dims=2))))
                𝛍 = vec(mean(X, dims=2))
                return normalise(X, 𝛍, C, standardise, 2)
            end
    return f
end
export reZero

# ------------------------------------------------------------------------------------------------ #
#                            Combine a high and low dimensional baseline                           #
# ------------------------------------------------------------------------------------------------ #
function reScale(x::AbstractVector, f::Function=identity)
    σ = std(x)
    σ′ = f(x)
    if σ′ == σ == 0.0
        σ′ = σ = 1.0 # Catch the limit
    end
    return σ′.*x./σ
end
function reScale(F::AbstractArray, f::Vector)
    F′ = deepcopy(F)
    for r = 1:size(F, 1)
        F′[r, :] = reScale(F′[r, :] , f[r])
    end
   return F′
end
function reScale(F::AbstractFeatureArray, f::AbstractFeatureVector)
    (Fᵣ, fᵣ) = intersectFeatures(F, f) # Assumes f has all of the features of F
    reScale(Fᵣ, vec(fᵣ))
end

function NonstationaryProcesses.rampInterval(Fₗ, Fₕ, F)
    𝛔ₗ, 𝛔ₕ = std(Array(Fₗ), dims=2), std(Array(Fₕ), dims=2)
    𝛔ₕ[𝛔ₕ .< 𝛔ₗ] .= Inf
    𝐟 = [𝐟 -> rampInterval(0.0, 1.0, 𝛔ₗ[i], 𝛔ₕ[i])(std(𝐟)) for i ∈ 1:size(F, 1)]
end
function NonstationaryProcesses.rampOn(Fₗ, Fₕ, F)
    𝛔ₗ, 𝛔ₕ = std(Array(Fₗ), dims=2), std(Array(Fₕ), dims=2)
    𝛔ₕ[𝛔ₕ .< 𝛔ₗ] .= Inf
    𝐟 = [𝐟 -> rampOn(0.0, 1.0, 𝛔ₗ[i], 𝛔ₕ[i])(std(𝐟)) for i ∈ 1:size(F, 1)]
end

function intervalscale(Fₗ::Array{Float64, 2}, Fₕ::Array{Float64, 2},
                    interval::Function=rampInterval)
    if size(Fₗ, 1) != size(Fₕ, 1)
        error("High and low dimensional baselines do not have the same number of features")
    end
    return F -> reScale(F, interval(Fₗ, Fₕ, F))
end
function intervalscale(Fₗ::AbstractFeatureMatrix, Fₕ::AbstractFeatureMatrix,
    interval::Function=rampInterval)
    if any(Catch22.featureDims(Fₗ) .!= Catch22.featureDims(Fₕ))
        error("High and low dimensional baselines do not have the same features")
    end
    return F -> reScale(F,  Catch22.featureVector(interval(Fₗ, Fₕ, F), Catch22.featureDims(F)))
end
function intervalscale(Fₗ::Array{Float64, 2}, Fₕ::AbstractFeatureMatrix,
    interval::Function=rampInterval)
    if size(Fₗ, 1) .!= size(Fₕ, 1)
        error("High and low dimensional baselines do not have the same features")
    end
    return F -> reScale(F,  Catch22.featureVector(interval(Fₗ, Fₕ, F), Catch22.featureDims(F)))
end
export intervalscale

# ------------------------------------------------------------------------------------------------ #
#                                    Scale a baseline using PCA                                    #
# ------------------------------------------------------------------------------------------------ #
function orthogonalise(F::AbstractArray, dimensionalityReduction=principalcomponents)
    M = dimensionalityReduction(F)
    F̂ = embed(M, F)
    return (F̂, M)
end
function orthogonalise(F::AbstractFeatureMatrix, dimensionalityReduction=principalcomponents)
    F̂, M = orthogonalise(Array(F), dimensionalityReduction)
    F̂ = Catch22.featureMatrix(F̂, [Symbol("PC$x") for x ∈ 1:size(F̂, 1)])
    return F̂, M
end
export orthogonalise

function orthogonalBaseline(F::AbstractArray, dimensionalityReduction=principalcomponents)
    function 𝑏(F_test)
        @assert size(F, 1) == size(F_test, 1)
        F̂, M = orthogonalise(F, dimensionalityReduction)
        F_out = embed(M, Array(F_test))
    end
    return 𝑏
end
function orthogonalBaseline(F::AbstractFeatureMatrix, dimensionalityReduction=principalcomponents)
    function 𝑏(F_test)
        F̂_test, F̂ = intersectFeatures(F_test, F)
        F̂, M = orthogonalise(F̂, dimensionalityReduction)
        F_out = embed(M, Array(F̂_test))
        F_out = Catch22.featureMatrix(F_out, [Symbol("PC$x") for x ∈ 1:size(F_out, 1)])
    end
    return 𝑏
end
export orthogonalBaseline

function orthonormalHiloBaseline(F::AbstractFeatureArray, ℱₗ::AbstractFeatureArray, ℱₕ::AbstractFeatureArray; interval::Function=(x, y) -> NonstationaryProcesses.rampInterval(0, 1, x, y))
    F, ℱₗ, ℱₕ = intersectFeatures(F, ℱₗ, ℱₕ) # Intersects to the feature set of F
    ℱₕ′, M = orthogonalise(Array(ℱₕ))
    ℱ′ = Catch22.featureMatrix(ℱₕ′, [Symbol("PC$x") for x ∈ 1:size(ℱₕ′, 1)])
    F′ = embed(M, F)
    ℱ′ₗ = embed(M, ℱₗ)
    𝑏′ = intervalscale(Array(ℱ′ₗ), Array(ℱₕ′), interval)
    return 𝑏′(F′)
end
orthonormalHiloBaseline(ℱₗ::AbstractFeatureArray, ℱₕ::AbstractFeatureArray; kwargs...) = F -> orthonormalHiloBaseline(F, ℱₗ, ℱₕ; kwargs...)
export orthonormalHiloBaseline

# ------------------------------------------------------------------------------------------------ #
#                               Filter features using correlations/MI                              #
# ------------------------------------------------------------------------------------------------ #
function meanDependence(f, F; metric=StatsBase.corspearman)
    ρ = [metric(f, f′) for f′ ∈ eachrow(F)]
    ρ̄ = mean(abs.(ρ))
end

function meanDependence(F; metric=StatsBase.corspearman)
    𝐟 = zeros(size(F, 1))
    𝐟 = [meanDependence(F[x, :], F[setdiff(1:size(F, 1), [x]), :]; metric) for x ∈ 1:size(F, 1)]
end
export meanDependence

function dependencyFilter(F, threshold=0.3; metric=StatsBase.corspearman, iterations=Inf, direction=:min)
    # Iteratively filter features according to algorithm in Ben's thesis, page 214
    iteration = 0
    if direction == :min
        direction = 1
    elseif direction == :max
        direction = -1
    end
    Ī = zeros(size(F, 1)) .+ (-direction + 1)/2
    while any(direction.*Ī .< direction.*threshold) && iteration < iterations # Yes Ī keeps changing size
        Ī = meanDependence(F; metric)
        I, f = findmin(direction.*Ī)
        F = F[setdiff(1:lastindex(F, 1), [f]), :]
        iteration += 1
    end
    return F
end
export dependencyFilter

# * We have four baselines associated with rescaling variances, plus the PCA whitening rotation that can be performed before all four rescalings

# ------------------------------------------- Rescaling ------------------------------------------ #
"""
Standardise the test features. Alternatively, just use the normalisation field of an Inference.
"""
standardbaseline(F::AbstractArray) = standardise(F, 2)
standardbaseline() = F -> standardbaseline(F)
export standardbaseline

"""
Scale the features so that the variance of a constant baseline is zero. Do this by setting mapping variances less than σₗ to 0, but keeping a gradient of 1.0 afterwards. Zscore features, zscore baseline from features and then map
"""
function lowbaseline(Fₗ::AbstractArray)
    interval = x -> NonstationaryProcesses.rampOn(0.0, 1.0, x, x+1.0)
    function lowscale(F::AbstractArray)
        F̂, F̂ₗ = intersectFeatures(F, Fₗ)
        𝛔 = StatsBase.std(F̂, dims=2)
        𝛍 = StatsBase.mean(F̂, dims=2)
        F̂ₗ = normalise(F̂ₗ, 𝛍, 𝛔, standardise, 2)
        F̂ = normalise(F̂, 𝛍, 𝛔, standardise, 2)
        𝛔 = StatsBase.std(F̂, dims=2)
        𝛍 = StatsBase.mean(F̂, dims=2)
        𝛔ₗ = StatsBase.std(F̂ₗ, dims=2)
        𝛍ₗ = StatsBase.mean(F̂ₗ, dims=2)
        𝐟 = [x -> interval(σ)(StatsBase.std(x)) for σ ∈ 𝛔ₗ][:]
        if F isa AbstractFeatureArray
            𝐟 = Catch22.FeatureVector(𝐟, Catch22.featureDims(F̂ₗ))
        end
        return reScale(F̂, 𝐟)
    end
    return lowscale
end
export lowbaseline

"""
Scale the features so that the variance of a high dimensional baselines is unity. Do this as an inverval with a constant variance of 0.0
"""
function highbaseline(Fₕ::AbstractArray)
    Fₗ = zeros(size(Fₕ))
    intervalscale(Fₗ, Fₕ)
end
function highbaseline(Fₕ::AbstractFeatureArray)
    Fₕ, Fₗ = intersectFeatures(Fₕ, zeros(size(Fₕ)))
    intervalscale(Fₗ, Fₕ)
end
export highbaseline

"""
Scale the features so that their variances map to a rampInterval between the low (0) and high dim (1) baselines.
"""
function intervalbaseline(Fₗ::AbstractArray, Fₕ::AbstractArray)
    intervalscale(Fₗ, Fₕ)
end
export intervalbaseline

# # ---------------------------------- High dim orthonormalisation --------------------------------- #
# """
# Add this to any baseline variables and the Inference normalisation to transform into the high dim. whitened space
# e.g. infer(S, var; parameters, features, baseline=intervalbaseline(𝑜(Fₗ), 𝑜(Fₕ)), normalisation=𝑜) # Note normalisation occurs before baseline
# """
# orthonormaliseto(Fₕ::AbstractArray, dimensionalityReduction=principalcomponents) = orthogonalBaseline(Fₕ, dimensionalityReduction)
# export orthonormaliseto

# ---------------------------------- High dim orthogonalisation --------------------------------- #
"""
Add this to a baseline after constructing a scaling baseline and using that as normalisation
e.g. 𝑏 = intervalbaseline(Fₗ, Fₕ)
infer(S, var; parameters, features, baseline=orthogonaliseto(𝑏(Fₕ)), normalisation=𝑏) # Note normalisation occurs before baseline
"""
orthogonaliseto(Fₕ::AbstractArray, dimensionalityReduction=principalcomponents) = orthogonalBaseline(Fₕ, dimensionalityReduction)
export orthogonaliseto

"""
Interval scaling informed by distributions. The idea is to decrease the scale of a feature if we are unsure of its location between the high and zero dim distributions
"""
function significanceinterval(Fₗ, Fₕ, F)
    𝐟 = rampInterval(Fₗ, Fₕ, F)
    𝑝₀ = [f -> pvalue(VarianceFTest(Fₗ[i, :], f), tail=:right) for i ∈ 1:size(Fₗ, 1)]
    𝑝ₕ = [f -> pvalue(VarianceFTest(Fₕ[i, :], f), tail=:left) for i ∈ 1:size(Fₕ, 1)]
    function out(f, i)
        𝐟[i](f[:]).*𝑝₀[i](f[:]).*𝑝ₕ[i](f[:])
    end
    return [f -> out(f, i) for i ∈ 1:size(𝐟, 1)]
end
export significanceinterval

"""
errorintervalscaling
"""
function errorintervalscaling(Fₗ, Fₕ, F)
    𝐟 = rampInterval(Fₗ, Fₕ, F)
    function out(f, i)
        Δ = rampInterval(0.0, 1.0, 0.0, 1.0)(bootstrapSEσ(f[:])/(std(Fₕ[i, :]) - std(Fₗ[i, :]))) # So it saturates at extremes
        𝐟[i](f[:])*(1 - Δ)
    end
    return [f -> out(f, i) for i ∈ 1:size(𝐟, 1)]
end
export errorintervalscaling

"""
Scale the features so that their variances map to a rampInterval between the low (0) and high dim (1) baselines.
"""
function errorintervalbaseline(Fₗ::AbstractArray, Fₕ::AbstractArray)
    intervalscale(Fₗ, Fₕ, errorintervalscaling)
end
export errorintervalbaseline

"""
What's this?
"""
# Total covariance correction
function dependencyscalingnorotation(𝑏, Fₕ)
    Fₕ′ = 𝑏(Fₕ)
    function g(F)
        F′, Fₕ′ = intersectFeatures(noconstantrows(F), Fₕ′)
        # Fₕ′, F′ = intersectFeatures(noconstantrows(Fₕ′), F′)
        Σₕ² = StatsBase.cov(Array(Fₕ′), dims=2)
        𝑜 = orthogonaliseto(Fₕ′, principalcomponents)
        𝐧 = sum(abs.(Array(Σₕ²)), dims=2)
        𝐧[𝐧 .== 0] .= Inf
        N⁻¹ = FeatureMatrix(inv(sqrt(Diagonal(𝐧[:]))), getnames(Fₕ′))
        return FeatureMatrix(N⁻¹*𝑏(F), getnames(F′))
    end
    return g
end
export dependencyscalingnorotation


text: # ispos, isneg, isnormal, ntneg, ntnan, ntnormal, ispow2


text:
"""

```
writexyz(x::Vector{T}, file::String) where T <: AbstractVector
```

Print test xyz file.

"""
function writexyz(x::Vector{T}, file::String) where {T<:AbstractVector}
    f = open(file, "w")
    nx = length(x)
    println(f, nx)
    println(f, "title")
    for i = 1:nx
        println(f, "H $(x[i][1]) $(x[i][2]) $(x[i][3])")
    end
    close(f)
    nothing
end


text:
using Statistics: mean

struct ReachableDiversityLogger <: Logger
    original_organisms_alive::Vector{UInt64}
    child_genotypes::Vector

    function ReachableDiversityLogger(snapshot)
        new(get_organism_ids(snapshot), [])
    end
end

function get_reachable_diversity(snapshot, rel_tolerance, min_samples, max_samples)
    reachable_diversity = estimate(rel_tolerance, min_samples, max_samples) do
        logger = ReachableDiversityLogger(snapshot)
        simulate_snapshot!(should_terminate, snapshot, logger)

        while length(logger.child_genotypes) < 2
            logger = ReachableDiversityLogger(snapshot)
            simulate_snapshot!(should_terminate, snapshot, logger)
        end

        phenotype_similarities = []
        num_samples_per_run = 5

        for _ in 1:num_samples_per_run
            genotype_1 = rand(logger.child_genotypes)
            genotype_2 = rand(logger.child_genotypes)

            snapshot_to_test = sample_snapshot_id() |> get_snapshot
            sample_id_to_replace = get_id(snapshot_to_test, get_organisms(snapshot_to_test) |> rand)

            phenotype_similarity = (get_fitness(snapshot_to_test, sample_id_to_replace, genotype_1) - get_fitness(snapshot_to_test, sample_id_to_replace, genotype_2))^2

            push!(phenotype_similarities, phenotype_similarity)
        end

        return phenotype_similarities
    end

    return reachable_diversity
end

should_terminate(logger::ReachableDiversityLogger, snapshot) = length(logger.original_organisms_alive) == 0

function log_step(::ReachableDiversityLogger, model) end
function save_log(::ReachableDiversityLogger) end

function log_birth(logger::ReachableDiversityLogger, model, child, parents=nothing)
    if parents === nothing return end

    child_genotype = get_genotype(model, child)
    parent_genotype = get_genotype(model, rand(parents))

    if child_genotype == parent_genotype return end

    push!(logger.child_genotypes, child_genotype)
end

function log_death(logger::ReachableDiversityLogger, model, organism)
    id = get_id(model, organism)

    if id in logger.original_organisms_alive
        delete!(logger.original_organisms_alive, id)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module GlobalSearchRegressionGUI
using GlobalSearchRegression, HTTP, WebSockets, DataStructures, Mux, JSON, CSV, Pkg, Distributed, UUIDs, Base64

const SERVER_BASE_DIR = "../front/dist"
const INSTALLED_PACKAGES = Pkg.installed()
const GSREG_VERSION = ( haskey(INSTALLED_PACKAGES, "GlobalSearchRegression") ) ? INSTALLED_PACKAGES["GlobalSearchRegression"] : v"1.0.3"

mutable struct GSRegJob
    file # tempfile of data
    hash # hash id of user
    options # options for calculation
    id # unique identifier for this job
    time_enqueued # time enqueued
    time_started
    time_finished
    res::GlobalSearchRegression.GSRegResult # results array
    GSRegJob(file, hash, options) = new(file, hash, options, string(uuid4()), time())
end

"""
    Enqueue the job and notify worker
"""
function enqueue_job(job::GSRegJob)
    global job_queue
    global job_queue_cond
    enqueue!(job_queue, job)
    notify(job_queue_cond)
    job
end

function sendMessage(id, data)
    global connections
    if haskey(connections, string(id))
        ws = connections[id]
        sendMessage(ws, data)
    end
end

function sendMessage(ws::WebSocket, data)
    write(ws, JSON.json(data))
end

function gsreg(job::GSRegJob)
    global jobs_finished
    try
        sendMessage(job.hash, Dict("operation_id" => job.id, "message" => "Reading data"))
        data = CSV.read(job.file)
        sendMessage(job.hash, Dict("operation_id" => job.id, "message" => "Executing GSReg"))

        opt = job.options["options"]
        opt[:parallel] = job.options["paraprocs"]

        job.time_started = time()
        job.res = GlobalSearchRegression.gsreg(job.options["depvar"] * " " * join(job.options["expvars"], " "), data; opt...,
                    onmessage = message -> sendMessage(job.hash, Dict("operation_id" => job.id, "message" => message)) )
        job.time_finished = time()
        push!(jobs_finished, Pair(job.id, job))

        sendMessage(job.hash, Dict("operation_id" => job.id, "done" => true, "message" => "Successful operation", "result" => GlobalSearchRegression.to_dict(job.res)))
    catch e
        io = IOBuffer()
        showerror(io, e)
        sendMessage(job.hash, Dict("operation_id" => job.id, "done" => false, "message" => String(take!(io))))
    end
end

"""
    Log request for analytics and debugging
"""
function logRequest(app, req)
#     log = string(
#         get(req[:headers], "Host", ""), " ",
#         get(req[:headers], "Origin", ""), " ",
#         Libc.strftime("%d/%b/%Y:%H:%M:%S %z", time()), " ",
#         "\"", req[:method], " ", req[:resource], "\"", " ",
#         size(req[:data],1), "B",
#         "\n")
#
#     open(joinpath(dirname(@__FILE__), "..", "access.log"), "a") do fp
#         write(fp, log)
#     end
    app(req)
end

"""
    Try to execute controller and if there is any exception, show exception message
"""
function errorCatch(app, req)
    try
        app(req)
    catch e
        io = IOBuffer()
        showerror(io, e)
        err_text = String(take!(io))
        toJsonWithCors(Dict("message" => err_text, "error" => true), req)
    end
end

"""
    Setting CORS headers and parsing it to JSON
"""
function toJsonWithCors(res, req)
    headers = []

    push!(headers, "Server" => "Julia/$VERSION")
    push!(headers, "Content-Type" => "text/html; charset=utf-8")

    if( req[:method] != "OPTIONS" )
        push!(headers, "Content-Type" => "application/json; charset=utf-8")
    end
    push!(headers, "Access-Control-Allow-Headers" => "X-User-Token, Content-Type")
    push!(headers, "Access-Control-Allow-Origin" => "*")

    Dict(
        :headers => headers,
        :body => codeunits((req[:method] == "OPTIONS") ? "" : JSON.json(res))
    )
end

"""
    Get Auth from header, and reply an error when it's not present
"""
function authHeader(app, req)
    headers = Dict(req[:headers])
    if(haskey(headers, "X-User-Token"))
        req[:token] = get(headers, "X-User-Token", "")
    end
    app(req)
end

function constructCommand(options)
    command = """Pkg.add("GlobalSearchRegression")
using GlobalSearchRegression, CSV
data = CSV.read("yourcsvfile.csv")
res = GlobalSearchRegression.gsreg("", data; opts...)"""
end

"""
    TODO: doc and d better validation, more secure for cloud environments
"""
function validateInput!(opt)
    options = Dict{Symbol,Any}()
    opt_types = Dict(
        "intercept" => Bool,
        "time" => String,
        "residualtest" => Bool,
        "ttest" => Bool,
        "orderresults" => Bool,
        "modelavg" => Bool,
        "outsample" => Int,
        "method" => String,
        "criteria" => Array{Symbol}
    )
    for (name, value) in opt["options"]
        if value != nothing && name == "criteria"
            push!(options, Pair(:criteria, map(Symbol, value)))
        elseif value != nothing && name == "time"
            push!(options, Pair(:time, Symbol(value)))
        elseif value != nothing && name != "csv" && name != "resultscsv"
            push!(options, Pair(Symbol(name), value))
        end
    end
    opt["options"] = options
    opt
end

"""
    Receives base64 encoded data in path with regression
    options and a content body with CSV file for processing

    req[:data] -> Array{UInt8,1}
    req[:params][:b64] -> String
"""
function upload(req)
    """
    Save the file to tmp folder, if there is any exception, it should be returned to the user for
    reuse. It must be available until the user explicit deletes it, or its deleted from
    operating system. If it is deleted, the file must be uploaded again.
    """
    tempfile = try
        temp = tempname()
        write(temp, IOBuffer(req[:data]))
        temp
    catch
        error("We can't save that file, try again later.")
    end

    """
    Check csv parsing & get data
    """
    data = try
        CSV.read(IOBuffer(req[:data]))
    catch
        error("The file must be a valid CSV")
    end

    global files_dict
    id = string(uuid4())
    push!(files_dict, Pair(id, tempfile))

    Dict(
        "filename" => id,
        "datanames" => names(data),
        "nobs" => size(data, 1)
    )
end

function server_info(req)
    global job_queue
    Dict(
        "ncores" => Sys.CPU_THREADS,
        "nworkers" => nworkers(),
        "gsreg_version" => string(GSREG_VERSION),
        "julia_version" => string(VERSION),
        "job_queue" => Dict(
            "length" => length(job_queue)
        )
    )
end

"""
    Enqueue the execution of regressions, expecting a confirmation. If any there is any error, should be a detailed report.
"""
function solve(req)
    """
    Try to get the filename from params
    """
    global files_dict

    if haskey(files_dict, req[:params][:hash])
        tempfile = files_dict[req[:params][:hash]]
        if (!isfile(tempfile))
            error("File was deleted")
        end
    else
        error("Filekey inexistent")
    end

    """
    Input options should be a valid base64:json encoded String.
    """
    options = try
        b64 = convert(String, req[:params][:options])
        JSON.parse(String(base64decode(b64)))
    catch
        error("Bad format in options")
    end

    # validate correct use of options
    validateInput!(options)

    # Enqueue the job
    job = GSRegJob(tempfile, req[:token], options)
    enqueue_job(job)

    Dict(
        "ok" => true,
        "operation_id" => job.id,
        "in_queue" => length(job_queue)
    )
end

function validpath(path; dirs = true)
    return true # fallback index.html will manage this
    (isfile(path) || (dirs && isdir(path))) || isfile(joinpath(path,"index.html"))
end

ormatch(r::RegexMatch, x) = r.match
ormatch(nothing, x) = x

extension(f) = ormatch(match(r"(?<=\.)[^\.\\/]*$", f), "")

mimetypes = Dict{AbstractString,AbstractString}([
    ("jpg", "image/jpeg"),
    ("png", "image/png"),
    ("html", "text/html"),
    ("csv", "text/csv"),
    ("css", "text/css"),
    ("js", "application/javascript"),
    ("ttf", "application/x-font-ttf"),
    ("svg", "image/svg+xml"),
    ("mp4", "video/mp4"),
    ("ico", "image/x-icon"),
    ("map", "text/plain")
])

fileheaders(f) = Dict("Content-Type" => get(mimetypes, extension(f), "application/octet-stream"))

fileresponse(f) = Dict(
                    :body => read(f),
                    :headers => fileheaders(f)
                    )

fresp(f) = isfile(f) ? fileresponse(f) : fileresponse(joinpath(dirname(@__FILE__), SERVER_BASE_DIR,"index.html"))

"""
    operation_id -> csv
"""
function result_file(req)
    global jobs_finished
    id = req[:params][:id]
    if haskey(jobs_finished, id)
        csv = IOBuffer()
        job = jobs_finished[id]
        export_csv(csv, job.res)
        delete!(jobs_finished, id)
        if job.options["csv"] != nothing
            filename = endswith(job.options["csv"], ".csv") ? job.options["csv"] : job.options["csv"] * ".csv"
        else
            filename = "gsreg" * Libc.strftime("%d/%b/%Y%H%M%S", time()) * ".csv"
        end
        Dict(
            :body => String(take!(csv)),
            :headers => Dict("Content-Type" => "application/octet-stream", "Content-Disposition" => "attachment; filename=$filename")
        )
    else
        Dict( :status => 404, :body => "Not found" )
    end
end

files_dict = Dict{String,String}()
jobs_finished = Dict{String,GSRegJob}()

job_queue = Queue{GSRegJob}()

job_queue_cond = Condition()

# Array of connections based in user token
connections = Dict{String,WebSocket}()

noop(app, req) = app(req)

function gui(;openbrowser=true, port=45872, cloud=false, log=false)
    global job_queue_cond
    global job_queue

    #Consume job_queue after notification
    function consume_job_queue()
        while true
            wait(job_queue_cond)
            while !(isempty(job_queue))
#                 try
                    gsreg(dequeue!(job_queue))
#                 catch
#                     println("Error executing gsreg job")
#                 end
            end
        end
    end

    schedule(Task(consume_job_queue))

    @app app = (
        Mux.stack(Mux.todict, logRequest, Mux.splitquery, errorCatch, authHeader, Mux.toresponse),
        page("/info", respond("
info
")),
        page("/upload", req -> toJsonWithCors(upload(req), req)),
        page("/server-info", req -> toJsonWithCors(server_info(req), req)),
        page("/result/:id", req -> result_file(req)),
        page("/solve/:hash/:options", req -> toJsonWithCors(solve(req), req)),
        branch( req -> validpath(joinpath(SERVER_BASE_DIR, req[:path]...)), req -> fresp(joinpath(dirname(@__FILE__), SERVER_BASE_DIR, req[:path]...))),
        Mux.notfound
    )

    """
     This WebSocket handler mantain a collection of ws opened with the user id.
    """
    function wshandler(ws)
        global connections
        while isopen(ws)
            msg, = readguarded(ws)
            try
                msg = JSON.parse(String(copy(msg)))
                id = msg["user-token"]
                connections[id] = ws
                sendMessage(ws, Dict("ok" => true, "message" => "Waiting in queue"))
            catch
                sendMessage(ws, Dict("ok" => false, "message" => "The next message must be in JSON format"))
            end
        end
        # TODO: find some way for delete connection with any reference
        # delete!(connections, id)
        # close(ws)
    end

    wsapp = Mux.App(Mux.mux(
        Mux.wdefaults,
        Mux.route("/ws", req -> wshandler(req[:socket])),
        Mux.wclose,
        Mux.notfound(),
    ))

    @async Mux.serve(app, wsapp, port)

    if(openbrowser)
        url = "http://127.0.0.1:" * string(port)
        sleep(3)
        if Sys.iswindows()
            run(`cmd /c start $url`)
        else
            start = ( Sys.isapple() ? "open" : "xdg-open" )
            run(Cmd([start,url]))
        end
    end
end

export gui

end # module GlobalSearchRegressionGUI


text:
mutable struct tcoeffs{T}
       mu::Array{T,2}
       hc::Array{T,1}
	   hb::Array{T,1}
	   nu::Array{T,2}
	   alpha::Array{T,2}
end

mutable struct tcache{uType,elTypeu,uLowType,low_prec_type}
	U::Array{uType,1}
	Uz::Array{uType,1}
	L::Array{uType,1}
	Lz::Array{uType,1}
	F::Array{uType,1}
	Dmin::Array{uType,1}
	Eval::Array{Bool,1}
	DY::Array{elTypeu,1}
	rejects::Array{Int64,1}
	nfcn::Array{Int64, 1}
	lambdas::Array{elTypeu,1}
end

mutable struct tcacheMix{uType,elTypeu,uLowType,low_prec_type}
	U::Array{uType,1}
	Uz::Array{uType,1}
	L::Array{uType,1}
	Lz::Array{uType,1}
	F::Array{uType,1}
	Dmin::Array{uType,1}
	Eval::Array{Bool,1}
	DY::Array{elTypeu,1}
	rejects::Array{Int64,1}
	nfcn::Array{Int64, 1}
	lambdas::Array{elTypeu,1}
	Ulow:: Array{uLowType,1}
	DU::Array{uLowType,1}
	DF::Array{uLowType,1}
	DL::Array{uLowType,1}
	delta::Array{uLowType,1}
	Fa::Array{uLowType,1}
	Fb::Array{uLowType,1}
	normU::Array{low_prec_type,1}
	lhb::Array{low_prec_type,1}
	lmu::Array{low_prec_type,2}
end

abstract type IRKAlgorithm  <: OrdinaryDiffEqAlgorithm end
struct IRKGL16 <: IRKAlgorithm end

function DiffEqBase.__solve(prob::DiffEqBase.AbstractODEProblem{uType,tType,isinplace},
     alg::IRKAlgorithm,args...;
#     dt::tType=0.,
     dt=0.,
     saveat=1,
     maxiters=100,
	 maxtrials=3,            # maximum of unsuccessful trials
     save_everystep=true,
     initial_interp=true,
	 adaptive=true,
	 reltol=1e-6,
	 abstol=1e-6,
	 myoutputs=false,
	 threading=false,
     mixed_precision=false,
	 low_prec_type = Float64,
     kwargs...) where{uType,tType,isinplace}

	 s = 8
     destats = DiffEqBase.DEStats(0)

    if (typeof(prob.f)<:DynamicalODEFunction)
	    @unpack tspan,p=prob
        f1=prob.f.f1
        f2=prob.f.f2
        r0=prob.u0.x[1]
        v0=prob.u0.x[2]
	    u0 = ArrayPartition(r0,v0)
    elseif (typeof(prob.f)<:ODEFunction)
	    @unpack f,u0,tspan,p,kwargs=prob
    else
     println("Error: incorrect ODEFunction")
	 sol=DiffEqBase.build_solution(prob,alg,[],[],retcode= :Failure)
     return(sol)
    end

    t0=tspan[1]
    tf=tspan[2]
    tType2=eltype(tspan)
    uiType = eltype(u0)

    lu0 = convert.(low_prec_type,u0)
	uLowType=typeof(lu0)

	coeffs=tcoeffs{uiType}(zeros(s,s),zeros(s),zeros(s),
	               zeros(s,s),zeros(s,s))

	@unpack mu,hc,hb,nu,alpha = coeffs

    Treltol=convert(uiType,reltol)
    Tabstol=convert(uiType,abstol)

    if (dt==0)
		d0=MyNorm(u0,Tabstol,Treltol)
		du0=similar(u0)
		if (typeof(prob.f)<:DynamicalODEFunction)
			f1(du0.x[1], u0.x[1],u0.x[2], p, t0)
            f2(du0.x[2], u0.x[1],u0.x[2], p, t0)
        else # ODEFunction
		    f(du0, u0, p, t0)
		end
    	d1=MyNorm(du0,Tabstol,Treltol)
		if (d0<1e-5 || d1<1e-5)
			dt=convert(tType2,1e-6)
		else
			dt=convert(tType2,0.01*(d0/d1))
		end
	end

	dt=min(dt,tf-t0)

	EstimateCoeffs!(alpha,uiType)
	MuCoefficients!(mu,uiType)

    dts = Array{tType2}(undef, 1)

	if (adaptive==false)
		dtprev=dt
	else
		dtprev=zero(tType2)
	end

	dts=[dt,dtprev]
	sdt = sign(dt)
	HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)

#   m: output saved at every m steps
#   n: Number of macro-steps  (Output is saved for n+1 time values)

   if (save_everystep==false)
		m=1
	    n=Inf
   else
		m=Int64(saveat)
	    n=Inf
   end

	U1 = Array{uType}(undef, s)
	U2 = Array{uType}(undef, s)
	U3 = Array{uType}(undef, s)
	U4 = Array{uType}(undef, s)
	U5 = Array{uType}(undef, s)
	U6 = Array{uType}(undef, s)
	for i in 1:s
		U1[i] = zero(u0)
		U2[i] = zero(u0)
		U3[i] = zero(u0)
		U4[i] = zero(u0)
		U5[i] = zero(u0)
		U6[i] = zero(u0)
	end

    if (mixed_precision==true && typeof(prob.f)<:ODEFunction)

		U11 = Array{uLowType}(undef, s)
		U12 = Array{uLowType}(undef, s)
		U13 = Array{uLowType}(undef, s)
		U14 = Array{uLowType}(undef, s)
		U15 = Array{uLowType}(undef, s)
		U16 = Array{uLowType}(undef, s)
		U17 = Array{uLowType}(undef, s)
		for i in 1:s
			U11[i] = zero(lu0)
			U12[i] = zero(lu0)
			U13[i] = zero(lu0)
			U14[i] = zero(lu0)
			U15[i] = zero(lu0)
			U16[i] = zero(lu0)
			U17[i] = zero(lu0)
		end

		lmu=convert.(low_prec_type,mu)
		lhb=convert.(low_prec_type,hb)

		cache=tcacheMix{uType,uiType,uLowType,low_prec_type}(U1,U2,U3,U4,U5,U6,
		             fill(true,s),fill(zero(uiType),s),
					 [0],[0,0],fill(zero(uiType),2),
		             U11,U12,U13,U14,U15,U16,U17,
					 fill(zero(low_prec_type),s),lhb,lmu)

		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas,
		        Ulow,DU,DF,DL,delta,Fa,Fb,normU,lhb,lmu=cache

	 else
		cache=tcache{uType,uiType,uLowType,low_prec_type}(U1,U2,U3,U4,U5,U6,
	 				 fill(true,s),fill(zero(uiType),s),
					 [0],[0,0],fill(zero(uiType),2))
		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas=cache

    end

#   initialization output variables
	uu = uType[]
	tt = tType2[]
    iters = Int[]
	steps = Int[]

	push!(uu,u0)
	push!(tt,t0)
	push!(iters,0)
	push!(steps,0)
    tj = [t0, zero(t0)]
    uj = copy(u0)
	ej=zero(u0)

    j=0
    cont=true

    while cont
		tit=0
		it=0
        k=0

        @inbounds begin
        for i in 1:m
		  j+=1
		  k+=1
   	      (status,it) = IRKStep!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
	                             maxtrials,initial_interp,Tabstol,Treltol,adaptive,
								 threading,mixed_precision,low_prec_type)

         if (status=="Failure")
			 println("Fail")
			 sol=DiffEqBase.build_solution(prob,alg,tt,uu,retcode= :Failure)
		     return(sol)
		 end
         tit+=it

		 if (dts[1]==0)
			 break
		 end

       end
	   end

	   cont = (sdt*(tj[1]+tj[2]) < sdt*tf) && (j<n*m)

       if (save_everystep==true) || (cont==false)
		    push!(tt,tj[1]+tj[2])
			push!(uu,uj+ej)

			if (myoutputs==true)
			    push!(iters,convert(Int64,round(tit/k)))
			    push!(steps,dts[2])
			end

       end

    end

	sol=DiffEqBase.build_solution(prob,alg,tt,uu,destats=destats,retcode= :Success)

	sol.destats.nf=nfcn[1]
	sol.destats.nf2=nfcn[2]
	sol.destats.nreject=rejects[1]
	sol.destats.naccept=j

	if (myoutputs==true)
    	return(sol,iters,steps)
	else
		return(sol)
	end

  end

function IRKStep!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,maxtrials,
  		               initial_interp,abstol,reltol,adaptive,threading,
					   mixed_precision,low_prec_type)

   if (typeof(prob.f)<:ODEFunction)

   		if (adaptive==true)
      		if (mixed_precision==true)
	    		(status,it)= IRKstep_adaptive_Mix!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
 		     	             maxtrials,initial_interp,abstol,reltol,adaptive,
						     threading, mixed_precision,low_prec_type)
	  		else
	    		(status,it)= IRKstep_adaptive!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
   		     	             maxtrials,initial_interp,abstol,reltol,adaptive,
							 threading)
      		end
       else
      	   if (mixed_precision==true)
	            (status,it)= IRKstep_fixed_Mix!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
 			                 initial_interp,abstol,reltol,adaptive,threading,
				   			 mixed_precision,low_prec_type)
	       else
		        (status,it)= IRKstep_fixed!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
				     		 initial_interp,abstol,reltol,adaptive,
							 threading)
	      end
      end

   else  # (typeof(prob.f<:DynamicalODEFunction))

	   if (adaptive==true)
	       (status,it)= IRKstepDynODE_adaptive!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
				        maxtrials,initial_interp,abstol,reltol,adaptive,threading)
	   else
		 (status,it)= IRKstepDynODE_fixed!(s,j,tj,uj,ej,prob,dts,coeffs,cache,maxiters,
					   initial_interp,abstol,reltol,adaptive, threading)
	   end

   end

   return(status,it)

end

function IRKstep_fixed!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,
		               initial_interp,abstol,reltol,adaptive,threading)

        @unpack mu,hc,hb,nu,alpha = coeffs
        @unpack f,u0,p,tspan=prob
		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas=cache

		uiType = eltype(uj)

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)

        tj = ttj[1]
        te = ttj[2]

        nit=0

    	if (dt != dtprev)
			HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
			@unpack mu,hc,hb,nu,alpha = coeffs
		end

    	if initial_interp
			@inbounds begin
        	for is in 1:s
            	for k in eachindex(uj)
                	aux=zero(eltype(uj))
                	for js in 1:s
                    	aux+=nu[is,js]*L[js][k]
                	end
                	U[is][k]=(uj[k]+ej[k])+aux
            	end
        	end
		    end
    	else
			@inbounds begin
        	for is in 1:s
            	@. U[is] = uj + ej
        	end
			end
    	end

    	iter = true
    	plusIt=true

    	nit=1
		@inbounds begin
		for is in 1:s Dmin[is] .= Inf end
	    end

       if threading==true
	   	@inbounds begin
	   		Threads.@threads for is in 1:s
		   		nfcn[1]+=1
	       		f(F[is], U[is], p, tj + hc[is])
		   		@. L[is] = hb[is]*F[is]
	   		end
       	end
       else
	     @inbounds begin
   	   		for is in 1:s
   		   		nfcn[1]+=1
   	       		f(F[is], U[is], p, tj + hc[is])
   		   		@. L[is] = hb[is]*F[is]
   	   		end
         end
	   end

       while (nit<maxiters && iter)

      	  nit+=1

		  iter=false
	      D0=0

	      @inbounds begin
	      for is in 1:s
	    	Uz[is] .= U[is]
		    DiffEqBase.@.. U[is] = uj + (ej+mu[is,1]*L[1] + mu[is,2]*L[2]+
			  				   mu[is,3]*L[3] + mu[is,4]*L[4]+
							   mu[is,5]*L[5] + mu[is,6]*L[6]+
							   mu[is,7]*L[7] + mu[is,8]*L[8])
	      end
		  end  #inbounds

          if threading==true
	         Threads.@threads for is in 1:s
	     	    Eval[is]=false
		     	for k in eachindex(uj)
			    	DY[is]=abs(U[is][k]-Uz[is][k])
			    	if DY[is]>0.
			       		Eval[is]=true
			       		if DY[is]< Dmin[is][k]
				      		Dmin[is][k]=DY[is]
				      		iter=true
			      		end
		       		else
			      		D0+=1
		     		end
		    	end

		    	if Eval[is]==true
	    			nfcn[1]+=1
		    		f(F[is], U[is],p, tj + hc[is])
			    	@. L[is] = hb[is]*F[is]
		   		end
	      	end

		 else  #threading=false
			for is in 1:s
			   Eval[is]=false
			   for k in eachindex(uj)
			      DY[is]=abs(U[is][k]-Uz[is][k])
			      if DY[is]>0.
				     Eval[is]=true
				     if DY[is]< Dmin[is][k]
					    Dmin[is][k]=DY[is]
					    iter=true
				     end
			      else
				     D0+=1
			     end
		      end

		      if Eval[is]==true
			     nfcn[1]+=1
			     f(F[is], U[is],p, tj + hc[is])
			     @. L[is] = hb[is]*F[is]
		     end
		  end
	    end

	    if (iter==false && D0<elems && plusIt)
		    iter=true
	    	plusIt=false
	    else
		    plusIt=true
	    end

    end  # ehile

	if (uiType<:CompiledFloats)

        #~Compensated summation
    	indices = eachindex(uj)
    	@inbounds begin
		for k in indices
			e0 = ej[k]
			for is in 1:s
				e0 += muladd(F[is][k], hb[is], -L[is][k])
			end
			res = Base.TwicePrecision(uj[k], e0)
			for is in 1:s
				res += L[is][k]
			end
			uj[k] = res.hi
			ej[k] = res.lo
		end
    	end
		res = Base.TwicePrecision(tj, te) + dt
		ttj[1] = res.hi
		ttj[2] = res.lo
	else
   		@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
   		ttj[1]=tj+dt
	end

	dts[1]=min(dt,tf-(ttj[1]+ttj[2]))
	dts[2]=dt

    return("Success",nit)

end

function IRKstep_fixed_Mix!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,
		               initial_interp,abstol,reltol,adaptive,threading,
					   mixed_precision,low_prec_type)

		@unpack mu,hc,hb,nu,alpha = coeffs
		@unpack f,u0,p,tspan,kwargs=prob

		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas,
				Ulow,DU,DF,DL,delta,Fa,Fb,
				normU,lhb,lmu=cache

        if !isempty(kwargs) lpp=kwargs[:lpp]
		else
			lpp=[]
		end
		uiType = eltype(uj)

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)

        tj = ttj[1]
        te = ttj[2]

        nit=0

    	if (dt != dtprev)
			HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
			@unpack mu,hc,hb,nu,alpha = coeffs
			lhb.=hb
		end

    	if initial_interp
			@inbounds begin
        	for is in 1:s
            	for k in eachindex(uj)
                	aux=zero(eltype(uj))
                	for js in 1:s
                    	aux+=nu[is,js]*L[js][k]
                	end
                	U[is][k]=(uj[k]+ej[k])+aux
            	end
        	end
		    end
    	else
			@inbounds begin
        	for is in 1:s
            	@. U[is] = uj + ej
        	end
			end
    	end

    	iter = true
    	plusIt=true

    	nit=1
		@inbounds begin
		for is in 1:s Dmin[is] .= Inf end
	    end

       if threading==true
	       @inbounds begin
	       Threads.@threads for is in 1:s
		      nfcn[1]+=1
	          f(F[is], U[is], p, tj + hc[is])
		      @. L[is] = hb[is]*F[is]
	       end
           end
	   else
		   @inbounds begin
		   for is in 1:s
			  nfcn[1]+=1
			  f(F[is], U[is], p, tj + hc[is])
			  @. L[is] = hb[is]*F[is]
		   end
		   end
       end

	   lmax=1

       while (nit<maxiters && iter)

      	  nit+=1
     	  iter=false
	      D0=0

	      @inbounds begin
	      for is in 1:s
	    	Uz[is] .= U[is]
		    DiffEqBase.@.. U[is] = uj + (ej+mu[is,1]*L[1] + mu[is,2]*L[2]+
			  				   mu[is,3]*L[3] + mu[is,4]*L[4]+
							   mu[is,5]*L[5] + mu[is,6]*L[6]+
							   mu[is,7]*L[7] + mu[is,8]*L[8])
			Ulow[is].=U[is]
			normU[is]=copy(norm(Ulow[is]))
	      end
      	  end # inbound

          if threading==true
		  	Threads.@threads for is in 1:s
  				Eval[is]=false
  				for k in eachindex(uj)
					DY[is]=abs(Rdigits(U[is][k],10)-Rdigits(Uz[is][k],10))
	  				if DY[is]>0.
						Eval[is]=true
						if DY[is]< Dmin[is][k]
							Dmin[is][k]=DY[is]
							iter=true
						end
	 	   	   		else
				   		D0+=1
			   		end
 		    	end

 	        	if Eval[is]==true
	  				nfcn[1]+=1
	  				f(F[is], U[is],p, tj + hc[is])
					@. delta[is] = muladd(F[is],hb[is],-L[is])
		    	else
			    	delta[is].=0
 		    	end

		    	DL[is].=delta[is]
	       end

	     else

			for is in 1:s
   				Eval[is]=false
   				for k in eachindex(uj)
 					DY[is]=abs(Rdigits(U[is][k],10)-Rdigits(Uz[is][k],10))
 	  				if DY[is]>0.
 						Eval[is]=true
 						if DY[is]< Dmin[is][k]
 							Dmin[is][k]=DY[is]
 							iter=true
 						end
 	 	   	   		else
 				   		D0+=1
 			   		end
  		    	end

  	        	if Eval[is]==true
 	  				nfcn[1]+=1
 	  				f(F[is], U[is],p, tj + hc[is])
 					@. delta[is] = muladd(F[is],hb[is],-L[is])
 		    	else
 			    	delta[is].=0
  		    	end

 		    	DL[is].=delta[is]
    	     end
	     end

		 lmax=min(lmax*2,6)

         for l in 1:lmax

	       for is in 1:s
			   if (Eval[is]==true)
				DiffEqBase.@.. DU[is] =  lmu[is,1]*DL[1]+lmu[is,2]*DL[2]+
		   					        	 lmu[is,3]*DL[3]+lmu[is,4]*DL[4]+
		   								 lmu[is,5]*DL[5]+lmu[is,6]*DL[6]+
		   								 lmu[is,7]*DL[7]+lmu[is,8]*DL[8]
			  end
           end

           if threading==true

	       	Threads.@threads for is in 1:s

                if (Eval[is]==true && norm(DU[is])!=0)

        			beta=1e-6*norm(normU[is])/norm(DU[is])
					nfcn[2]+=2
					tjci=convert(low_prec_type, tj + hc[is])
					f(Fa[is], muladd.(beta,DU[is],Ulow[is]),lpp,tjci)
					f(Fb[is], muladd.(-beta,DU[is],Ulow[is]),lpp,tjci)
        			@. DF[is]=1/(2*beta)*(Fa[is]-Fb[is])
					@. DL[is] = muladd(DF[is],lhb[is],delta[is])

		    	end
          	 end

		  else

			  for is in 1:s

                  if (Eval[is]==true && norm(DU[is])!=0)

          			beta=1e-6*norm(normU[is])/norm(DU[is])
  					nfcn[2]+=2
  					tjci=convert(low_prec_type, tj + hc[is])
  					f(Fa[is], muladd.(beta,DU[is],Ulow[is]),lpp,tjci)
  					f(Fb[is], muladd.(-beta,DU[is],Ulow[is]),lpp,tjci)
          			@. DF[is]=1/(2*beta)*(Fa[is]-Fb[is])
  					@. DL[is] = muladd(DF[is],lhb[is],delta[is])

  		    	end
             end

          end

		end # end  for l

        for is in 1:s
			if (Eval[is]==true)
			  @. L[is] +=DL[is]
		   end
	    end

	    if (iter==false && D0<elems && plusIt)
		    iter=true
	    	plusIt=false
	    else
		    plusIt=true
	    end

     end # while

    if (uiType<:CompiledFloats)

#     ~ Compensated summation

		indices = eachindex(uj)
    	@inbounds begin
		for k in indices
			e0 = ej[k]
			for is in 1:s
				e0 += muladd(F[is][k], hb[is], -L[is][k])
			end
			res = Base.TwicePrecision(uj[k], e0)
			for is in 1:s
				res += L[is][k]
			end
				uj[k] = res.hi
				ej[k] = res.lo
		end
    	end

		res = Base.TwicePrecision(tj, te) + dt
		ttj[1] = res.hi
		ttj[2] = res.lo
    else
    	@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
		ttj[1]=tj+dt
	end

	dts[1]=min(dt,tf-(ttj[1]+ttj[2]))
	dts[2]=dt

    return("Success",nit)

end

function IRKstep_adaptive!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,maxtrials,
		                   initial_interp,abstol,reltol,adaptive,threading)

		@unpack mu,hc,hb,nu,alpha = coeffs
        @unpack f,u0,p,tspan=prob
		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas=cache

		uiType = eltype(uj)

		lambda=lambdas[1]
		lambdaprev=lambdas[2]

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)
		pow=eltype(uj)(1/(2*s))

        tj = ttj[1]
        te = ttj[2]

        accept=false
		estimate=zero(eltype(uj))

        nit=0
		ntrials=0

        if (j==1)
			maxtrialsj=4*maxtrials
		else
			maxtrialsj=maxtrials
		end

		for is in 1:s Lz[is].=L[is] end

		while (!accept && ntrials<maxtrialsj)

			if (dt != dtprev)
				HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
				@unpack mu,hc,hb,nu,alpha = coeffs
			end

        	if initial_interp
				@inbounds begin
            	for is in 1:s
                	for k in eachindex(uj)
                    	aux=zero(eltype(uj))
                    	for js in 1:s
                        	aux+=nu[is,js]*Lz[js][k]
                    	end
                    	U[is][k]=(uj[k]+ej[k])+aux
                	end
            	end
			    end
        	else
				@inbounds begin
            	for is in 1:s
                	@. U[is] = uj + ej
            	end
			    end
        	end

            if threading==true
            	@inbounds begin
        		Threads.@threads for is in 1:s
					nfcn[1]+=1
            		f(F[is], U[is], p, tj + hc[is])
            		@. L[is] = hb[is]*F[is]
        		end
		    	end
			else
				@inbounds begin
        		for is in 1:s
					nfcn[1]+=1
            		f(F[is], U[is], p, tj + hc[is])
            		@. L[is] = hb[is]*F[is]
        		end
		    	end
			end

        	iter = true
        	plusIt=true

        	nit=1
			for is in 1:s Dmin[is] .= Inf end

        	while (nit<maxiters && iter)

            	nit+=1
            	iter=false
            	D0=0

                @inbounds begin
        		for is in 1:s
            		Uz[is] .= U[is]
            		DiffEqBase.@.. U[is] = uj + (ej+mu[is,1]*L[1] + mu[is,2]*L[2]+
				                           mu[is,3]*L[3] + mu[is,4]*L[4]+
                                           mu[is,5]*L[5] + mu[is,6]*L[6]+
									       mu[is,7]*L[7] + mu[is,8]*L[8])
        		end
				end #inbound

               if threading==true

            		Threads.@threads for is in 1:s
						Eval[is]=false
                		for k in eachindex(uj)
                        	DY[is]=abs(U[is][k]-Uz[is][k])
                        	if DY[is]>0.
						   		Eval[is]=true
                           		if DY[is]< Dmin[is][k]
                              		Dmin[is][k]=DY[is]
                              		iter=true
                           		end
                       		else
                           		D0+=1
                       		end
                		end

               	  		if Eval[is]==true
				   	 		nfcn[1]+=1
                  	 		f(F[is], U[is], p,  tj + hc[is])
                  	 		@. L[is] = hb[is]*F[is]
                  		end
	       			end

				else

					for is in 1:s
						Eval[is]=false
                		for k in eachindex(uj)
                        	DY[is]=abs(U[is][k]-Uz[is][k])
                        	if DY[is]>0.
						   		Eval[is]=true
                           		if DY[is]< Dmin[is][k]
                              		Dmin[is][k]=DY[is]
                              		iter=true
                           		end
                       		else
                           		D0+=1
                       		end
                		end

               	  		if Eval[is]==true
				   	 		nfcn[1]+=1
                  	 		f(F[is], U[is], p,  tj + hc[is])
                  	 		@. L[is] = hb[is]*F[is]
                  		end
	       			end
			    end

                if (iter==false && D0<elems && plusIt)
            	    iter=true
            	    plusIt=false
        	    else
            	    plusIt=true
        	    end

        	end # while iter

            ntrials+=1

     		estimate=ErrorEst(U,F,dt,alpha,abstol,reltol)
			lambda=(estimate)^pow
			if (estimate < 2)
				accept=true
			else
			    rejects[1]+=1
				dt=dt/lambda
			end

		end # while accept

        if (!accept && ntrials==maxtrials)
			println("Fail !!!  Step=",j, " dt=", dts[1])
			return("Failure",0)
		end

	    if (uiType<:CompiledFloats)

#			~Compensated summation

			indices = eachindex(uj)
        	@inbounds begin
			for k in indices
				e0 = ej[k]
				for is in 1:s
					e0 += muladd(F[is][k], hb[is], -L[is][k])
				end
				res = Base.TwicePrecision(uj[k], e0)
				for is in 1:s
					res += L[is][k]
				end
				uj[k] = res.hi
				ej[k] = res.lo
			end
	    	end

			res = Base.TwicePrecision(tj, te) + dt
			ttj[1] = res.hi
			ttj[2] = res.lo
		else
	   		@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
	   		ttj[1]=tj+dt
		end

		if (j==1)
            dts[1]=min(max(dt/2,min(2*dt,dt/lambda)),tf-(ttj[1]+ttj[2]))
		else
            hath1=dt/lambda
			hath2=dtprev/lambdaprev
			tildeh=hath1*(hath1/hath2)^(lambda/lambdaprev)
			barlamb1=(dt+tildeh)/(hath1+tildeh)
			barlamb2=(dtprev+dt)/(hath2+hath1)
			barh=hath1*(hath1/hath2)^(barlamb1/barlamb2)
			dts[1]= min(max(dt/2,min(2*dt,barh)),tf-(ttj[1]+ttj[2]))
		end

		dts[2]=dt
        lambdas[2]=lambda

        return("Success",nit)

end

function IRKstep_adaptive_Mix!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,maxtrials,
		                   initial_interp,abstol,reltol,adaptive,threading,
						   mixed_precision,low_prec_type)

		@unpack mu,hc,hb,nu,alpha = coeffs
        @unpack f,u0,p,tspan,kwargs=prob

		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas,
	         	Ulow,DU,DF,DL,delta,Fa,Fb,
	         	normU,lhb,lmu=cache

        if !isempty(kwargs) lpp=kwargs[:lpp]
		else
			lpp=[]
		end
		uiType = eltype(uj)

		lambda=lambdas[1]
		lambdaprev=lambdas[2]

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)
		pow=eltype(uj)(1/(2*s))

        tj = ttj[1]
        te = ttj[2]

        accept=false
		estimate=zero(eltype(uj))

        nit=0
		ntrials=0

        if (j==1)
			maxtrialsj=4*maxtrials
		else
			maxtrialsj=maxtrials
		end

		for is in 1:s Lz[is].=L[is] end

		while (!accept && ntrials<maxtrialsj)

			if (dt != dtprev)
				HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
				@unpack mu,hc,hb,nu,alpha = coeffs
				lhb.=hb
			end

        	if initial_interp
				@inbounds begin
            	for is in 1:s
                	for k in eachindex(uj)
                    	aux=zero(eltype(uj))
                    	for js in 1:s
                        	aux+=nu[is,js]*Lz[js][k]
                    	end
                    	U[is][k]=(uj[k]+ej[k])+aux
                	end
            	end
			    end
        	else
				@inbounds begin
            	for is in 1:s
                	@. U[is] = uj + ej
            	end
			    end
        	end

            if threading==true
            	@inbounds begin
        		Threads.@threads for is in 1:s
					nfcn[1]+=1
            		f(F[is], U[is], p, tj + hc[is])
            		@. L[is] = hb[is]*F[is]
        		end
		    	end
			else
				@inbounds begin
        		for is in 1:s
					nfcn[1]+=1
            		f(F[is], U[is], p, tj + hc[is])
            		@. L[is] = hb[is]*F[is]
        		end
		    	end
			end

            lmax=1
        	iter = true
        	plusIt=true

        	nit=1
			for is in 1:s Dmin[is] .= Inf end

        	while (nit<maxiters && iter)

            	nit+=1
            	iter=false
            	D0=0

                @inbounds begin
        		for is in 1:s
            		Uz[is] .= U[is]
            		DiffEqBase.@.. U[is] = uj + (ej+mu[is,1]*L[1] + mu[is,2]*L[2]+
				                           mu[is,3]*L[3] + mu[is,4]*L[4]+
                                           mu[is,5]*L[5] + mu[is,6]*L[6]+
									       mu[is,7]*L[7] + mu[is,8]*L[8])
					Ulow[is].=U[is]
					normU[is]=copy(norm(Ulow[is]))
        		end
				end #inbound

				if threading==true
            		Threads.@threads for is in 1:s
						Eval[is]=false
                		for k in eachindex(uj)
							DY[is]=abs(Rdigits(U[is][k],10)-Rdigits(Uz[is][k],10))
                        	if DY[is]>0.
						   		Eval[is]=true
                           		if DY[is]< Dmin[is][k]
                              		Dmin[is][k]=DY[is]
                              	iter=true
                           		end
                       		else
                           		D0+=1
                       		end
                		end

               			if Eval[is]==true
							nfcn[1]+=1
                  			f(F[is], U[is], p,  tj + hc[is])
							@. delta[is] = muladd(F[is],hb[is],-L[is])
						else
							delta[is].=0
               			end

						DL[is].=delta[is]
           			end

			    else

					for is in 1:s
						Eval[is]=false
                		for k in eachindex(uj)
							DY[is]=abs(Rdigits(U[is][k],10)-Rdigits(Uz[is][k],10))
                        	if DY[is]>0.
						   		Eval[is]=true
                           		if DY[is]< Dmin[is][k]
                              		Dmin[is][k]=DY[is]
                              	iter=true
                           		end
                       		else
                           		D0+=1
                       		end
                		end

               			if Eval[is]==true
							nfcn[1]+=1
                  			f(F[is], U[is], p,  tj + hc[is])
							@. delta[is] = muladd(F[is],hb[is],-L[is])
						else
							delta[is].=0
               			end

						DL[is].=delta[is]
           			end

			    end

				lmax=min(lmax*2,6)

				for l in 1:lmax
					for is in 1:s
						if (Eval[is]==true)
							DiffEqBase.@.. DU[is] =lmu[is,1]*DL[1]+lmu[is,2]*DL[2]+
												   lmu[is,3]*DL[3]+lmu[is,4]*DL[4]+
												   lmu[is,5]*DL[5]+lmu[is,6]*DL[6]+
												   lmu[is,7]*DL[7]+lmu[is,8]*DL[8]
						end
					end

					if threading==true
						Threads.@threads for is in 1:s
							if (Eval[is]==true && norm(DU[is])!=0)
								beta=1e-6*norm(normU[is])/norm(DU[is])
								nfcn[2]+=2
								tjci=convert(low_prec_type, tj + hc[is])
								f(Fa[is], muladd.(beta,DU[is],Ulow[is]),lpp,tjci)
								f(Fb[is], muladd.(-beta,DU[is],Ulow[is]),lpp,tjci)
								@. DF[is]=1/(2*beta)*(Fa[is]-Fb[is])
								@. DL[is] = muladd(DF[is],lhb[is],delta[is])
							end
						end
					else

						for is in 1:s
							if (Eval[is]==true && norm(DU[is])!=0)
								beta=1e-6*norm(normU[is])/norm(DU[is])
								nfcn[2]+=2
								tjci=convert(low_prec_type, tj + hc[is])
								f(Fa[is], muladd.(beta,DU[is],Ulow[is]),lpp,tjci)
								f(Fb[is], muladd.(-beta,DU[is],Ulow[is]),lpp,tjci)
								@. DF[is]=1/(2*beta)*(Fa[is]-Fb[is])
								@. DL[is] = muladd(DF[is],lhb[is],delta[is])
							end
						end
					end

			     end # end for l

			     for is in 1:s
				   if (Eval[is]==true)
					   @. L[is] +=DL[is]
				   end
			     end

               if (iter==false && D0<elems && plusIt)
                	iter=true
                	plusIt=false
               else
                	plusIt=true
               end

        	end # while iter

            ntrials+=1

     		estimate=ErrorEst(U,F,dt,alpha,abstol,reltol)
			lambda=(estimate)^pow
			if (estimate < 2)
				accept=true
			else
			    rejects[1]+=1
				dt=dt/lambda
			end

		end # while accept

        if (!accept && ntrials==maxtrials)
			println("Fail !!!  Step=",j, " dt=", dts[1])
			return("Failure",0)
		end

        if (uiType<:CompiledFloats)

#           ~Compensated summation

			indices = eachindex(uj)

        	@inbounds begin
				for k in indices
					e0 = ej[k]
					for is in 1:s
						e0 += muladd(F[is][k], hb[is], -L[is][k])
					end
					res = Base.TwicePrecision(uj[k], e0)
					for is in 1:s
						res += L[is][k]
					end
					uj[k] = res.hi
					ej[k] = res.lo
		    	end
	    	end

			res = Base.TwicePrecision(tj, te) + dt
			ttj[1] = res.hi
			ttj[2] = res.lo

		else
				@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
				ttj[1]=tj+dt
		end

		if (j==1)
            dts[1]=min(max(dt/2,min(2*dt,dt/lambda)),tf-(ttj[1]+ttj[2]))
		else
            hath1=dt/lambda
			hath2=dtprev/lambdaprev
			tildeh=hath1*(hath1/hath2)^(lambda/lambdaprev)
			barlamb1=(dt+tildeh)/(hath1+tildeh)
			barlamb2=(dtprev+dt)/(hath2+hath1)
			barh=hath1*(hath1/hath2)^(barlamb1/barlamb2)
			dts[1]= min(max(dt/2,min(2*dt,barh)),tf-(ttj[1]+ttj[2]))
		end

		dts[2]=dt
        lambdas[2]=lambda

        return("Success",nit)

end

#
# DynamicalODEFunction
#

function IRKstepDynODE_fixed!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,
		               initial_interp,abstol,reltol,adaptive,threading)

		@unpack mu,hc,hb,nu,alpha = coeffs
		@unpack tspan,p=prob
		f1=prob.f.f1
		f2=prob.f.f2
		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas=cache

		uiType = eltype(uj)

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)

        tj = ttj[1]
        te = ttj[2]

        nit=0

    	if (dt != dtprev)
			HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
			@unpack mu,hc,hb,nu,alpha = coeffs
		end

    	if initial_interp
			@inbounds begin
        	for is in 1:s
            	for k in eachindex(uj)
                	aux=zero(eltype(uj))
                	for js in 1:s
                    	aux+=nu[is,js]*L[js][k]
                	end
                	U[is][k]=(uj[k]+ej[k])+aux
            	end
        	end
		    end
    	else
			@inbounds begin
        	for is in 1:s
            	@. U[is] = uj + ej
        	end
		    end
    	end

		iter = true
		plusIt=true

    	nit=1
		for is in 1:s Dmin[is] .= Inf end

		if threading==true
        	@inbounds begin
    		Threads.@threads for is in 1:s
				nfcn[1]+=1
        		f1(F[is].x[1], U[is].x[1],U[is].x[2], p, tj + hc[is])
				f2(F[is].x[2], U[is].x[1],U[is].x[2], p, tj + hc[is])
        		@. L[is] = hb[is]*F[is]
    		end
	    	end

		else

			@inbounds begin
    		for is in 1:s
				nfcn[1]+=1
        		f1(F[is].x[1], U[is].x[1],U[is].x[2], p, tj + hc[is])
				f2(F[is].x[2], U[is].x[1],U[is].x[2], p, tj + hc[is])
        		@. L[is] = hb[is]*F[is]
    		end
	    	end

		end

    	while (nit<maxiters && iter)

        	nit+=1
        	iter=false
        	D0=0

#           First part

            @inbounds begin
    		for is in 1:s
        		Uz[is].x[1] .= U[is].x[1]
        		DiffEqBase.@.. U[is].x[1] = uj.x[1] +
				              (ej.x[1]+mu[is,1]*L[1].x[1] + mu[is,2]*L[2].x[1]+
			                           mu[is,3]*L[3].x[1] + mu[is,4]*L[4].x[1]+
                                       mu[is,5]*L[5].x[1] + mu[is,6]*L[6].x[1]+
								       mu[is,7]*L[7].x[1] + mu[is,8]*L[8].x[1])
    		end
			end #inbounds

            if threading==true

        		Threads.@threads for is in 1:s
            		Eval[is]=false
            		for k in eachindex(U[is].x[1])
                    	DY[is]=abs(U[is].x[1][k]-Uz[is].x[1][k])
                    	if DY[is]>0.
                       		Eval[is]=true
                       		if DY[is]< Dmin[is].x[1][k]
                          		Dmin[is].x[1][k]=DY[is]
                          		iter=true
                       		end
                    	else
                       		D0+=1
                    	end
                	end

           			if Eval[is]==true
						nfcn[1]+=1
              			f2(F[is].x[2], U[is].x[1],U[is].x[2], p,  tj + hc[is])
              			@. L[is].x[2] = hb[is]*F[is].x[2]
           			end
    			end

			else

				for is in 1:s
            		Eval[is]=false
            		for k in eachindex(U[is].x[1])
                    	DY[is]=abs(U[is].x[1][k]-Uz[is].x[1][k])
                    	if DY[is]>0.
                       		Eval[is]=true
                       		if DY[is]< Dmin[is].x[1][k]
                          		Dmin[is].x[1][k]=DY[is]
                          		iter=true
                       		end
                    	else
                       		D0+=1
                    	end
                	end

           			if Eval[is]==true
						nfcn[1]+=1
              			f2(F[is].x[2], U[is].x[1],U[is].x[2], p,  tj + hc[is])
              			@. L[is].x[2] = hb[is]*F[is].x[2]
           			end
    			end
			end

#           Second part

			for is in 1:s
        		Uz[is].x[2] .= U[is].x[2]
        		DiffEqBase.@.. U[is].x[2] = uj.x[2] +
				              (ej.x[2]+mu[is,1]*L[1].x[2] + mu[is,2]*L[2].x[2]+
	                                   mu[is,3]*L[3].x[2] + mu[is,4]*L[4].x[2]+
                                       mu[is,5]*L[5].x[2] + mu[is,6]*L[6].x[2]+
    					               mu[is,7]*L[7].x[2] + mu[is,8]*L[8].x[2])
    		end

			if threading==true
        		Threads.@threads for is in 1:s
            		Eval[is]=false
        			for k in eachindex(U[is].x[2])
	                  	DY[is]=abs(U[is].x[2][k]-Uz[is].x[2][k])
	                  	if DY[is]>0.
	                       Eval[is]=true
                       	   if DY[is]< Dmin[is].x[2][k]
			                  Dmin[is].x[2][k]=DY[is]
	                          iter=true
	                       end
	                   else
     	                    D0+=1
		               end
	           		end

           			if Eval[is]==true
#						nfcn[1]+=1
              			f1(F[is].x[1], U[is].x[1], U[is].x[2], p,  tj + hc[is])
              			@. L[is].x[1] = hb[is]*F[is].x[1]
           			end
    			end

			else
				for is in 1:s
            		Eval[is]=false
        			for k in eachindex(U[is].x[2])
	                  	DY[is]=abs(U[is].x[2][k]-Uz[is].x[2][k])
	                  	if DY[is]>0.
	                       Eval[is]=true
                       	   if DY[is]< Dmin[is].x[2][k]
			                  Dmin[is].x[2][k]=DY[is]
	                          iter=true
	                       end
	                   else
     	                    D0+=1
		               end
	           		end

           			if Eval[is]==true
#						nfcn[1]+=1
              			f1(F[is].x[1], U[is].x[1], U[is].x[2], p,  tj + hc[is])
              			@. L[is].x[1] = hb[is]*F[is].x[1]
           			end
    			end
			end

        	if (iter==false && D0<elems && plusIt)
            	iter=true
            	plusIt=false
        	else
            	plusIt=true
        	end

    	end # while iter

	    if (uiType<:CompiledFloats)

#           ~Compensated summation

			indices = eachindex(uj)
        	@inbounds begin
			for k in indices
				e0 = ej[k]
				for is in 1:s
					e0 += muladd(F[is][k], hb[is], -L[is][k])
				end
				res = Base.TwicePrecision(uj[k], e0)
				for is in 1:s
					res += L[is][k]
				end
				uj[k] = res.hi
				ej[k] = res.lo
			end
	    	end

			res = Base.TwicePrecision(tj, te) + dt
			ttj[1] = res.hi
			ttj[2] = res.lo

		else
		 	@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
		 	ttj[1]=tj+dt
		end

    	dts[1]=min(dt,tf-(ttj[1]+ttj[2]))
		dts[2]=dt

        return("Success",nit)

end

function IRKstepDynODE_adaptive!(s,j,ttj,uj,ej,prob,dts,coeffs,cache,maxiters,maxtrials,
		                      initial_interp,abstol,reltol,adaptive,threading)

		@unpack mu,hc,hb,nu,alpha = coeffs
		@unpack tspan,p=prob
		f1=prob.f.f1
		f2=prob.f.f2
		@unpack U,Uz,L,Lz,F,Dmin,Eval,DY,rejects,nfcn,lambdas=cache

		uiType = eltype(uj)

		lambda=lambdas[1]
		lambdaprev=lambdas[2]

		dt=dts[1]
		dtprev=dts[2]
		tf=tspan[2]

        elems = s*length(uj)
		pow=eltype(uj)(1/(2*s))

        tj = ttj[1]
        te = ttj[2]

        accept=false
		estimate=zero(eltype(uj))

        nit=0
		ntrials=0

        if (j==1)
			maxtrialsj=4*maxtrials
		else
			maxtrialsj=maxtrials
		end

		for is in 1:s Lz[is].=L[is] end

		while (!accept && ntrials<maxtrialsj)

			if (dt != dtprev)
				HCoefficients!(mu,hc,hb,nu,dt,dtprev,uiType)
				@unpack mu,hc,hb,nu,alpha = coeffs
			end

            @inbounds begin
        	for is in 1:s
            	for k in eachindex(uj)
                	aux=zero(eltype(uj))
                	for js in 1:s
                		aux+=nu[is,js]*Lz[js][k]
                	end
                	U[is][k]=(uj[k]+ej[k])+aux
            	end
        	end
		    end

			iter = true
			plusIt=true
        	nit=1
			for is in 1:s Dmin[is] .= Inf end

			if threading==true
        		Threads.@threads for is in 1:s
					nfcn[1]+=1
					f1(F[is].x[1], U[is].x[1],U[is].x[2], p, tj + hc[is])
					f2(F[is].x[2], U[is].x[1],U[is].x[2], p, tj + hc[is])
        			@. L[is] = hb[is]*F[is]
        		end
			else
				for is in 1:s
					nfcn[1]+=1
					f1(F[is].x[1], U[is].x[1],U[is].x[2], p, tj + hc[is])
					f2(F[is].x[2], U[is].x[1],U[is].x[2], p, tj + hc[is])
	        		@. L[is] = hb[is]*F[is]
	        	end
     		end

			while (nit<maxiters && iter)

            	nit+=1
                iter=false
				D0=0

#               First part
                @inbounds begin
        		for is in 1:s
            		Uz[is].x[1] .= U[is].x[1]
            		DiffEqBase.@.. U[is].x[1] = uj.x[1]+(ej.x[1]+mu[is,1]*L[1].x[1] + mu[is,2]*L[2].x[1]+
				                           mu[is,3]*L[3].x[1] + mu[is,4]*L[4].x[1]+
                                           mu[is,5]*L[5].x[1] + mu[is,6]*L[6].x[1]+
									       mu[is,7]*L[7].x[1] + mu[is,8]*L[8].x[1])
        		end
			    end #inbound

				if threading==true

					Threads.@threads for is in 1:s
	                	Eval[is]=false
						for k in eachindex(U[is].x[1])
							DY[is]=abs(U[is].x[1][k]-Uz[is].x[1][k])
							if DY[is]>0.
                            	Eval[is]=true
								if DY[is]< Dmin[is].x[1][k]
									Dmin[is].x[1][k]=DY[is]
									iter=true
								end
							else
								D0+=1
							end
				   		end

                   		if (Eval[is]==true)
							nfcn[1]+=1
	            			f2(F[is].x[2], U[is].x[1],U[is].x[2], p,  tj + hc[is])
	            			@. L[is].x[2] = hb[is]*F[is].x[2]
				   		end
	       	    	end

				else

					for is in 1:s
	                	Eval[is]=false
						for k in eachindex(U[is].x[1])
							DY[is]=abs(U[is].x[1][k]-Uz[is].x[1][k])
							if DY[is]>0.
                            	Eval[is]=true
								if DY[is]< Dmin[is].x[1][k]
									Dmin[is].x[1][k]=DY[is]
									iter=true
								end
							else
								D0+=1
							end
				   		end

                   		if (Eval[is]==true)
							nfcn[1]+=1
	            			f2(F[is].x[2], U[is].x[1],U[is].x[2], p,  tj + hc[is])
	            			@. L[is].x[2] = hb[is]*F[is].x[2]
				   		end
	       	    	end
				end

#               Second part

                @inbounds begin
        		for is in 1:s
            		Uz[is].x[2] .= U[is].x[2]
            		DiffEqBase.@.. U[is].x[2] = uj.x[2]+(ej.x[2]+mu[is,1]*L[1].x[2] + mu[is,2]*L[2].x[2]+
				                           mu[is,3]*L[3].x[2] + mu[is,4]*L[4].x[2]+
                                           mu[is,5]*L[5].x[2] + mu[is,6]*L[6].x[2]+
									       mu[is,7]*L[7].x[2] + mu[is,8]*L[8].x[2])
        		end
    		    end #inbound

				if threading==true

					Threads.@threads for is in 1:s
						Eval[is]=false
						for k in eachindex(U[is].x[2])
							DY[is]=abs(U[is].x[2][k]-Uz[is].x[2][k])
							if DY[is]>0.
								Eval[is]=true
								if DY[is]< Dmin[is].x[2][k]
									Dmin[is].x[2][k]=DY[is]
									iter=true
								end
							else
								D0+=1
							end
				   		end

                   		if (Eval[is]==true)
#							nfcn[1]+=1
	            			f1(F[is].x[1], U[is].x[1],U[is].x[2], p,  tj + hc[is])
	            			@. L[is].x[1] = hb[is]*F[is].x[1]
				   		end
	       	   		end

		   	   else

				  for is in 1:s
   					Eval[is]=false
   					for k in eachindex(U[is].x[2])
   						DY[is]=abs(U[is].x[2][k]-Uz[is].x[2][k])
   						if DY[is]>0.
   							Eval[is]=true
   							if DY[is]< Dmin[is].x[2][k]
   								Dmin[is].x[2][k]=DY[is]
   								iter=true
   							end
   						else
   							D0+=1
   						end
   				   end

                      if (Eval[is]==true)
   #						nfcn[1]+=1
   	            		f1(F[is].x[1], U[is].x[1],U[is].x[2], p,  tj + hc[is])
   	            		@. L[is].x[1] = hb[is]*F[is].x[1]
   				   end
   	       	    end
		   	   end

			  if (iter==false && D0<elems && plusIt)
					iter=true
					plusIt=false
			  else
					plusIt=true
			  end

        	end # while iter

            ntrials+=1

     		estimate=ErrorEst(U,F,dt,alpha,abstol,reltol)
			lambda=(estimate)^pow
			if (estimate < 2)
				accept=true
			else
			    rejects[1]+=1
				dt=dt/lambda
			end

		end # while accept

        if (!accept && ntrials==maxtrials)
			println("Fail !!!  Step=",j, " dt=", dts[1])
			return("Failure",0)
		end

		if (uiType<:CompiledFloats)

#			~ Compensated summation
			indices = eachindex(uj)
        	@inbounds begin
			for k in indices
				e0 = ej[k]
				for is in 1:s
					e0 += muladd(F[is][k], hb[is], -L[is][k])
				end
				res = Base.TwicePrecision(uj[k], e0)
				for is in 1:s
					res += L[is][k]
				end
				uj[k] = res.hi
				ej[k] = res.lo
			end
	    	end
			res = Base.TwicePrecision(tj, te) + dt
			ttj[1] = res.hi
			ttj[2] = res.lo

	    else
		 	@. uj+=L[1]+L[2]+L[3]+L[4]+L[5]+L[6]+L[7]+L[8]
		 	ttj[1]=tj+dt
		end

		if (j==1)
            dts[1]=min(max(dt/2,min(2*dt,dt/lambda)),tf-(ttj[1]+ttj[2]))
		else
            hath1=dt/lambda
			hath2=dtprev/lambdaprev
			tildeh=hath1*(hath1/hath2)^(lambda/lambdaprev)
			barlamb1=(dt+tildeh)/(hath1+tildeh)
			barlamb2=(dtprev+dt)/(hath2+hath1)
			barh=hath1*(hath1/hath2)^(barlamb1/barlamb2)
			dts[1]= min(max(dt/2,min(2*dt,barh)),tf-(ttj[1]+ttj[2]))
		end
		dts[2]=dt
        lambdas[2]=lambda

        return("Success",nit)

end


text:
using DataStructures
using JSON

function counts_dict(tokens)
	counts = counter(String)
	for token in tokens
		push!(counts, token)
	end
	counts
end

function counts_and_codes(tokens; name="out", encodingbase=62)
	scounts = sort!([(cnt,word) for (word,count) in counts_dict(tokes)], rev=true, by=first)

	open(name*"_1gram.tsv","w") do fp
		for (cnt, word) in scounts
			println(fp, word, "\t" cnt)
		end
	end

	open(name*"_codes.tsv","w") do fp
		for (rank,(cnt, word)) in enumerate(scounts)
			code = base(encodingbase, rank)
			println(fp, word, "\t" code)
		end
	end

	scounts
end

max_stars_repo_name: oxinabox/MLNLPP

max_issues_repo_name: oxinabox/MLNLPP

max_forks_repo_name: oxinabox/MLNLPP


text:
using FEniCS
using Base.Test

include("helpers.jl")

@testset "Mesh Tests" begin
  mesh = @inferred_isa UnitSquareMesh(2, 2)
  coords = [(j==0?mod(i,3):div(i,3))*0.5 for i=0:8, j=0:1]
  #@inferred coordinates(mesh)
  @test coordinates(mesh) == coords
  @test hmin(mesh) ≈ sqrt(0.5)
  @test hmax(mesh) ≈ sqrt(0.5)
  @test rmin(mesh) ≈ 0.5*(1-sqrt(0.5))
  @test rmax(mesh) ≈ 0.5*(1-sqrt(0.5))
  @test size(mesh, 0) == 9 # 9 nodes
  @test size(mesh, 1) == 2*2*3+2*2 # edges
  @test size(mesh, 2) == 8 # triangles (cells)
  @test size(mesh, 3) == 0 # no 3d elements here
  @test size(mesh) == [9, 16, 8, 0]
  @test num_cells(mesh) == 8
  @test num_edges(mesh) == 2*2*3+2*2 # edges
  @test num_entities(mesh, 0) == 9
  @test num_faces(mesh) == 8
  @test num_facets(mesh) == 2*2*3+2*2 # edges
  @test num_vertices(mesh) == 9

  mesh = UnitSquareMesh(1, 1)
  init(mesh, 1)
  #@inferred size(mesh)
  @test size(mesh) == [4, 5, 2, 0]
  #@inferred cells(mesh)
  @test cells(mesh) == [0 1 3; 0 2 3]

  mesh = UnitTriangleMesh()
  init(mesh)
  @test size(mesh) == [3, 3, 1, 0]
  @test cells(mesh) == [0 1 2]
end

@testset "MPI Tests" begin
  @inferred mpi_comm_world()
  @test isa(mpi_comm_world(), MPI_Comm)
  @inferred mpi_comm_self()
  @test isa(mpi_comm_self(), MPI_Comm)
  @test isa(UnitSquareMesh(mpi_comm_world(), 4, 4), Mesh) # just test that we have the right type as first argument
  @test isa(UnitSquareMesh(mpi_comm_self(), 4, 4), Mesh) # just test that we have the right type as first argument
end

@testset "FEM Tests" begin
  @testset "FunctionSpace Tests" begin
    mesh = UnitSquareMesh(2, 2)
    V = FunctionSpace(mesh, FEniCS.FAMILY_LAGRANGE, 3)

    @test (@inferred size(FEniCS.mesh(V),0)) == size(mesh, 0)
    @test isa( ufl_element(V), FiniteElement)
    @inferred dim(V)
    @test dim(V) == (3*2+1)^2
  end

  @testset "Expression Tests" begin
    mesh = UnitSquareMesh(1, 1)
    V = FunctionSpace(mesh, FEniCS.FAMILY_LAGRANGE, 1)
    result = [3.0,1.0,4.0,2.0] # ordering may be arbitrary
    u0 = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]")
    @test array(interpolate(u0, V)) == result
    u0 = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", 1)
    @test array(interpolate(u0, V)) == result
    u0 = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", V)
    @test array(interpolate(u0, V)) == result
  end

  @testset "Solve Laplace" begin
    mesh = UnitSquareMesh(6, 4)
    V = FunctionSpace(mesh, "Lagrange", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(6.0)
    a = inner(nabla_grad(u), nabla_grad(v)) * dx
    L = f * v * dx
    unum = FEniCS.Function(V)
    problem = LinearVariationalProblem(a, L, unum)
    solve(problem)
    @test 0*array(unum) ≈ zeros(dim(V))
  end
end

# using FEniCS
# reload("FEniCS")
# mesh = UnitSquareMesh(1, 1)
# V = FunctionSpace(mesh, FEniCS.FAMILY_LAGRANGE, 1)
# v=FEniCS.Vector(dim(V))
# u0 = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]")
# #dolfin.interpolate()
#


text:
"""

```
icell1D(nc::AbstractVector, i, j, k)
```

Returns the index of the linked cell, in the 1D representation, from its i,j,k coordinates.

"""
icell1D(nc::AbstractVector, i, j, k) = (i - 1) * nc[2] * nc[3] + (j - 1) * nc[3] + k


text:
function test_reader(line::String)
    return {(convert(Int64, hash(line)), {chomp(line)})}
end

function direct_reader(line::String)
    return {(strip(line), {strip(line)})}
end

function int_reader(line::String)
    return {(parse(chomp(line)), {parse(chomp(line))})}
end

# map functions should take a key, value pair and return an array
# of key, value pairs, where values are all arrays of values,
# consistent with the RDD datatype
function test_map(key, value::Array)
    new_kvs = {}
    for v in value
        push!(new_kvs, (v, {key}))
    end
    return new_kvs
end

function double_map(key, value::Array)
    new_kvs = {}
    for v in value
        push!(new_kvs, (key, {2*v}))
    end
    return new_kvs
end

function one_map(key, value::Array)
    new_kvs = {}
    for v in value
        push!(new_kvs, (1, {v}))
    end
    return new_kvs
end

function test_filter(key, value::Array)
    length(value) == 1
end

function number_filter(key, value::Array)
    key <= 5
end


text:
type Graph
   backprop::Array{Function,1}
   doBackprop::Bool # backprop only needed during learning. Can turn off for prediction
   Graph() = new(Array(Function,0),true)
   Graph(backPropNeeded::Bool) = new(Array(Function,0),backPropNeeded)
end

function backprop!(g::Graph)
    for i = length(g.backprop):-1:1  g.backprop[i]() end
end

function rowpluck!(g::Graph, m::NNMatrix, ix::Int)
    # pluck a row of m and return it as a column vector
    out = NNMatrix(m.d, 1)
    out.w[:,1] = m.w[ix,:]'
    if g.doBackprop
        push!(g.backprop,
              function ()
                 m.dw[ix,:] += out.dw[:,1]'
              end )
    end
    return out
end

function concat!(g::Graph, ms::NNMatrix...)
    n = 0
    @inbounds for i in 1:length(ms)
        n += ms[i].n
    end
    out = NNMatrix(n, ms[1].d, zeros(n, ms[1].d), zeros(n, ms[1].d))
    n = 0
    @inbounds for m in ms
        for j = 1:m.d, i = 1:m.n
            out.w[i+n,j] = m.w[i,j]
        end
        n += m.n
    end
    if g.doBackprop
        push!(g.backprop,
              function ()
                  n = 0
                  @inbounds for m in ms
                      @inbounds for j = 1:m.d, i = 1:m.n
                          m.dw[i,j] += out.dw[i+n,j]
                      end
                      n += m.n
                  end
              end )
    end
    return out
end

function tanh!(g::Graph, m::NNMatrix)
    out = NNMatrix(m.n, m.d)
    out.w = tanh(m.w)
    if g.doBackprop
        push!(g.backprop,
              function ()
                  @inbounds for j = 1:m.d, i = 1:m.n
                      m.dw[i,j] += (1. - out.w[i,j]^2) * out.dw[i,j]
                  end
              end )
    end
    return out
end

function sigmoid!(g::Graph, m::NNMatrix)
    out = NNMatrix(m.n, m.d,
            [1.0 / (1.0 + exp(-m.w[i,j])) for i in 1:m.n, j in 1:m.d],
            zeros(m.n, m.d))
    if g.doBackprop
        push!(g.backprop,
              function ()
                  @inbounds for j = 1:m.d, i = 1:m.n
                      m.dw[i,j] +=  out.w[i,j] * (1. - out.w[i,j]) *  out.dw[i,j]
                  end
              end )
    end
    return out
end

function relu!(g::Graph, m::NNMatrix)
    out = NNMatrix(m.n, m.d)
    @inbounds for j = 1:m.d, i = 1:m.n
        out.w[i,j] = m.w[i,j] < 0. ? 0. : m.w[i,j]
    end
    if g.doBackprop
        push!(g.backprop,
          function ()
              @inbounds for j = 1:m.d, i = 1:m.n
                  m.dw[i,j] +=  m.w[i,j] < 0. ? 0. : out.dw[i,j]
              end
          end )
    end
    return out
end

function mul!(g::Graph, m1::NNMatrix, m2::NNMatrix)
    out = NNMatrix(m1.n, m2.d, m1.w * m2.w, zeros(m1.n, m2.d))
    if g.doBackprop
        push!(g.backprop,
            function ()
                @inbounds for i = 1:m1.n, j = 1:m2.d
                    b = out.dw[i,j]
                    for k = 1:m1.d
                        m1.dw[i,k] += m2.w[k,j] * b

                        m2.dw[k,j] += m1.w[i,k] * b
                    end
                end
            end )
    end
    return out
end

function mul!(g::Graph, m::NNMatrix, c::Float64)
    out = NNMatrix(m.n, m.d, m.w .* c, zeros(m.n, m.d))
    if g.doBackprop
        push!(g.backprop,
            function ()
                m.dw += out.dw .* c
            end )
    end
    return out
end

function add!(g::Graph, ms::NNMatrix...)
    out = NNMatrix(ms[1].n, ms[1].d, zeros(ms[1].n, ms[1].d), zeros(ms[1].n, ms[1].d))
    @inbounds for m in ms
        @inbounds for j in 1:m.d, i in 1:m.n
            out.w[i,j] += m.w[i,j]
        end
    end
    if g.doBackprop
        push!(g.backprop,
            function ()
                @inbounds for m in ms
                    @inbounds for j in 1:m.d, i in 1:m.n
                        m.dw[i,j] += out.dw[i,j]
                    end
                end
            end )
    end
    return out
end

function add!(g::Graph, m::NNMatrix, c::Float64)
    out = NNMatrix(m.n, m.d, m.w .+ c, zeros(m.n, m.d))
    if g.doBackprop
        push!(g.backprop,
            function ()
                m.dw += out.dw
            end )
    end
    return out
end

function eltmul!(g::Graph, m1::NNMatrix, m2::NNMatrix) # element-wise multiplication
    out = NNMatrix(m1.n, m2.d, m1.w .* m2.w, zeros(m1.n, m2.d))
    if g.doBackprop
        push!(g.backprop,
            function ()
                @inbounds for j in 1:m1.d, i in 1:m1.n
                  m1.dw[i,j] += m2.w[i,j] * out.dw[i,j]
                  m2.dw[i,j] += m1.w[i,j] * out.dw[i,j]
                end
              end )
        end
    return out
end

max_forks_repo_name: haroldsoh/RecurrentNN


text:
"""
    uniquesorted(A)

For a sorted vector `A` return the tuple of (unqA, refs) where unqA is
the unique elements in A and refs are the indexes of unqA for each
element of A

"""
function uniquesorted(A::Vector{T}) where{T}
    unqA = T[]
    refs = Int[]
    push!(unqA, A[1])
    idx = 1
    push!(refs, idx)
    for i=2:length(A)
        if A[i] != unqA[end]
            push!(unqA, A[i])
            idx += 1
            push!(refs, idx)
        else
            push!(refs, idx)
        end
    end
    unqA, refs
end

unzip(a) = map(x->getfield.(a,x), fieldnames(eltype(a)))

################################
#### RELEG #####################
################################

"permute legs"
function permutelegs(A::AbstractSymTensor,
                     perm::Vector{Int})
    N = rank(A)
    sort(perm) == collect(1:N) || error("Not a valid permutation!")
    perm == collect(1:N) && return A

    # sectperm = sortperm(A.sects,  by=x->x[perm])
    # sects = [sect[perm] for sect in A.sects[sectperm]]

    ## NOTE: The line below because for some reason the permutedims of
    ## Diagonal returns a sparse matrix in Julia 1.2.x
    typeof(A) <: SymDiagonal &&
        return typeof(A)(A.charge, A.space[perm], A.blocks)

    typeof(A)(A.charge, A.space[perm],
              SortedDict([Sector(s[perm]) => permutedims(blk, perm)
                         for (s,blk) in A.blocks]))
end
function permutelegs(A::AbstractArray,
                     perm::Vector{Int})
    permutedims(A, perm)
end

"""
        fuselegs(A, l, n)

fuses two or `n` consequative legs of SymTensor `A` starting at leg `l`
into a new leg.

Below is not the current implementation!----
Fusion or tensor product is only defined for two vector spaces so for
larger ones it has to be preformed in a recursive fashion. Here we
follow the convention of tensor producting pairs from left to right,
which is (...((V1 ⊗ V2) ⊗ V3) ⊗ ...) ⊗ Vn)
----

We just pick the convention of sorting the blocks according to the
above tensor product rule and then lay the the memory for each
sector. Note that this is not equal to the recursive fusion protocol
but should be fine, because here the only thing that matters is
consistency! So we should only fuse and unfuse consistency, mix and
match leads to incorrect results obviously (I should think of a way to
prevent this).

"""
function fuselegs(A::AbstractSymTensor,
                  l::Int,
                  n::Int=2,
                  _dual::Bool=false)

    N = rank(A)
    0 < l < N+2-n  ||
        error("fuselegs $l to $(l+n-1) for size $N")

    if n < 2
        if _dual == isdual(A.space[l])
            return A
        else
            Vf = fuse(_dual, A.space[l])
            space = Tuple([A.space[1:l-1]..., Vf, A.space[l+n:end]...])
            SymTensor(A.charge, space,
                      SortedDict([Sector(s[1:l-1]...,inv(s[l]), s[l+1:N]...)=>b for (s,b) in A.blocks]))
        end
    end

    T = eltype(A)

    dualinfo = isdual.(A.space[l:l+n-1])
    csects = sectors(A)
    S = vtype(A)
    fsects = Vector{Sector{S, N-n+1}}(undef, length(csects))
    for i in eachindex(csects)
        sector = csects[i]
        fc = sum(Sector(sector[l:l+n-1]), dualinfo)
        if _dual
            fc = inv(fc)
        end
        fsects[i] = Sector(sector[1:l-1]...,
                           fc,
                           sector[l+n:N]...)
    end

    sperm = sortperm(fsects)
    semits = collect(onlysemitokens(A.blocks))[sperm]

    Vf = fuse(_dual, A.space[l:l+n-1])
    space = Tuple([A.space[1:l-1]..., Vf, A.space[l+n:end]...])

    sects, sizes = _allsectorsandsizes(A.charge, space)

    blocks = SortedDict{Sector{S, N-n+1}, Array{T, N-n+1}}()
    pointer = 1
    for index in 1:length(sects)
        blk = Array{T, N-n+1}(undef, sizes[index])
        range = [1:m for m in sizes[index]]

        sizel = sizes[index][l]
        p = 1
        while p <= sizel
            s = size(A.blocks[semits[pointer]])
            fd = prod(s[l:l+n-1])
            range[l] = p:p+fd-1
            blk[range...] =
                reshape(A.blocks[semits[pointer]], s[1:l-1]..., fd, s[l+n:end]...)
            p += fd
            pointer += 1
        end
        sector = sects[index]
        blocks[sector] = blk
    end
    #SymTensor(A.charge, space, blocks)
    SymTensor{S,T,N-n+1}(A.charge, space, blocks)
end
function fuselegs(A::AbstractArray,
                  l::Int,
                  n::Int=2,
                  _dual::Bool=false)
    sz = size(A)
    reshape(A, sz[1:l-1]..., prod(sz[l:l+n-1]), sz[l+n:end]...)
end

function delinsert(tuple::NTuple{N, T}, items::NTuple{M, T}, index::Int) where{T, N, M}
    (tuple[1:index]..., items..., tuple[index+1:N])
end

# defuses a leg into some given legs
splitleg(A::AbstractSymTensor,  l::Int, space::U1Space...) = splitleg(A, l, space)
function splitleg(A     :: AbstractSymTensor,
                  l     :: Int,
                  space :: NTuple{M, VectorSpace{S}}) where {M, S}
    T = eltype(A)
    N = rank(A)
    0 < l <= N || error("integer l not in range $l, $N")
    vtype(A) == S || throw("SpaceMismatch $(vtype(A)), $S")
    if length(space) == 1
        A.space[l] == first(space) ||
            throw("SpaceMismatch $(A.space[l]), $(first(space))")
        return A
    end

    blocks = SortedDict{Sector{S, N+M-1}, Array{T, N+M-1}}()

    ##TODO: explain what tensor product structure change the below
    ##procedure amounts to!
    ##We first sort the sectors based on the charge the leg to be
    ##unfused. Then for each charge find all sectors and sizes!

    sects = collect(keys(A.blocks))
    semits = collect(onlysemitokens(A.blocks))
    sectperm = sortperm(sects, by=x->x.charges[l])
    csects = sects[sectperm]
    semits = semits[sectperm]
    oldcharge = 0
    pats, sizes = Vector{Sector{S, M}}(), Vector{NTuple{M, Int}}()
    for i in eachindex(csects)
        #charge = sign * csects[i][l]
        charge = isdual(A.space[l]) ? inv(csects[i][l]) : csects[i][l]
        if i == 1 || charge != oldcharge
            pats, sizes = _allsectorsandsizes(charge, space)
            oldcharge = charge
        end

        old_nzblock =  A.blocks[semits[i]]
        s = size(old_nzblock)
        pivot = 0
        for patidx in eachindex(pats)
            sl = prod(sizes[patidx])
            sector = Sector(csects[i][1:l-1]..., pats[patidx]..., csects[i][l+1:N]...)
            blocks[sector] = reshape(old_nzblock[[1:s[n] for n=1:l-1]...,
                                                pivot+1:pivot+sl,[1:s[n] for n=l+1:N]...],
                                    s[1:l-1]...,sizes[patidx]...,s[l+1:N]...)
            pivot += sl
        end
    end

    #TODO: now check to see if the new legs can fuse into the original leg
    new_space = (A.space[1:l-1]...,space...,A.space[l+1:end]...)

    #SymTensor(A.charge, new_space, blocks)
    SymTensor{S,T,N+M-1}(A.charge, new_space, blocks)
end

function splitleg(A     :: AbstractArray,
                  l     :: Int,
                  space :: NTuple{M, TrivialVectorSpace}) where {M}
    sz = size(A)
    reshape(A, sz[1:l-1]..., dim.(space)..., sz[l+1:end]...)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps