text:
"""
Example:
'''
using FileToRigidTransform; FileToRigidTransform.subscribe()
'''
Follow directions, or update configuration files in ~/.julia_hid/ConfigFo/
"""
module FileToRigidTransform
import FileWatching.watch_file
import Base.getindex

const FOLDER = joinpath(homedir(), ".julia_hid")
const CONFIGFO = joinpath(FOLDER, "ConfigFo")
const TRANSFORMFO = joinpath(FOLDER, "TransformFo")
const TIMEOUT = 30.0 # seconds
const RETRY_TIME = 0.5 # seconds
const TIMEZERO = time()

"String descriptors for a device from file name"
struct DevConfig
    filename::String
    specbynamedict::Dict{String, NamedTuple}
end
DevConfig()= DevConfig("", Dict{String, Int}())

"""A structure for 'show' dispatching. Could
otherwise use Base.CodeUnits <: DenseArray, but we don't
need the functionality.
"""
struct Bytevector
    data::Vector{UInt8}
end
Bytevector() = Bytevector(UInt8[])
Bytevector(v::Array{Int8,1 }) = Bytevector(Array{UInt8,1}(v))
Base.getindex(bv::Bytevector, i::Int) = getindex(bv.data, i)
Base.length(bv::Bytevector) = length(bv.data)
struct DevState
    values::Bytevector
    timestamp::Float64
    proceed::Bool
end
DevState() = DevState(Bytevector(), localtime(), true)
DevState(v::Array, t, p) = DevState(Bytevector(v), t, p)

include("show_devstate.jl")
include("extract_from_rawdata.jl")
include("loggers.jl")
include("file_operations.jl")

"""
    subscribe(timeout = TIMEOUT, logger = log_by_channel) -> Vector{Task}

Start reacting to changes in device(s) state as given and updated in files.
The files which are present in a folder structure defines how the state is interpreted.

Default behaviour is logging the interpreted device state to a file and to screen, see keyword arguments.

Tip: To see error messages from failed tasks, display them individually, e.g.
    julia> ans[1]
"""
function subscribe(;timeout = TIMEOUT, logger = log_by_channel)
    devconfvec = devices_configuration_vector()
    monitor(devconfvec, timeout, logger)
end

"""
    devices_configuration_vector() -> Vector{DevConfig}

Use folder structure to locate:
- devices monintored (subscribed), i.e. another process is updating a file with device state
- device configuration, i.e. a text file where the user can describe the conversion from device state to input variables for transformations

If a configuration file is missing, a template file is generated based on the number of bytes in the device state file.

Give warnings when folder structure is missing. Creating folders programmatically might cause difficulties concerning folder ownership.
"""
function devices_configuration_vector()
    dev_conf_vec = Vector{DevConfig}()
    if ispath(FOLDER)
        filenames = filter(isfile, readdir(FOLDER, join=true))
        shortfilenames = filenames .|> splitpath .|> last
        if length(filenames) != 0
            for (fi, shortfi) in zip(filenames, shortfilenames)
                if ispath(CONFIGFO)
                    configfi = joinpath(CONFIGFO, shortfi)
                    if isfile(configfi)
                        # Found both a device state and a device configuration file.
                        chdi = specbyname(configfi)
                        if length(chdi) > 0
                            devconfig = DevConfig(shortfi, chdi)
                            push!(dev_conf_vec, devconfig)
                            if !ispath(TRANSFORMFO)
                                @warn("Please create folder $TRANSFORMFO, then rerun")
                                #return dev_conf_vec
                            else
                                @info configfi
                            end
                        else
                            @warn "Empty device configuration $shortfi"
                        end
                    else
                        @info("Creating template configuration file\n\t$configfi")
                        sourcefile = joinpath(@__DIR__, "..", "example", "config_template.txt")
                        @assert isfile(sourcefile) joinpath(pwd(), sourcefile)
                        cp(sourcefile, configfi, force= false)
                        @warn "Retry with new device configuration by rerunning FileToRigidTransform.subscribe()"
                        return dev_conf_vec
                    end
                else
                    @warn("Please create folder $CONFIGFO, then rerun")
                    return dev_conf_vec
                end
            end
            @info("Configured $(length(dev_conf_vec)) usb pipelines")
        else
            @info "Could not find a *.txt file in $fo .\n\tGenerate one using WinControllerToFile."
        end
    else
        @info("You need to create folder $fo")
        return dev_conf_vec
    end
    dev_conf_vec
end

"""
    devstate_updated(timeout, filename) -> DevState

Read state from file when file is updated.
If the file is (momentarily) unreadable, retry during RETRY_TIME.
"""
function devstate_updated(timeout, filename)
    prevstate = devstate(filename)
    # Yield to other tasks while waiting for file change
    fileevent = watch_file(filename, timeout)  # does not exit by itself after timeout actually.
    if fileevent.renamed || fileevent.timedout
        return DevState(prevstate.values, prevstate.timestamp, false)
    end

    ds = devstate(filename)
    t0 = localtime()
    while length(ds.values) == 0 && localtime() - t0 < RETRY_TIME
        # May occur if we are reading before the file is closed for writing. Retry twice before accepting.
        sleep(0.01)
        ds = devstate(filename)
    end
    if localtime() - t0 >= RETRY_TIME + 0.01   ##?
        @warn "No state read from \n\t$filename \n\tin $RETRY_TIME s. If WinControllerToFile.subscribe() is running, you may have to populate the file by giving input through the device"
    end
    ds
end

"""
    monitor(devconfvec, timeout, logger) -> Vector{Task}

Monitor vector of devices. The 'logger' argument is typically used for
for logging.

Arguments to 'func' are (ios::IOstream, d::DevConfig, dsprev::DevState, ds::DevState)
"""
function monitor(devconfvec, timeout, logger)
    # Create file monitors
    monitors = Vector{Task}()
    for d in devconfvec
        fina = joinpath(FOLDER, d.filename)
        monitor =  @async monitor_file(fina, d, timeout, logger)
        push!(monitors, monitor)
    end
    monitors
end

"""
Monitor a single device. The given function argument is called at every update.
Arguments to 'logger' are (ios::IOstream, d::DevConfig, dsprev::DevState, ds::DevState)
"""
function monitor_file(filename, d::DevConfig, timeout, logger)
    _, shfina = splitdir(filename)
    logfile = joinpath(TRANSFORMFO, shfina)
    t0 = localtime()
    accumulated = Vector{Float64}()
    open(logfile, write = true) do ios
        ds = DevState()
        dsprev = ds
        tpassed = localtime() - t0
        while true
            tpassed = localtime() - t0
            ds = devstate_updated(timeout-tpassed, filename)
            !dsprev.proceed && break
            accumulated = logger(ios, d, accumulated, dsprev, ds)
            dsprev = ds
            if (tpassed > timeout)
                @info "Exit monitor_file due to timeout $timeout s"
                break
            end
            flush(ios)
        end
        ds
    end
    @info "Exit logging to \n\t$logfile \n\tafter $(floor(localtime()-t0)) s"
end

localtime() = time() - TIMEZERO

end # module


text:
## Feacalc.jl Full implementation of feature calculation for speech files.
## (c) 2013-2016 David A. van Leeuwen

## Feacalc.  Feature calculation as used for speaker and language recognition.

nrow(x) = size(x,1)
ncol(x) = size(x,2)

## compute features according to standard settingsm directly from a wav file.
## this does channel at the time.
function feacalc(wavfile::AbstractString; method=:wav, kwargs...)
    if method == :wav
        (x, sr) = wavread(wavfile)
    elseif method == :sox
        (x, sr) = soxread(wavfile)
    elseif method == :sphere
        (x, sr) = sphread(wavfile)
    end
    sr = Float64(sr)       # more reasonable sr
    feacalc(x; sr=sr, source=wavfile, kwargs...)
end

## assume we have an array already
function feacalc(x::Array; augtype=:ddelta, normtype=:warp, sadtype=:energy, dynrange::Real=30., nwarp::Int=399, chan=:mono, sr::AbstractFloat=8000.0, source=":array", defaults=:nbspeaker, mfccargs...)
    if ndims(x) > 1
        nsamples, nchan = size(x)
        if chan == :mono
            x = vec(mean(x, 2))            # averave multiple channels for now
        elseif in(chan, [:a, :b])
            channum = findin([:a, :b], [chan])
            x = vec(x[:,channum])
        elseif isa(chan, Integer)
            if !(chan in 1:nchan)
                error("Bad channel specification: ", chan)
            end
            x = vec(x[:,chan])
            chan=[:a, :b][chan]
        else
            error("Unknown channel specification: ", chan)
        end
    else
        nsamples, nchan = length(x), 1
    end
    ## save some metadata
    meta = Dict("nsamples" => nsamples, "sr" => sr, "source" => source, "nchan" => nchan,
            "chan" => chan)
    preemph = 0.97
    preemph ^= 16000. / sr

    ## basic features
    (m, pspec, params) = mfcc(x, sr, defaults; preemph=preemph, mfccargs...)
    meta["totnframes"] = nrow(m)

    ## augment features
    if augtype==:delta || augtype==:ddelta
        d = deltas(m)
        if augtype==:ddelta
            dd = deltas(d)
            m = hcat(m, d, dd)
        else
            m = hcat(m, d)
        end
    elseif augtype==:sdc
        m = sdc(m)
    end
    meta["augtype"] = augtype

    if nrow(m)>0
        if sadtype==:energy
            ## integrate power
            deltaf = size(pspec,2) / (sr/2)
            minfreqi = round(Int, 300deltaf)
            maxfreqi = round(Int, 4000deltaf)
            power = 10log10(sum(pspec[:,minfreqi:maxfreqi], 2))

            maxpow = maximum(power)
            speech = find(power .> maxpow - dynrange)
            params["dynrange"] = dynrange
        elseif sadtype==:none
            speech = collect(1:nrow(m))
        end
    else
        speech=Int[]
    end
    meta["sadtype"] = sadtype
    ## perform SAD
    m = m[speech,:]
    meta["speech"] = map(UInt32, speech)
    meta["nframes"] , meta["nfea"] = size(m)

    ## normalization
    if nrow(m)>0
        if normtype==:warp
            m = warp(m, nwarp)
            params["warp"] = nwarp          # the default
        elseif normtype==:mvn
            if nrow(m)>1
                znorm!(m,1)
            else
                fill!(m, 0)
            end
        end
        meta["normtype"] = normtype
    end
    return(map(Float32, m), meta, params)
end

## When called with a specific application in mind, call with two arguments
function feacalc(wavfile::AbstractString, application::Symbol; kwargs...)
    if (application in [:speaker, :nbspeaker])
        feacalc(wavfile; defaults=:nbspeaker, kwargs...)
    elseif application==:wbspeaker
        feacalc(wavfile; defaults=:wbspeaker, kwargs...)
    elseif (application==:language)
        feacalc(wavfile; defaults=:rasta, nwarp=299, augtype=:sdc, kwargs...)
    elseif (application==:diarization)
        feacalc(wavfile; defaults=:rasta, sadtype=:none, normtype=:mvn, augtype=:none, kwargs...)
    else
        error("Unknown application ", application)
    end
end

function sad(pspec::Matrix{Float64}, sr::Float64, method=:energy; dynrange::Float64=30.)
    deltaf = size(pspec,2) / (sr/2)
    minfreqi = round(Int, 300deltaf)
    maxfreqi = round(Int, 4000deltaf)
    power = 10log10(sum(pspec[:,minfreqi:maxfreqi], 2))
    maxpow = maximum(power)
    speech = find(power .> maxpow - dynrange)
end

## listen to SAD
function sad(wavfile::AbstractString, speechout::AbstractString, silout::AbstractString; dynrange::Float64=30.)
    (x, sr, nbits) = wavread(wavfile)
    sr = float64(sr)               # more reasonable sr
    x = vec(mean(x, 2))            # averave multiple channels for now
    (m, pspec, meta) = mfcc(x, sr; preemph=0)
    sp = sad(pspec, sr, dynrange=dynrange)
    sl = round(Int, meta["steptime"] * sr)
    xi = falses(size(x))
    for i in sp
        xi[(i-1)*sl+(1:sl)] = true
    end
    y = x[find(xi)]
    wavwrite(y, speechout, Fs=sr, nbits=nbits, compression=WAVE_FORMAT_PCM)
    y = x[find(.!xi)]
    wavwrite(y, silout, Fs=sr, nbits=nbits, compression=WAVE_FORMAT_PCM)
