text:
struct OrbitGrid{T}
    energy::AbstractVector{T}
    pitch::AbstractVector{T}
    r::AbstractVector{T}
    counts::Vector{Int}
    orbit_index::Array{Int,3}
    class::Array{Symbol,3}
    tau_p::Array{T,3}
    tau_t::Array{T,3}
end

function Base.show(io::IO, og::OrbitGrid)
    print(io, "OrbitGrid: $(length(og.energy))×$(length(og.pitch))×$(length(og.r)):$(length(og.counts))")
end

function orbit_grid(M::AbstractEquilibrium, eo::AbstractVector, po::AbstractVector, ro::AbstractVector;
                    q = 1, amu = H2_amu, kwargs...)

    nenergy = length(eo)
    npitch = length(po)
    nr = length(ro)

    orbit_index = zeros(Int,nenergy,npitch,nr)
    class = fill(:incomplete,(nenergy,npitch,nr))
    tau_t = zeros(Float64,nenergy,npitch,nr)
    tau_p = zeros(Float64,nenergy,npitch,nr)

    norbs = nenergy*npitch*nr
    subs = CartesianIndices((nenergy,npitch,nr))

    p = Progress(norbs)
    channel = RemoteChannel(()->Channel{Bool}(norbs), 1)
    orbs = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            orbs = @distributed (vcat) for i=1:norbs
                ie,ip,ir = Tuple(subs[i])
                c = EPRCoordinate(M,eo[ie],po[ip],ro[ir],q=q,amu=amu)
                try
                    o = get_orbit(M, c; kwargs...)
                catch
                    o = Orbit(EPRCoordinate(;q=q,amu=amu),:incomplete)
                end

                if o.class in (:incomplete,:invalid,:lost)
                    o = Orbit(o.coordinate,:incomplete)
                end
                put!(channel, true)
                o
            end
            put!(channel, false)
            orbs
        end
    end)

    for i=1:norbs
        class[subs[i]] = orbs[i].class
        tau_p[subs[i]] = orbs[i].tau_p
        tau_t[subs[i]] = orbs[i].tau_t
    end

    grid_index = filter(i -> orbs[i].class != :incomplete, 1:norbs)
    orbs = filter(x -> x.class != :incomplete, orbs)
    norbs = length(orbs)
    orbit_index[grid_index] = 1:norbs

    return orbs, OrbitGrid(eo,po,ro,fill(1,norbs),orbit_index,class,tau_p,tau_t)

end

function segment_orbit_grid(M::AbstractEquilibrium, orbit_grid::OrbitGrid, orbits::Vector;
                        norbits=1000, combine=(length(orbits[1].path) != 0),
                        q = 1, amu = H2_amu, kwargs...)

    eo = orbit_grid.energy
    po = orbit_grid.pitch
    ro = orbit_grid.r

    nenergy = length(eo)
    npitch = length(po)
    nr = length(ro)
    norbs = length(orbits)

    e_range = extrema(eo)
    p_range = extrema(po)
    r_range = extrema(ro)
    orbs_index = zeros(Int,norbs)
    for i = 1:length(orbit_grid.orbit_index)
        ii = orbit_grid.orbit_index[i]
        ii == 0 && continue
        orbs_index[ii] != 0 && continue
        orbs_index[ii] = i
    end

    orbit_index = zeros(Int,nenergy,npitch,nr)

    norm = abs.([-(e_range...), -(p_range...), -(r_range...)])
    mins = [e_range[1],p_range[1],r_range[1]]
    oclasses = [:potato, :stagnation, :trapped, :co_passing, :ctr_passing]

    orbit_counts = Dict{Symbol,Int}(o=>count(x -> x.class == o, orbits)
                                    for o in oclasses)

    nclusters = 0
    orbs = eltype(orbits)[]
    orbit_num = 0
    for oc in oclasses
        nk = max(Int(ceil(norbits*orbit_counts[oc]/norbs)),1)
        if (nclusters + nk) > norbits
            nk = norbits - nclusters
        end
        nk == 0 && continue

        inds_oc = findall([o.class == oc for o in orbits])
        coords = hcat((([o.coordinate.energy, o.coordinate.pitch, o.coordinate.r] .- mins)./norm
                       for o in orbits if o.class == oc)...)

        if nk == 1
            if !combine
                c = coords .* norm .+ mins
                cc = EPRCoordinate(M,mean(c,dims=2)...;q = q, amu=amu)
                try
                    o = get_orbit(M, GCParticle(cc.energy,cc.pitch,cc.r,cc.z,cc.m,cc.q); kwargs...)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                catch
                    o = Orbit(cc)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                end
            else
                o = combine_orbits(orbits[inds_oc])
                push!(orbs,o)
                orbit_num = orbit_num + 1
            end
            orbit_index[orbs_index[inds_oc]] .= orbit_num
            nclusters = nclusters + nk
            continue
        end

        k = kmeans(coords,nk)
        if !combine
            coords = k.centers.*norm .+ mins
            for i=1:size(coords,2)
                w = k.assignments .== i
                sum(w) == 0 && continue
                cc = EPRCoordinate(M,coords[1,i],coords[2,i],coords[3,i], q=q, amu=amu)
                try
                    o = get_orbit(M, GCParticle(cc.energy,cc.pitch,cc.r,cc.z,cc.m,cc.q); kwargs...)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                catch
                    o = Orbit(cc)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                end
                orbit_index[orbs_index[inds_oc[w]]] .= orbit_num
            end
        else
            for i=1:nk
                w = k.assignments .== i
                sum(w) == 0 && continue
                o = combine_orbits(orbits[inds_oc[w]])
                push!(orbs,o)
                orbit_num = orbit_num + 1
                orbit_index[orbs_index[inds_oc[w]]] .= orbit_num
            end
        end
        nclusters = nclusters + nk
    end

    counts = [count(x -> x == i, orbit_index) for i=1:length(orbs)]
    return orbs, OrbitGrid(eo,po,ro,counts,orbit_index,orbit_grid.class,orbit_grid.tau_p,orbit_grid.tau_t)

end

"""
    get_orbel_volume(myOrbitGrid, equidistant)

Get the orbit element volume of this specific orbit-grid.
If equidistant=false, then return a 3D array with all the orbit volumes.
Otherwise, just a scalar.
"""
function get_orbel_volume(og::OrbitGrid, equidistant::Bool)

    if equidistant
        dRm = abs(og.r[2]-og.r[1])
        dE = abs(og.energy[2]-og.energy[1])
        dpm = abs(og.pitch[2]-og.pitch[1])
        dO = dE*dpm*dRm
        return dO
    else
        return get3DVols(og.energy, og.pitch, og.r)
    end
end

"""
    get3DDiffs(E, pm, Rm)

Calculate and return 3D-array of the diffs of orbit-space coordinates. Assume edge diff to be same as next-to-edge diff.
"""
function get3DDiffs(E, pm, Rm)
    dE = vcat(abs.(diff(E)),abs(E[end]-E[end-1]))
    dE = reshape(dE,length(dE),1,1)
    dE3D = repeat(dE,1,length(pm),length(Rm))
    dpm = vcat(abs.(diff(pm)),abs(pm[end]-pm[end-1]))
    dpm = reshape(dpm,1,length(dpm),1)
    dpm3D = repeat(dpm,length(E),1,length(Rm))
    dRm = vcat(abs.(diff(Rm)),abs(Rm[end]-Rm[end-1]))
    dRm = reshape(dRm,1,1,length(dRm))
    dRm3D = repeat(dRm,length(E),length(pm),1)

    return dE3D, dpm3D, dRm3D
end

"""
    get3DVols(E, pm, Rm)

This function will calculate the volume of all the hyper-voxels pertaining to the 3D grid. It assumes the hyper-voxels pertaining to the
upper-end (edge) grid-points have the same volumes as the hyper-voxels just inside of them. It return a 3D array, containing all the hyper-voxel
volumes. The 3D array will have size()=(length(E), length(pm), length(Rm)). The function assumes a rectangular 3D grid.
"""
function get3DVols(E, pm, Rm)

    # Safety-check to ensure vectors
    if !(1==length(size(E))==length(size(pm))==length(size(Rm)))
        throw(ArgumentError("E, pm, Rm inputs are not all vectors. Please correct and re-try."))
    end

    dE3D, dpm3D, dRm3D = get3DDiffs(E, pm, Rm)
    dvols = dE3D .*dpm3D .*dRm3D
    return dvols
end

"""
    get4DDiffs(E, p, R, Z)

Calculate and return 4D-arrays of the diffs of particle-space coordinates. Assume edge diff to be same as next-to-edge diff.
"""
function get4DDiffs(E, p, R, Z)
    dR = vcat(abs.(diff(R)),abs(R[end]-R[end-1]))
    dR = reshape(dR,1,1,length(dR),1)
    dR4D = repeat(dR,length(E),length(p),1,length(Z))
    dZ = vcat(abs.(diff(Z)),abs(Z[end]-Z[end-1]))
    dZ = reshape(dZ,1,1,1,length(dZ))
    dZ4D = repeat(dZ,length(E),length(p),length(R),1)
    dE = vcat(abs.(diff(E)),abs(E[end]-E[end-1]))
    dE = reshape(dE,length(dE),1,1,1)
    dE4D = repeat(dE,1,length(p),length(R),length(Z))
    dp = vcat(abs.(diff(p)),abs(p[end]-p[end-1]))
    dp = reshape(dp,1,length(dp),1,1)
    dp4D = repeat(dp,length(E),1,length(R),length(Z))

    return dE4D, dp4D, dR4D, dZ4D
end

"""
    get4DVols(E, p, R, Z)

This function will calculate the volume of all the hyper-voxels pertaining to the 4D grid. It assumes the hyper-voxels pertaining to the
upper-end (edge) grid-points have the same volumes as the hyper-voxels just inside of them. It return a 4D array, containing all the hyper-voxel
volumes. The 4D array will have size()=(length(E), length(p), length(R), length(Z)). The function assumes a rectangular 4D grid.
"""
function get4DVols(E, p, R, Z)

    # Safety-check to ensure vectors
    if !(1==length(size(E))==length(size(p))==length(size(R))==length(size(Z)))
        throw(ArgumentError("E, p, R, Z inputs are not all vectors. Please correct and re-try."))
    end

    dE4D, dp4D, dR4D, dZ4D = get4DDiffs(E, p, R, Z)
    dvols = dE4D .*dp4D .*dR4D .*dZ4D
    return dvols
end

function orbit_matrix(M::AbstractEquilibrium, grid::OrbitGrid, energy, pitch, r, z; kwargs...)
    nenergy = length(energy)
    npitch = length(pitch)
    nr = length(r)
    nz = length(z)
    norbits = length(grid.counts)
    subs = CartesianIndices((nenergy,npitch,nr,nz))
    nsubs = length(subs)

    p = Progress(nsubs)
    channel = RemoteChannel(()->Channel{Bool}(nsubs), 1)
    R = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            R = @distributed (hcat) for i=1:nsubs
                ie,ip,ir,iz = Tuple(subs[i])
                gcp = GCParticle(energy[ie],pitch[ip],r[ir],z[iz])
                o = get_orbit(M,gcp;store_path=false,kwargs...)
                Rcol = spzeros(norbits)
                if !(o.class in (:lost,:incomplete)) && o.coordinate.r > magnetic_axis(M)[1]
                    oi = orbit_index(grid,o.coordinate)
                    (oi > 0) && (Rcol[oi] = 1.0)
                end
                put!(channel,true)
                Rcol
            end
            put!(channel,false)
            R
        end
    end)
    return R
end

function write_orbit_grid(grid::OrbitGrid;filename="orbit_grid.h5")
    h5open(filename,"w") do file
        file["energy"] = collect(grid.energy)
        file["pitch"] = collect(grid.pitch)
        file["r"] = collect(grid.r)
        file["counts"] = grid.counts
        file["orbit_index"] = grid.orbit_index
        file["class"] = String.(grid.class)
        file["tau_p"] = grid.tau_p
        file["tau_t"] = grid.tau_p
    end
    nothing
end

function read_orbit_grid(filename)
    isfile(filename) || error("File does not exist")

    f = h5open(filename)
    energy = read(f["energy"])
    pitch = read(f["pitch"])
    r = read(f["r"])
    counts = read(f["counts"])
    orbit_index = read(f["orbit_index"])
    class = Symbol.(read(f["class"]))
    tau_p = read(f["tau_p"])
    tau_t = read(f["tau_t"])
    close(f)

    return OrbitGrid(energy,pitch,r,counts,orbit_index,class,tau_p,tau_t)
end

function map_orbits(grid::OrbitGrid, f::Vector)
    if length(grid.counts) != length(f)
        throw(ArgumentError("Incompatible sizes"))
    end
    dorb = abs((grid.r[2]-grid.r[1])*(grid.energy[2]-grid.energy[1])*(grid.pitch[2]-grid.pitch[1]))
    return [i == 0 ? zero(f[1]) : f[i]/(grid.counts[i]*dorb) for i in grid.orbit_index]
end

"""
    map_orbits(og, f, equidistant)

Unlike regular map_orbits, take non-equidistant 3D grid-points into consideration.
"""
function map_orbits(grid::OrbitGrid, f::Vector, os_equidistant::Bool)
    if length(grid.counts) != length(f)
        throw(ArgumentError("Incompatible sizes"))
    end
    if os_equidistant
        dorb = abs((grid.r[2]-grid.r[1])*(grid.energy[2]-grid.energy[1])*(grid.pitch[2]-grid.pitch[1])) # Scalar
    else
        dorb = get3DVols(og.energy, og.pitch, og.r) # 3D matrix
    end

    if os_equidistant
        return [i == 0 ? zero(f[1]) : f[i]/(grid.counts[i]*dorb) for i in grid.orbit_index]
    else
        return [i == 0 ? zero(f[1]) : f[i]/(grid.counts[i]*dorb[i]) for i in grid.orbit_index]
    end
end

function orbit_index(grid::OrbitGrid, o::EPRCoordinate; nearest=false)

    if !nearest
        if (grid.energy[1] <= o.energy <= grid.energy[end]) &&
           (grid.pitch[1] <= o.pitch <= grid.pitch[end]) &&
           (grid.r[1] <= o.r <= grid.r[end])

            i = argmin(abs.(o.energy .- grid.energy))
            j = argmin(abs.(o.pitch .- grid.pitch))
            k = argmin(abs.(o.r .- grid.r))
            ind = grid.orbit_index[i,j,k]
        else
            ind = 0
        end
    else
        inds = filter(x->grid.orbit_index[x] != 0, CartesianIndices(grid.orbit_index))
        data = hcat( ([grid.energy[I[1]], grid.pitch[I[2]], grid.r[I[3]]] for I in inds)...)
        tree = KDTree(data)
        idxs, dists = knn(tree,[o.energy,o.pitch,o.r],1,false)
        ind = grid.orbit_index[inds[idxs[1]]]
    end

    return ind
end

function bin_orbits(grid::OrbitGrid, orbits; weights::Union{Nothing,Vector},nearest=false)

    if weights != nothing
        length(weights) == length(orbits) || error("Incompatible weight vector size")
        w = weights
    else
        w = fill(1.0,length(orbits))
    end

    f = zeros(length(grid.counts))
    if !nearest
        for (io,o) in enumerate(orbits)
            i = argmin(abs.(o.energy .- grid.energy))
            j = argmin(abs.(o.pitch .- grid.pitch))
            k = argmin(abs.(o.r .- grid.r))
            ind = grid.orbit_index[i,j,k]
            if ind != 0
                f[ind] += w[io]
            end
        end
    else
        inds = filter(x->grid.orbit_index[x] != 0, CartesianIndices(grid.orbit_index))
        data = hcat( ([grid.energy[I[1]], grid.pitch[I[2]], grid.r[I[3]]] for I in inds)...)
        tree = KDTree(data)
        for (io, o) in enumerate(orbits)
            idxs, dists = knn(tree,[o.energy,o.pitch,o.r],1,false)
            f[grid.orbit_index[inds[idxs[1]]]] += w[io]
        end
    end

    return f
end

function combine_orbits(orbits)
    norbits = length(orbits)
    norbits == 1 && return orbits[1]

    o = orbits[1]
    r = o.path.r
    z = o.path.z
    phi = o.path.phi
    pitch = o.path.pitch
    energy = o.path.energy
    dt = o.path.dt
    #dl = o.path.dl

    c = o.coordinate
    isa(c, EPRCoordinate) || error("Wrong orbit coordinate. Expected EPRCoordinate")
    ec = c.energy
    pc = c.pitch
    rc = c.r
    zc = c.z
    tau_p = o.tau_p
    tau_t = o.tau_t

    for i=2:norbits
        oo = orbits[i]
        ec = ec + oo.coordinate.energy
        pc = pc + oo.coordinate.pitch
        rc = rc + oo.coordinate.r
        zc = zc + oo.coordinate.z
        tau_t = tau_t + oo.tau_t
        tau_p = tau_p + oo.tau_p
        append!(r, oo.path.r)
        append!(z, oo.path.z)
        append!(phi, oo.path.phi)
        append!(pitch, oo.path.pitch)
        append!(energy, oo.path.energy)
        append!(dt, oo.path.dt)
        #append!(dl, oo.path.dl)
    end
    ec = ec/norbits
    pc = pc/norbits
    rc = rc/norbits
    zc = zc/norbits
    tau_p = tau_p/norbits
    tau_t = tau_p/norbits

    cc = EPRCoordinate(ec,pc,rc,zc,zero(zc),c.m,c.q)
    path = OrbitPath(o.path.vacuum, o.path.drift, energy,pitch,r,z,phi,dt)#,dl)

    if all(x -> x.class == orbits[1].class, orbits)
        class = orbits[1].class
    else
        class = :meta
    end

    return Orbit(cc, class, tau_p, tau_t, path)
end

function mc2orbit(M::AbstractEquilibrium, d::FIDASIMGuidingCenterParticles, GCP::T;  kwargs...) where T <: Function
    p = Progress(d.npart)
    channel = RemoteChannel(()->Channel{Bool}(d.npart),1)

    t = @sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end

        @async begin
            orbs = @distributed (vcat) for i=1:d.npart
                o = get_orbit(M,GCP(d.energy[i],B0Ip_sign(M)*d.pitch[i],d.r[i]/100,d.z[i]/100); kwargs...,store_path=false)
                put!(channel,true)
                o.coordinate
            end
            put!(channel,false)
            orbs
        end
    end
    return fetch(t)
end

function fbm2orbit(M::AbstractEquilibrium,d::FIDASIMGuidingCenterFunction; GCP=GCDeuteron, n=1_000_000, kwargs...)
    dmc = fbm2mc(d,n=n)
    return mc2orbit(M, dmc, GCP; kwargs...)
end

"""
    removeBadIndices!(dmc, energy, pitch, R, Z)
    removeBadIndices!(-||-, verbose=true)

Remove bad samples from the dmc (FIDASIMGuidingCenterParticles). Samples that might have ended up outside of the original
(E,p,R,Z)-ranges.
"""
function removeBadIndices!(dmc, energy, pitch, R, Z; verbose=false)
    if verbose
        println("Removing bad samples... ")
    end
    bad_indices_E = findall(x-> x<=minimum(energy) || x>=maximum(energy),dmc.energy)
    bad_indices_p = findall(x-> x<=minimum(pitch) || x>=maximum(pitch),dmc.pitch)
    bad_indices_R = findall(x-> x<=minimum(R) || x>=maximum(R),dmc.r)
    bad_indices_Z = findall(x-> x<=minimum(Z) || x>=maximum(Z),dmc.z)
    bad_indices_tot = sort!(unique(vcat(bad_indices_E,bad_indices_p,bad_indices_R,bad_indices_Z)))
    deleteat!(dmc.energy,bad_indices_tot)
    deleteat!(dmc.pitch,bad_indices_tot)
    deleteat!(dmc.r,bad_indices_tot)
    deleteat!(dmc.z,bad_indices_tot)
    deleteat!(dmc.weight,bad_indices_tot)
    deleteat!(dmc.class,bad_indices_tot)
    dmc = FIDASIMGuidingCenterParticles(length(dmc.energy), length(unique(dmc.class)), dmc.class, dmc.weight, dmc.r, dmc.z, dmc.energy, dmc.pitch, sum(dmc.weight))
end

struct OrbitSpline{T<:Function}
    n::Int
    itp::T
end

@inline (os::OrbitSpline)(x) = os.itp(x)
@inline Base.length(os::OrbitSpline) = os.n

function OrbitSpline(p::OrbitPath, t)
    length(p) == 0 && return OrbitSpline(0, x -> S4(zeros(4)))
    eprz = hcat(p.energy,p.pitch,p.r,p.z)
    oi = scale(interpolate(eprz, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp())), t, 1:4)
    return OrbitSpline(length(t), x -> S4(oi.(x,1:4)))
end

function OrbitSpline(o::Orbit, t)
    return OrbitSpline(o.path, t)
end

function OrbitSpline(o::Orbit)
    t = range(0.0, 1.0, length=length(o))
    return OrbitSpline(o.path, t)
end

function OrbitSpline(p::OrbitPath)
    t = range(0.0, 1.0, length=length(p))
    return OrbitSpline(p, t)
end

OrbitSpline(o::OrbitSpline) = o
OrbitSpline(o::OrbitSpline, t) = o


text:
# Copyright (c) 2014, Dreki Þórgísl <dreki@billo.systems>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#  * Redistributions of source code must retain the above copyright notice,
#    this list of conditions and the following disclaimer.
#  * Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions and the following disclaimer in the documentation
#    and/or other materials provided with the distribution.
#  * Neither the name of the copyright holders nor the names of its
#    contributors may be used to endorse or promote products derived from this
#    software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
module Encode

export encode, encodeterm,
charintpack, charint4pack, charint2pack, charsignedint4pack

include("Tags.jl")
include("Util.jl")

function encode(term; compressed=false)
    encoded = encodeterm(term)
    if compressed
        encoded = vcat(compressedtag, compressterm(encoded, compressed))
    end
    vcat(version, encoded)
end

function encodeboolorsym(term)
    str = string(term)
    vcat(atomtag, charint2pack(length(str)), convert(Array{Uint8}, str))
end

function encodeterm(term::Symbol)
    encodeboolorsym(term)
end

function encodeterm(term::Nothing)
    encodeboolorsym(:nothing)
end

function encodeterm(term::Bool)
    encodeboolorsym(term)
end

function encodeterm(term::Array{Uint8,1})
    len = length(term)
    if len == 0
        return niltag
    elseif len <= typemax(Uint16)
        return vcat(stringtag, charint2pack(len), term)
    elseif len > typemax(Uint32)
        throw(InvalidListLength(len))
    return vcat(listtag, charint4pack(len), map(encodeterm, term), niltag)
    end
end

function encodeterm(term::UTF8String)
    encodeterm(convert(Array{Uint8}, term))
end