end

## this should probably be called soxread...
function soxread(file)
    nch = parse(Int, readall(`soxi -c $file`))
    sr = parse(i=Int, readall(`soxi -r $file`))
    sox = `sox $file -t raw -e signed -r $sr -b 16 -`
    fd, proc = open(sox, "r")
    x = Int16[]
    while !eof(fd)
        push!(x, read(fd, Int16))
    end
    ns = length(x) ÷ nch
    reshape(x, nch, ns)' / (1<<15), sr
end

## similarly using sphere tools
function sphread(file)
    nch = parse(Int, readstring(`h_read -n -F channel_count $file`))
    sr = parse(Int, readstring(`h_read -n -F sample_rate $file`))
    open(pipeline(`w_decode -o pcm $file -`, `h_strip - - `), "r") do fd
        x = Int16[]
        while !eof(fd)
            push!(x, read(fd, Int16))
        end
        ns = length(x) ÷ nch
        return reshape(x, nch, ns)' / (1<<15), sr
    end
end

ext: jl


text:
@doc raw"""
Holds the tableau of a stochastic explicit Runge-Kutta method.

Order of the tableau is not included, because unlike in the deterministic
setting, it depends on the properties of the noise (e.g., the dimension of
the Wiener process and the commutativity properties of the diffusion matrix).

Orders stored in `qdrift`, `qdiff` and `qdiff2` are understood as the classical orders of these methods.
"""
struct TableauSERK{T} <: AbstractTableauERK{T}
    name::Symbol
    s::Int

    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}
    qdiff2::CoefficientsRK{T}

    function TableauSERK{T}(name, qdrift, qdiff, qdiff2) where {T}
        @assert qdrift.s == qdiff.s == qdiff2.s
        @assert qdrift.c[1] == qdiff.c[1] == qdiff2.c[1] == 0.
        @assert istrilstrict(qdrift.a)
        @assert istrilstrict(qdiff.a)
        @assert istrilstrict(qdiff2.a)
        @assert !(qdrift.s==1. && qdrift.a[1,1] ≠ 0.)
        @assert !(qdiff.s==1. && qdiff.a[1,1] ≠ 0.)
        @assert !(qdiff2.s==1. && qdiff2.a[1,1] ≠ 0.)
        new(name, qdrift.s, qdrift, qdiff, qdiff2)
    end
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}, qdiff2::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, qdiff2)
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, CoefficientsRK(T, :NULL, 0, zero(qdrift.a), zero(qdrift.b), zero(qdrift.c)) )
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T},
                                    order_diff2::Int, a_diff2::Matrix{T}, b_diff2::Vector{T}, c_diff2::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(name, order_diff2, a_diff2, b_diff2, c_diff2))
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(:NULL, 0, zero(a_drift), zero(b_drift), zero(c_drift)))
end

"Parameters for stochastic explicit Runge-Kutta methods."
struct ParametersSERK{DT, TT, D, M, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::TableauSERK{TT}
    Δt::TT

    function ParametersSERK{DT,D,M}(equs::ET, tab::TableauSERK{TT}, Δt::TT) where {DT, TT, D, M, ET <: NamedTuple}
        new{DT, TT, D, M, tab.s, ET}(equs, tab, Δt)
    end
end

"Stochastic Explicit Runge-Kutta integrator cache."
struct IntegratorCacheSERK{DT,D,M,S} <: SDEIntegratorCache{DT,D,M}
    Δy::Vector{DT}

    Q::Vector{Vector{DT}}     # Q[j][k] - the k-th component of the j-th internal stage
    V::Vector{Vector{DT}}     # V[j][k] - the k-th component of v(Q[j])
    B::Vector{Matrix{DT}}     # B[j]    - the diffusion matrix B(Q[j])

    function IntegratorCacheSERK{DT,D,M,S}() where {DT,D,M,S}
        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        B = create_internal_stage_matrix(DT, D, M, S)

        new(zeros(DT,M), Q, V, B)
    end
end

function Integrators.IntegratorCache{ST}(params::ParametersSERK{DT,TT,D,M,S}; kwargs...) where {ST,DT,TT,D,M,S}
    IntegratorCacheSERK{ST,D,M,S}(; kwargs...)
end

@inline Integrators.CacheType(ST, params::ParametersSERK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = IntegratorCacheSERK{ST,D,M,S}

"Stochastic Explicit Runge-Kutta integrator."
struct IntegratorSERK{DT, TT, D, M, S, ET <: NamedTuple} <: StochasticIntegratorRK{DT,TT,D,M,S}
    params::ParametersSERK{DT,TT,D,M,S,ET}
    caches::CacheDict{ParametersSERK{DT,TT,D,M,S,ET}}

    function IntegratorSERK(params::ParametersSERK{DT,TT,D,M,S,ET}, caches) where {DT,TT,D,M,S,ET}
        new{DT, TT, D, M, S, ET}(params, caches)
    end

    function IntegratorSERK{DT,D,M}(equations::ET, tableau::TableauSERK{TT}, Δt::TT) where {DT, TT, D, M, ET <: NamedTuple}
        # get number of stages
        S = tableau.s

        # create params
        params = ParametersSERK{DT,D,M}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create integrator
        IntegratorSERK(params, caches)
    end

    function IntegratorSERK(equation::SDE{DT,TT}, tableau::TableauSERK{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorSERK{DT, ndims(equation), equation.m}(get_function_tuple(equation), tableau, Δt; kwargs...)
    end
end

"Integrate SDE with explicit Runge-Kutta integrator."
function Integrators.integrate_step!(int::IntegratorSERK{DT,TT}, sol::AtomicSolutionSDE{DT,TT},
                                     cache::IntegratorCacheSERK{DT}=int.caches[DT]) where {DT,TT}
    local tᵢ::TT
    local ydrift::DT

    # reset cache
    reset!(sol, timestep(int))

    for i in eachstage(int)
        for k in eachindex(cache.Q[i])
            # contribution from the drift part
            ydrift = 0
            for j = 1:i-1
                ydrift += tableau(int).qdrift.a[i,j] * cache.V[j][k]
            end

            # ΔW contribution from the diffusion part
            cache.Δy .= 0
            for j = 1:i-1
                for l in eachnoise(int)
                    cache.Δy[l] += tableau(int).qdiff.a[i,j] * cache.B[j][k,l]
                end
            end

            cache.Q[i][k] = sol.q̅[k] + timestep(int) * ydrift + dot(cache.Δy, sol.ΔW)

            # ΔZ contribution from the diffusion part
            if tableau(int).qdiff2.name ≠ :NULL
                cache.Δy .= 0
                for j = 1:i-1
                    for l in eachnoise(int)
                        cache.Δy[l] += tableau(int).qdiff2.a[i,j] * cache.B[j][k,l]
                    end
                end

                cache.Q[i][k] += dot(cache.Δy, sol.ΔZ)/timestep(int)
            end

        end

        tᵢ = sol.t̅ + timestep(int) * tableau(int).qdrift.c[i]
        equation(int, :v)(tᵢ, cache.Q[i], cache.V[i])
        equation(int, :B)(tᵢ, cache.Q[i], cache.B[i])
    end

    # compute final update
    if tableau(int).qdiff2.name == :NULL
        update_solution!(sol, cache.V, cache.B, tableau(int).qdrift.b, tableau(int).qdiff.b, timestep(int), sol.ΔW, cache.Δy)
        # update_solution!(sol, cache.V, cache.B, tableau(int).qdrift.b̂, tableau(int).qdiff.b̂, timestep(int), int.ΔW, cache.Δy)
    else
        update_solution!(sol, cache.V, cache.B, tableau(int).qdrift.b, tableau(int).qdiff.b, tableau(int).qdiff2.b, timestep(int), sol.ΔW, sol.ΔZ, cache.Δy)
        # update_solution!(sol, cache.V, cache.B, tableau(int).qdrift.b̂, tableau(int).qdiff.b̂, tableau(int).qdiff2.b̂, timestep(int), sol.ΔW, sol.ΔZ, cache.Δy)
    end
end


text:
function GPSBlockIIF(arg0::jdouble, arg1::jdouble, arg2::AbsoluteDate, arg3::AbsoluteDate, arg4::ExtendedPVCoordinatesProvider, arg5::Frame)
    return GPSBlockIIF((jdouble, jdouble, AbsoluteDate, AbsoluteDate, ExtendedPVCoordinatesProvider, Frame), arg0, arg1, arg2, arg3, arg4, arg5)
end


text:
using Colors

@testset "din99" begin
    # Test data from the DIN 6176 specification
    testdata = [
        ((50,  10,  10), (61.43,   9.70,   3.76)),
        ((50,  50,  50), (61.43,  28.64,  11.11)),
        ((50, -10,  10), (61.43,  -5.57,   7.03)),
        ((50, -50,  50), (61.43, -17.22,  21.75)),
        ((50, -10, -10), (61.43,  -9.70,  -3.76)),
        ((50, -50, -50), (61.43, -28.64, -11.11)),
        ((50,  10, -10), (61.43,   5.57,  -7.03)),
        ((50,  50, -50), (61.43,  17.22, -21.75)),
        (( 0,   0,   0), ( 0,      0,      0)),
        ((100,  0,   0), (100,     0,      0))]

    # A high error threshold has been chosen because converting from DIN99
    # to CIELAB with only two decimal places of accuracy yields fairly inaccurate
    # results due to the exponentiation.
    conveps = 0.05
    diffeps = 0.01
    let a, b, metric
        for (i, (a, b)) in enumerate(testdata)
            converted = convert(DIN99, Lab(a...))
            test = DIN99(b...)

            @test (abs(converted.l - test.l) < conveps)
            @test (abs(converted.a - test.a) < conveps)
            @test (abs(converted.b - test.b) < conveps)

            converted = convert(Lab, DIN99(b...))
            test = Lab(a...)

            @test (abs(converted.l - test.l) < conveps)
            @test (abs(converted.a - test.a) < conveps)
            @test (abs(converted.b - test.b) < conveps)

            # This is not a real test of the color difference metric, but at least
            # makes sure it isn't doing anything really crazy.
            metric = DE_DIN99()
            @test (abs(colordiff(convert(DIN99, Lab(a...)), DIN99(b...), metric)) < diffeps)
        end
    end
end  # @testset


text:
using ITensors

include("runtests.jl")
include("itensor.jl")
include("dmrg.jl")


text:
#=
    Vector valued constraints
=#
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Nonnegatives}) = cmap(m).nonnegatives
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Nonpositives}) = cmap(m).nonpositives
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Zeros})        = cmap(m).zeros

function MOI.addconstraint!(m::LinQuadOptimizer, func::VecLin, set::S) where S <: VecLinSets
    @assert MOI.dimension(set) == length(func.constant)

    nrows = get_number_linear_constraints(m)
    addlinearconstraint!(m, func, backend_type(m,set))
    nrows2 = get_number_linear_constraints(m)

    m.last_constraint_reference += 1
    ref = VLCI{S}(m.last_constraint_reference)

    dict = constrdict(m, ref)
    dict[ref] = collect(nrows+1:nrows2)
    for i in 1:MOI.dimension(set)
        push!(m.constraint_primal_solution, NaN)
        push!(m.constraint_dual_solution, NaN)
        push!(m.constraint_constant, func.constant[i])
    end
    ref
end

function addlinearconstraint!(m::LinQuadOptimizer, func::VecLin, sense::Cchar)
    @assert length(func.outputindex) == length(func.variables) == length(func.coefficients)
    # get list of unique rows
    rows = unique(func.outputindex)
    @assert length(rows) == length(func.constant)
    # sort into row order
    pidx = sortperm(func.outputindex)
    cols = getcol.(m, func.variables)[pidx]
    vals = func.coefficients[pidx]
    # loop through to gte starting position of each row
    rowbegins = Vector{Int}(length(rows))
    rowbegins[1] = 1
    cnt = 1
    for i in 2:length(pidx)
        if func.outputindex[pidx[i]] != func.outputindex[pidx[i-1]]
            cnt += 1
            rowbegins[cnt] = i
        end
    end
    add_linear_constraints!(m, rowbegins, cols, vals, fill(sense, length(rows)), -func.constant)
end

MOI.canmodifyconstraint(m::LinQuadOptimizer, ::VLCI{<: VecLinSets}, ::Type{MOI.VectorConstantChange{Float64}}) = true
function MOI.modifyconstraint!(m::LinQuadOptimizer, ref::VLCI{<: VecLinSets}, chg::MOI.VectorConstantChange{Float64})
    @assert length(chg.new_constant) == length(m[ref])
    for (r, v) in zip(m[ref], chg.new_constant)
        change_coefficient!(m, r, 0, -v)
        m.constraint_constant[r] = v
    end