function encodeterm(term::Tuple)
    len = length(term)
    if len <= typemax(Uint8)
        header = vcat(smalltupletag, convert(Uint8, len))
    elseif arity <= typemax(Uint32)
        header = charint4pack(arity)
    else
        throw(InvalidTupleArity(arity))
    end
    vcat(header, map(encodeterm, term))
end

function encodeterm(term::Integer)
    if 0 <= term <= typemax(Uint8)
        return vcat(smallinttag, uint8(term))
    elseif typemin(Int32) <= term <= typemax(Int32)
        return vcat(inttag, charsignedint4pack(term))
    end
    if term >= 0
        sign = 0
    else
        sign = 1
        term = -term
    end
    bytes = charintpack(term)
    len = length(bytes)
    if len <= typemax(Uint8)
        return vcat(smallbiginttag, len, sign, bytes)
    elseif len <= typemax(Uint32)
        return vcat(largebiginttag, charint4pack(len), sign, bytes)
    end
    msg = "Got length: $len"
    throw(InvalidIntLength(msg))
end

function encodeterm(term::Float64)
end

end

ext: jl


text:
module TransformUtils

import Base: convert, promote_rule, *, transpose

export
  Quaternion,
  AngleAxis,
  AxisAngle,
  SO3,
  so3,
  SE3,
  Euler,
  skew,
  vee!,
  vee,
  *,
  transpose,
  matrix,
  inverse,
  compare,
  normalize!,
  normalize,
  q_conj,
  q_conj!,
  convert,
  rotate!,
  rotate,
  wrapRad,
  logmap,
  rightJacExmap,
  rightJacExmapinv,
  # Should be good to go
  veeAngleAxis,
  veeQuaternion,
  veeEuler,

  # type aliases
  FloatInt,
  VectorFloatInt,

  # basic Taylor exported to try against Pade version in expm
  expmOwn,
  expmOwn1,
  expmOwnT,
  expmOwn2,
  expmOwn3,
  expmOwn4,

  # needs refactoring
  SE2,
  se2vee,
  se2vee!

  function skew(v::Array{Float64,1})
      S = zeros(3,3)
      S[1,:] = [0., -v[3], v[2]]
      S[2,:] = [v[3],0.,-v[1]]
      S[3,:] = [-v[2],v[1],0.]
      return S
  end

typealias FloatInt Union{Float64,Int}
typealias VectorFloatInt Union{Vector{Float64},Vector{Int}}

type Quaternion
    s::Float64
    v::Array{Float64,1}
    Quaternion() = new()
    Quaternion(s::FloatInt) = new(1.0,zeros(3))
    Quaternion(s::FloatInt,v::VectorFloatInt) = new(s,v)
end

type AngleAxis
    theta::Float64
    ax::Array{Float64,1}
    AngleAxis() = new()
    AngleAxis(s::FloatInt) = new(0,[1;0;0])
    AngleAxis(s::FloatInt,v::VectorFloatInt) = new(s,v)
end

typealias AxisAngle AngleAxis

type SO3
    R::Array{Float64,2}
    SO3() = new()
    SO3(dummy::FloatInt) = new(eye(3))
    SO3(r::Array{Float64,2}) = new(r)
end

type so3
    S::Array{Float64,2}
    so3() = new()
    so3(s::FloatInt) = new(zeros(3,3))
    so3(v::VectorFloatInt) = new(skew(v))
    so3(S::Array{Float64,2}) = new(S)
end

type Euler
    R::Float64
    P::Float64
    Y::Float64
    Euler() = new()
    Euler(s::FloatInt) = new(0.0,0.0,0.0)
    Euler(r::FloatInt,p::FloatInt,y::FloatInt) = new(r,p,y)
    Euler(v::VectorFloatInt) = new(v[1],v[2],v[3])
end

type SE3
  R::SO3
  t::Vector{Float64}
  SE3() = new()
  SE3(dummy::FloatInt) = new(SO3(0.0), zeros(3))
  SE3(t::VectorFloatInt, r::SO3) = new(r,t)
  SE3(v::VectorFloatInt, E::Euler) = new(convert(SO3,E), v[1:3] )
  SE3(v::VectorFloatInt, aa::AngleAxis) = new(convert(SO3,aa), v[1:3] ) # maybe replace with promote via dispatch
  SE3(v::VectorFloatInt, q::Quaternion) = new(convert(SO3,q), v[1:3] ) # maybe replace with promote via dispatch
  SE3(M::Array{Float64,2}) = new(SO3(M[1:3,1:3]), vec(M[1:3,4]) )
end

function normalize!(q::Quaternion, tol=0.00001)
    mag2 = sum(q.v.^2) + q.s^2
    if abs(mag2 - 1.0) > tol
        mag = sqrt(mag2)
        q.v = q.v ./ mag
        q.s = q.s / mag
    end
    nothing
end
function normalize(q::Quaternion, tol=0.00001)
  qq = deepcopy(q)
  normalize!(qq,tol)
  return qq
end

function normalize(v::Array{Float64,1})
  return v / norm(v)
end

function matrix(a::SE3)
  T = eye(4)
  T[1:3,1:3] = a.R.R
  T[1:3,4] = a.t
  return T
end

function q_conj!(q::Quaternion)
    normalize!(q)
    q.v = -q.v
    nothing
end
function q_conj(q::Quaternion)
    qq = deepcopy(q)
    q_conj!(qq)
    return qq
end