end

MOI.candelete(m::LinQuadOptimizer, c::VLCI{<:VecLinSets}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::VLCI{<:VecLinSets})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    # we delete rows from largest to smallest here so that we don't have
    # to worry about updating references in a greater numbered row, only to
    # modify it later.
    for row in sort(dict[c], rev=true)
        delete_linear_constraints!(m, row, row)
        deleteat!(m.constraint_primal_solution, row)
        deleteat!(m.constraint_dual_solution, row)
        deleteat!(m.constraint_constant, row)
        # shift all the other references
        shift_references_after_delete_affine!(m, row)
    end
    delete!(dict, c)
end

#=
    Constraint set of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{VLCI{S}}) where S <: VecLinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::VLCI{S}) where S <: VecLinSets
    constraint_indices = m[c]
    S(length(constraint_indices))
end

#=
    Constraint function of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{VLCI{S}}) where S <: VecLinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::VLCI{<: VecLinSets})
    ctrs = m[c]
    n = length(ctrs)
    out = MOI.VectorAffineFunction(Int[],VarInd[],Float64[],Float64[])
    for i in 1:n
        rhs = get_rhs(m, ctrs[i])
        push!(out.constant, -rhs)

        # TODO more efficiently
        colidx, coefs = get_linear_constraint(m, ctrs[i])
        append!(out.variables, m.variable_references[colidx+1])
        append!(out.coefficients, coefs)
        append!(out.outputindex, i*ones(Int,length(coefs)))
    end
    return out
end


text:
## abstract types ##

abstract type Pi0Estimator end

abstract type Pi0Fit end

abstract type PValueAdjustment end

abstract type PValueCombination end

## concrete types ##

## PValues

struct PValues{T <: AbstractFloat} <: AbstractVector{T}
    values::Vector{T}
    min::T
    max::T

    function (::Type{PValues})(values::AbstractVector{T}) where T
        min, max = extrema(values)
        if min < zero(T) || max > one(T)
            throw(DomainError("p-values must all be in [0, 1]"))
        end
        return new{T}(copy(values), min, max)
    end
end

Base.convert(::Type{PValues}, x::AbstractVector{T}) where T <: AbstractFloat = PValues(x)

Base.size(pv::PValues) = size(pv.values)
Base.getindex(pv::PValues, i::Integer) = pv.values[i]
Base.setindex!(pv::PValues, x::AbstractFloat, i::Integer) =
    throw(ErrorException("Modification of values is not permitted"))

Base.minimum(pv::PValues) = pv.min
Base.maximum(pv::PValues) = pv.max
Base.extrema(pv::PValues) = (minimum(pv), maximum(pv))
Base.length(pv::PValues) = length(pv.values)


text:
using HetaSimulator
using Plots

p = load_platform("./src")

m = p.models[:nameless]

Scenario(m, tspan = (0, 50)) |> sim |> plot

### Add scenarios

scenarios_df = read_scenarios("./data/scenarios.csv")
add_scenarios!(p, scenarios_df)

### Add measurements
data1 = read_measurements("./data/NoNP_67.csv")
data2 = read_measurements("./data/U16c_61.csv")
data3 = read_measurements("./data/U16p_61.csv")
data4 = read_measurements("./data/T_67.csv")
add_measurements!(p, data1)
add_measurements!(p, data2)
add_measurements!(p, data3)
add_measurements!(p, data4)

# measure = get_measuments(p, scenarios=[:T_67_low])
# measure1 = get_measuments(scenarios(p)[:T_67_low])

# measure2 = get_measurements(res)
# simulate2 = DataFrame(res)

res = sim(p)
plot(res, yscale=:log10, ylims=(1e-1,1e3))
#plotd = plot(res, yscale=:log10, ylims=(1e-1,1e3))
#savefig(plotd, "sim1.png")

### fitting

# to_fit = read_parameters("./julia/parameters.csv")

to_fit = [
    :sigma_K => 0.1,
    :sigma_L => 0.1,
    :sigma_S => 0.1,
    :sigma_T => 0.1,
    :sigma_U => 0.1,
    :sigma_P => 0.1,

    :Kp_K_D => 5.562383e+01,
    :Kp_R_D => 5.562383e+01,
    :Kp_L_D => 4.196478e+01,
    :Kp_S_D => 6.475789e+01,
    :Kp_T_D => 6.322612e+01,
    :Kp_U_D => 7.643254e+01,
    :PS_K_061 => 9.271306e-04,
    :PS_L_061 => 2.142868e-02,
    :PS_R_061 => 5.792745e-02,
    :PS_S_061 => 2.022008e-03,
    :PS_T_061 => 1.305856e-03,
    :PS_U_061 => 1.217994e-04,
    :Q_T_g => 2.062306e+00,
    :bto_K_m => 8.579112e-02,
    :bto_L_m => 1.408861e-01,
    :bto_S_m => 1.529400e-01,
    :bto_T_m => 9.959412e-03,
    :bto_U_m => 2.230023e-01,
    :deg_resp => 2.455650e+00,
    :kel_L_D_m => 4.749195e+00,
    :koef_resp_K => 7.934936e+01,
    :koef_resp_L => 2.091350e+02,
    :koef_resp_R => 2.456379e+02,
    :tar_S => 9.044673e-01,
    :w_B_m => 8.138074e-02
]
# OFV: 2584
fit_res = fit(p, to_fit)

optim(fit_res)
# save_as_heta("julia/fitted.heta", fit_res)

res_optim = sim(p, parameters_upd = optim(fit_res))
plot(res_optim, yscale=:log10, ylims=(1e-4,1e3))

## plot fitted results in files
plot1 = plot(res_optim[1:3], yscale=:log10, ylims=(1e-1,1e3))
plot2 = plot(res_optim[4:6], yscale=:log10, ylims=(1e-1,1e3))
plot3 = plot(res_optim[7:9], yscale=:log10, ylims=(1e-1,1e3))
plot4 = plot(res_optim[10:11], yscale=:log10, ylims=(1e-4,1e3))

savefig(plot1, "diagnostics/heta-simulator/U16c_61.png")
savefig(plot2, "diagnostics/heta-simulator/U16p_61.png")
savefig(plot3, "diagnostics/heta-simulator/T_67.png")
savefig(plot4, "diagnostics/heta-simulator/NoNP_67.png")

max_stars_repo_name: insysbio/drug-loaded-nanoparticles

max_issues_repo_name: insysbio/drug-loaded-nanoparticles

max_forks_repo_name: insysbio/drug-loaded-nanoparticles


text:
"""
This module provides a framework for defining systems of units.
"""
module units
import Base: +, -, *, ^, /
import Base: convert, show, string, promote_rule

export BaseUnit, Unit, Unitful, SI, reciprocal

"""
Represents quantities that have units.
"""
abstract Unitful
*(x::Unitful, y::Number) = *(promote(x, y)...)
*(x::Number, y::Unitful) = *(promote(x, y)...)
*(x::Unitful, y::Unitful) = *(promote(x, y)...)
/(x::Unitful, y::Number) = /(promote(x, y)...)
/(x::Number, y::Unitful) = /(promote(x, y)...)
/(x::Unitful, y::Unitful) = /(promote(x, y)...)

"""
Represents the unit (e.g. seconds, meters) of a quantity.
"""
abstract BaseUnit{n} <: Unitful

"""
Represents a whole system of units.
"""
type UnitSystem
    baseunits::Array{AbstractString}
    abbreviations::Array{AbstractString}
    function UnitSystem()
        s = new([], [])
        push!(s.baseunits, "Dimensionless")
        push!(s.abbreviations, "#")
        s
    end
end

# The SI system of units
SI = UnitSystem()
push!(SI.baseunits, "Meter")
push!(SI.baseunits, "Second")
push!(SI.baseunits, "Kilogram")
push!(SI.baseunits, "Radian")
push!(SI.abbreviations, "m")
push!(SI.abbreviations, "s")
push!(SI.abbreviations, "kg")
push!(SI.abbreviations, "rad")

# This loop defines all the standard behaviors for base units
for unitabbr in zip(SI.baseunits, SI.abbreviations)
    unit = symbol(unitabbr[1])
    abbrtext = unitabbr[2]
    abbr = symbol(abbrtext)

    @eval begin
        export $unit, $abbr
        """
        The unit.
        """
        type $unit{n} <: BaseUnit{n}
        end
        """
        Constructor for parameter-less type defaults to parameter n=1.
        """
        function $unit()
            $unit{1}()
        end
        """
        Extract the base parametric type (i.e. Meter for Meter{2}).
        """
        function unit(x::$unit)
            $unit
        end
        function unit{T<:$unit}(::Type{T})
            $unit
        end
        """
        Get the power that the unit is raised to (e.g 2 for m^2).
        """
        pow{n}(quantity::$unit{n}) = n
        """
        Get the power that the unit is raised to (e.g 2 for m^2).
        """
        pow{n}(::Type{$unit{n}}) = n
        """
        Get a power of 1 for parameter-less type.
        """
        pow(::Type{$unit}) = 1
        $abbr = $unit()
        function string(x::$unit{1})
            string($abbrtext)
        end
        function string{n}(x::$unit{n})
            string($abbrtext, "^$n")
        end
        function show(io::IO, quantity::$unit)
            print(io, string(quantity))
        end
        function *{n, m}(x::$unit{n}, y::$unit{m})
            $unit{n+m}()
        end
        function /{n, m}(x::$unit{n}, y::$unit{m})
            $unit{n-m}()
        end
        function ^{n}(x::$unit{n}, y::Integer)
            $unit{n*y}()
        end
        function ^{n}(x::$unit{n}, y::Number)
            $unit{n*y}()
        end
        reciprocal{n}(x::$unit{n}) = $unit{-n}()
        reciprocal{n}(::Type{$unit{n}}) = $unit{-n}
    end
end

"""
Represents a general combination of different base units.

Unit objects have a type invariant which is enforced through all of
the operations defined below. The invariant is that each BaseUnit
type appears at most once in the type parameter. This is what you'd
expect—although technicaly units like seconds per second exist, I don't
want them to show up here.
"""
typealias BaseUnitTuple Tuple{Vararg{BaseUnit}}
type Unit{T<:BaseUnitTuple} <: Unitful
end

# Constructor out of BaseUnit type
Unit(x::BaseUnit) = convert(Unit, x)

function types{T<:BaseUnitTuple}(x::Unit{T})
    T.types
end

function show{T<:BaseUnitTuple}(io::IO, x::Unit{T})
    unitsstring = mapreduce(
        S->string(S()),
        (a, b)->string(a, " ", b),
        "", T.types)[2:end]
    print(io, unitsstring)
end

function promote_rule{T<:Unitful, S<:Unitful}(::Type{T}, ::Type{S})
    Unit
end
function convert{T<:BaseUnit}(::Type{Unit}, x::T)
    Unit{Tuple{T}}()
end

function *{T<:BaseUnit, S<:BaseUnit}(x::T, y::S)
    tupletypes = (string(x) < string(y) ? (T, S)
        : (S, T))
    Unit{Tuple{tupletypes...}}()
end
function *{T<:BaseUnitTuple}(x::BaseUnit, y::Unit{T})
    xsearch = unit(x)
    ylist = types(y)
    ylist = map(x->unit(x), ylist)
    xiny = findin(ylist, (xsearch,))
    if length(xiny) == 0
        # Find the last index for which typeof(x) comes later in the
        # alphabet than typeof(T.types[index])
        index = indexin([true], map(t->string(t) < string(typeof(x)),
            T.types))[1] + 1
        tupletype = Tuple{T.types[1:index-1]...,typeof(x),
            T.types[index:end]...}
    elseif length(xiny) == 1
        index = xiny[1]
        newtype = unit(x){pow(T.types[index]) + 1}
        tupletype = Tuple{T.types[1:index-1]..., newtype,
        T.types[index+1:end]...}
    else
        error("Corrupted Unit object: $y")
    end
    result = Unit{tupletype}()
end

function *{T<:BaseUnitTuple, S<:BaseUnitTuple}(x::Unit{T}, y::Unit{S})
    result = Unit{T}()
    yfactors = map(x->x(), types(y))
    for factor in yfactors
        result = factor * result
    end
    result
end

/(x::Unit, y::Unit) = x * reciprocal(y)

function reciprocal{T<:BaseUnitTuple}(x::Unit{T})
    newtypes = map(reciprocal, T.types)
    Unit{Tuple{newtypes...}}()
end

end


text:
using AstroBase
using LinearAlgebra: norm
using OrdinaryDiffEq: ODEProblem, solve, Vern9
using Test

function rhs(du, u, p, t)
    du[1:3] = u[4:6]
    r3 = norm(u[1:3])^3
    du[4:6] = -p .* u[1:3] ./ r3
end

r = [
        6068279.27,
        -1692843.94,
        -2516619.18,
    ]/1000

v = [
        -660.415582,
        5495.938726,
        -5303.093233,
    ]/1000

rv = [r; v]

ep = UTCEpoch("2016-05-30T12:00:00.000")

ele = keplerian(r, v, grav_param(earth))
t1 = period(ele[1], grav_param(earth))
s0_exp = State(ep, r, v)
k0_exp = KeplerianState(s0_exp)

@testset "Coords" begin
    @testset "States" begin
        eph = VSOP87()

        tdb = TDBEpoch(ep)

        iss = State(ep, r, v)
        rot = Rotation(icrf, iau_earth, ep)
        pos_mars, vel_mars = state(eph, ep, earth, mars)
        r_mars = r .- pos_mars
        v_mars = v .- vel_mars
        mars_rot = State(tdb, r, v, frame=iau_mars, body=mars)

        @testset "Operators" begin
            @test iss == State(ep, r, v, frame=icrf, body=earth)
            @test iss ≈ State(ep, r .+ eps(), v .+ eps(), frame=icrf, body=earth)
        end
        @testset "Change frame" begin
            @test State(iss, eph, frame=iau_earth) ≈ State(ep, rot(r, v)..., frame=iau_earth)
            @test iss ≈ State(State(iss, eph, frame=iau_earth), eph, frame=icrf)
        end
        @testset "Change time scale" begin
            @test State(iss, eph, scale=TDB) == State(tdb, r, v)
        end
        @testset "Change central body" begin
            @test State(iss, eph, body=mars) ==
                State(ep, r_mars, v_mars, frame=icrf, body=mars)
            @test State(State(mars_rot, eph, body=earth), eph, body=mars) ≈ mars_rot
        end
        @testset "Change scale and body" begin
            @test State(iss, eph, scale=TDB, body=mars) ==
                State(tdb, r_mars, v_mars, frame=icrf, body=mars)
        end
        @testset "Change frame and body" begin
            @test State(iss, eph, frame=iau_earth, body=mars) ≈
                State(ep, rot(r_mars, v_mars)..., frame=iau_earth, body=mars)
            @test State(State(mars_rot, eph, frame=itrf, body=earth), eph, frame=iau_mars,
                body=mars) ≈ mars_rot
        end
        @testset "Change frame and scale" begin
            @test State(iss, eph, frame=iau_earth, scale=TDB) ≈
                State(tdb, rot(r, v)..., frame=iau_earth)
        end
        @testset "Change scale, frame, and body" begin
            @test State(State(mars_rot, eph, frame=itrf, scale=UT1, body=earth), eph,
                frame=iau_mars, scale=TDB, body=mars) ≈ mars_rot
        end
    end
    @testset "Time Series" begin
        t = 1seconds:4seconds
        u = collect(1.0:4.0)
        ep = UTCEpoch(2000, 1, 1)
        ts = TimeSeries(ep, t, u)
        @test collect(ts) == [1.0, 2.0, 3.0, 4.0]
        @test ts[1] == 1.0
        @test ts(1.5seconds) ≈ 1.5
        @test ts(ep + 1.5seconds) ≈ 1.5
    end
    @testset "Trajectories" begin
        @testset "Transpose - Array of Arrays" begin
            a = [[1.0, 2.0, 3.0], [1.0, 2.0, 3.0]]
            b = [[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]]
            @test AstroBase.Coords.transpose_aoa(a) == b
        end
        @testset "Propagation" begin
            prob = ODEProblem(rhs, rv, (0.0, value(t1)), p=grav_param(earth))
            sol = solve(prob, Vern9())
            tra = Trajectory(ep, sol.t * seconds, sol.u)
            rv0 = tra(0.0seconds)
            ele0 = keplerian(tra, 0.0seconds)
            s0 = State(tra, 0.0seconds)
            k0 = KeplerianState(tra, 0.0seconds)
            @test rv0 ≈ rv
            @test collect(ele0) ≈ collect(ele)
            @test s0_exp ≈ s0
            @test k0_exp ≈ k0
            rv1 = tra(t1)
            ele1 = keplerian(tra, t1)
            s1 = State(tra, t1)
            k1 = KeplerianState(tra, t1)
            @test rv1 ≈ rv atol=0.1
            @test collect(ele1) ≈ collect(ele) atol=0.1
            @test s0_exp ≈ s1 atol=0.1
            @test k0_exp ≈ k1 atol=0.1
        end
    end
end


text:
using Roots
include("TypeDefine.jl")
include("PressureGuess.jl")

#function (4.6) and (4.7) from Toro
f(p, status::HydroStatus, ::Shock) = (p - status.p) * sqrt(status.A / (p + status.B))
f(p, status::HydroStatus, ::Rarefaction) = 2. * status.c / (status.gamma - 1.) * ((p / status.p)^((status.gamma - 1.) / 2. /status.gamma) - 1.)
function f(p, status::HydroStatus)
	wave_type = p > status.p ? Shock() : Rarefaction()
	return f(p, status, wave_type)
end

#solver for p_star based on the first guess of the pressure(see PressureGuess.jl)
#uses secant method to find solve the equation
function p_star_calc(left::HydroStatus, right::HydroStatus, guess_scheme::T, TOL::Float64) where {T<:FirstGuessScheme}
	p0 = max(TOL, p_guess(left, right, guess_scheme))

	#we define p = exp(u) so inside the interation it won't get negative value
	u0 = log(p0) #p = exp(u) => u = log(p)
	f_star(u) = f(exp(u), left) + f(exp(u), right) + right.u - left.u #function (4.5) from Toro
	u = find_zero(f_star, u0, Order1(), rtol = TOL) #Secant method
	return exp(u) #p = exp(u)
end

#equation (4.9) from Toro
u_star_calc(p_star, left::HydroStatus, right::HydroStatus) = 0.5 * (left.u + right.u) +
															 0.5 * (f(p_star, right) - f(p_star, left))

#equations (4.23) and (4.32) from Toro
rho_star_calc(p_star, status::HydroStatus, ::Rarefaction) = status.rho * (p_star / status.p) ^ (1. / status.gamma)

#equation (4.25) and (4.34) from Toro
c_star_calc(p_star, status::HydroStatus, ::Rarefaction) = status.c * (p_star / status.p) ^ ((status.gamma - 1.) / 2. / status.gamma)

#equations (4.50) and (4.57) from Toro
function rho_star_calc(p_star, status::HydroStatus, ::Shock)
	gamma_ratio = (status.gamma - 1.) / (status.gamma + 1.)
	pressure_ratio = p_star / status.p
	return status.rho * (pressure_ratio + gamma_ratio) / (pressure_ratio * gamma_ratio + 1.)
end

#equation (4.55) and (4.62) from Toro
head_speed_calc(p_star, u_star, status::HydroStatus,
				side::T, ::Rarefaction) where {T<:Side} = status.u + side_factor(side) * status.c

#equation (4.55) and (4.62) from Toro
function tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Rarefaction) where {T<:Side}
	c_star = c_star_calc(p_star, status, Rarefaction())
	return u_star + side_factor(side) * c_star
end

#equations (4.52) and (4.59) from Toro
function head_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Shock) where {T<:Side}
	gamma = status.gamma
	gamma_ratio1 = (gamma + 1.) / 2. / gamma
	gamma_ratio2 = (gamma - 1.) / 2. / gamma
	return status.u + side_factor(side) * status.c * sqrt(gamma_ratio1 * p_star / status.p + gamma_ratio2)
end

#for shock the head and the tail are the same
tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Shock) where {T<:Side} = head_speed_calc(p_star, u_star, status, side, Shock())

#for vaccum generation the density in the middle is zero (= vaccum)
rho_star_calc(p_star, status::HydroStatus, ::Vaccum) = 0.

#those are the velocities of the head and the tail of the rarefaciton wave when a vaccum is generated
#see (4.76), (4.77), (4.79) and (4.80) from Toro
head_speed_calc(p_star, u_star, status::HydroStatus,
				side::T, ::Vaccum) where {T<:Side} = status.u + side_factor(side) * status.c
tail_speed_calc(p_star, u_star, status::HydroStatus, side::T, ::Vaccum) where {T<:Side} = status.u - side_factor(side) * 2 * status.c / (status.gamma - 1.)

#this is the density, velocity and pressure profile of a rarefaction wave in some specific x/t
#see (4.56) and (4.63) from Toro
function rarefaction_profile(x_over_t, status::HydroStatus, side::T) where {T<:Side}
	gamma = status.gamma
	gamma_ratio = (gamma - 1.) / (gamma + 1.)
	gamma_plus = 2. / (gamma + 1.)
	gamma_minus = 2. / (gamma - 1.)

	rarefaction_factor = (gamma_plus - side_factor(side) * gamma_ratio / status.c * (status.u - x_over_t)) ^ gamma_minus

	rho = status.rho * rarefaction_factor
	u = gamma_plus * (-side_factor(side)*status.c + (1. /gamma_minus) * status.u + x_over_t)
	p = status.p * rarefaction_factor ^ gamma

	return HydroStatusSimple(rho, u, p)
end

#sample the riemann problem in the situation where both sides are not vaccum (but vaccum can be generated through the dynamics)
function sample_riemann_regular(x, t::Float64, left::HydroStatus, right::HydroStatus,
								  guess_scheme::T, TOL::Float64) where {T<:FirstGuessScheme}
	@assert left.rho != 0.
	@assert right.rho != 0.

	profile = similar(x, HydroStatusSimple) #we return the values of rho, u and p in each point x in space

	#the status "far" from the center = like the starting condition
	status_left = HydroStatusSimple(left)
	status_right = HydroStatusSimple(right)

	x_over_t = x ./ t

	if 2. / (left.gamma - 1.) * (left.c + right.c) <= right.u - left.u #vaccum generation
		#if vaccum is generated - the density in the middle will be zero
		p_star = 0.
		u_star = 0.
		wave_type_left = Vaccum()
		wave_type_right = Vaccum()
	else
		#if vaccum is not generated we can calculate the star profile in the middle
		p_star = p_star_calc(left, right, guess_scheme, TOL)
		u_star = u_star_calc(p_star, left, right)
		#check what kind of waves are in every direction (shock or rarefaction)
		wave_type_left = p_star > left.p ? Shock() : Rarefaction()
		wave_type_right = p_star > right.p ? Shock() : Rarefaction()
	end

	rho_star_left = rho_star_calc(p_star, left, wave_type_left)
	status_left_star = HydroStatusSimple(rho_star_left, u_star, p_star) #the profile in the left near the contact discontinuity
	head_speed_left = head_speed_calc(p_star, u_star, left, LeftSide(), wave_type_left)
	tail_speed_left = tail_speed_calc(p_star, u_star, left, LeftSide(), wave_type_left)

	rho_star_right = rho_star_calc(p_star, right, wave_type_right)
	status_right_star = HydroStatusSimple(rho_star_right, u_star, p_star) #the profile in the right near the contact discontinuity
	head_speed_right = head_speed_calc(p_star, u_star, right, RightSide(), wave_type_right)
	tail_speed_right = tail_speed_calc(p_star, u_star, right, RightSide(), wave_type_right)

	for i = 1:length(x)
		S = x_over_t[i] #this is like the S which Toro use in Section 4.5

		#see Figure 4.14 in Toro for the flow of the following lines
		side = S < u_star ? LeftSide() : RightSide()
		status = choose_side(left, right, side)
		head_speed = choose_side(head_speed_left, head_speed_right, side)
		tail_speed = choose_side(tail_speed_left, tail_speed_right, side)

		#this is used to flip the direction of the inequallity in the branching between the left and the right side
		#the xor which will be in the following lines uses this boolean like some "controlled not"
		#when right_condition is "true" the inequallity will be flipped
		#when right_condition is "false" the inequallity will be stay the same
		right_condition = isa(side, RightSide)

		if xor(S < head_speed, right_condition)
			profile[i] = choose_side(status_left, status_right, side)
		elseif xor(S < tail_speed, right_condition) #can only happen in Rarefaction, because in Shock  head_speed == tail_speed
			profile[i] = rarefaction_profile(S, status, side)
		else
			profile[i] = choose_side(status_left_star, status_right_star, side)
		end
	end

	return profile