transpose(a::SO3) = SO3(a.R')
inverse(a::SO3) = transpose(a)
# backslash inverse -- con't remember if it is transpose for SE3 also TODO
inverse(a::SE3) = SE3( matrix(a) \ eye(4) )
A_invB(a::SE3, b::SE3) = SE3( ( matrix(b)' \ (matrix(a)') )' )

function *(a::SO3, b::SO3)
  return SO3(a.R*b.R)
end

function *(a::SO3, b::so3)
  return SO3(a*convert(SO3,b))
end
function *(a::so3, b::SO3)
  return SO3(convert(SO3,a)*b)
end

function *(a::SE3, b::SE3)
  return SE3(vec(a.R.R*b.t + a.t), a.R*b.R)
end

function *(q1::Quaternion, q2::Quaternion)
  ee  = [q1.s; q1.v] * [q2.s; q2.v]'
  w = ee[1,1] - ee[2,2] - ee[3,3] - ee[4,4]
  x = ee[1,2] + ee[2,1] + ee[3,4] - ee[4,3]
  y = ee[1,3] - ee[2,4] + ee[3,1] + ee[4,2]
  z = ee[1,4] + ee[2,3] - ee[3,2] + ee[4,1]
  if (w < 0.0)
    return Quaternion(-w, -[x; y; z])
  end
  return Quaternion(w, [x; y; z])
end

*(a::AngleAxis, b::AngleAxis) = convert(AngleAxis,convert(Quaternion,a)*convert(Quaternion,b))
#mangled type products, return first type or nearest Group type (doesn't return an Algebra)
*(a::so3,b::so3) = convert(SO3,a)*convert(SO3,b)
*(a::SO3, bq::Quaternion) = a*convert(SO3,bq)
*(aq::Quaternion, b::SO3) = aq*convert(Quaternion,b)
*(a::SO3, b::so3) = a*convert(SO3,b)
*(a::so3, b::SO3) = convert(SO3,a)*b
*(aq::Quaternion, b::so3) = aq*convert(Quaternion, convert(SO3,b) )
*(a::so3, bq::Quaternion) = convert(Quaternion, convert(SO3,a) )*bq
*(aq::Quaternion, baa::AngleAxis) = aq*convert(Quaternion, baa)
*(a::SO3, baa::AngleAxis) = a*convert(SO3,baa)
*(a::AngleAxis, bq::Quaternion) = convert(AngleAxis, convert(Quaternion,a)*bq)
*(a::AngleAxis, b::SO3) = convert(AngleAxis, convert(Quaternion,a)*b )
*(aa::AngleAxis, b::so3) = aa*convert(SO3,b)
*(a::so3, b::AngleAxis) = convert(AngleAxis, convert(SO3,a))*b

# comparison functions

compare(a::SO3, b::SO3; tol::Float64=1e-14) = norm((a*transpose(b)).R-eye(3)) < tol

function compare(a::SE3, b::SE3; tol::Float64=1e-14)
  norm(a.t-b.t) < tol ? nothing : return false
  return compare(a.R,b.R, tol=tol)
end
function compare(a::Quaternion, b::Quaternion; tol::Float64=1e-14)
  qiq = a*q_conj(b)
  return tol <= qiq.s <= 1.0+tol && norm(qiq.v) < tol
end
function compare(a::AngleAxis,b::AngleAxis; tol::Float64=1e-14)
  aTb = q_conj(convert(Quaternion,a))*b
  return compare(Quaternion(0),aTb, tol=tol)
end

# convert functions

function convert(::Type{Quaternion}, v::VectorFloatInt)
    return Quaternion(v[1],v[2:4])
end

function convert(::Type{Quaternion}, aa::AngleAxis)
    v = normalize(aa.ax)
    theta = aa.theta/2
    w = cos(theta)
    x = v[1] * sin(theta)
    y = v[2] * sin(theta)
    z = v[3] * sin(theta)
    return Quaternion(w, [x, y, z])
end

function convert(::Type{AngleAxis}, q::Quaternion)
    theta = acos(q.s) * 2.0
    if norm(q.v)==0
      return AngleAxis(theta, [1.,0,0])
    end
    return AngleAxis(theta, normalize(q.v))
end

function convert(::Type{SO3}, q::Quaternion)
    w = q.s;
    x = q.v[1];
    y = q.v[2];
    z = q.v[3];
    if w < 0.0
        w,x,y,z = -w, -x, -y, -z
    end
    nrm = sqrt(w*w+x*x+y*y+z*z)
    if (nrm < 0.999)
        println("q2C -- not a unit quaternion nrm = $(nrm)")
        R = eye(3)
    else
        nrm = 1.0/nrm
        w = w*nrm
        x = x*nrm
        y = y*nrm
        z = z*nrm
        x2 = x*x
        y2 = y*y
        z2 = z*z
        w2 = w*w
        xy = 2.0*x*y
        xz = 2.0*x*z
        yz = 2.0*y*z
        wx = 2.0*w*x
        wy = 2.0*w*y
        wz = 2.0*w*z
        #build matrix
        R = zeros(3,3)
        R[1,1] = w2+x2-y2-z2
        R[1,2] = xy-wz
        R[1,3] = xz+wy
        R[2,1] = xy+wz
        R[2,2] = w2-x2+y2-z2
        R[2,3] = yz-wx
        R[3,1] = xz-wy
        R[3,2] = yz+wx
        R[3,3] = w2-x2-y2+z2
    end
    return SO3(R)
end

function convert(::Type{Quaternion}, S::SO3)
  rot = S.R
  tr = rot[1,1] + rot[2,2] + rot[3,3];

  if (tr > 0.0)
    S = sqrt(tr+1.0) * 2.0; # S=4*qw
    qw = 0.25 * S;
    qx = (rot[3,2] - rot[2,3]) / S;
    qy = (rot[1,3] - rot[3,1]) / S;
    qz = (rot[2,1] - rot[1,2]) / S;
  else
    if ((rot[1,1] > rot[2,2]) && (rot[1,1] > rot[3,3]))
      S = sqrt(1.0 + rot[1,1] - rot[2,2] - rot[3,3]) * 2.0; # S=4*qx
      qw = (rot[3,2] - rot[2,3]) / S;
      qx = 0.25 * S;
      qy = (rot[1,2] + rot[2,1]) / S;
      qz = (rot[1,3] + rot[3,1]) / S;
    else
      if (rot[2,2] > rot[3,3])
        S = sqrt(1.0 + rot[2,2] - rot[1,1] - rot[3,3]) * 2.0 # S=4*qy
        qw = (rot[1,3] - rot[3,1]) / S;
        qx = (rot[1,2] + rot[2,1]) / S;
        qy = 0.25 * S;
        qz = (rot[2,3] + rot[3,2]) / S;
      else
        S = sqrt(1.0 + rot[3,3] - rot[1,1] - rot[2,2]) * 2.0 # S=4*qz
        qw = (rot[2,1] - rot[1,2]) / S
        qx = (rot[1,3] + rot[3,1]) / S
        qy = (rot[2,3] + rot[3,2]) / S
        qz = 0.25 * S;
      end
    end
  end

  if (qw >= 0)
    q = Quaternion(qw,[qx;qy;qz])
  else
    q = Quaternion(-qw,-[qx;qy;qz])
  end
  return q
end

function convert(::Type{SO3}, aa::AngleAxis)
  convert(SO3,convert(Quaternion, aa))
end

function convert(::Type{AngleAxis}, s::SO3)
  convert(AngleAxis, convert(Quaternion, s))
end

function convert(::Type{Euler}, q::Quaternion)
  #Function to convert quaternion to Euler angles, from Titterton and Weston
  #04, following quaternion convention of Titterton
  #UJ UAV Research 11/05/09

  a = q.s
  b = q.v[1]
  c = q.v[2]
  d = q.v[3]

  phi = -atan2(2.0*(c*d - b*a),1.0-2.0*(b^2+c^2)); # -atan2(2.0*(c*d - b*a), a*a - b*b - c*c + d*d) #numerical more stable
  theta = -asin(-2.0*(c*a - b*d));
  psi = -atan2(2.0*(b*c - d*a),1.0-2.0*(c^2+d^2)); # -atan2(2.0*(b*c - d*a), a*a + b*b - c*c - d*d)

  return  Euler(phi,theta,psi)
end

function convert(::Type{so3}, G::SO3)
  so3(logmap(G))
end

function convert(::Type{Euler}, R::SO3)
  convert(Euler, convert(Quaternion, R))
end

function convert(::Type{Quaternion}, E::Euler)
  # Using fixed frame rotation scheme, as used in MIT libbot
  q = zeros(4)

  halfroll = E.R/2.0;
  halfpitch = E.P/2.0;
  halfyaw = E.Y/2.0;

  sin_r2 = sin(halfroll);
  sin_p2 = sin(halfpitch);
  sin_y2 = sin(halfyaw);

  cos_r2 = cos(halfroll);
  cos_p2 = cos(halfpitch);
  cos_y2 = cos(halfyaw);

  q[1] = cos_r2 * cos_p2 * cos_y2 + sin_r2 * sin_p2 * sin_y2;
  q[2] = sin_r2 * cos_p2 * cos_y2 - cos_r2 * sin_p2 * sin_y2;
  q[3] = cos_r2 * sin_p2 * cos_y2 + sin_r2 * cos_p2 * sin_y2;
  q[4] = cos_r2 * cos_p2 * sin_y2 - sin_r2 * sin_p2 * cos_y2;

  q = q./norm(q);

  # Enforce positive scalar quaternions following conversion from Euler angles
  if (q[1]<0)
    q = -q;
  end
  return Quaternion(q[1],q[2:4])
end

function convert(::Type{SO3},E::Euler)
  return convert(SO3,convert(Quaternion, E))
end

function rotate!(q1::Quaternion, v1::Array{Float64,1})
    #v = (q1*Quaternion(0.0,v1)*q_conj(q1)).v
    R = convert(SO3, q1);
    v = R.R * v1
    for i = 1:3 v1[i] = v[i] end
    nothing
end
function rotate(q1::Quaternion, v1::Array{Float64,1})
    vv1 = deepcopy(v1)
    rotate!(q1, vv1)
    return vv1
end

function expmOwn(alg::so3)
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  I = eye(3)
  if (v_norm>1e-6)
    #1E-6 is chosen because double numerical LSB is around 1E-18 for nominal values [-pi..pi] and (1E-7)^2 is at 1E-14, but 1E-7 rad/s is 0.02 deg/h
    R = I + sin(v_norm)/v_norm*alg.S + (1.0-cos(v_norm))/(v_norm^2)*(alg.S^2)
  else
    R = I + alg.S + 0.5*(alg.S^2)
  end
  return R
end

function logmap(grp::SO3)
  R = grp.R
  tr = trace(R)
  if abs(tr+1.0) < 1e-10
    if abs(R[3,3]+1.0) > 1e-10
      ret = (pi / sqrt(2.0+2.0*R[3,3] )) * [R[1,3]; R[2,3]; (1.0+R[3,3])]
    else
      if abs(grp.R[1,1]+1.0) > 1e-10
        ret = (pi / sqrt(2.0+2.0*R[2,2])) * [R[1,2]; (1.0+R[2,2]); R[3,2]]
      else
        ret = (pi / sqrt(2.0+2.0*R[1,1])) * [(1.0+R[1,1]); R[2,1]; R[3,1]]
      end
    end
  else
    tr_3 = tr-3.0
    if (tr_3<-1e-7)
      theta = acos((tr-1.0)/2.0)
      magnitude = theta/(2.0*sin(theta))
    else
      magnitude = 0.5 - tr_3*tr_3/12.0
    end
    ret = magnitude*[ (R[3,2]-R[2,3]); (R[1,3]-R[3,1]); (R[2,1]-R[1,2])]
  end
  return ret
end

# Chikjain's book p??
function rightJacExmap(alg::so3)
  Gam = alg.S
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  I = eye(3)
  if (v_norm>1e-7)
    Jr = I + (v_norm-sin(v_norm))/(v_norm^3)*(Gam^2) - (1.0-cos(v_norm))/(v_norm^2+0.0)*Gam
  else
    Jr = I
  end
  return Jr
end

function rightJacExmapinv(alg::so3)
  Gam = alg.S
  v_norm = sqrt(alg.S[1,2]^2 + alg.S[1,3]^2 + alg.S[2,3]^2)
  I = eye(3)
  if (v_norm>1e-7)
    brace = (  1.0/(v_norm^2) - (1.0+cos(v_norm))/(2.0*v_norm*sin(v_norm))  )*Gam
    Jr_inv = I + 0.5*Gam + brace*Gam
  else
    Jr_inv = I
  end
  return Jr_inv
end

# function computeRatesFromRPY(Q, prevQ, dt)
#   if (dt<1):
#     # Assume oneside rectangular derivative, we could use a trapezoidal
#     # approach to obtain a more accurate estimate of the rotation rate
#     # TODO: Improve derivative scheme
#     S = (    np.dot(q2R(Q),np.transpose(q2R(prevQ)))  - np.identity(3)    ) / dt;
#     bw = np.matrix([[-S[2,1].item()],[-S[0,2].item()],[-S[1,0].item()]])
#   else
#     bw = zeros(3,1);
#   end
#   return bw
# end

# See Julia's implementation the matrix exponential function -- expm!()
# https://github.com/acroy/julia/blob/0bce8951f19e8f47d361e73b5b5bd6283926c01b/base/linalg/dense.jl
expmOwnT(M::Array{Float64,2}) = (eye(size(M,1)) + 0.5*M)*((eye(size(M,1)) - 0.5*M) \ eye(size(M,1)));
expmOwn1(M::Array{Float64,2}) = eye(size(M,1)) + M;
expmOwn2(M::Array{Float64,2}) = eye(size(M,1)) + M + 0.5*(M^2);
function expmOwn3(M::Array{Float64,2})
  M052 = 0.5*M^2;
  return eye(size(M,1)) + M + M052 + M052*M/3.0;
end
function expmOwn4(M::Array{Float64,2})
  M052 = 0.5*M^2;
  return eye(size(M,1)) + M + M052 + M052*M/3.0 + M052*M052/6.0
end

function convert(::Type{SO3}, alg::so3)
  v = vee(alg)
  nv = norm(v)
  if nv < 1e-3
    return SO3(expm(alg.S))
  else
    invnv = 1.0/nv
    return SO3(eye(3) + invnv*(sin(nv)*alg.S + invnv*(1.0-cos(nv))*(alg.S^2) ) )
  end
end

function wrapRad(th::Float64)
  if th >= pi
    th -= 2.0*pi
  end
  if th < -pi
    th += 2.0*pi
  end
  return th
end

R(th::Float64) = [[cos(th);-sin(th)]';[sin(th);cos(th)]'];

function vee!(rv::Vector{Float64}, alg::so3)
  rv[1] = -alg.S[2,3]
  rv[2] = alg.S[1,3]
  rv[3] = -alg.S[1,2]
  nothing
end

function vee(alg::so3)
   rv = zeros(3)
   vee!(rv, alg)#[-alg.S[2,3], alg.S[1,3], -alg.S[1,2]]
   return rv
end

# vectorize SE3 group to translation and AngleAxis numbers
# dim7 out: [xyz, theta, axyz]
function veeAngleAxis(G::SE3)
  aa = convert(AngleAxis, G.R)
  v = zeros(7)
  v[1:3] = G.t
  v[4] = aa.theta
  v[5:7] = aa.ax[:]
  return v
end

# vectorize SE3 group to translation and Quaternion numbers
# dim7 out: [xyz, cos(th/2), sin(th/2)xyz]
function veeQuaternion(G::SE3)
  q = convert(Quaternion, G.R)
  v = zeros(7)
  v[1:3] = G.t
  v[4] = q.s
  v[5:7] = q.v[:]
  return v
end

# vectorize SE3 group to translation and Quaternion numbers
# dim7 out: [xyz, cos(th/2), sin(th/2)xyz]
function veeEuler(G::SE3)
  E = convert(Euler, G.R)
  v = zeros(6)
  v[1:3] = G.t
  v[4],v[5],v[6] = E.R,E.P,E.Y
  return v
end

# TODO -- Change to type and overload the operators
# TODO -- uncomment
# function vee!(rv::Vector{Float64,1},T::SE2)
#   rv[1] = T.t[1]
#   rv[2] = T.t[2]
#   rv[3] = wrapRad(atan2(-T.R.R[1,2], T.R.R[1,1]))
#   nothing
# end
# function *(a::SE2, b::SE2)
#   return SE2(R.R*b.R, vec(a.R.R*b.t + a.t))
# end

function SE2(X::Array{Float64,1})
    T = eye(3)
    T[1:2,1:2] = R(X[3])
    T[1,3] = X[1]
    T[2,3] = X[2]
    return T
end

function se2vee!(retval::Array{Float64,1}, T::Array{Float64,2})
    retval[1] = T[1,3]
    retval[2] = T[2,3]
    retval[3] = wrapRad(atan2(-T[1,2], T[1,1]))
    nothing
end

function se2vee(T::Array{Float64,2})
    retval = zeros(3)
    se2vee!(retval, T)
    return retval
end

end #module


text:
module TestReaderWriter

using TSML
using TSMLextra
using Test

function test_csv()
    gcdims = (8761,2)
    ssum = 97564.0
    resdf=DataFrame()
    datapath=joinpath(dirname(pathof(TSMLextra)),"../data")
    outputfname = joinpath(tempdir(),"testdateval.csv")
    basefilename = "testdateval"
    fname = joinpath(datapath,basefilename*".csv")
    lcsv=DataReader(Dict(:filename=>fname))
    fit!(lcsv)
    dateval=transform!(lcsv)
    @test sum(size(dateval) .== gcdims ) == 2
    @test sum(dateval.Value) |> round == ssum
    csvname = replace(outputfname,"test"=>"out")
    wcsv = DataWriter(Dict(:filename=>csvname))
    fit!(wcsv)
    transform!(wcsv,dateval)
    pcsv = DataReader(Dict(:filename=>csvname))
    fit!(pcsv)
    resdf=transform!(pcsv)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf.Value) |> round == ssum
    rm(csvname,force=true)