end

#sample the riemann problem in the situation where the gas in "side" is actually a vaccum
function sample_riemann_side_vaccum(x, t::Float64, status::HydroStatus, side::T) where {T<:Side}

	profile = similar(x, HydroStatusSimple) #we return the values of rho, u and p in each point x in space
	vaccum = HydroStatusSimple(0., 0., 0.) #vaccum
	simple = HydroStatusSimple(status)

	x_over_t = x ./ t

	tail_speed = tail_speed_calc(0., 0., status, side, Vaccum()) #the head of the rarefaciton wave
	head_speed = head_speed_calc(0., 0., status, side, Vaccum()) #the tail of the rarefaciton wave (=the boundary with the vaccum)

	#see "sample_riemann_regular" for explnation about this boolean
	right_condition = isa(side, RightSide)

	for i = 1:length(x)
		S = x_over_t[i]

		if xor(S < head_speed, right_condition)
			profile[i] = simple
		elseif xor(S < tail_speed, right_condition)
			profile[i] = rarefaction_profile(S, status, side)
		else
			profile[i] = vaccum
		end
	end

	return profile
end

#this is the main function that samples the Riemann problem and returns the density, pressure and velocity in the space points
#specified by the array x. It gets the following parameters:
# x - an array type which saves all the points in space in which we want to sample the problem (pay attention - the contact discontinuity begins at x=0.)
# t - a floating number which tells us in what time we sample to problem
# left + right - struct from type HydroStatus (see "TypeDefine.jl") to get the initial state of the system
# guess_scheme - the type of guess which will be used to solve the p_star equation (see "PressureGuess.jl")
# TOL - the tolerance for the convergance of the p_star finding algorithm
function sample_riemann(x, t::Float64, left::HydroStatus =  HydroStatus(1.,    0., 1.,  7. / 5.),
									   right::HydroStatus = HydroStatus(0.125, 0., 0.1, 7. / 5.),
									   guess_scheme::T = PrimitiveValue(), TOL::Float64 = 1.e-6) where {T<:FirstGuessScheme}
	@assert left.gamma == right.gamma

	if left.rho == 0 #material is on the right, vaccum on the left
		return sample_riemann_side_vaccum(x, t, right, RightSide())
	elseif right.rho == 0 #material is on the left, vaccum on the right
		return sample_riemann_side_vaccum(x, t, left, LeftSide())
	else
		return sample_riemann_regular(x, t, left, right, guess_scheme, TOL)
	end

end

max_stars_repo_name: ryarazi/ExactRiemannProblemSolver

max_issues_repo_name: ryarazi/ExactRiemannProblemSolver

max_forks_repo_name: ryarazi/ExactRiemannProblemSolver


text:
using Documenter, GRASP

makedocs(
    format=:html,
    sitename = "GRASP.jl",
    modules = [GRASP]
)


text:
using Documenter, GeoData, Plots, Logging, Statistics, Dates

using GeoData.LookupArrays, GeoData.Dimensions

ENV["GKSwstype"] = "100"

# Plots warnings are brWarn doctests. They dont warn the second time.
# Downloads also show op in doctests. So download everything first.
function flush_info_and_warnings()
    # GeoStack(AWAP, (:tmin, :tmax); date=DateTime(2001, 1, 1))
    # GeoStack(EarthEnv{HabitatHeterogeneity}, (:evenness, :range, :contrast, :correlation))
    # GeoStack(WorldClim{BioClim}, (1, 3, 5, 7, 12))
    st = GeoStack(WorldClim{Climate}; month=1);
    plot(st)
end
flush_info_and_warnings()

Logging.disable_logging(Logging.Warn)

# Make the docs, without running the tests again
makedocs(
    modules = [GeoData],
    sitename = "GeoData.jl",
    strict = true,
    clean = false,
)

# Enable logging to console again
Logging.disable_logging(Logging.BelowMinLevel)

deploydocs(
    repo = "github.com/rafaqz/GeoData.jl.git",
)


text:
using Test
using CompScienceMeshes

p = simplex(
    point(0,0,0),
    point(1,0,0),
    point(0,1,0))

q1 = simplex(
    point(0.6, 0.6, 0.0),
    point(1.6, 0.6, 0.0),
    point(0.6, 1.6, 0.0),
)

q2 = simplex(
    point(0.4, 0.4, 0.0),
    point(1.4, 0.4, 0.0),
    point(0.4, 1.4, 0.0),
)

@test overlap(p, q1) == false
@test overlap(p, q2) == true

q1 = simplex(
    point(0.949726718617726,-0.278864925637586,0.142314838273285),
    point(0.989821441880933,0.0,-0.142314838273285),
    point(0.989821441880933,0.0,0.142314838273285),
)

q2 = simplex(
    point(0.949726718617726,-0.278864925637586,-0.142314838273285),
    point(0.989821441880933,0.0,-0.142314838273285),
    point(0.949726718617726,-0.278864925637586,0.142314838273285),
)

@test overlap(q1, q2) == false

## make sure the submesh function work for 1D meshes

l1 = meshsegment(1.0,1/2)
vt = skeleton(l1,0)
bd = boundary(l1)

overlaps = overlap_gpredicate(bd)
pred1 = c -> overlaps(simplex(vertices(vt,c)))
@test pred1(SVector(1))
@test !pred1(SVector(2))
@test pred1(SVector(3))

# test a case where the segments are:
#   not of unit length
#   colinear and opposite
#   meet in a common point
ch1 = simplex(point(1/3,0,0), point(1/3,1/3,0))
ch2 = simplex(point(1/3,1/3,0), point(1/3,2/3,0))
@test !overlap(ch1, ch2)


text:
module TestsForBayesianNonparametricStatistics

if v"0.6" ≤ VERSION < v"0.7-"
	include("testsforversion0dot6.jl")
elseif VERSION ≥ v"0.7-"
	include("testsforversion0dot7and1dot0.jl")
else
	error("BayesianNonparametricStatistics is only available for Julia version 0.6 and higher, no test functions available.")
end

end # module.


text:
const AllowedXY = Union{Integer,Colon,AbstractRange}
const AllowedBand = Union{Integer,Colon,AbstractArray}

# AbstractRasterBand indexing

Base.size(band::AbstractRasterBand) = width(band), height(band)
Base.firstindex(band::AbstractRasterBand, d) = 1
Base.lastindex(band::AbstractRasterBand, d) = size(band)[d]

Base.getindex(band::AbstractRasterBand, x::AllowedXY, y::AllowedXY) = begin
    I = Base.to_indices(band, (x, y))
    buffer = Array{pixeltype(band),2}(undef, length.(I)...)
    # Calculate `read!` args and read
    xoffset, yoffset = first.(I) .- 1
    xsize, ysize = length.(I)
    read!(band, buffer, xoffset, yoffset, xsize, ysize)
    # Drop dimensions of integer indices
    dimsize = map(length, _dropint(I...))
    if dimsize == ()
        return buffer[1]
    else
        return reshape(buffer, dimsize...)
    end
end

Base.setindex!(band::AbstractRasterBand, value, x::AllowedXY, y::AllowedXY) = begin
    I = Base.to_indices(band, (x, y))
    if value isa AbstractArray
        Base.setindex_shape_check(value, length.(I)...)
        value = reshape(value, length.(I))
    else
        value = reshape([value], 1, 1)
    end
    # Calculate `write!` args and write
    xoffset, yoffset = first.(I) .- 1
    xsize, ysize = length.(I)
    write!(band, value, xoffset, yoffset, xsize, ysize)
end

# AbstractDataset indexing

Base.size(dataset::AbstractDataset) = (size(getband(dataset, 1))..., nraster(dataset))
Base.firstindex(dataset::AbstractDataset, d) = 1
Base.lastindex(dataset::AbstractDataset, d) = size(dataset)[d]

Base.getindex(dataset::AbstractDataset, x::AllowedXY, y::AllowedXY, bands::AllowedBand=1) = begin
    I = x, y, bands = Base.to_indices(dataset, (x, y, bands))
    buffer = Array{pixeltype(getband(dataset, 1)),3}(undef, length.(I)...)
    # Calculate `read!` args and read
    xoffset, yoffset = first.((x, y)) .- 1
    xsize, ysize = length.((x, y))
    indices = _bandindices(bands)
    read!(dataset, buffer, indices, xoffset, yoffset, xsize, ysize)
    # Drop dimensions of integer indices
    dimsize = map(length, _dropint(I...))
    if dimsize == ()
        return buffer[1]
    else
        return reshape(buffer, dimsize...)
    end
end

Base.setindex!(dataset::AbstractDataset, value, x::AllowedXY, y::AllowedXY, bands::AllowedBand=1) = begin
    I = x, y, bands = Base.to_indices(dataset, (x, y, bands))
    if value isa AbstractArray
        Base.setindex_shape_check(value, length.((x, y, bands))...)
        value = reshape(value, length.(I))
    else
        value = reshape([value], 1, 1, 1)
    end
    # Calculate `write!` args and write
    xoffset, yoffset = first.((x, y)) .- 1
    xsize, ysize = length.((x, y))
    indices = _bandindices(bands)
    write!(dataset, value, indices, xoffset, yoffset, xsize, ysize)
end

# Index conversion utilities

_dropint(x::Integer, args...) = _dropint(args...)
_dropint(x::AbstractRange, args...) = (x, _dropint(args...)...)
_dropint() = ()

_bandindices(band::Integer) = Cint[band]
_bandindices(bands::AbstractArray) = Cint[b for b in bands]

Array(dataset::AbstractDataset) = read(dataset)


text:
function getinfo(s::Solver, m::DRO)
    return "Benchmark of the $(typeof(s).name) solver for the $(typeof(m).name) model with $(name(m.d)) \n"
end

function getinfo(s::Solver, m::Simplex)
    return "Benchmark of the $(typeof(s).name) solver for the $(typeof(m).name) model \n"
end

function statstodataframe(stats::Stats)
    fields = fieldnames(Stats)
    vals   = [[getfield(stats, field)] for field in fields]
    return DataFrames.DataFrame(vals, [fields...])
end

function benchmark(s::Solver, m::Model; maxevals::Int = 10, verbose::Bool = true)
    verbose && printstyled(getinfo(s, m); bold = true, color = :green)
    solve(s, m)

    progress = ProgressMeter.Progress(maxevals, 0.1)

    table = map(1:maxevals) do i
        solve(s, m)
        row = statstodataframe(stats)
        verbose && ProgressMeter.next!(progress)
        return row
    end |> rows -> reduce(vcat, rows)

    return DataFrames.aggregate(table, [:model, :constraint, :solver, :optimizer], [Statistics.mean, Statistics.std])
end

function benchmark(s::Solver, getmodel::Function, N::AbstractVector{<:Int}; maxevals::Int = 10, verbose::Bool = true)
    m = getmodel(10)
    verbose && printstyled(getinfo(s, m); bold = true, color = :green)

    ## Compiling ...
    benchmark(s, m; maxevals = 10, verbose = false)

    progress = ProgressMeter.Progress(length(N), 0.1)

    table = map(N) do n
        m   = getmodel(n)
        row = benchmark(s, m; maxevals = maxevals, verbose = false)
        verbose && ProgressMeter.next!(progress)
        return row
    end |> rows -> reduce(vcat, rows)
    table.dimension = N

    return table
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
Base.@kwdef struct Cell
    char::Char
    style::Crayons.Crayon
end

Cell(char) = Cell(char, Crayons.Crayon())

convert(::Type{Cell}, char::Char) = Cell(char)

struct Buffer
    area::Rect
    content::Matrix{Cell}
end

Buffer(rect::Rect) = Buffer(rect, Cell(' '))

function Buffer(rect::Rect, cell::Cell)
    content = Matrix{Cell}(undef, rect.height, rect.width)
    for i in eachindex(content)
        content[i] = deepcopy(cell)
    end
    return Buffer(rect, content)
end

function Buffer(lines::Vector{String})
    height = length(lines)
    width = maximum(length(l) for l in lines)
    buffer = Buffer(
                Rect(
                    1,
                    1,
                    width,
                    height,
                ),
            )
    for (y, line) in enumerate(lines)
        set(buffer, 1, y, line)
    end
    return buffer