end
@testset "Data Readers/Writers: csv" begin
    test_csv()
end

function test_hdf5()
    gcdims = (8761,2)
    ssum = 97564.0
    resdf=DataFrame()
    datapath=joinpath(dirname(pathof(TSMLextra)),"../data")
    outputfname = joinpath(tempdir(),"testdateval.csv")
    basefilename = "testdateval"
    fname = joinpath(datapath,basefilename*".csv")
    lcsv=DataReader(Dict(:filename=>fname))
    fit!(lcsv)
    dateval=transform!(lcsv)
    # check hdf5
    hdf5name = replace(outputfname,"csv"=>"h5")
    lhdf5 = DataWriter(Dict(:filename=>hdf5name))
    fit!(lhdf5)
    transform!(lhdf5,dateval)
    whdf5 = DataReader(Dict(:filename=>hdf5name))
    fit!(whdf5)
    resdf = transform!(whdf5)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf.Value) |> round == ssum
    rm(hdf5name,force=true)
end
@testset "Data Readers/Writers: hdf5" begin
    test_hdf5()
end

function test_feather()
    gcdims = (8761,2)
    ssum = 97564.0
    resdf=DataFrame()
    datapath=joinpath(dirname(pathof(TSMLextra)),"../data")
    outputfname = joinpath(tempdir(),"testdateval.csv")
    basefilename = "testdateval"
    fname = joinpath(datapath,basefilename*".csv")
    lcsv=DataReader(Dict(:filename=>fname))
    fit!(lcsv)
    dateval=transform!(lcsv)
    # check feather
    feathername = replace(outputfname,"csv"=>"feather")
    lfeather = DataWriter(Dict(:filename=>feathername))
    fit!(lfeather)
    transform!(lfeather,dateval)
    wfeather = DataReader(Dict(:filename=>feathername))
    fit!(wfeather)
    resdf = transform!(wfeather)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf.Value) |> round == ssum
    rm(feathername,force=true)
end
@testset "Data Readers/Writers: feather" begin
    test_feather()
end

function test_jld()
    gcdims = (8761,2)
    ssum = 97564.0
    resdf=DataFrame()
    datapath=joinpath(dirname(pathof(TSMLextra)),"../data")
    outputfname = joinpath(tempdir(),"testdateval.csv")
    basefilename = "testdateval"
    fname = joinpath(datapath,basefilename*".csv")
    lcsv=DataReader(Dict(:filename=>fname))
    fit!(lcsv)
    dateval=transform!(lcsv)
    # check jld
    jldname = replace(outputfname,"csv"=>"jld")
    ljld = DataWriter(Dict(:filename=>jldname))
    fit!(ljld)
    transform!(ljld,dateval)
    wjld = DataReader(Dict(:filename=>jldname))
    fit!(wjld)
    resdf = transform!(wjld)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf.Value) |> round == ssum
    rm(jldname,force=true)
end
@testset "Data Readers/Writers: jld" begin
    test_jld()
end

end

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MiniFB"
version = v"0.2.0"

# Collection of sources required to complete build
sources = [
	    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/11.3/MacOSX10.13.sdk.tar.xz",
                  "1d2984acab2900c73d076fbd40750035359ee1abe1a6c61eafcd218f68923a5a"),
	    GitSource("https://github.com/emoon/minifb.git", "5066489cd81b23b0c79952f7d6f464b20c54867c") #master as of 24Feb2022
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/minifb/
sed -i -e 's/add_library(minifb STATIC/add_library(minifb SHARED/' \
    -e 's/-Wall/-I$ENV{includedir} -Wall/' \
    -e 's/Opengl32/opengl32/' \
    CMakeLists.txt
sed -i -e 's?<gl/gl.h>?<GL/gl.h>?' src/gl/MiniFB_GL.c
mkdir build && cd build
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    # Workaround for https://github.com/emoon/minifb/issues/88
    export MACOSX_DEPLOYMENT_TARGET=10.13
    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd
fi
if [[ "$target}" == *-musl* || "${target}" == *-freebsd* ]]; then
    export CFLAGS="-I${includedir}"
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DMINIFB_BUILD_EXAMPLES=OFF \
    -DMINIFB_AVOID_CPP_HEADERS=ON \
    ..
make -j${nproc}
install -Dm 755 "libminifb.${dlext}" "${libdir}/libminifb.${dlext}"
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=p->arch(p)=="armv6l"||BinaryBuilder.os(p)=="freebsd")

# The products that we will ensure are always built
products = [
    LibraryProduct("libminifb", :libminifb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"); platforms=filter(p->Sys.islinux(p) || Sys.isfreebsd(p), platforms))
    Dependency(PackageSpec(name="xkbcommon_jll", uuid="d8fb68d0-12a3-5cfd-a85a-d49703b185fd"); platforms=filter(Sys.islinux, platforms))
    Dependency(PackageSpec(name="Libglvnd_jll", uuid="7e76a0d4-f3c7-5321-8279-8d96eeed0f29"); platforms=filter(p->Sys.islinux(p) || Sys.isfreebsd(p), platforms))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"); platforms=filter(p->Sys.islinux(p) || Sys.isfreebsd(p), platforms))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0", julia_compat="1.6")


text:
function toh03()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    k = [(1:20)...]
    z = 10/11 .- exp2.(-k)
    p = poly(z)
    z = [z ones(20)]

    p, PolyZeros(z)
end


text:
using Parameters, Distributions, StatsPlots

@with_kw struct Sensor
    noise::Distribution = Normal(0, 0.2)
    # κ1, κ2 = cupula_dynamics(Δt)
end

@with_kw struct NigmatullinaTrial
    D::String # direction
    onset::Real = 0.0
end

@with_kw struct HeadTurn
    A::Real = 20 # amplitude
    D::String = "right" # direction
    onset::Real = 1
    duration::Real = 1
end

@with_kw struct PlannedMovement
    A::Real = 20 # amplitude
    D::String = "right" # direction
    onset::Real = 1
    duration::Real = 1
end

@with_kw struct ExternalDisturbance
    A::Real = 10 # amplitude
    D::String = "right" # direction
    onset::Real = 1
    duration::Real = 2
end

@with_kw struct Translation
    A::Real = 20 # amplitude
    D::String = "right" # direction
    onset::Real = 1
    duration::Real = 1
end

function cupula_dynamics(Δt::Real; τ::Real = 4)
    # τ = 4 same as Laurens & Angelaki (2017)
    κ₁ = τ / (τ + Δt)
    κ₂ = Δt / (τ + Δt)
    return κ₁, κ₂
end

function velocity(D::Number, A::Number, f::Number, t::Number, start::Number)
    D * A * 1/(2*π*f) * (1-cos.(2*π*f*(t-start)))
end

function acceleration(D::Number, A::Number, f::Number, t::Number, start::Number;
    noisy::Bool = false)
    noise = noisy == 0 ? 0 : rand(Normal(0, 2.0))
    D * A * sin.(2*π*f*(t-start)) + noise
end

function displacement(A::Number; T::Number = 1)
    A * T^2 /(2π)
end

function amplitude(θ::Number; T::Number = 1)
    (2π/T^2) * θ
end

function peak_velocity(A::Number; T::Number = 1)
    (T/π) * A
end

function simulate(m::NigmatullinaTrial;
    onset::Real = 0, trial_duration::Real = 15)
    #= An alternative “limit-finding” approach was adopted with each trial applying
    progressively increasing acceleration in the yaw plane, increasing by 0.5 deg/s^2
    every three seconds, reaching 82 deg/sec

    Cutfield, N. J., Cousins, S., Seemungal, B. M., Gresty, M. A., &
    Bronstein, A. M. (2011). Vestibular perceptual thresholds to angular
    rotation in acute unilateral vestibular paresis and with galvanic
    stimulation. Annals of the New York Academy of Sciences, 1233(1), 256–262.
    =#

    @assert onset >= 0
    @assert trial_duration > 0

    m.D ∈ ["left", "right"] || error("Direction not specified.")
    D = m.D == "left" ? -1 : 1
    noise = m.sensor
    Δt = m.Δt
    total_duration = onset + trial_duration

    timesteps = range(0, stop = total_duration, step = Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(Δt)

    # either use undef intialiser, or intialise with zeros
    # ω = Array{Float64}(undef, nsteps) # velocity
    # θ = Array{Float64}(undef, nsteps) # position
    # c = Array{Float64}(undef, nsteps) # cupula dynamics
    # y = Array{Float64}(undef, nsteps) # sensory observations
    α = zeros(T)
    ω = zeros(T)
    θ = zeros(T)
    c = zeros(T)
    y = zeros(T)

    increment = D * 0.3
    tstep = 3

    # start until motion onset
    j = 1
    t = round(timesteps[j]; digits = 3)

    while t < onset
        y[j] = ω[j] - c[j] + rand(noise)
        j += 1
        t = round(timesteps[j]; digits = 3)
    end

    k = 1 # counts the number of times that α is incremented
    # start moving at time j
    α[j] = α[j] + increment

    for i ∈ (j+1):T
        t = round(timesteps[i]; digits = 3)

        if t >= (onset + k * tstep)
            α[i] = α[i-1] + increment
            k += 1
        else
            α[i] = α[i-1]
        end
        ω[i] = ω[i-1] + Δt * α[i-1]
        θ[i] = θ[i-1] + Δt * ω[i-1] + 1/2 * Δt^2 * α[i-1]
        c[i] = κ1 * c[i-1] + κ2 * ω[i]
        y[i] = ω[i] - c[i] + rand(noise)
    end

    # return a named tuple
    out = (timesteps = collect(timesteps),
            y = y, α = α, ω = ω, θ = θ, c = c,
            Δt = Δt, onset = onset)
    return out
end

function simulate(mᵤ::HeadTurn,
    mₑ::ExternalDisturbance,
    sensor::Sensor;
    Δt = 0.01,
    duration::Real = 5,
    stochastic::Bool = false)

    mᵤ = @isdefined(mᵤ) ? mᵤ : HeadTurn(A = 0)
    mₑ = @isdefined(mₑ) ? mₑ : ExternalDisturbance(A = 0)

    onsetᵤ = mᵤ.onset
    onsetₑ = mₑ.onset
    @assert onsetᵤ >= 0
    @assert onsetₑ >= 0
    @assert duration > 0

    # onset = m.onset
    duration1 = mᵤ.duration
    duration2 = mₑ.duration
    total_duration = duration

    onsetᵤ + duration1 <= total_duration || error("Head turn extends beyond simulation event.")
    onsetₑ + duration2 <= total_duration || error("Disturbance extends beyond simulation event.")

    mᵤ.D ∈ ["left", "right"] || error("Direction not specified.")
    mₑ.D ∈ ["left", "right"] || error("Direction not specified.")

    D1 = mᵤ.D == "left" ? -1 : 1
    D2 = mₑ.D == "left" ? -1 : 1

    # motion amplitude is noisy if stochastic simulation
    A1 = stochastic == 0 ? mᵤ.A : rand(Normal(mᵤ.A, 2.0))
    A2 = stochastic == 0 ? mₑ.A : rand(Normal(mₑ.A, 2.0))

    noise = sensor.noise
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    f2 = 1/(duration2)

    timesteps = range(0, stop = total_duration, step = Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(Δt)

    α1 = zeros(T)
    α2 = zeros(T)
    α = zeros(T)
    ω1 = zeros(T)
    ω2 = zeros(T)
    ω = zeros(T) # ω = Array{Float64}(undef, T) # velocity
    θ = zeros(T)
    c = zeros(T)
    y = zeros(T)

    for i ∈ Iterators.drop(1:T, 1)

        t = round(timesteps[i]; digits = 3)

        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0
        α2[i] = (t > onsetₑ) & (t < onsetₑ + duration2) ? acceleration(D2, A2, f2, t, onsetₑ, noisy = stochastic) : 0
        α[i] = α1[i] + α2[i]
        ω1[i] = ω1[i-1] + Δt * α1[i]
        ω2[i] = ω2[i-1] + Δt * α2[i]

        ω[i] =  ω1[i] + ω2[i]
        θ[i] = θ[i-1] + Δt * ω[i] + 1/2 * Δt^2 * α[i]
        c[i] = κ1 * c[i-1] + κ2 * ω[i]
        y[i] = ω[i] - c[i] + rand(noise)
    end

    out = (timesteps = collect(timesteps),
            y = y,
            αᵤ = α1, αₑ = α2,
            α = α,
            ωᵤ = ω1, ωₑ = ω2,
            ω = ω,
            θ = θ, c = c,
            Δt = Δt, onsetᵤ = onsetᵤ, onsetₑ = onsetₑ,
            sensor = noise, stochastic = stochastic)
    return out

end

simulate(mᵤ::HeadTurn,
    sensor::Sensor;
    Δt = 0.01,
    duration::Real = 5,
    stochastic::Bool = false) = simulate(mᵤ, ExternalDisturbance(A = 0, duration = 0.1), sensor,
        Δt = Δt, duration = duration, stochastic = stochastic)

simulate(mₑ::ExternalDisturbance,
    sensor::Sensor;
    Δt = 0.01,
    duration::Real = 5,
    stochastic::Bool = false) = simulate(HeadTurn(A = 0, duration = 0.1), mₑ, sensor,
        Δt = Δt, duration = duration, stochastic = stochastic)


text:
using Test, Example
@test domath(2) == 4


text:
## This function prints some information on what should be the contents the user has to give in the different sections of the package.
# To select from which section you want the information from, just type the keywords:
#        Model, Simulation/s, Pseudo-Data/Pseudodata, MLE, Inference/Stan

function infoAll(woo=[])

    if typeof(woo) != String
        println("Sorry, but the only allowed input to this function is a string :(");

        print(string("

Please, if you want information about one of the package structures, type:
    1) ","\"","model","\""," for the model generation section
    2) ","\"","simulation","\"",", ","\"","simulation","\""," or ","\"","simul","\""," for the model simulation section
    3) ","\"","pseudo-data","\"",", ","\"","pseudodata","\"",", ","\"","pseudo data","\""," for the model pseudo-data generation section
    4) ","\"","mle","\"",", ","\"","likelihood","\""," for the maximum likelihood estimation section
    5) ","\"","inference","\"",", ","\"","stan","\"",", ","\"","stan inference","\"",", ","\"","staninference","\""," for the Bayesian parameter inference section
    5b) ","\"","turing","\"",", ","\"","turing inference","\"",", ","\"","turinginference","\""," for the Bayesian parameter inference section using Turing.jl
    6) ","\"","oedms","\"",", ","\"","model selection","\"",", ","\"","modelselection","\"",", ","\"","oed model selection","\""," for the optimal experimental design for model selection section
    7) ","\"","oedmc","\"",", ","\"","model calibration","\"",", ","\"","modelcalibration","\"",", ","\"","oed model calibration","\""," for the optimal experimental design for model calibration section

        "))

        return
    end

    woohoo = lowercase(woo);

    if woohoo == "model"
        print(string(raw"

                CALL defModStruct()

model_def[","\"","NameF","\"","] = [];
        # String containing the name of the model. Scripts and results will be stored using this name

model_def[","\"","nStat","\"","] = [];
        # Integer indicating the total number of steps of the model.

model_def[","\"","nPar","\"","] = [];
        # Integer indicating the total number of parameters of the model.

model_def[","\"","nInp","\"","] = [];
        # Integer indicating the total number of stimuli (inducers) of the model. If the model has no inputs, set it to 0.

model_def[","\"","stName","\"","] = [];
        # Vector of strings indicating the name of all the states of the model (without a d letter in front).

model_def[","\"","parName","\"","] = [];
        # Vector of strings indicating the name of all the parameters of the model.

model_def[","\"","inpName","\"","] = [];
        # Vector of strings indicating the name of all the stimuli (inducers) of the model. If the model has no inputs
        # just give an empty vector.

model_def[","\"","eqns","\"","] = [];
        # Vector of strings containing all the equations for the model (left and right-hand sides).
        # If an equation represents a state, the left-hand side has to be one of the strings contained in
        # stName but with a d in front (example: Prot -> dProt = ...).
        # Equations that are not states of ODEs are also allowed. Same as Julia expressions (println, for, if, etc.)
        # If you want to include a condition (if) for a state variable, each one of the if statement has to be
        # written in a separate string. Be careful with this, since then your Stan model will not work. To make it
        # work you have to write it in Stan language, but then the Julia code will not work. We recommend that if so,
        # first generate the stan code and then apply all the necessary modifications there. However, if you do not care
        # for the Julia code and just want the Stan code go on. Just know that in this case, if there is a condition in
        # one of the states, in stan you need to type the full if statement in one same string.
        # Note that models without external inputs are also supported except for the optimal experimental design sections.

model_def[","\"","Y0eqs","\"","] = [];
        # Vector of strings containing the steady-state equations of the model if desired (if not, just leave
        # it as an empty vector). These equations will be used to compute y0 assuming steady-state reached
        # before the experiment.
        # All the states must appear (identified in the left-hand side, but this time without the d in front),
        # however other equations are allowed.
        # If some  element of the equation requires an experimental value for the calculation,
        # please add exp at the beginning of the state (example: Cmrna -> expCmrna).
        # Please, do not use the name alp for anything, since this is reserved.

model_def[","\"","Y0Sim","\"","] = [];
        # If analytical solution of the model at steady state is not accurate enough and you want to add
        # an Over-Night simulation before the actual experiment simulation.
        # Allowed vales are true, false, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\"",".
        # Default value is false.
        # Time-scale for the simulation is assumed in minutes (1440 min). If this wants to be changed, what should be
        # introduced here is a number for the time conversion (e.g. 1/60 if to convert to days or 60 if to convert to seconds)

model_def[","\"","tols","\"","] = [];
        # Vector of 2 floats containing the relative and absolute tolerances for the solver (in this order).
        # If left empty, 1e-6 will be assumed for both.

model_def[","\"","solver","\"","] = [];
        # IVP solver to solve the ODEs. If nothing specified, the default will be Tsit5().
        # For more info check https://diffeq.sciml.ai/v2.0/tutorials/ode_example.html

"));

    elseif woohoo == "simulation" || woohoo == "simulations" || woohoo == "simul"
        print(string(raw"

                CALL defSimulStruct()
        MAIN STRUCTURE

simul_def[","\"","Nexp","\"","] = [];
        # Integer indicating the number of experiments to be simulated

simul_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

simul_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

simul_def[","\"","y0","\"","] = [];
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

simul_def[","\"","preInd","\"","] = [];
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

simul_def[","\"","uInd","\"","] = [];
        # Array containing the values for the stimuli at each step for each experiment. In each sub-array, columns
        # indicate a step and rows indicate an inducer (if multiple ones are considered).

simul_def[","\"","theta","\"","] = [];
        # Vector/Matrix with the parameter samples or directory and file location of CSV file with them.

simul_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

simul_def[","\"","plot","\"","] = [];
        # Boolean or yes/no string to save the resulting simulations in the results directory (false will be
        # considered as default).

simul_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated scripts and result files so they are not overwritten.
        # If empty, nothing will be added.

"))

print(string(raw"

            CALL defSimulStructFiles()

        STRUCTURE IF CSV FILES ARE USED

simul_def[","\"","ObservablesFile","\"","] = [];
        # Vector of strings containing the name of the files that have the information about the sampling
        # times and y0 values.

simul_def[","\"","EventInputsFile","\"","] = [];
        # Vector of strings containing the name of all the files that have the information about the stimuli
        # and events of the experiment.

simul_def[","\"","theta","\"","] = [];
        # Vector/Matrix with the parameter samples or directory and file location of CSV file with them.

simul_def[","\"","MainDir","\"","] = [];
        # Main directory path for the files in ObervablesFile and EventInputsFile.

simul_def[","\"","plot","\"","] = [];
        # Boolean or yes/no string to save the resulting simulations in the results directory (false will be
        # considered as default).

simul_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated scripts and result files so they are not overwritten.
        # If empty, nothing will be added.

"))

    elseif woohoo == "pseudo-data" || woohoo == "pseudodata" || woohoo == "pseudo data"
        print(string(raw"

            CALL defPseudoDatStruct()

        MAIN STRUCTURE

pseudo_def[","\"","Nexp","\"","] = [];
        # Integer indicating the number of experiments to be simulated

pseudo_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

pseudo_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

pseudo_def[","\"","y0","\"","] = [];
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

pseudo_def[","\"","preInd","\"","] = [];
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

pseudo_def[","\"","uInd","\"","] = [];
        # Array containing the values for the stimuli at each step for each experiment. In each sub-array, columns
        # indicate a step and rows indicate an inducer (if multiple ones are considered).

pseudo_def[","\"","theta","\"","] = [];
        # Vector/Matrix with the parameter samples or directory and file location of CSV file with them.

pseudo_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

pseudo_def[","\"","plot","\"","] = [];
        # Boolean or yes/no string to save the resulting simulations in the results directory (false will be
        # considered as default).

pseudo_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated scripts and result files so they are not overwritten.
        # If empty, nothing will be added.

pseudo_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is either a vector of strings, a vector of integers
        # indicating which entries from model_def[","\"","stName","\"","] are observables. If a vector of
        # strings is given, these could also be an expression combining states (Only +,-,*,/ and ^ will be
        # considered).

pseudo_def[","\"","NoiseType","\"","] = [];
        # String indicating if the desired noise is homoscedastic or heteroscedastic. The allowed strings are
        # homo, homoscedastic or hetero, heteroscedastic (case independent). If empty, default will be heteroscedastic.
        # Note that for high homoscedastic noise data could become negative (which to me does not make much sense), for
        # now this will be set to a value close to 0.

pseudo_def[","\"","Noise","\"","] = [];
        # Percentage of heteroscedastic noise (introduced as value from 0 to 1). If empty 10% heteroscedastic
        # noise will be assumed. In the homoscedastic case, if empty a standard deviation of 1 will be assumed.
        # This has to be a vector of noise values for each observable.

"))

        print(string(raw"

            CALL defPseudoDatStructFiles()

        STRUCTURE IF CSV FILES ARE USED

pseudo_def[","\"","ObservablesFile","\"","] = [];
        # Vector of strings containing the name of the files that have the information about the sampling
        # times and y0 values.

pseudo_def[","\"","EventInputsFile","\"","] = [];
        # Vector of strings containing the name of all the files that have the information about the stimuli
        # and events of the experiment.

pseudo_def[","\"","theta","\"","] = [];
        # Vector/Matrix with the parameter samples or directory and file location of CSV file with them.

pseudo_def[","\"","MainDir","\"","] = [];
        # Main directory path for the files in ObervablesFile and EventInputsFile.

pseudo_def[","\"","plot","\"","] = [];
        # Boolean or yes/no string to save the resulting simulations in the results directory (false will be
        # considered as default).

pseudo_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated scripts and result files so they are not overwritten.
        # If empty, nothing will be added.

pseudo_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is either a vector of strings, a vector of integers
        # indicating which entries from model_def[","\"","stName","\"","] are observables. If a vector of
        # strings is given, these could also be an expression combining states (Only +,-,*,/ and ^ will be
        # considered).

pseudo_def[","\"","NoiseType","\"","] = [];
        # String indicating if the desired noise is homoscedastic or heteroscedastic. The allowed strings are
        # homo, homoscedastic or hetero, heteroscedastic (case independent). If empty, default will be heteroscedastic.
        # Note that for high homoscedastic noise data could become negative (which to me does not make much sense), for
        # now this will be set to a value close to 0.

pseudo_def[","\"","Noise","\"","] = [];
        # Percentage of heteroscedastic noise (introduced as value from 0 to 1). If empty 10% heteroscedastic
        # noise will be assumed. In the homoscedastic case, if empty a standard deviation of 1 will be assumed.
        # This has to be a vector of noise values for each observable.

"))

    elseif woohoo == "mle" || woohoo == "likelihood"
        print(string(raw"

                CALL defMLEStruct()

mle_def[","\"","Nexp","\"","] = [];
        # Integer indicating the number of experiments to be simulated

mle_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

mle_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

mle_def[","\"","y0","\"","] = [];
        # Array (single simulation) of Y0s for the simulations for the
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

mle_def[","\"","preInd","\"","] = [];
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

mle_def[","\"","uInd","\"","] = [];
        # Array containing the values for the stimuli at each step for each experiment. In each sub-array, columns
        # indicate a step and rows indicate an inducer (if multiple ones are considered).

mle_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

mle_def[","\"","plot","\"","] = [];
        # Boolean or yes/no string to save the resulting simulations in the results directory (false will be
        # considered as default).

mle_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated scripts and result files so they are not overwritten.
        # If empty, nothing will be added.

mle_def[","\"","thetaMAX","\"","] = [];
        # Vector containing the maximum bounds for theta (no files can be introduced)

mle_def[","\"","thetaMIN","\"","] = [];
        # Vector containing the minimum bounds for theta (no files can be introduced)

mle_def[","\"","runs","\"","] = [];
        # Integer indicating how many runs of Optimisation will be done. You will get as many theta vectors
        # as runs selected.

mle_def[","\"","parallel","\"","] = [];
        # Boolean or yes/no string indicating if the different runs want to be done in parallel (true) or
        # series (false). Default is false.

    # For the two following fields, you can introduce a string pointing to the observable files (same strings in)
    #     both fields having the same structure as the ones generated in the PseudoData section. If multiple theta are
    #     considered in the file, then the covariance matrix will be taken.
mle_def[","\"","DataMean","\"","] = [];
        # Array containing the vector of means for each experiment.

mle_def[","\"","DataError","\"","] = [];
        # Array containing the vector or matrices (covariance included) of errors for the data for each experiment.

                                        IMPORTANT!!!!

        # Whilst each entry (experiment) of DataMean can be a matrix where each column is an observable of the
        # system, for DataError this is not the case. Each entry of the array (experiment) will have as many
        # entries as observables, where it would be a vector of errors or a matrix. This is done this way to
        # generalise the presence of both options.

mle_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is either a vector of strings, a vector of integers
        # indicating which entries from model_def[","\"","stName","\"","] are observables. If a vector of
        # strings is given, these could also be an expression combining states (Only +,-,*,/ and ^ will be
        # considered).

mle_def[","\"","OPTsolver","\"","] = [];
        # For now we only use the package BlackBoxOptim, so any of their options can be used. The default
        # is adaptive_de_rand_1_bin_radiuslimited. Please, introduce it as a string.

mle_def[","\"","MaxTime","\"","] = [];
        # Integer indicating the maximum number of time allowed for the optimisation as a stop criterion.
        # If this is selected MaxFuncEvals has to be empty (only 1 stop criterion). If both are empty,
        # the default will be to do 1000 function evaluations

mle_def[","\"","MaxFuncEvals","\"","] = [];
        # Integer indicating the maximum number of function evaluations allowed for the optimisation as a
        # stop criterion. If this is selected MaxTime has to be empty (only 1 stop criterion). If both are
        # empty, the default will be to do 1000 function evaluations

"))

    elseif woohoo == "inference" || woohoo == "stan" || woohoo == "stan inference" || woohoo == "staninference"
        print(string(raw"

            CALL defBayInfStruct()

        MAIN STRUCTURE

bayinf_def[","\"","Priors","\"","] = []; # Five options for this:
        # 1) A 2*N array containing the bounds for the parameters. Order of parameters will be assumed as
            # the one introduced in model_def[","\"","parName","\"","]. As a prior a truncated Normal covering 2
            # standard deviations in the bounds given will be generated as prior.
        # 2) Path to a CSV file containing samples for the parameters. Fitting of the samples to
            # different type of distributions will be done to generate the priors. Order of parameters
            # will be assumed as the one introduced in model_def[","\"","parName","\"","]. You can also
            # introduce a 2D array of floats with the samples (Array{Float64,2}).
        # 3) Dictionary with fields pars, transpars and pridis defining the
            # parameters, subsequent desired transformations and prior distributions.
            # In the field pars, parameters
            # have to be defined with the same names and order as in model_def[","\"","parName","\"","],
            # otherwise the script will not proceed.
        # 4) An empty array. If this is the case, the stan model will be generated with nothing in
            # the parameters and transformed parameters section. The path to the stan file will
            # be given so the user can fill these sections.
        # 5) Path to a Stan model file if you already have one. For example, if first you introduce an
            # empty array, then you can fill the parameters section and run from that script instead
            # of having to copy the things in here.

bayinf_def[","\"","Data","\"","] = []; # Two options:
        # 1) A dictionary containing 3 fields. Observables, for the path to the files containing the
            # experimental data. Inputs, for the path to the files containing the input profiles for the
            # experiments. Obs, a string vector containing the observables of the experiments.
            # The format of the files has to be the same one as the ones generated in the
            # pseudo-data section. For each of the 2 entries, more than one file can be given if a
            # multi-experimental inference wants to be done. Obs file also needs to be given.
            # A y0 in a field with the same name has to be given for each experiment. The field is
            # compulsory, even if not used.
            # WARNING: For now, this option does not include extraction of covariance matrix for data.
        # 2) A dictionary containing the same structure as the simul_def plus the fields DataMean and
            # DataError containing the experimental data. You can call the function defBayInfDataStruct()
            # to obtain the empty structure of the dictionary.

bayinf_def[","\"","StanSettings","\"","] = []; # Two options:
        # 1) A dictionary with the basic fields from a stan run. The structure of the dictionary can be
            # extracted calling the function defBasicStanSettingsStruct().
        # 2) An empty array. If this is the case, no inference will be done after calling the main
            # function. Instead, you will be given a StanModel file and data structure and you will
            # have to use this to run a call of Stan by yourself. An example will be provided.

bayinf_def[","\"","flag","\"","] = [];
            # String to attach a unique flag to the generated files so it is not overwritten.
            # If empty, nothing will be added.

bayinf_def[","\"","plot","\"","] = [];
            # true, flase, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\""," or []
            # indicating if plots with the results will be generated. Default is false.

bayinf_def[","\"","runInf","\"","] = [];
            # true, flase, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\""," or []
            # indicating if inference wants to be performed or only stan structure is given.
            # Default is true (if all the necessary elements are present)

bayinf_def[","\"","MultiNormFit","\"","] = [];
            # true, false, ","\"","yes","\"",", ","\"","no","\"",", ","\"","Yes","\"",", ","\"","No","\""," or []
            # indicating that if samples are given to fit, the prior will be fitted as a MultiNormal distribution
            # or not. the default is false. If some parameter is better fit with a Log-Normal distribution, this
            # will be used with a distribution reparameterisation to begin to be included in the multinormal.
            # If for some parameter a Uniform distribution is better, this parameter will be excluded from
            # the multinormal and defined as a separate parameter. This can be used as an example of how to set
            # your Stan Model to use multi_normal priors.

"))

print(string(raw"

            CALL defBayInfDataStruct()

        DATA DICTIONARY STRUCTURE

data_def[","\"","Nexp","\"","] = [];
        # Integer indicating the number of experiments to be simulated

data_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

data_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

data_def[","\"","y0","\"","] = [];
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

data_def[","\"","preInd","\"","] = [];
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

data_def[","\"","uInd","\"","] = [];
        # Array containing the values for the stimuli at each step for each experiment. In each sub-array, columns
        # indicate a step and rows indicate an inducer (if multiple ones are considered).

data_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

data_def[","\"","DataMean","\"","] = [];
        # Array containing the vector of means for each experiment.

data_def[","\"","DataError","\"","] = [];
        # Array containing the vector or matrices (covariance included) of errors for the data for each experiment.

                                        IMPORTANT!!!!

        # Whilst each entry (experiment) of DataMean can be a matrix where each column is an observable of the
        # system, for DataError this is not the case. Each entry of the array (experiment) will have as many
        # entries as observables, where it would be a vector of errors or a matrix. This is done this way to
        # generalise the presence of both options.

data_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is either a vector of strings, a vector of integers
        # indicating which entries from model_def[","\"","stName","\"","] are observables. If a vector of
        # strings is given, these could also be an expression combining states (Only +,-,*,/ and ^ will be
        # considered).

"))

print(string(raw"

            CALL defBayInfDataFromFilesStruct()

        DATA DICTIONARY STRUCTURE IF FILES ARE GIVEN

data_def[","\"","Observables","\"","] = [];
        # Vector of strings containing the name of the files that have the information about the sampling
        # times and y0 values.

data_def[","\"","Inputs","\"","] = [];
        # Vector of strings containing the name of all the files that have the information about the stimuli
        # and events of the experiment.

data_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is either a vector of strings, a vector of integers
        # indicating which entries from model_def[","\"","stName","\"","] are observables. If a vector of
        # strings is given, these could also be an expression combining states (Only +,-,*,/ and ^ will be
        # considered).

data_def[","\"","y0","\"","] = [];
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

"))

print(string(raw"

            CALL defBasicStanSettingsStruct()

        STAN SETTINGS DICTIONARY STRUCTURE

stan_def[","\"","cmdstan_home","\"","] = [];
        # String containing the full path to the installation directory for cmdstan

stan_def[","\"","nchains","\"","] = [];
        # Number of chains for the inference (integer)

stan_def[","\"","nsamples","\"","] = [];
        # Number of post-warmup samples for each chain (integer)

stan_def[","\"","nwarmup","\"","] = [];
        # Number of warm-up samples for each chain (integer)

stan_def[","\"","printsummary","\"","] = [];
        # Print sumary of inference at the end. This can be either true or false. Default will be true.

stan_def[","\"","init","\"","] = [];
        # Initial point for the parameters in each chain. See output of MLE results. This field can be empty.
        # Please introduce the parameter values in their true parameter range.

stan_def[","\"","maxdepth","\"","] = [];
        #  Maximum tree-depth. This field can be empty. Check Stan documentation for more information.

stan_def[","\"","adaptdelta","\"","] = [];
        # Delta value between 0 and 1. This field can be empty.  Check Stan documentation for more information.

stan_def[","\"","jitter","\"","] = [];
        # Jitter value between 0 and 1. This field can be empty.  Check Stan documentation for more information.

"));

elseif woohoo == "turing" || woohoo == "turinginference"
    print(string(raw"

                CALL defTurInfStruct()

            MAIN STRUCTURE

    turinf_def[","\"","Priors","\"","] = []; # Three options for this:
            # 1) A 2*N array containing the bounds for the parameters. Order of parameters will be assumed as
                # the one introduced in model_def[","\"","parName","\"","]. As a prior a truncated Normal covering 2
                # standard deviations in the bounds given will be generated as prior.
            # 2) Path to a CSV file containing samples for the parameters. Fitting of the samples to
                # different type of distributions will be done to generate the priors. Order of parameters
                # will be assumed as the one introduced in model_def[","\"","parName","\"","]. You can also
                # introduce a 2D array of floats with the samples (Array{Float64,2}).
            # 3) Array containing any distribution from the Julia package Distributions.jl with their
                # corresponding parameters. Order of parameters will be assumed as
                # the one introduced in model_def[","\"","parName","\"","].

    turinf_def[","\"","Data","\"","] = []; # Two options:
            # 1) A dictionary containing 3 fields. Observables, for the path to the files containing the
                # experimental data. Inputs, for the path to the files containing the input profiles for the
                # experiments. Obs, a string vector containing the observables of the experiments.
                # The format of the files has to be the same one as the ones generated in the
                # pseudo-data section. For each of the 2 entries, more than one file can be given if a
                # multi-experimental inference wants to be done. Obs file also needs to be given.
                # A y0 in a field with the same name has to be given for each experiment. The field is
                # compulsory, even if not used.
                # WARNING: For now, this option does not include extraction of covariance matrix for data.
            # 2) A dictionary containing the same structure as the simul_def plus the fields DataMean and
                # DataError containing the experimental data. You can call the function defBayInfDataStruct()
                # to obtain the empty structure of the dictionary.

    turinf_def[","\"","flag","\"","] = [];
                # String to attach a unique flag to the generated files so it is not overwritten.
                # If empty, nothing will be added.

    turinf_def[","\"","plot","\"","] = [];
                # true, flase, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\""," or []
                # indicating if plots with the results will be generated. Default is false.

    turinf_def[","\"","MultiNormFit","\"","] = [];
                # true, false, ","\"","yes","\"",", ","\"","no","\"",", ","\"","Yes","\"",", ","\"","No","\""," or []
                # indicating that if samples are given to fit, the prior will be fitted as a set of Normal distributions
                # or instead the best type of distributions that can be fit with the samples.

    turinf_def[","\"","Trunc","\"","] = [];
                # true, false, ","\"","yes","\"",", ","\"","no","\"",", ","\"","Yes","\"",", ","\"","No","\"",", [] or
                # 2*N array indicating upper and lower bounds for the priors. If true then the truncations will be
                # the smaller sample minus 10% its value (never allowing this go bellow 0) and the bigest sample plus
                # 10% its value. If false no truncation will be applaied to any prior. In all cases, order of parameters
                # will be assumed as the one introduced in model_def[","\"","parName","\"","]. Default will be true.

    "))

    elseif woohoo == "oedms" || woohoo == "model selection" || woohoo == "modelselection" || replace("woohoo", " "=>"") == "oedmodelselection"
print(string(raw"

                CALL defODEModelSelectStruct()

oedms_def[","\"","Model_1","\"","] = [];
        # Model structure for Model 1. Dict. See Model Generation Section.
        # Note that the order of the inputs is the one defined in this model. If there is a chance that
        # Model_2 has an input that does not exist in this model, these will be appended at the end
        # of the ones in Model_1.

oedms_def[","\"","Model_2","\"","] = [];
        # Model structure for Model 2. Dict. See Model Generation Section.

oedms_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is a vector of strings. These could also be an
        # expression combining states (Only +,-,*,/ and ^ will be considered).
        # Note that due to the use of the covariance matrix in some computations an addition of a small
        # number (0.1) is done in the diagonal elements (for now) to ensure this to be positive definite.
        # If the observable(s) of your model are normalised or in a range between 0 and 1 we recommend to
        # re-scale this to a larger range (0 to 100 for example) so the variances considered are as close
        # as possible to the real ones. This can be done by just multiplying the observable by a constant
        # (Obs*100 for example) in each entry of the vector.

oedms_def[","\"","Theta_M1","\"","] = [];
        # Theta vector (frequentist OED or model with 1 parameter) or matrix (Bayesian OED) for model 1.
        # Path to file with samples is also allowed.

oedms_def[","\"","Theta_M2","\"","] = [];
        # Theta vector (frequentist OED or model with 2 parameter) or matrix (Bayesian OED) for model 2.
        # Path to file with samples is also allowed.

oedms_def[","\"","y0_M1","\"","] = [];
        # For Model 1:
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

oedms_def[","\"","y0_M2","\"","] = [];
        # For Model 2:
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

oedms_def[","\"","preInd_M1","\"","] = [];
        # For Model 1:
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

oedms_def[","\"","preInd_M2","\"","] = [];
        # For Model 2:
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

oedms_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

oedms_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

oedms_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

oedms_def[","\"","fixedInp","\"","] = [];
        # If more than 1 inducer for the system exists, but only 1 input can be dynamic, give a vector
        # of strings indicating which inputs are going to be optimised but as a constant input instead
        # of dynamic. If none, just give an empty vector.
        # If there is only 1 stimuli in the model, this field will be ignored.

oedms_def[","\"","fixedStep","\"","] = [];
        # If you want any of the steps to be fixed to a value. This has to be an empty array if none is fixed
        # or an array of tuples, where each tuple is a step to be fixed. The first entry of the tuple is the
        # index of the step (as an Integer), and the second and array of values for each inducer. Note that
        # the fixed inputs will be ignored, so do not take them into account here.
        # The example type should be Array{Array{Int,1},1} or an empty array ([]) if not used.

oedms_def[","\"","equalStep","\"","] = [];
        # If you want a series of steps to have the same optimised value (for example if you want to design a
        # pulse experiment) you can introduce inside this array different arrays with the indexes of the steps
        # that will have the same value. The values introduced in each array need to be integers.

oedms_def[","\"","plot","\"","] = [];
        # true, flase, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\""," or []
        # indicating if plots with the results will be generated. Default is false.

oedms_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated files so it is not overwritten.
        # If empty, nothing will be added.

            # The order of uUpper and uLower will be taken as the order of inducers defined in Model 1.
            # If Model 2 has an aditional input, this will be added after, so consider this when introducing
            # the bounds for them.
oedms_def[","\"","uUpper","\"","] = [];
        # Vector indicating the upper bounds for the inducers

oedms_def[","\"","uLower","\"","] = [];
        # Vector indicating the lower bounds for the inducers

oedms_def[","\"","maxiter","\"","] = [];
        # Maximum number of iterations for the Bayesian Optimisation. If nothing is introduced a default
        # of 100 iterations will be taken

oedms_def[","\"","maxtime","\"","] = [];
        # Maximum time allowed for the Bayesian Optimisation. If nothing is introduced, only the
        # maximum number of iterations will be taken into account.

"))

    elseif woohoo == "oedmc" || woohoo == "model calibration" || woohoo == "modelcalibration" || replace("woohoo", " "=>"") == "oedmodelcalibration"
print(string(raw"

                CALL defODEModelCalibrStruct()

oedmc_def[","\"","Model","\"","] = [];
        # Dict with Model. See Model Generation Section.

oedmc_def[","\"","Obs","\"","] = [];
        # States of the model that are observables. This is a vector of strings.These could also be an
        # expression combining states (Only +,-,*,/ and ^ will be considered).

oedmc_def[","\"","Theta","\"","] = [];
        # Theta matrix (Bayesian OED) for the model. No single vectors will be allowed.
        # Path to file is also allowed.

oedmc_def[","\"","y0","\"","] = [];
        # Array (single simulation) or matrix (multiple simulations) of Y0s for the simulations for each
        # experiment. If you are computing the steady-state this vector might not be used, however, you
        # still need to introduce it with some random numbers.

oedmc_def[","\"","preInd","\"","] = [];
        # Vector of numbers with the values for the stimuli (inducer) in the over-night. It might be the case
        # that this entry is not required since only the y0 vector is considered for the initial point of
        # the simulation. However, you still need to introduce a random value for it to avoid future issues.

oedmc_def[","\"","finalTime","\"","] = [];
        # Vector of final times for each simulation (initial time will always be assumed as 0,
        # so please consider that).

oedmc_def[","\"","switchT","\"","] = [];
        # Array with the switching times of the inducer in the simulation (time 0 and final
        # time need to be considered)

oedmc_def[","\"","tsamps","\"","] = [];
        # Array of sampling time vectors for the experiments.

oedmc_def[","\"","fixedInp","\"","] = [];
        # If more than 1 inducer for the system exists, but only 1 input can be dynamic, give a vector
        # of strings indicating which inputs are going to be optimised but as a constant input instead
        # of dynamic. If none, just give an empty vector.
        # If there is only 1 stimuli in the model, this field will be ignored.

oedmc_def[","\"","fixedStep","\"","] = [];
        # If you want any of the steps to be fixed to a value. This has to be an empty array if none is fixed
        # or an array of tuples, where each tuple is a step to be fixed. The first entry of the tuple is the
        # index of the step (as an Integer), and the second and array of values for each inducer. Note that
        # the fixed inputs will be ignored, so do not take them into account here.

oedmc_def[","\"","equalStep","\"","] = [];
        # If you want a series of steps to have the same optimised value (for example if you want to design a
        # pulse experiment) you can introduce inside this array different arrays with the indexes of the steps
        # that will have the same value. The values introduced in each array need to be integers.

oedmc_def[","\"","plot","\"","] = [];
        # true, flase, ","\"","Yes","\"",", ","\"","yes","\"",", ","\"","No","\"",", ","\"","no","\""," or []
        # indicating if plots with the results will be generated. Default is false.

oedmc_def[","\"","flag","\"","] = [];
        # String to attach a unique flag to the generated files so it is not overwritten.
        # If empty, nothing will be added.

oedmc_def[","\"","uUpper","\"","] = [];
        # Vector indicating the upper bounds for the inducers

oedmc_def[","\"","uLower","\"","] = [];
        # Vector indicating the lower bounds for the inducers

oedmc_def[","\"","maxiter","\"","] = [];
        # Maximum number of iterations for the Bayesian Optimisation. If nothing is introduced a default
        # of 100 iterations will be taken

oedmc_def[","\"","maxtime","\"","] = [];
        # Maximum time allowed for the Bayesian Optimisation. If nothing is introduced, only the
        # maximum number of iterations will be taken into account.

oedmc_def[","\"","util","\"","] = [];
        # String indicating entropy or perc (or percentile) as the core of the utility function to
        # compute the uncertainty of the model simulations. The default will be to use percentiles.

"))

    else
        print(string("

Please, if you want information about one of the package structures, type:
    1) ","\"","model","\""," for the model generation section
    2) ","\"","simulation","\"",", ","\"","simulation","\""," or ","\"","simul","\""," for the model simulation section
    3) ","\"","pseudo-data","\"",", ","\"","pseudodata","\"",", ","\"","pseudo data","\""," for the model pseudo-data generation section
    4) ","\"","mle","\"",", ","\"","likelihood","\""," for the maximum likelihood estimation section
    5) ","\"","inference","\"",", ","\"","stan","\"",", ","\"","stan inference","\"",", ","\"","staninference","\""," for the Bayesian parameter inference section
    6) ","\"","oedms","\"",", ","\"","model selection","\"",", ","\"","modelselection","\"",", ","\"","oed model selection","\""," for the optimal experimental design for model selection section
    7) ","\"","oedmc","\"",", ","\"","model calibration","\"",", ","\"","modelcalibration","\"",", ","\"","oed model calibration","\""," for the optimal experimental design for model calibration section

        "))
    end

end