end

function set(buffer::Buffer, col::Integer, row::Integer, style::Crayons.Crayon)
    row = min(row, size(buffer.content, 1))
    col = min(col, size(buffer.content, 2))
    buffer.content[row, col] = Cell(buffer.content[row, col].char, style)
end

function set(buffer::Buffer, col::Integer, row::Integer, lines::Vector{T}) where T <: AbstractString
    for line in lines
        set(buffer, col, row, string(line))
        row += 1
    end
end

function set(buffer::Buffer, col::Integer, row::Integer, line::String)
    for (_col, c) in enumerate(line)
        set(buffer, col + _col - 1, row, c)
    end
end

function set(buffer::Buffer, col::Integer, row::Integer, line::String, style::Crayons.Crayon)
    for (_col, c) in enumerate(line)
        set(buffer, col + _col - 1, row, c, style)
    end
end

function set(buffer::Buffer, col::Integer, row::Integer, c::Char)
    set(buffer, col, row, c, buffer.content[row, col].style)
end

function set(buffer::Buffer, col::Integer, row::Integer, c::Char, style::Crayons.Crayon)
    set(buffer, col, row, Cell(c, style))
end

function set(buffer::Buffer, col::Integer, row::Integer, c::Cell)
    row = min(row, size(buffer.content, 1))
    col = min(col, size(buffer.content, 2))
    buffer.content[row, col] = c
end

reset(buffer::Buffer) = fill!(buffer.content, Cell(' '))

function diff(buffer1::Buffer, buffer2::Buffer)
    arr = Tuple{Tuple{Int, Int}, Cell}[]
    R, C = size(buffer1.content)
    for c in 1:C
        for r in 1:R
            if buffer1.content[r, c] != buffer2.content[r, c]
                push!(arr, ((r, c), buffer2.content[r, c]))
            end
        end
    end
    return arr
end

function background(buf::Buffer, rect::Rect, color::Crayons.Crayon)
    for y in top(rect):bottom(rect)
        for x in left(rect):right(rect)
            set(buf, x, y, color);
        end
    end
end

draw(t, widget, r::Rect) = draw(widget, r, current_buffer(t))


max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.7 P(x) (d+e x^2)^q (a+b x^2+c x^4)^p.jl


text:
using Test

@testset "MCMCChains" begin
    # run plotting tests
    include("plot_test.jl")

    # run function tests
    include("diagnostic_tests.jl")

    # run tests for missing values
    include("missing_tests.jl")

    # run tests for missing values
    include("sections_tests.jl")

    # run tests for missing values
    include("serialization_tests.jl")

    # run tests for sampling api
    include("sampling_tests.jl")

    # run tests for array constructor
    include("arrayconstructor_tests.jl")

    # run tests for sampoling api
    include("sampling_tests.jl")

    # run tests for array constructor
    include("arrayconstructor_tests.jl")

    # run tests for dataframe constructor
    include("dfconstructor_tests.jl")

    # run tests for dataframe constructor
    include("summarize_tests.jl")

    # run tests for posterior stats
    include("modelstats_test.jl")
end


text:
function addition(;x, y)
    result = x+y
    return (;result)  # implicit field name
end

function subtraction(;x, y)
    result = x-y
    return (;result)
end

function multiplication(;x, y)
    result = x*y
    return (;result)
end

function division(;x, y)
    result = x/y
    return (;result)
end

function conquer(;kwargs...) # the function MUST accept all keyword arguments passed to it, but is not obliged to process them
    result = +Inf
    return (;result)
end

fns = (;addition, subtraction, multiplication, division, conquer)


text:
module RandomBasedArrays

using Random

export RandomBasedArray

const _default_rng = if VERSION < v"1.3"
    () -> Random.GLOBAL_RNG
else
    Random.default_rng
end

struct RandomBasedArray{T,N,P<:AbstractArray,R<:AbstractRNG} <: AbstractArray{T,N}
    parent::P
    rng::R
end
RandomBasedArray(p::P, rng::R = _default_rng()) where {T,N,P<:AbstractArray{T,N},R<:AbstractRNG} =
    RandomBasedArray{T,N,P,R}(p, rng)

index(A::RandomBasedArray) = rand(A.rng, eachindex(parent(A)))

Base.parent(A::RandomBasedArray) = A.parent
Base.size(A::RandomBasedArray) = size(parent(A))
Base.getindex(A::RandomBasedArray, ::Int) =
    @inbounds getindex(parent(A), index(A))
Base.getindex(A::RandomBasedArray{T,N}, ::Vararg{Int,N}) where {T,N} =
    getindex(parent(A), Base._ind2sub(parent(A), index(A))...)
Base.setindex!(A::RandomBasedArray, v, ::Int) =
    @inbounds setindex!(parent(A), v, index(A))
Base.setindex!(A::RandomBasedArray{T,N}, v, ::Vararg{Int,N}) where {T,N} =
    setindex!(parent(A), v, Base._ind2sub(parent(A), index(A))...)
Base.show(io::IO, m::MIME"text/plain", A::RandomBasedArray) = show(io, m, copy(A))

end # module

max_stars_repo_name: giordano/RandomBasedArrays

max_issues_repo_name: giordano/RandomBasedArrays

max_forks_repo_name: giordano/RandomBasedArrays


text:
@testset "Kriging" begin
  # -----------------
  # 1D PROBLEM (ALL)
  # -----------------

  data1D = georef((z=[0.0,0.1,0.2,0.3,0.4,0.5,0.4,0.3,0.2,0.1,0.0],), collect(0.0:10.0:100.0)')
  grid1D = CartesianGrid(100)
  problem = EstimationProblem(data1D, grid1D, :z)

  global_kriging = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.), maxneighbors=nothing))
  nearest_kriging = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.), maxneighbors=3))
  local_kriging = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.),
                                 maxneighbors=3, neighborhood=NormBall(100.)))

  Random.seed!(2021)
  solvers   = [global_kriging, nearest_kriging, local_kriging]
  solnames  = ["global", "nearest", "local"]
  solutions = [solve(problem, solver) for solver in solvers]

  if visualtests
    for i in 1:3
      solution, sname = solutions[i], solnames[i]
      @test_reference "data/krig-1D-$(sname).png" plot(solution, size=(600,200))
    end
  end

  # --------------------
  # 2D PROBLEM (GLOBAL)
  # --------------------

  data2D = georef((z=[1.,0.,1.],), [(25.,25.), (50.,75.), (75.,50.)])
  grid2D = CartesianGrid((100,100), (0.5,0.5), (1.0,1.0))

  problem = EstimationProblem(data2D, grid2D, :z)

  solver = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.), maxneighbors=nothing))

  Random.seed!(2021)
  solution = solve(problem, solver)
  Z = asarray(solution, :z)
  @test isapprox(Z[25,25], 1., atol=1e-3)
  @test isapprox(Z[50,75], 0., atol=1e-3)
  @test isapprox(Z[75,50], 1., atol=1e-3)

  if visualtests
    @test_reference "data/krig-2D-global.png" contourf(solution, size=(800,400))
  end

  # ---------------------
  # 2D PROBLEM (NEAREST)
  # ---------------------

  problem = EstimationProblem(data2D, grid2D, :z)

  solver = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.),))

  Random.seed!(2021)
  solution = solve(problem, solver)
  Z = asarray(solution, :z)
  @test isapprox(Z[25,25], 1., atol=1e-3)
  @test isapprox(Z[50,75], 0., atol=1e-3)
  @test isapprox(Z[75,50], 1., atol=1e-3)

  if visualtests
    @test_reference "data/krig-2D-nearest.png" contourf(solution, size=(800,400))
  end

  # -------------------
  # 2D PROBLEM (LOCAL)
  # -------------------

  problem = EstimationProblem(data2D, grid2D, :z)

  solver = Kriging(:z => (variogram=GaussianVariogram(range=35.,nugget=0.),
                          maxneighbors=3, neighborhood=NormBall(100.)))

  Random.seed!(2021)
  solution = solve(problem, solver)

  # basic checks
  inds = LinearIndices(size(grid2D))
  S = solution[:z]
  @test isapprox(S[inds[25,25]], 1., atol=1e-3)
  @test isapprox(S[inds[50,75]], 0., atol=1e-3)
  @test isapprox(S[inds[75,50]], 1., atol=1e-3)

  if visualtests
    @test_reference "data/krig-2D-local.png" contourf(solution, size=(800,400))
  end
end


text:
t = parse(Int, readline())
h = div(t, 3600)
m = div(t, 60) % 60
s = t % 60
println("$h:$m:$s")


text:
@testset "call" begin
    @testset "basic" begin
        @system SCall(Controller) begin
            fa(; x) => x ~ call
            a(fa) => fa(1) ~ track
            fb(i; x) => i + x ~ call
            b(fb) => fb(1) ~ track
            i => 1 ~ preserve
        end
        s = instance(SCall)
        @test s.a' === 1.0
        @test s.b' === 2.0
    end

    @testset "unit" begin
        @system SCallUnit(Controller) begin
            fa(; x(u"m")) => x ~ call(u"m")
            a(fa) => fa(1u"m") ~ track(u"m")
            fb(i; x(u"m")) => i + x ~ call(u"m")
            b(fb) => fb(1u"m") ~ track(u"m")
            i => 1 ~ preserve(u"m")
        end
        s = instance(SCallUnit)
        @test s.a' === 1.0u"m"
        @test s.b' === 2.0u"m"
    end

    @testset "type and unit" begin
        @system SCallTypeUnit(Controller) begin
            fa(; x::Int(u"m")) => x ~ call::int(u"m")
            a(fa) => fa(1u"m") ~ track::int(u"m")
            fb(i; x::Int(u"m")) => i + x ~ call::int(u"m")
            b(fb) => fb(1u"m") ~ track::int(u"m")
            i => 1 ~ preserve::int(u"m")
        end
        s = instance(SCallTypeUnit)
        @test s.a' === 1u"m"
        @test s.b' === 2u"m"
        @test s.a' |> Cropbox.deunitfy isa Int
        @test s.b' |> Cropbox.deunitfy isa Int
    end
end


text:
#=
Copyright 2018 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    laodSimulationFiles(input_filename :: String)

Load, check, and return `.yml` input file content.
"""
function loadSimulationFiles(input_filename :: String)
    data = loadYAMLFile(input_filename)

    checkInputFile(data)

    return data
end

"""
    loadYAMLFile(filename :: String)

Check existence and open a `.yml` file and return the content as `Dict{Any,Any}`.
"""
function loadYAMLFile(filename :: String)
    if ~isfile(filename)
        error("missing file $filename")
    end

    return YAML.load(open(filename))
end

"""
    checkInputFile(input_filename :: String)

Check YAML input file. Raise errors if any of the important variables is ill
defined.
"""
function checkInputFile(data :: Dict{Any,Any})
    checkSections(data)
    checkNetwork(data["network"])
end

"""
    checkSections(data :: Dict{Any,Any})

Look for the four sections in the input data. Run integrity checks for `blood` and
`solver` sections.
"""
function checkSections(data :: Dict{Any,Any})
    keys = ["project name", "network", "blood", "solver"]
    for key in keys
        if ~haskey(data, key)
            error("missing section $key in YAML input file")
        end
    end

    checkSection(data, "blood", ["mu", "rho"])
    checkSection(data, "solver", ["Ccfl", "cycles", "convergence tolerance"])

    if ~haskey(data["solver"], "jump")
        data["solver"]["jump"] = 100
    end
end

"""
    checkSection(data :: Dict{Any,Any}, section :: String, keys :: Array{String,1})

Look for a list of keys in the given data section.
"""
function checkSection(data :: Dict{Any,Any}, section :: String, keys :: Array{String,1})
    for key in keys
        if ~haskey(data[section], key)
            error("missing $key in $section section")
        end
    end
end

"""
    checkNetwork(network :: Array{Dict{Any,Any},1})

Loop trough the network and run checks on each single vessel. Check also if at least one
inlet and one oulet has been defined.
"""
function checkNetwork(network :: Array{Dict{Any,Any},1})
    has_inlet = false
    has_outlet = false
    for i = 1:length(network)
        checkVessel(i, network[i])

        if haskey(network[i], "inlet")
            has_inlet = true
        end
        if haskey(network[i], "outlet")
            has_outlet = true
        end
    end

    if ~has_inlet
        error("missing inlet(s) definition")
    end

    if ~has_outlet
        error("missing outlet(s) definition")
    end
end

"""
    checkVessel(i :: Int, vessel :: Dict{Any,Any})

Check if all the important parameters are defined for the current vessel. If the vessel
has been indicated to be an inlet/outlet segment, check also for the definition of
boundary condition parameters.
"""
function checkVessel(i :: Int, vessel :: Dict{Any,Any})
    keys = ["label", "sn", "tn", "L", "E"]
    for key in keys
        if ~haskey(vessel, key)
            error("vessel $i is missing $key value")
        end
    end

    if ~haskey(vessel, "R0")
        if ~haskey(vessel, "Rp") && ~haskey(vessel, "Rd")
            error("vessel $i is missing lumen radius value(s)")
        end
    end

    if haskey(vessel, "inlet")
        if ~haskey(vessel, "inlet file")
            error("inlet vessel $i is missing the inlet file path")
        elseif ~isfile(vessel["inlet file"])
            file_path = vessel["inlet file"]
            error("vessel $i inlet file $file_path not found")
        end

        if ~haskey(vessel, "inlet number")
            error("inlet vessel $i is missing the inlet number")
        end
    end

    if haskey(vessel, "outlet")
        outlet = vessel["outlet"]
        if outlet == "wk3"
            if ~haskey(vessel, "R1") || ~haskey(vessel, "Cc")
                error("outlet vessel $i is missing three-element windkessel values")
            end
        elseif outlet == "wk2"
            if ~haskey(vessel, "R1") || ~haskey(vessel, "Cc")
                error("outlet vessel $i is missing two-element windkessel values")
            end
        elseif outlet == "reflection"
            if ~haskey(vessel, "Rt")
                error("outlet vessel $i is missing reflection coefficient value")
            end
        end
    end
end

"""
    makeResultsFolder(data :: Dict{Any,Any})

Create results folder and cd in.
"""
function makeResultsFolder(data :: Dict{Any,Any})
    project_name = data["project name"]
    r_folder = join([project_name, "_results"])

    if isdir(r_folder) == false
      mkdir(r_folder)
    end

    cd(r_folder)
end

"""
    buildBlood(project_constants :: Dict{Any,Any})

Return an instance of type `::Blood`.
"""
function buildBlood(blood_data :: Dict{Any,Any})
    mu = blood_data["mu"]
    rho = blood_data["rho"]
    rho_inv = 1.0/rho

    return Blood(mu, rho, rho_inv)
end

"""
    buildArterialNetwork(network :: Array{Dict{Any,Any},1}, blood :: Blood, jump :: Int)

Populate `vessels` and `edges` lists containing mechanical properties and topology og the
system, respectively.
"""
function buildArterialNetwork(network :: Array{Dict{Any,Any},1}, blood :: Blood,
                              jump :: Int)
    vessels = [buildVessel(1, network[1], blood, jump)]
    edges = zeros(Int, length(network), 3)
    @inbounds edges[1,1] = vessels[1].ID
    @inbounds edges[1,2] = vessels[1].sn
    @inbounds edges[1,3] = vessels[1].tn

    for i = 2:length(network)
        push!(vessels, buildVessel(i, network[i], blood, jump))
        @inbounds edges[i,1] = vessels[i].ID
        @inbounds edges[i,2] = vessels[i].sn
        @inbounds edges[i,3] = vessels[i].tn
    end

    return vessels, edges
end

"""
    buildVessel(vessel_data :: Dict{Any,Any}, blood :: Blood, jump :: Int)

Allocate arrays, set initial conditions, and create a new instance of `Vessel` with data
from the `.yml`.
"""
function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, jump :: Int)
    vessel_name = vessel_data["label"]
    sn = vessel_data["sn"]
    tn = vessel_data["tn"]
    L = vessel_data["L"]
    E = vessel_data["E"]

    Rp, Rd = computeRadii(vessel_data)
    Pext = getPext(vessel_data)
    M, dx, invDx, halfDx, invDxSq = meshVessel(vessel_data, L)
    h0 = initialiseThickness(vessel_data, M)
    outlet, Rt, R1, R2, Cc = addOutlet(vessel_data)
    viscT = computeViscousTerm(vessel_data, blood)
    inlet, heart = buildHeart(vessel_data)
    phi = getPhi(vessel_data)

    Q = zeros(Float64, M)
    P = zeros(Float64, M)
    A = zeros(Float64, M)
    u = zeros(Float64, M)
    c = zeros(Float64, M)
    A0 = zeros(Float64, M)
    R0 = zeros(Float64, M)
    s_A0 = zeros(Float64, M)
    beta = zeros(Float64, M)
    vA = zeros(Float64, M+2)
    vQ = zeros(Float64, M+2)
    Al = zeros(Float64, M+2)
    Ar = zeros(Float64, M+2)
    Ql = zeros(Float64, M+2)
    Qr = zeros(Float64, M+2)
    wallE = zeros(Float64, M)
    gamma = zeros(Float64, M)
    slope = zeros(Float64, M)
    wallVa = zeros(Float64, M)
    wallVb = zeros(Float64, M)
    inv_A0 = zeros(Float64, M)
    dU = zeros(Float64, 2, M+2)
    Fl = zeros(Float64, 2, M+2)
    Fr = zeros(Float64, 2, M+2)
    s_inv_A0 = zeros(Float64, M)
    slopesA = zeros(Float64, M+2)
    slopesQ = zeros(Float64, M+2)
    Q_t = zeros(Float64, jump, 6)
    P_t = zeros(Float64, jump, 6)
    A_t = zeros(Float64, jump, 6)
    u_t = zeros(Float64, jump, 6)
    c_t = zeros(Float64, jump, 6)
    Q_l = zeros(Float64, jump, 6)
    P_l = zeros(Float64, jump, 6)
    A_l = zeros(Float64, jump, 6)
    u_l = zeros(Float64, jump, 6)
    c_l = zeros(Float64, jump, 6)
    flux  = zeros(Float64, 2, M+2)
    uStar = zeros(Float64, 2, M+2)
    s_15_gamma = zeros(Float64, M)
    gamma_ghost = zeros(Float64, M+2)

    s_pi = sqrt(pi)
    s_pi_E_over_sigma_squared = s_pi*E/0.75
    one_over_rho_s_p = 1.0/(3.0*blood.rho*s_pi)
    radius_slope = computeRadiusSlope(Rp, Rd, L)

    ah = 0.2802
    bh = -5.053e2
    ch = 0.1324
    dh = -0.1114e2

    if h0 == 0.0
        Rmean = 0.5*(Rp + Rd)
        h0 = computeThickness(Rmean, ah, bh, ch, dh)
    end
    Cv = 0.5*s_pi*phi*h0/(blood.rho*0.75)

    @fastmath @inbounds for i = 1:M
      R0[i] = radius_slope*(i - 1)*dx + Rp
      A0[i] = pi*R0[i]*R0[i]
      s_A0[i] = sqrt(A0[i])
      inv_A0[i] = 1.0/A0[i]
      s_inv_A0[i] = sqrt(inv_A0[i])
      A[i] = A0[i]
      beta[i] = s_inv_A0[i]*h0*s_pi_E_over_sigma_squared
      gamma[i] = beta[i]*one_over_rho_s_p/R0[i]
      s_15_gamma[i] = sqrt(1.5*gamma[i])
      gamma_ghost[i+1] = gamma[i]
      P[i] = pressure(1.0, beta[i], Pext)
      c[i] = waveSpeed(A[i], gamma[i])
      wallE[i] = 3.0*beta[i]*radius_slope*inv_A0[i]*s_pi*blood.rho_inv
      if phi != 0.0
          wallVb[i] = Cv*s_inv_A0[i]*invDxSq
          wallVa[i] = 0.5*wallVb[i]
      end
    end

    gamma_ghost[1] = gamma[1]
    gamma_ghost[end] = gamma[end]

    if outlet == "wk2"
        R1, R2 = computeWindkesselInletImpedance(R2, blood, A0, gamma)
        outlet = "wk3"
    end

    U00A = A0[1]
    U01A = A0[2]
    UM1A = A0[M]
    UM2A = A0[M-1]

    U00Q = 0.0
    U01Q = 0.0
    UM1Q = 0.0
    UM2Q = 0.0

    W1M0 = u[end] - 4.0*c[end]
    W2M0 = u[end] + 4.0*c[end]

    node2 = convert(Int, floor(M*0.25))
    node3 = convert(Int, floor(M*0.5))
    node4 = convert(Int, floor(M*0.75))

    Pcn = 0.0

    last_A_name = join((vessel_name,"_A.last"))
    last_Q_name = join((vessel_name,"_Q.last"))
    last_u_name = join((vessel_name,"_u.last"))
    last_c_name = join((vessel_name,"_c.last"))
    last_P_name = join((vessel_name,"_P.last"))

    out_A_name = join((vessel_name,"_A.out"))
    out_Q_name = join((vessel_name,"_Q.out"))
    out_u_name = join((vessel_name,"_u.out"))
    out_c_name = join((vessel_name,"_c.out"))
    out_P_name = join((vessel_name,"_P.out"))

    return Vessel(vessel_name, ID, sn, tn, inlet, heart,
                  M, dx, invDx, halfDx,
                  beta, gamma, s_15_gamma, gamma_ghost,
                  A0, s_A0, inv_A0, s_inv_A0, Pext,
                  viscT, wallE, wallVa, wallVb,
                  A, Q, u, c, P,
                  A_t, Q_t, u_t, c_t, P_t,
                  A_l, Q_l, u_l, c_l, P_l,
                  W1M0, W2M0,
                  U00A, U00Q, U01A, U01Q, UM1A, UM1Q, UM2A, UM2Q,
                  last_P_name, last_Q_name, last_A_name,
                  last_c_name, last_u_name,
                  out_P_name, out_Q_name, out_A_name,
                  out_c_name, out_u_name,
                  node2, node3, node4,
                  Rt, R1, R2, Cc,
                  Pcn,
                  slope, flux, uStar, vA, vQ,
                  dU, slopesA, slopesQ,
                  Al, Ar, Ql, Qr, Fl, Fr,
                  outlet)
end

"""
    computeRadiusSlope(Rd :: Float64, Rp :: Float64, L :: Float64).

Calculate the slope for the lumen radius linear tapering as

(Rd - Rp)/L

"""
function computeRadiusSlope(Rp :: Float64, Rd :: Float64, L :: Float64)
    return (Rd - Rp)/L
end

"""
    computeThickness(R0i :: Float64)

Compute wall thickness based on local lumen radius as

h_{0i} = R_{0i} ( a_h e^{b_h R_{0i}} + c_h e^{d_h R_{0i}} )

where

- ah = 0.2802
- bh = -5.053e2 m^{-1}
- ch = 0.1324
- dh = -0.1114e2 m^{-1}

as reported in

> Avolio AP. Multi-branched model of the human arterial system. Medical and Biological Engineering and Computing. 1980 Nov 1;18(6):709-18.
"""
function computeThickness(R0i :: Float64,
                          ah :: Float64, bh :: Float64, ch :: Float64, dh :: Float64)
    return R0i*(ah*exp(bh*R0i) + ch*exp(dh*R0i))
end

"""
    computeRadii(vessel_data :: Dict{Any,Any})

If only a constant lumen radius is defined, return the same value for proximal and
distal variables, `Rp` and `Rd`, respectively.
"""
function computeRadii(vessel :: Dict{Any,Any})
    if ~haskey(vessel, "R0")
        Rp = vessel["Rp"]
        Rd = vessel["Rd"]
        return Rp, Rd
    else
        R0 = vessel["R0"]
        return R0, R0
    end
end

"""
    getPext(vessel :: Dict{Any,Any})

Extract Pext value for current vessels; return default `Pext = 0.0` if no value is
specified.
"""
function getPext(vessel :: Dict{Any,Any})
    if ~haskey(vessel, "Pext")
        return 0.0
    else
        return vessel["Pext"]
    end
end

"""
    getPhi(vessel :: Dict{Any,Any})

Extract `phi` value for current vessels; return default `phi = 0.0` if no value is
specified.
"""
function getPhi(vessel :: Dict{Any,Any})
    if ~haskey(vessel, "phi")
        return 0.0
    else
        return vessel["phi"]
    end
end

"""
    meshVessel(vessel :: Dict{Any,Any}, L :: Float64)

Pre-compute Delta x, frac{1}{Delta x} and frac{Delta x}{2} for the current
vessel. The Delta x is computed as

Delta x = frac{L}{M}

where `M` is the maximum value between `5` (the minimum needed by the solver), the value
defined in the `.yml`, and `ceil(L*1e3)` (which would make Delta x = 1mm).
"""
function meshVessel(vessel :: Dict{Any,Any}, L :: Float64)
    if ~haskey(vessel, "M")
        M = maximum([5, convert(Int, ceil(L*1e3))])
    else
        M = vessel["M"]
        M = maximum([5, M, convert(Int, ceil(L*1e3))])
    end

    dx = L/M
    invDx = M/L
    halfDx = 0.5*dx
    invDxSq = invDx*invDx

    return M, dx, invDx, halfDx, invDxSq
end

"""
    initialiseThickness(vessel :: Dict{Any,Any}, M :: Int)

If vessel thickness is not specified in the `.yml` file, return a `zeroes` array to be
filed with radius dependent values. Otherwise, return an array with the specified `h0`.
"""
function initialiseThickness(vessel :: Dict{Any,Any}, M :: Int)
    if ~haskey(vessel, "h0")
        return 0.0
    else
        return vessel["h0"]
    end
end

"""
    addOutlet(vessel :: Dict{Any,Any})

Parse outlet information for the current vessel and return windkessel and reflection
coeffiecient values.
"""
function addOutlet(vessel :: Dict{Any,Any})
    if haskey(vessel, "outlet")
        outlet = vessel["outlet"]
        if outlet == "wk3"
            Rt = 0.0
            R1 = vessel["R1"]
            R2 = vessel["R2"]
            Cc = vessel["Cc"]
        elseif outlet == "wk2"
            Rt = 0.0
            R1 = 0.0
            R2 = vessel["R1"]
            Cc = vessel["Cc"]
        elseif outlet == "reflection"
            Rt = vessel["Rt"]
            R1 = 0.0
            R2 = 0.0
            Cc = 0.0
        end
    else
        outlet = "none"
        Rt = 0.0
        R1 = 0.0
        R2 = 0.0
        Cc = 0.0
    end

    return outlet, Rt, R1, R2, Cc
end

"""
    computeViscousTerm(vessel_data :: Dict{Any,Any}, blood :: Blood)

Return

2(gamma_v + 2)pi*mu

where gamma_v (`gamma_profile`) is either specified in the vessel definition or
assumed equal to `9` (plug-flow).
"""
function computeViscousTerm(vessel_data :: Dict{Any,Any}, blood :: Blood)
    if haskey(vessel_data, "gamma_profile")
        gamma_profile = vessel_data["gamma_profile"]
    else
        gamma_profile = 9
    end
    return 2*(gamma_profile + 2)*pi*blood.mu*blood.rho_inv
end

"""
    buildHeart(vessel_data)

If the current vessel is an inlet vessel, return `true` flag and an `Heart` struct.
"""
function buildHeart(vessel :: Dict{Any,Any})
    if haskey(vessel, "inlet")
        inlet_type = vessel["inlet"]
        input_data = loadInletData(vessel["inlet file"])
        cardiac_period = input_data[end, 1]
        inlet_number = vessel["inlet number"]

        return true, Heart(inlet_type, cardiac_period, input_data, inlet_number)
    else
        return false, Heart("none", 0.0, zeros(1,2), 0)
    end
end

"""
    loadInletData(inlet_file :: String)

Read discretised inlet data from inlet file. Return an `Array{Float64, 2}` whose first
columun contains time variable and second column contains the inlet time function.
"""
function loadInletData(inlet_file :: String)
    return readdlm(inlet_file)
end

"""
    computeWindkesselInletImpedance(R2 :: Float64, blood :: Blood, A0 :: Array{Float64,1},
                                    gamma :: Array{Float64,1})

In case only one peripheral resistance is defined (two-element windkessel), the second
one is set as equal to the outlet vessel impedance.
"""
function computeWindkesselInletImpedance(R2 :: Float64, blood :: Blood,
    A0 :: Array{Float64,1}, gamma :: Array{Float64,1})

    R1 = blood.rho*waveSpeed(A0[end], gamma[end])/A0[end]
    R2 -= R1

    return R1, R2
end

# http://carlobaldassi.github.io/ArgParse.jl/stable/index.html
function parseCommandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "input_filename"
            help = ".yml input file name"
            required = true
        "--verbose", "-v"
            help = "Print STDOUT - default false"
            action = :store_true
        "--out_files", "-f"
            help = "Save complete results story rather than only the last cardiac cycle"
            action = :store_true
    end

    return parse_args(s)
end


text:
getpath(dataset::DatasetName) = @datadep_str name(dataset)

"""
    getfilename(path)

Get name of the file without filename extension.
`path` must include filename.
"""
function getfilename(path)
    filename = basename(path)
    return splitext(filename)[1]
end

"""
    df_or_array(df::DataFrame, toarray::Bool; cols::Int=1)

Based on `toarray` return either a `DataFrame` (false) or an `Array` (true).
"""
function df_or_array(df::DataFrame, toarray::Bool; cols::Int=1)
    return toarray ? df_to_array(df, cols=cols) : df
end

"""
    df_to_array(df::DataFrame; cols::Int=1)

Converts `DataFrame` to `Array` (1:end-1 columns) and `Vector` (last column).
"""
function df_to_array(df::DataFrame; cols::Int=1)
    if cols >= Base.size(df,2) || cols <= 0
        error("cols has to be bigger than 1 and smaller than number of all columns of DataFrame.")
    end
        return Array(df[:,1:end-cols]), Array(df[:,end-cols+1:end])
end

"""
    load(dataset::DatasetName, type::Symbol; kwargs...)

Loads dataset of given type.

- `type::Type{<:Split}=Train`: other possible type is `Valid` or `Test`.
"""
load(dataset::DatasetName, type::Type{<:Split}=Train) = getdata(dataset, type)

"""
    getdata(dataset::DatasetName, type::Type{<:Split})

Method not specified for type `DatasetName`. `type` is either `Train`, `Test` or `Valid`.
"""
getdata(dataset::DatasetName, type::Type{<:Split}) = error("Method getdata not specified for $(typeof(dataset)).")

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
# Integrating factor

intfact(x, a) = exp(-2a)besseli(x,2a)

intfact(x, y,a) = exp(-4a)besseli(x,2a)besseli(y,2a)

intfact(x, y, z, a) = exp(-6a)besseli(x,2a)besseli(y,2a)besseli(z,2a)

"""
    plan_intfact(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])

Constructor to set up an operator for evaluating the integrating factor with
real-valued parameter `a`. This can then be applied with the `*` operation on
data of the appropriate size.

The `dims` argument can be replaced with `w::Nodes` to specify the size of the
domain.

# Example

```jldoctest
julia> w = Nodes(Dual,(6,6));

julia> w[4,4] = 1.0;

julia> E = plan_intfact(1.0,(6,6))
Integrating factor with parameter 1.0 on a (nx = 6, ny = 6) grid

julia> E*w
Nodes{Dual,6,6,Float64} data
Printing in grid orientation (lower left is (1,1))
6×6 Array{Float64,2}:
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00888233   0.028765    0.0664124   0.0951774   0.0664124   0.028765
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.000828935  0.00268447  0.00619787  0.00888233  0.00619787  0.00268447
```
"""
function plan_intfact end

"""
    plan_intfact!(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])

Same as [`plan_intfact`](@ref), but the resulting operator performs an in-place
operation on data.
"""
function plan_intfact! end

struct IntFact{NX, NY, a, inplace}
    conv::Union{CircularConvolution{NX, NY},Nothing}
end

for (lf,inplace) in ((:plan_intfact,false),
                     (:plan_intfact!,true))

    @eval function $lf(a::Real,dims::Tuple{Int,Int};fftw_flags = FFTW.ESTIMATE)
        NX, NY = dims

        if a == 0
          return IntFact{NX, NY, 0.0, $inplace}(nothing)
        end

        #qtab = [intfact(x, y, a) for x in 0:NX-1, y in 0:NY-1]
        Nmax = 0
        while abs(intfact(Nmax,0,a)) > eps(Float64)
          Nmax += 1
        end
        qtab = [max(x,y) <= Nmax ? intfact(x, y, a) : 0.0 for x in 0:NX-1, y in 0:NY-1]
        #IntFact{NX, NY, a, $inplace}(Nullable(CircularConvolution(qtab, fftw_flags)))
        IntFact{NX, NY, a, $inplace}(CircularConvolution(qtab, fftw_flags))
      end

      @eval $lf(a::Real,nodes::Nodes{T,NX,NY}; fftw_flags = FFTW.ESTIMATE) where {T<:CellType,NX,NY} =
          $lf(a,node_inds(T,(NX,NY)), fftw_flags = fftw_flags)

end

function Base.show(io::IO, E::IntFact{NX, NY, a, inplace}) where {NX, NY, a, inplace}
    nodedims = "(nx = $NX, ny = $NY)"
    isinplace = inplace ? "In-place integrating factor" : "Integrating factor"
    print(io, "$isinplace with parameter $a on a $nodedims grid")
end

function mul!(out::Nodes{T,NX, NY},
                   E::IntFact{MX, MY, a, inplace},
                   s::Nodes{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, a, inplace}

    mul!(out.data, E.conv, s.data)
    out
end

function mul!(out::Nodes{T,NX, NY},
                   E::IntFact{MX, MY, 0.0, inplace},
                   s::Nodes{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}
    out .= deepcopy(s)
end

*(E::IntFact{MX,MY,a,false},s::Nodes{T,NX,NY}) where {MX,MY,a,T <: CellType, NX,NY} =
  mul!(Nodes(T,s), E, s)

*(E::IntFact{MX,MY,a,true},s::Nodes{T,NX,NY}) where {MX,MY,a,T <: CellType, NX,NY} =
    mul!(s, E, deepcopy(s))


text:
using Dates
using Base.Threads

"""
    struct KeyHeader

Helper struct to store header information of .key files.

# Fields

| Name                   | Meaning                                        |
| :---                   | :------                                        |
| `nkeys_file`           | Number of keys in this .key file               |
| `domain_corners`       | Corners of the domains defined by these keys   |
| `domain_fac`           | Factor needed for reconstructung int positions |
| `bits`                 | Size of PH keys it bits                        |
| `nkeys_total`          | Total number of keys in all files              |
| `nkeys_total_highword` | Total number of keys in all files              |
"""
struct KeyHeader
    nkeys_file::Vector{Int32}
    domain_corners::Vector{Float64}
    domain_fac::Float64
    bits::Int32
    nkeys_total::Vector{UInt32}
    nkeys_total_highword::Vector{UInt32}
end

"""
    read_keyheader(filename::String)

Reads the header of a .key file.
"""
function read_keyheader(filename::String)

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize == 8
        swap = 0
        snap_format = 2
    elseif blocksize == 256
        swap = 0
        snap_format = 1
    else
        blocksize = bswap(blocksize)
        if blocksize == 8
            swap = 1
            snap_format = 2
        elseif blocksize == 256
            swap = 1
            snap_format = 1
        else
            error("incorrect file format encountered when reading header of $filename")
        end
    end

    if snap_format == 2
        seek(f, 16)
        skip_line = read(f, Int32)
    end

    nkeys = read!(f, Array{Int32,1}(undef,6))
    domain_corners = read!(f, Array{Float64,1}(undef,3))
    domain_fac = read(f, Float64)

    read!(f, Array{Float64,1}(undef,2))
    read(f, Float64)
    read(f, Float64)
    read(f, Int32)

    bits = read(f, Int32)
    nkeys_total = read!(f, Array{UInt32,1}(undef,6))

    read(f, Int32)
    read(f, Int32)
    read(f, Float64)
    read(f, Float64)
    read(f, Float64)
    read(f, Float64)
    read(f, Int32)
    read(f, Int32)

    nkeys_total_highword = read!(f, Array{UInt32,1}(undef,6))

    close(f)

    return KeyHeader( nkeys, domain_corners, domain_fac, bits, nkeys_total, nkeys_total_highword )
end

"""
    read_key_index(file_key_index::String)

Reads the .key.index file.
"""
function read_key_index(file_key_index::String)

    finfo = stat(file_key_index)

    f = open(file_key_index)
    n = read(f, Int32)

    if finfo.size == 4+n*8*2+n*4
        int_type = UInt64
    else
        int_type = UInt32
    end

    low_list = read!(f, Array{int_type,1}(undef,n))
    high_list = read!(f, Array{int_type,1}(undef,n))
    file_list = read!(f, Array{UInt32,1}(undef,n))

    close(f)

    return low_list, high_list, file_list
end