text:
function contract!!(
  R::DenseTensor{<:Number,NR},
  labelsR::NTuple{NR},
  T1::DenseTensor{<:Number,N1},
  labelsT1::NTuple{N1},
  T2::DenseTensor{<:Number,N2},
  labelsT2::NTuple{N2},
  α::Number=1,
  β::Number=0,
) where {NR,N1,N2}
  if N1 == 0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    # TODO: replace with an add! function?
    # What about doing `R .= T1[] .* PermutedDimsArray(T2,perm)`?
    perm = getperm(labelsR, labelsT2)
    R = permutedims!!(R, T2, perm, (r, t2) -> T1[] * t2)
  elseif N2 == 0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    perm = getperm(labelsR, labelsT1)
    R = permutedims!!(R, T1, perm, (r, t1) -> T2[] * t1)
  elseif N1 + N2 == NR
    (α ≠ 1 || β ≠ 0) && error(
      "contract!! not yet implemented for outer product tensor contraction with non-trivial α and β",
    )
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R, T1, T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR, labelsRp)
    if !is_trivial_permutation(perm)
      R = permutedims!!(R, copy(R), perm)
    end
  else
    #if dim(T1) > 2^13 && dim(T2) > 2^13
    #    R = _big_contract!!(R,labelsR,T1,labelsT1,T2,labelsT2, α, β)
    #else
    if α ≠ 1 || β ≠ 0
      R = _contract!!(R, labelsR, T1, labelsT1, T2, labelsT2, α, β)
    else
      R = _contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)
    end
    #end
  end
  return R
end

function _big_contract!!(
  R::DenseTensor{<:Number,NR},
  labelsR,
  T1::DenseTensor{ElT1,N1},
  labelsT1,
  T2::DenseTensor{ElT2,N2},
  labelsT2,
  α::Number=1,
  β::Number=0,
) where {ElT1,ElT2,N1,N2,NR}
  props = ContractionProperties(labelsT1, labelsT2, labelsR)
  compute_contraction_properties!(props, T1, T2, R)
  _big_contract!(R, T1, T2, props, α, β)
  #val, t, _ = @timed _blasmg_contract!(R,T1,T2,props,α,β)
  return R
end

function _big_contract!(
  CT::DenseTensor{El,NC},
  AT::DenseTensor{El,NA},
  BT::DenseTensor{El,NB},
  props::ContractionProperties,
  α::Number=one(El),
  β::Number=zero(El),
) where {El,NC,NA,NB}
  Ainds = inds(AT)
  Adims = dims(Ainds)
  Binds = inds(BT)
  Bdims = dims(Binds)
  Cinds = inds(CT)
  Cdims = dims(Cinds)
  Adata = reshape(data(store(AT)), Adims)
  Bdata = reshape(data(store(BT)), Bdims)
  Cdata = reshape(data(store(CT)), Cdims)
  contracted = commoninds(Ainds, Binds)
  A_only = uniqueinds(Ainds, Binds)
  B_only = uniqueinds(Binds, Ainds)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(contracted)
    push!(ind_dict, i)
  end
  if length(A_only) > 0
    for (idx, i) in enumerate(A_only)
      push!(ind_dict, i)
    end
  end
  if length(B_only) > 0
    for (idx, i) in enumerate(B_only)
      push!(ind_dict, i)
    end
  end
  ctainds = zeros(Int, length(Ainds))
  ctbinds = zeros(Int, length(Binds))
  ctcinds = zeros(Int, length(Cinds))
  for (ii, ia) in enumerate(Ainds)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Binds)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  for (ii, ic) in enumerate(Cinds)
    ctcinds[ii] = findfirst(x -> x == ic, ind_dict)
  end
  id_op = CUDA.CUTENSOR.CUTENSOR_OP_IDENTITY
  dict_key = ""
  for cc in zip(ctcinds, Cdims)
    dict_key *= string(cc[1]) * "," * string(cc[2]) * ","
  end
  for aa in zip(ctainds, Adims)
    dict_key *= string(aa[1]) * "," * string(aa[2]) * ","
  end
  for bb in zip(ctbinds, Bdims)
    dict_key *= string(bb[1]) * "," * string(bb[2]) * ","
  end
  #=synchronize()
  if haskey(ENV, "CUTENSOR_AUTOTUNE") && tryparse(Int, ENV["CUTENSOR_AUTOTUNE"]) == 1
      if haskey(ContractionPlans, dict_key)
          dict_val = ContractionPlans[dict_key]
          algo  = dict_val
          Cdata = CUDA.CUTENSOR.contraction!(α, Adata, Vector{Char}(ctainds), id_op, Bdata, Vector{Char}(ctbinds), id_op, β, Cdata, Vector{Char}(ctcinds), id_op, id_op; algo=algo)
          synchronize()
      else
          # loop through all algos
          # pick the fastest one
          # store that plan!
          best_time = 1e6
          best_plan = nothing
          best_algo = nothing
          max_algos = Ref{Int32}(C_NULL)
          CUDA.CUTENSOR.cutensorContractionMaxAlgos(max_algos)
          # fix once the other options are documented
          #algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(max_algos[] - 1))
          algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(-1))
          for algo in reverse(algos)
              try
                  Cdata, this_time, bytes, gctime, memallocs = @timed CUDA.CUTENSOR.contraction!(α, Adata, Vector{Char}(ctainds), id_op, Bdata, Vector{Char}(ctbinds), id_op, β, Cdata, Vector{Char}(ctcinds), id_op, id_op; algo=CUDA.CUTENSOR.cutensorAlgo_t(algo))
                  synchronize()
                  if this_time < best_time
                      best_time = this_time
                      best_algo = CUDA.CUTENSOR.cutensorAlgo_t(algo)
                  end
              catch err
                  @warn "Algorithm $algo not supported"
              end
          end
          ContractionPlans[dict_key] = best_algo
      end
  else
  =#
  Cdata .= zero(eltype(Cdata))
  #@show size(Adata)
  #@show size(Bdata)
  #@show size(Cdata)
  @assert !any(isnan.(Adata))
  @assert !any(isnan.(Bdata))
  @assert !any(isnan.(Cdata))
  #@show ctainds
  #@show ctbinds
  #@show ctcinds
  #flush(stdout)
  CUDA.Mem.pin(Adata)
  CUDA.Mem.pin(Bdata)
  CUDA.Mem.pin(Cdata)
  synchronize()
  #AC = CuArray(Adata)
  #BC = CuArray(Bdata)
  #CC = CuArray(Cdata)
  @assert !any(isnan.(Adata))
  @assert !any(isnan.(Bdata))
  @assert !any(isnan.(Cdata))
  #@assert !any(isnan.(AC))
  #@assert !any(isnan.(BC))
  #@assert !any(isnan.(CC))
  #CC = CUDA.CUTENSOR.contraction!(α, AC, ctainds, id_op, BC, ctbinds, id_op, β, CC, ctcinds, id_op, id_op)
  #synchronize()
  #@assert !any(isnan.(AC))
  #@assert !any(isnan.(BC))
  #@assert !any(isnan.(CC))
  Cdata = CUDA.CUTENSOR.contraction!(
    α, Adata, ctainds, id_op, Bdata, ctbinds, id_op, β, Cdata, ctcinds, id_op, id_op
  )
  synchronize()
  #end
  #CCh = collect(CC)
  #@assert !any(isnan.(CCh))
  #Cdata .= CCh
  @assert !any(isnan.(Adata))
  @assert !any(isnan.(Bdata))
  @assert !any(isnan.(Cdata))
  return parent(Cdata)
end

function _blasmg_contract!(
  CT::DenseTensor{El,NC},
  AT::DenseTensor{El,NA},
  BT::DenseTensor{El,NB},
  props::ContractionProperties,
  α::Number=one(El),
  β::Number=zero(El),
) where {El,NC,NA,NB}
  # TODO: directly use Tensor instead of Array
  C = array(CT)
  A = array(AT)
  B = array(BT)

  tA = 'N'
  if props.permuteA
    pA = NTuple{NA,Int}(props.PA)
    @strided Ap = permutedims(A, pA)
    AM = reshape(Ap, props.dmid, props.dleft)
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(props)
      AM = reshape(A, props.dmid, props.dleft)
      tA = 'T'
    else
      AM = reshape(A, props.dleft, props.dmid)
    end
  end

  tB = 'N'
  if props.permuteB
    pB = NTuple{NB,Int}(props.PB)
    @strided Bp = permutedims(B, pB)
    BM = reshape(Bp, props.dmid, props.dright)
  else
    if Btrans(props)
      BM = reshape(B, props.dright, props.dmid)
      tB = 'T'
    else
      BM = reshape(B, props.dmid, props.dright)
    end
  end

  # TODO: this logic may be wrong
  if props.permuteC
    # Need to copy here since we will be permuting
    # into C later
    CM = reshape(copy(C), props.dleft, props.dright)
  else
    if Ctrans(props)
      CM = reshape(C, props.dleft, props.dright)
      (AM, BM) = (BM, AM)
      if tA == tB
        tA = tB = (tA == 'T' ? 'N' : 'T')
      end
    else
      CM = reshape(C, props.dleft, props.dright)
    end
  end

  if length(AM) > 4096 && length(BM) > 4096 && length(CM) > 4096
    CM = CUBLASMG.mg_gemm!(
      tA,
      tB,
      El(α),
      AM,
      BM,
      El(β),
      CM;
      devs=devs[],
      dev_rows=dev_rows[],
      dev_cols=dev_cols[],
    )
  else
    BLAS.gemm!(tA, tB, El(α), AM, BM, El(β), CM)
  end

  if props.permuteC
    pC = NTuple{NC,Int}(props.PC)
    Cr = reshape(CM, props.newCrange...)
    @strided C .= permutedims(Cr, pC)
  end
  return C
end


text:
using Test
using ClimateMachine
using ClimateMachine.VariableTemplates: @vars, Vars
using ClimateMachine.BalanceLaws
using ClimateMachine.DGMethods: AbstractCustomFilter, apply!
import ClimateMachine
import ClimateMachine.BalanceLaws:
    vars_state, init_state_auxiliary!, init_state_prognostic!
using MPI
using LinearAlgebra

struct CustomFilterTestModel <: BalanceLaw end
struct CustomTestFilter <: AbstractCustomFilter end

vars_state(::CustomFilterTestModel, ::Auxiliary, FT) = @vars()
vars_state(::CustomFilterTestModel, ::Prognostic, FT) where {N} = @vars(q::FT)

init_state_auxiliary!(::CustomFilterTestModel, _...) = nothing

function init_state_prognostic!(
    ::CustomFilterTestModel,
    state::Vars,
    aux::Vars,
    localgeo,
)
    coord = localgeo.coord
    state.q = hypot(coord[1], coord[2])
end

@testset "Test custom filter" begin
    let
        ClimateMachine.init()
        N = 4
        Ne = (2, 2, 2)

        function ClimateMachine.DGMethods.custom_filter!(
            ::CustomTestFilter,
            bl::CustomFilterTestModel,
            state::Vars,
            aux::Vars,
        )
            state.q = state.q^2
        end

        @testset for FT in (Float64, Float32)
            dim = 2
            brickrange =
                ntuple(j -> range(FT(-1); length = Ne[j] + 1, stop = 1), dim)
            topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                MPI.COMM_WORLD,
                brickrange,
                periodicity = ntuple(j -> true, dim),
            )

            grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                topl,
                FloatType = FT,
                DeviceArray = ClimateMachine.array_type(),
                polynomialorder = N,
            )

            model = CustomFilterTestModel()
            dg = ClimateMachine.DGMethods.DGModel(
                model,
                grid,
                nothing,
                nothing,
                nothing;
                state_gradient_flux = nothing,
            )

            Q = ClimateMachine.DGMethods.init_ode_state(dg)
            data = Array(Q.realdata)

            apply!(CustomTestFilter(), grid, model, Q, dg.state_auxiliary)

            @test all(Array(Q.realdata) .== data .^ 2)
        end
    end
end

ext: jl


text:
immutable Stockholm <: Format end

@inline function _fill_with_sequence_line!(line, IDS, SEQS)
    if line[1:1] != "#" && line[1:2] != "//"
        words = get_n_words(line,2)
        push!(IDS, words[1])
        push!(SEQS,words[2])
    end
end

function _fill_with_line!(line, IDS, SEQS, GF, GS, GC, GR)
    if line[1:4] == "#=GF"
        words = get_n_words(line,3)
        id = words[2]
        if id in keys(GF)
            GF[ id ] = GF[ id ] * "\n" * words[3]
        else
            GF[ id ] = words[3]
        end
    elseif line[1:4] == "#=GS"
        words = get_n_words(line,4)
        idtuple = (words[2],words[3])
        if idtuple in keys(GS)
            GS[ idtuple ] = GS[ idtuple ] * "\n" * words[4]
        else
            GS[ idtuple ] = words[4]
        end
    elseif line[1:4] == "#=GC"
        words = get_n_words(line,3)
        GC[words[2]] = words[3]
    elseif line[1:4] == "#=GR"
        words = get_n_words(line,4)
        GR[(words[2],words[3])] = words[4]
    else
        _fill_with_sequence_line!(line, IDS, SEQS)
    end
end

function _pre_readstockholm(io::Union{IO, AbstractString})
    IDS  = ASCIIString[]
    SEQS = ASCIIString[]
    GF = OrderedDict{ASCIIString,ByteString}()
    GC = Dict{ASCIIString,ASCIIString}()
    GS = Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}()
    GR = Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}()

    @inbounds for line in eachline(io)
        if length(line) >= 4
            _fill_with_line!(line, IDS, SEQS, GF, GS, GC, GR)
            if line[1:2] == "//"
               break
            end
        end
    end

    GF = sizehint!(GF, length(GF))
    GC = sizehint!(GC, length(GC))
    GS = sizehint!(GS, length(GS))
    GR = sizehint!(GR, length(GR))
    (IDS, SEQS, GF, GS, GC, GR)
end

function _pre_readstockholm_sequences(io::Union{IO, AbstractString})
    IDS  = ASCIIString[]
    SEQS = ASCIIString[]
    @inbounds for line in eachline(io)
        if length(line) >= 4
            _fill_with_sequence_line!(line, IDS, SEQS)
        end
    end
    (IDS, SEQS)
end

function parse(io::Union{IO, AbstractString}, format::Type{Stockholm},
               output::Type{AnnotatedMultipleSequenceAlignment}; generatemapping::Bool=false,
               useidcoordinates::Bool=false, deletefullgaps::Bool=true,
               checkalphabet::Bool=false, keepinserts::Bool=false)
    IDS, SEQS, GF, GS, GC, GR = _pre_readstockholm(io)
    annot = Annotations(GF, GS, GC, GR)
    if keepinserts
        _keepinserts!(SEQS, annot)
    end
    if generatemapping
        MSA, MAP = useidcoordinates && hascoordinates(IDS[1]) ? _to_msa_mapping(SEQS, IDS) : _to_msa_mapping(SEQS)
        setannotfile!(annot, "NCol", string(size(MSA,2)))
        setannotfile!(annot, "ColMap", join(vcat(1:size(MSA,2)), ','))
        for i in 1:length(IDS)
            setannotsequence!(annot, IDS[i], "SeqMap", MAP[i])
        end
    else
        MSA = convert(Matrix{Residue}, SEQS)
    end
    msa = AnnotatedMultipleSequenceAlignment(IndexedArray(IDS), MSA, annot)
    if checkalphabet
        deletenotalphabetsequences!(msa, SEQS)
    end
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

function parse(io::Union{IO, AbstractString}, format::Type{Stockholm}, output::Type{MultipleSequenceAlignment};
               deletefullgaps::Bool=true, checkalphabet::Bool=false)
    IDS, SEQS = _pre_readstockholm_sequences(io)
    msa = MultipleSequenceAlignment(IndexedArray(IDS), convert(Matrix{Residue}, SEQS))
    if checkalphabet
        deletenotalphabetsequences!(msa, SEQS)
    end
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

function parse(io::Union{IO,AbstractString}, format::Type{Stockholm}, output::Type{Matrix{Residue}};
               deletefullgaps::Bool=true, checkalphabet::Bool=false)
    IDS, SEQS = _pre_readstockholm_sequences(io)
    _strings_to_msa(SEQS, deletefullgaps, checkalphabet)
end

parse(io, format::Type{Stockholm};  generatemapping::Bool=false,
      useidcoordinates::Bool=false, deletefullgaps::Bool=true,
      checkalphabet::Bool=false, keepinserts::Bool=false) = parse(io, Stockholm, AnnotatedMultipleSequenceAlignment,
                                                                  generatemapping=generatemapping,
                                                                  useidcoordinates=useidcoordinates,
                                                                  deletefullgaps=deletefullgaps,
                                                                  checkalphabet=checkalphabet,
                                                                  keepinserts=keepinserts)

# Print Pfam
# ==========

function _to_sequence_dict(annotation::Dict{Tuple{ASCIIString,ASCIIString},ASCIIString})
    seq_dict = Dict{ASCIIString,Vector{ASCIIString}}()
    for (key, value) in annotation
        seq_id = key[1]
        if seq_id in keys(seq_dict)
            push!(seq_dict[seq_id], string(key[2], '\t', value))
        else
            seq_dict[seq_id] = [ string(key[2], '\t', value) ]
        end
    end
    sizehint!(seq_dict, length(seq_dict))
end

function print(io::IO, msa::AnnotatedMultipleSequenceAlignment, format::Type{Stockholm})
    _printfileannotations(io, msa.annotations)
    _printsequencesannotations(io, msa.annotations)
    res_annotations = _to_sequence_dict(msa.annotations.residues)
    for i in 1:nsequences(msa)
        id = msa.id[i]
        seq = asciisequence(msa, i)
        println(io, string(id, "\t\t", seq))
        if id in keys(res_annotations)
            for line in res_annotations[id]
                println(io, string("#=GR ", id, '\t', line))
            end
        end
    end
    _printcolumnsannotations(io, msa.annotations)
    println(io, "//")
end

function print(io::IO, msa::MultipleSequenceAlignment, format::Type{Stockholm})
    for i in 1:nsequences(msa)
        id = msa.id[i]
        seq = asciisequence(msa, i)
        println(io, string(id, "\t\t", seq))
    end
    println(io, "//")
end

print(msa::AnnotatedMultipleSequenceAlignment) = print(STDOUT, msa, Stockholm)


text:
rename(env, ex) = IRTools.prewalk(x -> x isa Variable ? env[x] : x, ex)
rename(::Nothing, ex) = ex

function xidentity()
  ir = IR()
  argument!(ir)
  return ir
end

function branches(ir, b)
  brs = IRTools.branches(b)
  meta = Dict()
  for br in brs[1:end-1]
    cond = push!(ir, xcall(Not(), br.condition))
    meta[br.block] = (cond, br.args)
  end
  meta[brs[end].block] = (true, brs[end].args)
  return meta
end

function extract(b)
  ir = IR(b)
  empty!(ir.blocks[1].branches)
  return ir
end

function reloop!(out, ir, cfg::IRTools.Multiple, brs)
  nbrs = length(brs)
  (cond, targs) = pop!(brs, cfg.inner[1].block)
  tfunc = push!(out, reloop_tupleargs(ir, cfg.inner[1]))
  if length(cfg.inner) == 2 && brs[cfg.inner[2].block][1] == true
    fargs = brs[cfg.inner[2].block][2]
    ffunc = push!(out, reloop_tupleargs(ir, cfg.inner[2]))
  elseif length(cfg.inner) == 1 && length(brs) == 1
    fargs = first(brs)[2][2]
    ffunc = push!(out, xidentity())
  else
    error("Multiple block: not implemented")
  end
  args = push!(out, xcall(Conditional(), cond,
                          xcall(XTuple(), targs...), tfunc,
                          xcall(XTuple(), fargs...), ffunc))
  reloop!(out, ir, cfg.next, args)
  return out
end

function reloop_next!(out, cfg, b, env = nothing)
  if cfg.next != nothing
    reloop!(out, b.ir, cfg.next, branches(out, b)) # TODO pass env
  elseif IRTools.isreturn(b)
    return
  else
    @assert length(IRTools.branches(b)) == 1 "Unstructured control flow not implemented"
    push!(out, xcall(XTuple(), rename.((env,), arguments(IRTools.branches(b)[1]))...))
  end
end

function reloop!(out, ir, cfg::IRTools.Simple, args)
  b = block(ir, cfg.block)
  env = Dict()
  for (i, k) in enumerate(arguments(b))
    env[k] = push!(out, xcall(GetTupleElement(i-1), args))
  end
  for (v, st) in b
    env[v] = push!(out, rename(env, st.expr))
  end
  reloop_next!(out, cfg, b, env)
  return out
end

function reloop_tupleargs(ir, cfg::IRTools.Simple)
  out = IR()
  args = argument!(out)
  reloop!(out, ir, cfg, args)
  return out
end

function reloop(ir, cfg::IRTools.Simple)
  b = block(ir, cfg.block)
  out = extract(b)
  reloop_next!(out, cfg, b)
  return out
end

function controlflow(ir::IR)
  length(blocks(ir)) == 1 && return ir
  ir = ir |> copy |> explicitbranch! |> merge_returns! |> expand!
  cfg = IRTools.reloop(CFG(ir))
  reloop(ir, cfg)
end


text:
if basename(pwd()) == "aoc"
    cd("2021/15")
end
using Pkg
Pkg.activate(".")
using Base.Iterators
using Graphs, SimpleWeightedGraphs

function loadfile(filename::AbstractString)
    transpose(reshape(
        [parse(Int, c) for line in eachline(filename) for c in line],
        length(readline(filename)),
        :
    ))
end

function loadbig(filename::AbstractString)
    A = loadfile(filename)
    B = hcat((A .+ n for n = 0:4)...)
    C = vcat((B .+ n for n = 0:4)...)
    (C .- 1) .% 9 .+ 1
end

function neighbourhood(A::AbstractMatrix, i::CartesianIndex)
    (
        n + i
        for n in CartesianIndex.([(1, 0), (0, 1), (-1, 0), (0, -1)])
        if n + i in CartesianIndices(A)
    )
end

tolinear(A, i::CartesianIndex) = size(A, 2) * (i[1] - 1) + i[2]

function todigraph(A::AbstractMatrix{Int})::SimpleWeightedDiGraph
    g = SimpleWeightedDiGraph{Int,Int}(length(A))
    for i in CartesianIndices(A)
        for n in neighbourhood(A, i)
            @assert add_edge!(g, tolinear(A, i), tolinear(A, n), A[n])
        end
    end
    g
end

minpath(g::AbstractGraph) = last(dijkstra_shortest_paths(g, 1).dists)
minpath(A::AbstractMatrix) = minpath(todigraph(A))

part1(filename) = minpath(loadfile(filename))
part2(filename) = minpath(loadbig(filename))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Reproject a Geometry" begin
    @testset "Method 1" begin
        source = AG.unsafe_importEPSG(2927); target = AG.unsafe_importEPSG(4326)
            transform = AG.unsafe_createcoordtrans(source, target)
                point = AG.unsafe_fromWKT("POINT (1120351.57 741921.42)")
                    @test AG.toWKT(point) == "POINT (1120351.57 741921.42)"
                    AG.transform!(point, transform)
                    @test AG.toWKT(point) == "POINT (-122.598135130878 47.3488013802885)"
                AG.destroy(point)
            AG.destroy(transform)
        AG.destroy(target); AG.destroy(source)
    end

    @testset "Method 2" begin
        AG.importEPSG(2927) do source; AG.importEPSG(4326) do target
            AG.createcoordtrans(source, target) do transform
                AG.fromWKT("POINT (1120351.57 741921.42)") do point
                    @test AG.toWKT(point) == "POINT (1120351.57 741921.42)"
                    AG.transform!(point, transform)
                    @test AG.toWKT(point) == "POINT (-122.598135130878 47.3488013802885)"
        end end end end
    end
end

@testset "Get Projection" begin
    AG.registerdrivers() do
        AG.read("ospy/data1/sites.shp") do dataset
            layer = AG.getlayer(dataset, 0)
            spatialref = AG.getspatialref(layer)
            @test AG.toWKT(spatialref) == """PROJCS["WGS_1984_UTM_Zone_12N",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-111.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG","32612"]]"""
            @test AG.toWKT(spatialref, false) == """
            PROJCS["WGS_1984_UTM_Zone_12N",
                GEOGCS["GCS_WGS_1984",
                    DATUM["WGS_1984",
                        SPHEROID["WGS_84",6378137.0,298.257223563]],
                    PRIMEM["Greenwich",0.0],
                    UNIT["Degree",0.0174532925199433],
                    AUTHORITY["EPSG","4326"]],
                PROJECTION["Transverse_Mercator"],
                PARAMETER["False_Easting",500000.0],
                PARAMETER["False_Northing",0.0],
                PARAMETER["Central_Meridian",-111.0],
                PARAMETER["Scale_Factor",0.9996],
                PARAMETER["Latitude_Of_Origin",0.0],
                UNIT["Meter",1.0],
                AUTHORITY["EPSG","32612"]]"""
            @test AG.toWKT(spatialref, true) == """
            PROJCS["WGS_1984_UTM_Zone_12N",
                GEOGCS["GCS_WGS_1984",
                    DATUM["WGS_1984",
                        SPHEROID["WGS_84",6378137.0,298.257223563]],
                    PRIMEM["Greenwich",0.0],
                    UNIT["Degree",0.0174532925199433]],
                PROJECTION["Transverse_Mercator"],
                PARAMETER["False_Easting",500000.0],
                PARAMETER["False_Northing",0.0],
                PARAMETER["Central_Meridian",-111.0],
                PARAMETER["Scale_Factor",0.9996],
                PARAMETER["Latitude_Of_Origin",0.0],
                UNIT["Meter",1.0]]"""
            @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs "
            @test AG.toXML(spatialref)[1:17] == "<gml:ProjectedCRS"
            @test AG.toMICoordSys(spatialref) == "Earth Projection 8, 104, \"m\", -111, 0, 0.9996, 500000, 0"
            AG.nextfeature(layer) do feature
                @test AG.toWKT(AG.getspatialref(AG.getgeom(feature))) == """PROJCS["WGS_1984_UTM_Zone_12N",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-111.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG","32612"]]"""
            end
        end
    end

    AG.importEPSG(26912) do spatialref
        @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
        @test AG.toWKT(spatialref) == """PROJCS["NAD83 / UTM zone 12N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26912"]]"""
        AG.morphtoESRI!(spatialref)
        @test AG.toPROJ4(spatialref) == "+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs "
        AG.toWKT(spatialref) == """PROJCS["NAD_1983_UTM_Zone_12N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]"""
    end
end


text:
using UnitfulAstro; using Unitful

"""
kep_deriv() function created by examples/create_keplerian_derivatives.jl.
Derivative of a Keplerian radial velocity signal using h and k instead of
e and ω

Parameters:
K (Unitful.Velocity): velocity semi-amplitude
P (Unitful.Time): period
M0 (real): initial mean anomaly
h (real): eccentricity * sin(argument of periastron)
k (real): eccentricity * cos(argument of periastron)
γ (Unitful.Velocity): velocity offset
t (Unitful.Time): time
dorder (vector of integers): A vector of how many partial derivatives you want
	to take with respect to each variable in the input order

Returns:
Unitful Quantity: The derivative specified with dorder
"""
function kep_deriv(
	K::Unitful.Velocity,
	P::Unitful.Time,
	M0::Real,
	h::Real,
	k::Real,
	γ::Unitful.Velocity,
	t::Unitful.Time,
	dorder::Vector{<:Integer})

	validate_kepler_dorder(dorder)

	esq = h * h + k * k
	e = sqrt(esq)
	jsq = 1 - esq
	j = sqrt(jsq)
	EAval = ecc_anomaly(t, P, M0, e)
	cosEA = cos(EAval)
	sinEA = sin(EAval)
	qmod = 1 - e * cosEA

    if dorder ==[0, 0, 0, 0, 0, 0]
        func = γ + K*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod
    end

    if dorder ==[1, 0, 0, 0, 0, 0]
        func = (-h*j*sinEA/e + k*jsq*cosEA/e)/qmod
    end

    if dorder ==[0, 1, 0, 0, 0, 0]
        func = K*(2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/qmod + 2*K*e*t*sinEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^2*qmod^3)
    end

    if dorder ==[0, 0, 1, 0, 0, 0]
        func = K*(h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))/qmod + K*e*sinEA*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^3
    end

    if dorder ==[0, 0, 0, 1, 0, 0]
        func = K*(-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))/qmod - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[0, 0, 0, 0, 1, 0]
        func = K*(jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))/qmod - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[0, 0, 0, 0, 0, 1]
        func = 1
    end

    if dorder ==[2, 0, 0, 0, 0, 0]
        func = 0
    end

    if dorder ==[1, 1, 0, 0, 0, 0]
        func = (2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/qmod + 2*e*t*sinEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^2*qmod^3)
    end

    if dorder ==[0, 2, 0, 0, 0, 0]
        func = K*(4*k*t^2*jsq*sinEA^2*π^2/(P^4*qmod^3) + 4*h*j*t^2*sinEA*π^2/(P^4*e*qmod^2) - 4*k*t^2*jsq*cosEA*π^2/(P^4*e*qmod^2) + 4*h*j*t^2*sinEA*cosEA*π^2/(P^4*qmod^3) - 4*h*j*t*cosEA*π/(P^3*e*qmod) - 4*k*t*jsq*sinEA*π/(P^3*e*qmod))/qmod - 4*K*e*t^2*cosEA*π^2*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^4*qmod^4) + 12*K*t^2*esq*sinEA^2*π^2*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^4*qmod^5) - 4*K*e*t*sinEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^3*qmod^3) + 4*K*e*t*sinEA*π*(2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/(P^2*qmod^3)
    end

    if dorder ==[1, 0, 1, 0, 0, 0]
        func = (h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))/qmod + e*sinEA*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^3
    end

    if dorder ==[0, 1, 1, 0, 0, 0]
        func = K*(2*k*t*jsq*sinEA^2*π/(P^2*qmod^3) + 2*h*j*t*sinEA*π/(P^2*e*qmod^2) - 2*k*t*jsq*cosEA*π/(P^2*e*qmod^2) + 2*h*j*t*sinEA*cosEA*π/(P^2*qmod^3))/qmod + K*e*sinEA*(2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/qmod^3 - 2*K*e*t*cosEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^2*qmod^4) + 2*K*e*t*sinEA*π*(h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))/(P^2*qmod^3) + 6*K*t*esq*sinEA^2*π*(-h*j*sinEA/e + k*jsq*cosEA/e)/(P^2*qmod^5)
    end

    if dorder ==[0, 0, 2, 0, 0, 0]
        func = K*(k*jsq*sinEA^2/qmod^3 + h*j*sinEA/(e*qmod^2) - k*jsq*cosEA/(e*qmod^2) + h*j*sinEA*cosEA/qmod^3)/qmod - K*e*cosEA*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^4 + 2*K*e*sinEA*(h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))/qmod^3 + 3*K*esq*sinEA^2*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^5
    end

    if dorder ==[1, 0, 0, 1, 0, 0]
        func = (-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))/qmod - (-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[0, 1, 0, 1, 0, 0]
        func = K*(2*j*t*cosEA*π/(P^2*e*qmod) - 2*h^2*j*t*cosEA*π/(P^2*e^3*qmod) + 2*h^2*j*t*cosEA^2*π/(P^2*e^2*qmod^2) - 2*h^2*j*t*sinEA^2*π/(P^2*e^2*qmod^2) - 4*h*k*t*sinEA*π/(P^2*e*qmod) - 2*h^2*t*cosEA*π/(P^2*e*j*qmod) - 2*h*k*t*jsq*sinEA*π/(P^2*e^3*qmod) - 2*h*k*t*jsq*sinEA^3*π/(P^2*e*qmod^3) - 2*h^2*j*t*sinEA^2*cosEA*π/(P^2*e*qmod^3) + 4*h*k*t*jsq*sinEA*cosEA*π/(P^2*e^2*qmod^2))/qmod - K*(-h*cosEA/e + h*sinEA^2/qmod)*(2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/qmod^2 - K*(-2*h*t*sinEA*π/(P^2*e*qmod) + 2*e*h*t*sinEA^3*π/(P^2*qmod^3) - 4*h*t*sinEA*cosEA*π/(P^2*qmod^2))*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^2 + 2*K*e*t*sinEA*π*(-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))/(P^2*qmod^3) - 4*K*e*t*sinEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*cosEA/e + h*sinEA^2/qmod)/(P^2*qmod^4)
    end

    if dorder ==[0, 0, 1, 1, 0, 0]
        func = K*(j*cosEA/(e*qmod) - h^2*j*cosEA/(e^3*qmod) + h^2*j*cosEA^2/(e^2*qmod^2) - h^2*j*sinEA^2/(e^2*qmod^2) - 2*h*k*sinEA/(e*qmod) - h^2*cosEA/(e*j*qmod) - h*k*jsq*sinEA/(e^3*qmod) - h*k*jsq*sinEA^3/(e*qmod^3) - h^2*j*sinEA^2*cosEA/(e*qmod^3) + 2*h*k*jsq*sinEA*cosEA/(e^2*qmod^2))/qmod - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*sinEA/(e*qmod) + e*h*sinEA^3/qmod^3 - 2*h*sinEA*cosEA/qmod^2)/qmod^2 - K*(h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2 + K*e*sinEA*(-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))/qmod^3 - 2*K*e*sinEA*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^4
    end

    if dorder ==[0, 0, 0, 2, 0, 0]
        func = K*(-2*k*cosEA/e - 3*h^3*j*sinEA/e^5 + 3*h*j*sinEA/e^3 + 4*h^2*k*cosEA/e^3 - 2*h^3*sinEA/(e^3*j) - k*jsq*cosEA/e^3 + 3*h*sinEA/(e*j) + h^3*sinEA/(e*j^3) + 3*h^2*k*jsq*cosEA/e^5 + h^3*j*sinEA^3/(e^3*qmod^2) + 4*h^2*k*sinEA^2/(e^2*qmod) - k*jsq*sinEA^2/(e^2*qmod) + 3*h^2*k*jsq*sinEA^2/(e^4*qmod) + 3*h^3*j*sinEA*cosEA/(e^4*qmod) - h^3*j*sinEA*cosEA^2/(e^3*qmod^2) - 3*h*j*sinEA*cosEA/(e^2*qmod) + 2*h^3*sinEA*cosEA/(e^2*j*qmod) - 2*h^2*k*jsq*sinEA^2*cosEA/(e^3*qmod^2) + h*k*jsq*sinEA^2*(-h*cosEA/e + h*sinEA^2/qmod)/(e^2*qmod^2) + h^2*j*sinEA*cosEA*(-h*cosEA/e + h*sinEA^2/qmod)/(e^2*qmod^2))/qmod + 2*K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-h*cosEA/e + h*sinEA^2/qmod)^2/qmod^3 - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-cosEA/e + sinEA^2/qmod + h^2*cosEA/e^3 + h^2*sinEA^2/(e^2*qmod) - h*sinEA^2*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2 + 2*h^2*sinEA^2*cosEA/(e*qmod^2))/qmod^2 - 2*K*(-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[1, 0, 0, 0, 1, 0]
        func = (jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))/qmod - (-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[0, 1, 0, 0, 1, 0]
        func = K*(-4*k^2*t*sinEA*π/(P^2*e*qmod) + 2*t*jsq*sinEA*π/(P^2*e*qmod) - 2*k^2*t*jsq*sinEA*π/(P^2*e^3*qmod) - 2*k^2*t*jsq*sinEA^3*π/(P^2*e*qmod^3) - 2*h*j*k*t*cosEA*π/(P^2*e^3*qmod) + 2*h*j*k*t*cosEA^2*π/(P^2*e^2*qmod^2) - 2*h*j*k*t*sinEA^2*π/(P^2*e^2*qmod^2) + 4*k^2*t*jsq*sinEA*cosEA*π/(P^2*e^2*qmod^2) - 2*h*k*t*cosEA*π/(P^2*e*j*qmod) - 2*h*j*k*t*sinEA^2*cosEA*π/(P^2*e*qmod^3))/qmod - K*(-k*cosEA/e + k*sinEA^2/qmod)*(2*h*j*t*cosEA*π/(P^2*e*qmod) + 2*k*t*jsq*sinEA*π/(P^2*e*qmod))/qmod^2 - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-2*k*t*sinEA*π/(P^2*e*qmod) + 2*e*k*t*sinEA^3*π/(P^2*qmod^3) - 4*k*t*sinEA*cosEA*π/(P^2*qmod^2))/qmod^2 + 2*K*e*t*sinEA*(jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))*π/(P^2*qmod^3) - 4*K*e*t*sinEA*π*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)/(P^2*qmod^4)
    end

    if dorder ==[0, 0, 1, 0, 1, 0]
        func = K*(jsq*sinEA/(e*qmod) - 2*k^2*sinEA/(e*qmod) - k^2*jsq*sinEA/(e^3*qmod) - k^2*jsq*sinEA^3/(e*qmod^3) - h*j*k*cosEA/(e^3*qmod) + h*j*k*cosEA^2/(e^2*qmod^2) - h*j*k*sinEA^2/(e^2*qmod^2) + 2*k^2*jsq*sinEA*cosEA/(e^2*qmod^2) - h*k*cosEA/(e*j*qmod) - h*j*k*sinEA^2*cosEA/(e*qmod^3))/qmod - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*sinEA/(e*qmod) + e*k*sinEA^3/qmod^3 - 2*k*sinEA*cosEA/qmod^2)/qmod^2 - K*(h*j*cosEA/(e*qmod) + k*jsq*sinEA/(e*qmod))*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2 + K*e*sinEA*(jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))/qmod^3 - 2*K*e*sinEA*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^4
    end

    if dorder ==[0, 0, 0, 1, 1, 0]
        func = K*(-2*h*cosEA/e - h*jsq*cosEA/e^3 + 4*h*k^2*cosEA/e^3 + j*k*sinEA/e^3 + k*sinEA/(e*j) + 3*h*k^2*jsq*cosEA/e^5 - 3*h^2*j*k*sinEA/e^5 - 2*h^2*k*sinEA/(e^3*j) - h*jsq*sinEA^2/(e^2*qmod) + 4*h*k^2*sinEA^2/(e^2*qmod) + h^2*k*sinEA/(e*j^3) + 3*h*k^2*jsq*sinEA^2/(e^4*qmod) + h^2*j*k*sinEA^3/(e^3*qmod^2) - j*k*sinEA*cosEA/(e^2*qmod) + k^2*jsq*sinEA^2*(-h*cosEA/e + h*sinEA^2/qmod)/(e^2*qmod^2) + 3*h^2*j*k*sinEA*cosEA/(e^4*qmod) - 2*h*k^2*jsq*sinEA^2*cosEA/(e^3*qmod^2) - h^2*j*k*sinEA*cosEA^2/(e^3*qmod^2) + 2*h^2*k*sinEA*cosEA/(e^2*j*qmod) + h*j*k*sinEA*cosEA*(-h*cosEA/e + h*sinEA^2/qmod)/(e^2*qmod^2))/qmod - K*(h*k*cosEA/e^3 - k*sinEA^2*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2 + h*k*sinEA^2/(e^2*qmod) + 2*h*k*sinEA^2*cosEA/(e*qmod^2))*(-h*j*sinEA/e + k*jsq*cosEA/e)/qmod^2 - K*(-j*sinEA/e + h^2*j*sinEA/e^3 - 2*h*k*cosEA/e + h^2*sinEA/(e*j) - h*k*jsq*cosEA/e^3 - h*k*jsq*sinEA^2/(e^2*qmod) - h^2*j*sinEA*cosEA/(e^2*qmod))*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2 - K*(jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^2 + 2*K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)*(-h*cosEA/e + h*sinEA^2/qmod)/qmod^3
    end

    if dorder ==[0, 0, 0, 0, 2, 0]
        func = K*(4*k^3*cosEA/e^3 - 6*k*cosEA/e + 3*k^3*jsq*cosEA/e^5 + h*j*sinEA/e^3 - 3*k*jsq*cosEA/e^3 + 4*k^3*sinEA^2/(e^2*qmod) + h*sinEA/(e*j) - 3*h*j*k^2*sinEA/e^5 + 3*k^3*jsq*sinEA^2/(e^4*qmod) - 2*h*k^2*sinEA/(e^3*j) - 3*k*jsq*sinEA^2/(e^2*qmod) + h*k^2*sinEA/(e*j^3) + h*j*k^2*sinEA^3/(e^3*qmod^2) - 2*k^3*jsq*sinEA^2*cosEA/(e^3*qmod^2) - h*j*sinEA*cosEA/(e^2*qmod) + k^2*jsq*sinEA^2*(-k*cosEA/e + k*sinEA^2/qmod)/(e^2*qmod^2) + 3*h*j*k^2*sinEA*cosEA/(e^4*qmod) - h*j*k^2*sinEA*cosEA^2/(e^3*qmod^2) + 2*h*k^2*sinEA*cosEA/(e^2*j*qmod) + h*j*k*sinEA*cosEA*(-k*cosEA/e + k*sinEA^2/qmod)/(e^2*qmod^2))/qmod + 2*K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-k*cosEA/e + k*sinEA^2/qmod)^2/qmod^3 - K*(-h*j*sinEA/e + k*jsq*cosEA/e)*(-cosEA/e + sinEA^2/qmod + k^2*cosEA/e^3 + k^2*sinEA^2/(e^2*qmod) - k*sinEA^2*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2 + 2*k^2*sinEA^2*cosEA/(e*qmod^2))/qmod^2 - 2*K*(jsq*cosEA/e - 2*k^2*cosEA/e - k^2*jsq*cosEA/e^3 + h*j*k*sinEA/e^3 - k^2*jsq*sinEA^2/(e^2*qmod) + h*k*sinEA/(e*j) - h*j*k*sinEA*cosEA/(e^2*qmod))*(-k*cosEA/e + k*sinEA^2/qmod)/qmod^2
    end

    if dorder ==[1, 0, 0, 0, 0, 1]
        func = 0
    end

    if dorder ==[0, 1, 0, 0, 0, 1]
        func = 0
    end

    if dorder ==[0, 0, 1, 0, 0, 1]
        func = 0
    end

    if dorder ==[0, 0, 0, 1, 0, 1]
        func = 0
    end

    if dorder ==[0, 0, 0, 0, 1, 1]
        func = 0
    end

    if dorder ==[0, 0, 0, 0, 0, 2]
        func = 0
    end

    return float(func)

end
kep_deriv(ks::kep_signal, t::Unitful.Time, dorder::Vector{<:Integer}) =
    kep_deriv(ks.K, ks.P, ks.M0, ks.h, ks.k, ks.γ, t, dorder)


text:
#the strategy is the following:
#== UnivariateFunction is a type hierarchy do dispatch on functions with ranges.
it's tighly coupled with the VariableProperty struct.
Calling a variable property will return a result depending on the type of UnivariateFunction
This gives the advantage to decouple the use of a dict in each call of a property.
==#
abstract type UnivariateFunction end
struct ConstantEq <: UnivariateFunction end
struct LinearEq <: UnivariateFunction end
struct CuadraticEq <: UnivariateFunction end
struct CubicPolyEq <: UnivariateFunction end
struct EqNo5 <: UnivariateFunction end
struct EqNo6 <: UnivariateFunction end
struct AntoineEq <: UnivariateFunction end
struct Poly4Eq <: UnivariateFunction end
struct DIPPREq <: UnivariateFunction end
struct EqNo102 <: UnivariateFunction end
struct EqNo103 <: UnivariateFunction end
struct EqNo104 <: UnivariateFunction end
struct EqNo105 <: UnivariateFunction end
struct EqNo106 <: UnivariateFunction end
struct EqNo107 <: UnivariateFunction end
struct EqNo114 <: UnivariateFunction end
struct EqNo117 <: UnivariateFunction end

struct EqNo11 <: UnivariateFunction end
struct EqNo12 <: UnivariateFunction end
struct EqNo13 <: UnivariateFunction end
struct EqNo14 <: UnivariateFunction end
struct EqNo15 <: UnivariateFunction end
struct EqNo16 <: UnivariateFunction end
struct EqNo17 <: UnivariateFunction end
struct EqNo18 <: UnivariateFunction end
struct EqNo19 <: UnivariateFunction end
struct EqNo45 <: UnivariateFunction end
struct EqNo75 <: UnivariateFunction end
struct EqNo118 <: UnivariateFunction end
struct EqNo120 <: UnivariateFunction end
struct EqNo121 <: UnivariateFunction end
struct EqNo122 <: UnivariateFunction end
struct EqNo123 <: UnivariateFunction end
struct EqNo124 <: UnivariateFunction end
struct EqNo125 <: UnivariateFunction end
struct EqNo130 <: UnivariateFunction end
struct EqNo131 <: UnivariateFunction end
struct EqNo150 <: UnivariateFunction end
struct EqNo200 <: UnivariateFunction end
struct EqNo201 <: UnivariateFunction end
struct EqNo206 <: UnivariateFunction end
struct EqNo207 <: UnivariateFunction end
struct EqNo208 <: UnivariateFunction end
struct EqNo209 <: UnivariateFunction end
struct EqNo210 <: UnivariateFunction end
struct EqNo211 <: UnivariateFunction end
struct EqNo212 <: UnivariateFunction end
struct EqNo213 <: UnivariateFunction end
struct EqNo221 <: UnivariateFunction end
struct EqNo230 <: UnivariateFunction end
struct EqNo231 <: UnivariateFunction end
#==
_EqNo transforms a number to the appropiate UnivariateFunction Type
This is basically the Dict, but at compile time.
Some names are used, to make the VariableProperty Struct more readable.
==#
_EqNo(x::Int) = _EqNo(Val(x))
_EqNo(x::Val{1}) =  ConstantEq() #there was a value 1, that i can interpret as a constant.
_EqNo(x::Val{2}) =  LinearEq()
_EqNo(x::Val{3}) =  CuadraticEq()
_EqNo(x::Val{4}) =  CubicPolyEq()
_EqNo(x::Val{5}) =  EqNo5()
_EqNo(x::Val{6}) =  EqNo6()
_EqNo(x::Val{10}) = AntoineEq()
_EqNo(x::Val{100}) = Poly4Eq()
_EqNo(x::Val{101}) = DIPPREq()
_EqNo(x::Val{102}) = EqNo102()
_EqNo(x::Val{103}) = EqNo103()
_EqNo(x::Val{104}) = EqNo104()
_EqNo(x::Val{105}) = EqNo105()
_EqNo(x::Val{106}) = EqNo106()
_EqNo(x::Val{107}) = EqNo107()
_EqNo(x::Val{114}) = EqNo114()
_EqNo(x::Val{117}) = EqNo117()

_EqNo(x::Val{11}) = EqNo11()
_EqNo(x::Val{12}) = EqNo12()
_EqNo(x::Val{13}) = EqNo13()
_EqNo(x::Val{14}) = EqNo14()
_EqNo(x::Val{15}) = EqNo15()
_EqNo(x::Val{16}) = EqNo16()
_EqNo(x::Val{17}) = EqNo17()
_EqNo(x::Val{18}) = EqNo18()
_EqNo(x::Val{19}) = EqNo19()
_EqNo(x::Val{45}) = EqNo45()
_EqNo(x::Val{75}) = EqNo75()
_EqNo(x::Val{118}) = EqNo118()
_EqNo(x::Val{120}) = EqNo120()
_EqNo(x::Val{121}) = EqNo121()
_EqNo(x::Val{122}) = EqNo122()
_EqNo(x::Val{123}) = EqNo123()
_EqNo(x::Val{124}) = EqNo124()
_EqNo(x::Val{125}) = EqNo125()
_EqNo(x::Val{130}) = EqNo130()
_EqNo(x::Val{131}) = EqNo131()
_EqNo(x::Val{150}) = EqNo150()
_EqNo(x::Val{200}) = EqNo200()
_EqNo(x::Val{201}) = EqNo201()
_EqNo(x::Val{206}) = EqNo206()
_EqNo(x::Val{207}) = EqNo207()
_EqNo(x::Val{208}) = EqNo208()
_EqNo(x::Val{209}) = EqNo209()
_EqNo(x::Val{210}) = EqNo210()
_EqNo(x::Val{211}) = EqNo211()
_EqNo(x::Val{212}) = EqNo212()
_EqNo(x::Val{213}) = EqNo213()
_EqNo(x::Val{221}) = EqNo221()
_EqNo(x::Val{230}) = EqNo230()
_EqNo(x::Val{231}) = EqNo231()
#max_length is used to check if some function hasn't enough parameters.
max_length(x :: T) where T<:UnivariateFunction= 0
max_length(x::ConstantEq) = 1
max_length(x::LinearEq) = 2
max_length(x::CuadraticEq) = 3
max_length(x::CubicPolyEq) = 4
max_length(x::AntoineEq) = 3
max_length(x::Poly4Eq) = 5
max_length(x::DIPPREq) = 5
max_length(x::EqNo102) = 4
max_length(x::EqNo103) = 4
max_length(x::EqNo104) = 5
max_length(x::EqNo105) = 4
max_length(x::EqNo106) = 5
max_length(x::EqNo107) = 5
max_length(x::EqNo114) = 4
max_length(x::EqNo117) = 5

#struct to hold the variable property,can hold parameters, or Tc where necessary.
#with those types, its is fully stack-allocated. if passed to a thermodynamic routine,
#the evaluation will be fast.
struct VariableProperty{F <: UnivariateFunction ,T <: Real,N}
    f::F
    p::NTuple{N,T}
    valid_range::NTuple{2,T}
end

#Variable property constructor.
#it does some analisis in some cases, where a function is used, but all parameters are zero,
#efectively making that function constant.
function VariableProperty(eqno::Int,p,valid_range) where _FF <: UnivariateFunction
    _T = promote_type(eltype(p),eltype(valid_range))
    p0 = first(p)
    p_nonzero = filter(!iszero,[i for i in p])
    _F = ConstantEq()
    _N = 0
    #check if property is just a constant,then apply constant function (f1)
    if (length(p_nonzero) == 1) && p0 == first(p_nonzero) && (eqno != 1) && (eqno in (2,3,4,100,103,104,107))
        _F = ConstantEq()
        _N = 1
    else
        _F = _EqNo(eqno)
        _N = length(p)
    end

    #check if length of respective type is correct
    if max_length(_F) > length(p)
        throw(ArgumentError("not enough parameters for $eqno property"))
    end
    P = NTuple{_N,_T}(convert(_T,i) for i in p)
    R = NTuple{2,_T}(convert(_T,i) for i in valid_range)
return VariableProperty{typeof(_F),_T,_N}(_F,P,R)
end

#@inline f11(T,p) = p[1] p[2] p[3] p[4] p[5]
@inline (prop::EqNo11)(T,p) = exp(p[1])
@inline (prop::EqNo12)(T,p) = exp(p[1]+p[2]*T)
@inline (prop::EqNo13)(T,p) = exp(p[1]+p[2]*T+p[3]*T^2)
@inline (prop::EqNo14)(T,p) = exp(p[1]+p[2]*T+p[3]*T^2+p[4]*T^3)
@inline (prop::EqNo15)(T,p) = exp(p[1]+p[2]*T+p[3]*T^2+p[4]*T^3+p[5]*T^4)
@inline (prop::EqNo16)(T,p) = p[1] + exp(p[2]/T + p[3] + p[4]*T +p[5]*T^2)
@inline (prop::EqNo17)(T,p) = p[1] + exp(p[2] + p[3]*T + p[4]*T^2 +p[5]*T^3)
@inline (prop::EqNo18)(T,p) = p[1] + p[2]*(1+log(T)*(1+p[3]/T))*exp(-p[3]/T)
@inline (prop::EqNo19)(T,p) = p[1] + p[2]*T*log(T)*exp(-p[3]/T)
@inline (prop::EqNo45)(T,p) = p[1]*T + (p[2]*T^2)/2 + (p[3]*T^3)/3 + (p[4]*T^4)/4 + (p[5]*T^5)/5
@inline (prop::EqNo75)(T,p) = p[2]*(2*p[3]*T + 3*p[4]*T^2 + 4*p[5]*T^3) #weird too, check
@inline (prop::EqNo118)(T,p) = exp(p[1] + p[2]/(T^p[5]) + p[3]*log(T) + p[4]*T^2) #weird too check
@inline (prop::EqNo120)(T,p) = p[1] - p[2]/(p[3] + T) #LogAntoine
@inline (prop::EqNo121)(T,p) =  p[1] + p[2]/T + p[3]*log(T) + p[4]*T^p[5]
@inline (prop::EqNo122)(T,p) =  p[1] + p[2]/T + p[3]*log(T) + p[4]*T^2 + p[5]/(T^2)
@inline (prop::EqNo123)(T,p) =  p[1] + p[2]/T + p[3]*T + p[4]*T^2 + p[5]*(T^3)
@inline (prop::EqNo124)(T,p) =  p[1] + p[2]/T + p[3]/(T^2) + p[4]*T + p[5]*(T^2)
@inline (prop::EqNo125)(T,p) =  exp(p[1] + p[2]/T + p[3]/(T^2) + p[4]*T + p[5]*(T^2))
@inline (prop::EqNo130)(T,p) =  exp(p[1] + p[2]/T + p[3]*log(T) + p[4]*T + p[5]/(T^2))
@inline (prop::EqNo131)(T,p) =  p[1] + p[2]*T + p[3]/(T-p[4])
@inline (prop::EqNo150)(T,p) = p[1]+p[2]*T+p[3]*T^2+p[4]*T^3 + p[5]/(T^2)

@inline function (prop::EqNo200)(Tr,p)
    t = 1-Tr
    return p[5]*exp((p[1]*t + p[2]*t^1.5 + p[3]*t^2.5 + p[4]*t^5)/Tr)
end

@inline function (prop::EqNo201)(Tr,p)
    t = 1-Tr
    return p[5]+ exp((p[1]*t + p[2]*t^1.5 + p[3]*t^2.5 + p[4]*t^5)/Tr)
end

@inline function (prop::EqNo206)(Tr,p)
    t =Tr
    return p[5]+ exp(p[1]+p[2]/(p[3]+t+p[4]*t^2 + p[5]*t^3))
end
@inline (prop::EqNo207)(T,p) = exp(p[1] - p[2]/(p[3] + T))
@inline (prop::EqNo208)(T,p) = exp10(p[1] - p[2]/(p[3] + T))
@inline (prop::EqNo209)(T,p) = exp10(p[1]*(1/T - 1/p[2]))
@inline (prop::EqNo210)(T,p) = exp10(p[1] + p[2]/T + p[3]*T + p[4]*T^2)
@inline (prop::EqNo211)(T,p) = p[1]*((p[2]-T)/(p[2]-p[3]))^p[4]
@inline function (prop::EqNo212)(T,p)
    a = (1-T/p[5])
    pol = p[1]*a + p[2]*a^1.5 + p[3]*a^3 + p[4]^6
    return exp((p[5]*T)*pol)
end

@inline function (prop::EqNo213)(T,p)
    t = T - 1
    pol = 1 + p[3]*t + p[4]*t^2 + p[5]*t^3
    return (p[1]-p[2])/pol
end
@inline (prop::EqNo221)(T,p) = -p[2]/(T^2) + p[3]/T + p[4]*p[5]*T^(p[5]-1)
@inline (prop::EqNo230)(T,p) = -p[2]/(T^2) + p[3]/T + p[4] - 2*p[5]/T^3
@inline (prop::EqNo231)(T,p) = p[2] - p[3]/(T-p[4])^2
#functor,calling the property like a function.
#==as these function calls are cheap,
the function returns for now the value and a bool indicating if thats a valid range.
an alternative would be adding a keyword to allow the evaluation of the property outside the range
==#
@inline function (prop::VariableProperty)(T)
    low,high = prop.valid_range
    if low <= T <= high
        return prop.f(T,prop.p),true
    else
        return prop.f(T,prop.p),false
    end
end

#dispatch to call the f functions from UnivariateFunction.
#technically you can write the implementation right here
@inline (prop::ConstantEq)(T,p) = p[1]
@inline (prop::LinearEq)(T,p) = p[1] + p[2]*T
@inline (prop::CuadraticEq)(T,p) = p[1] + p[2]*T + p[3]*T^2
@inline (prop::CubicPolyEq)(T,p) =  p[1] + p[2]*T + p[3]*T^2 + p[4]*T^3
@inline (prop::EqNo5)(T,p) = p[1] + p[2]*T + p[3]*T^2 + p[4]*T^3 + p[5]*T^4
@inline (prop::EqNo6)(T,p) =  p[1] + p[2]*T + p[3]*T^2 + p[4]*T^3 + p[5]/T^2
@inline (prop::AntoineEq)(T,p) =  exp(p[1] - p[2]/(p[3] + T))
@inline (prop::Poly4Eq)(T,p)  =  p[1] + p[2]*T + p[3]*T^2 + p[4]*T^3 + p[5]*T^4
@inline (prop::DIPPREq)(T,p) =  exp(p[1] + p[2]/T + p[3]*log(T) + p[4]*T^p[5])
@inline (prop::EqNo102)(T,p) =  p[1]*T^p[2]/(1 + p[3]/T + p[4]/T^2)
@inline (prop::EqNo103)(T,p) =  p[1] + p[2]*exp(-p[3]/(T^p[4]))
@inline (prop::EqNo104)(T,p) =   p[1] + p[2]/T + 1E6*p[3]/T^3 + 1E16*p[4]/T^8 + 1E18*p[5]/T^9
@inline (prop::EqNo105)(T,p) =  p[1]/(p[2]^(1 + (1 - T/p[3])^p[4]))
@inline (prop::EqNo106)(Tr,p) = p[1]*(1-Tr)^(p[2]+p[3]*Tr+p[4]*Tr^2+p[5]*Tr^3)
@inline (prop::EqNo107)(T,p) =  p[1] + p[2]*((p[3]/T)/sinh(p[3]/T))^2 + p[5]*((p[4]/T)*cosh(p[4]/T))^2
@inline (prop::EqNo114)(T,p) =  p[1]*T + p[2]*T^2/2 + p[3]*T^3/3 + p[4]*T^4/4
@inline (prop::EqNo117)(T,p) =  p[1]*T + p[2]*(p[3]/T)/tanh(p[3]/T) - p[4]*(p[5]/T)/tanh(p[5]/T)

# Generic function to call the right temperature dependency
#function calcvectorproperty(T, prop::VectorProperty)
#    tfunc = tfuncs[prop.EqNo]
#    tfunc(T, prop.A, prop.B, prop.C, prop.D, prop.E)
#end


text:
## Table with dynamic scope for iteration variable.
## using a[2] for iteration variable still does not work.
##
## We need to reimplement tricks like writing the value
## of the iteration variable into a copy of the expression (or restoring it).
## In this case, we still have to localize the variable with dynamic scope so
## it has the expected (possible) side-effects.
##
## Also, this is copied from older code. There are probably no inefficiencies introduced.
## But, it is more complicated than it needs to be

abstract type TableEval end

struct TableSetEval{T,V} <: TableEval
    expr::T
    sym::V
end

struct TableNoSetEval{T} <: TableEval
    expr::T
end

@mkapprule Table

@sjdoc Table """
    Table(expr,[imax])

return a `List` of `imax` copies of `expr`.

    Table(expr,[i,imax])

returns a `List` of `expr` evaluated `imax` times with `i` set successively to `1` through `imax`.

    Table(expr,iter)

`iter` can be any standard iterator.

    Table(expr,iter1,iter2,...)

is equivalent to `Table(Table(expr,iter2),iter1)...`

This example calls an anonymous Julia function
```
f = :( g(x) = x^2 )
Table( f(i), [i,10])
```
"""

set_attribute(:Table, :HoldAll)

@doap Table(x) = mx

@doap function Table(inexpr, iters...)
    if length(iters) > 1
        riters = reverse(iters)
        mxout = mxpr(:Table,inexpr,riters[1])
        for i in 2:length(riters)
            mxout = mxpr(:Table,mxout,riters[i])
        end
        return mxout
    end
    expr = is_Mxpr(inexpr,:Jxpr) ? eval(inexpr) : inexpr
    it = iters[1]
    sjiter = make_sjiter(it)
    dosaveval::Bool = true
    if typeof(sjiter) <: SJIter1
        dosaveval = false
    end
    local saveval
    if dosaveval
        saveval = symval(it[1])
    end
    local t = NullMxpr
    try          # at the moment it doesn't matter much because we can't interrupt anyway..., could be errors though
        t = table(expr, sjiter)
    catch
        symerror("error in Table")  # or interrupt
    finally
        if dosaveval
            setsymval(it[1], saveval)
        end
    end
    t
end

function table(expr, sjiter)
    mxpr(:List, intable(expr, sjiter))
end

function intable(expr, iter::SJIter1)
    tableloop1(expr, iter.imax)
end

function intable(expr, iter::SJIter2)
    (sym, nexpr) = (iter.i,expr)
    tableloop2(sym, nexpr, iter.imax)
end

function intable(expr, iter::SJIterList)
    (sym, nexpr) = (iter.i,expr)
    tableloop_list(sym, nexpr, iter.list)
end

function evaltable(model::TableSetEval, val)
    setandeval(model.expr, model.sym, val)
end

function evaltable(model::TableNoSetEval)
    nosetandeval(model.expr)
end

function setandeval(expr,sym,val)
    setsymval(sym,val)
    unsetfixed(expr)
    r = doeval(expr)
    setfixed(r) # trying to gain a little efficiency
    setcanon(r) # usually fails in some more or less edge case
    r
end

function nosetandeval(expr)
    unsetfixed(expr)
    r = doeval(expr)
    setfixed(r) # trying to gain a little efficiency
    setcanon(r) # usually fails in some more or less edge case
    r
end

function intable(expr0,iter::SJIter3{T,V}) where {T<:Real,V<:Real}
    (sym, expr) = (iter.i,expr0)
    imax = doeval(iter.imax) # maybe this should be done earlier. When iter is created ?
    imin = doeval(iter.imin)
    args = newargs(imax-imin+1)
    if iter.i == expr
        @inbounds for i in imin:imax
            args[i-imin+1] = i
        end
    else
        model = TableSetEval(expr,sym)
        _table3(model,imin,imax,args)
    end
    return args
end

function _table3(model,imin,imax,args)
    for i in imin:imax
        args[i-imin+1] = evaltable(model,i)
    end
end

function intable(expr0, iter::SJIter4{T,V,W}) where {T<:Real, V<:Real, W<:Real}
    (sym, expr) = (iter.i,expr0)
    imax = iter.imax  # why do we not evaluate here
    imin = iter.imin
    di = iter.di
    args = newargs(iter.num_iters)
    j::Int = 0
    if iter.i == expr
        @inbounds for i in imin:di:imax
            j += 1
            args[j] = i
        end
    else
        model = TableSetEval(expr,sym)
        @inbounds for i in imin:di:imax
            j += 1
            args[j] = evaltable(model,i)
        end
    end
    return args
end

# Symbolic value and increment
function intable(expr0, iter::SJIter4{T,V,W}) where {T, V, W}
    (sym, expr) = (iter.i,expr0)
    imax = iter.imax
    imin = iter.imin
    di = iter.di
    num_iters = iter.num_iters
    args = newargs(iter.num_iters)
    model = TableSetEval(expr,sym)
    val = deepcopy(imin)
    @inbounds for i in 1:num_iters
        args[i] = evaltable(model,val)
        val = doeval(mxpr(:Plus, val, di))
    end
    args
end

function tableloop1(expr, imax)
    args = newargs(imax)
    model = TableNoSetEval(expr)
    @inbounds  for i=1:imax
        args[i] = evaltable(model)
    end
    args
end

function tableloop2(sym, expr, imax)
    args = newargs(imax)
    model = TableSetEval(expr,sym)
    @inbounds  for i=1:imax
        args[i] = evaltable(model,i)
    end
    args
end

function table3(expr,iter::SJIter3{T,V}) where {T<:Real,V<:Real}
    imax = doeval(iter.imax) # maybe this should be done earlier. When iter is created ?
    imin = doeval(iter.imin)
    args = newargs(imax-imin+1)
    if iter.i == expr
        @inbounds for i in imin:imax
            args[i-imin+1] = i
        end
    else
        model = TableSetEval(expr,sym)
        @inbounds for i in imin:imax
            args[i-imin+1] =  evaltable(model,i)
        end
    end
    return args
end

function tableloop_list(sym, expr, list)
    len = length(list)
    args = newargs(len)
    model = TableSetEval(expr,sym)
    for i in 1:len
        args[i] = evaltable(model,list[i])
    end
    args
end


text:
"""Functions to run firing and hebbian plasticity updates
on the `AssemblyCalculus` types defined in `types.jl`

Author: DJ Passey
Email: djpasseyjr@gmail.com
"""

## SIMULATION WRAPPERS

function simulate!(
    stims::Array{Stimulus{T}, 1},
    assemblies::Array{Assembly{T}, 1},
    stop_criteria::StopCriteria;
    random_initial::Bool = true
) where T
    # TODO Validate stims and assemblies
    active_areas = [inp.area for inp in stims]
    # Consolidate stim current with current from assemblies
    currents, assem_attrib = simulation_currents(stims, assemblies)
    # Save initial stim currents
    init_currents = deepcopy(currents)
    # Assign random neurons to fire initially
    random_initial && map(random_firing!, active_areas)
    # Allocate firing and stim current arrays for each area
    firing_and_currents = deepcopy(currents)
    # stop_criteria checks for convergence
    while ! stop_criteria(active_areas, currents)
        # All neurons fire and modulate firing_and_currents
        attrib = fire!(active_areas, firing_and_currents)
        # Compute k winners based on firing and stimulus current vectors
        map(winners!, active_areas, firing_and_currents)
        # Hebbian update for synapses
        map(area -> hebb_update!(area, attrib), active_areas)
        # Hebbian update for currents
        map(hebb_update!, active_areas, currents)
        # Reset firing_and_currents to exclude firing data
        firing_and_currents = deepcopy(currents)
    end
    # Hebbian update to assembly synapses
    map(x -> hebb_update!(assem_attrib, x...), zip(stims, currents, init_currents))
    final_assemblies = map(x -> Assembly(x, stims, assemblies), active_areas)
    return final_assemblies
end

function simulate!(
    stims::Array{Stimulus{T}, 1},
    assemblies::Array{Assembly{T}, 1},
    timesteps::Int;
    random_initial::Bool = true
) where T
    stop_criteria = MaxIters(max_iters=timesteps, record=true)
    new_assemblies = simulate!(stims, assemblies, stop_criteria, random_initial=random_initial)
    return new_assemblies, stop_criteria.spikes
end

function simulate!(
    stims::Array{Stimulus{T}, 1},
    timesteps::Int;
    random_initial::Bool = true
) where T
    return simulate!(stims, Assembly{T}[], timesteps; random_initial=random_initial)
end

function simulate!(
    stims::Array{Stimulus{T}, 1},
    assemblies::Array{Assembly{T}, 1};
    record = true,
    tol::Int = 2,
    fire_past_convergence::Int = 0,
    max_iters::Int = 50,
    random_initial::Bool = true
) where T
    stop_criteria = Converge(
        tol=tol,
        fire_past_convergence=fire_past_convergence,
        max_iters=max_iters,
        record=record
    )
    new_assemblies = simulate!(stims, assemblies, stop_criteria, random_initial=random_initial)
    if record
        return new_assemblies, stop_criteria.spikes, stop_criteria.distances
    else
        return new_assemblies
    end
end

function simulate!(
    stims::Array{Stimulus{T}, 1};
    record = true,
    tol::Int = 2,
    fire_past_convergence::Int = 0,
    max_iters::Int = 50,
    random_initial::Bool = true
) where T
    return simulate!(
                stims,
                Assembly{T}[],
                record=record,
                tol=tol,
                fire_past_convergence=fire_past_convergence,
                max_iters=max_iters,
                random_initial=random_initial
            )
end
"""
    simulation_currents(
        stims::Array{Stimulus{T}, 1},
        assemblies::Array{Assembly{T}, 1}
    ) where T -> currents, attribution

Computes the ionic current flow into each area and creates an attribution
graph to be used for ion current and assembly plasticity updates.

**Parameters**

1. `stims` The list of IonCurrents to be used in the simulation
2. `assemblies` The list of active Assemblies for the simulation

"""
function simulation_currents(
    stims::Array{Stimulus{T}, 1},
    assemblies::Array{Assembly{T}, 1} = Assembly{T}[]
) where T
    # Allocate
    currents = [zeros(T, size(inp.currents)) for inp in stims]
    attrib = neuron_attrib_graph(T, typeof(stims[1].area))
    # Compute currents and attributions
    for (i, inp) in enumerate(stims)
        currents[i] += inp.currents
        target_area = inp.area
        # Fire each assembly into the target area
        for assem in assemblies
            fire!(assem, target_area, currents[i], attrib)
        end
    end
    return currents, attrib
end

## FIRE! FUNCTIONS

"""
    fire!(
        assembly::Assembly{T},
        area::BrainArea{T},
        currents::Array{T, 1},
        attrib::AttributionGraph{Neuron{T, U}}
    ) where {T, U}

Fires each neuron in `assembly` into the target area. Adds resulting ion
currents to corresponding neuron indexes in `current` array.
"""
function fire!(
    assembly::Assembly{T},
    area::BrainArea{T},
    currents::Array{T, 1},
    attrib::AttributionGraph{Neuron{T, U}}
) where {T, U}
    # Fire each neuron in the assembly into the target area
    for idx in assembly.neurons
        neuron = assembly.area.neurons[idx]
        fire!(neuron, area, currents, attrib)
    end
end

"""
    fire!(
        source_area::BrainArea{T},
        target_area::BrainArea{T},
        target_currents::Array{T, 1},
        attrib::AttributionGraph{Neuron{T, U}}
    ) where {T, U}

Fires each neuron in `source_area` into `target_area`. Adds resulting ion
currents to corresponding neuron indexes in `currents` array. Collects
neurons into the attribution graph.
"""
function fire!(
    source_area::BrainArea{T},
    target_area::BrainArea{T},
    target_currents::Array{T, 1},
    attrib::AttributionGraph{Neuron{T, U}}
) where {T, U}
    # Fire each neuron in the assembly into the target area
    for idx in source_area.firing
        neuron = source_area.neurons[idx]
        fire!(neuron, target_area, target_currents, attrib)
    end
end

"""
    fire!(
        neuron::Neuron{T, U},
        target_area::BrainArea{T},
        target_currents::Array{T, 1},
        attrib::AttributionGraph{Neuron{T, U}}
    ) where {T, U}

Fires neuron into the target area. Adds resulting ion
currents to corresponding neuron indexes in `currents` array. Adds neruon
to the `AttributionGraph`.
"""
function fire!(
    neuron::Neuron{T, U},
    target_area::BrainArea{T},
    target_currents::Array{T, 1},
    attrib::AttributionGraph{Neuron{T, U}}
) where {T, U}
    # For each synapse into the target area
    for (target_neuron, w) in get(neuron, target_area, Dict{Neuron{T}, T}())
        # Add ion current to the current array
        target_currents[target_neuron.idx] += w
        # Update attribution graph entry for target neuron
        contributors = get!(attrib, target_neuron)
        push!(contributors, neuron)
    end
end

"""
    fire!(
        areas::Array{<:BrainArea{T}, 1},
        area_currents::Array{Array{T, 1}}
    ) where T -> attrib

Fires all active neurons in each area. Collects ion currents in `currents`.
Returns an `AttributionGraph` of the firing.
"""

function fire!(areas::Array{U, 1}, area_currents::Array{Array{T, 1}}) where {T, U <:BrainArea{T}}
    attrib = neuron_attrib_graph(T, U)
    for source_area in areas
        for (target_area, target_currents) in zip(areas, area_currents)
            fire!(source_area, target_area, target_currents, attrib)
        end
    end
    return attrib
end

## HEBB_UPDATE! FUNCTIONS

"""
    hebb_update!(
        neuron::Neuron{T, U},
        attrib::AttributionGraph{Neuron{T, U}};
        pow::Int=1
    ) where {T, U}

Uses the attribution graph to find which synapses fired on to `neuron`
in the previous timestep. Increases these synapses multiplicatively according
to the plasticity of their parent neuron's area. The parameter `pow` controls
how many multiplicative updates to apply.
"""
function hebb_update!(
    neuron::Neuron{T, U},
    attrib::AttributionGraph{Neuron{T, U}};
    pow::T = T(1)
) where {T, U}
    source_neurons = get(attrib, neuron, Neuron{T, U}[])
    if length(source_neurons) != 0
        for source_neuron in source_neurons
            β = source_neuron.area.plasticity
            w = source_neuron[neuron.area][neuron]
            source_neuron[neuron.area][neuron] = w * (1 + β)^pow
        end
    end
end

"""
    hebb_update!(area::BrainArea{T}, attrib::AttributionGraph{Neuron{T, U}})

Uses the attribution graph to update the synapses pointing to all currently
firing neurons in `area`.
"""
function hebb_update!(area::BrainArea{T}, attrib::AttributionGraph{Neuron{T, U}}) where {T, U}
    map(n -> hebb_update!(n, attrib), area.neurons[area.firing])
end

"""
    hebb_update!(area::BrainArea{T}, current::Array{T, 1}) where T

Multiplicatively increases the current into currently firing neurons in `area`
by a factor of `(1 + area.plasticity)`.
"""
function hebb_update!(area::BrainArea{T}, current::Array{T, 1}) where T
    current[area.firing] .*= (1 + area.plasticity)
end

"""
    hebb_update!(
        attrib::AttributionGraph{Neuron{T, U}},
        stim::Stimulus{T},
        currents::Array{T, 1},
        init_currents::Array{T, 1}
    ) where {T, U}

Updates the assembly synapse weights all at once after the simulation is
finished.

Since the contribution of these synapses is included in the simulation current,
we can update the synapses once the simulation is finished based on changes in
the current array due to plasticity. This is faster than updating them each
timestep.
"""
function hebb_update!(
    attrib::AttributionGraph{Neuron{T, U}},
    ic::Stimulus{T},
    currents::Array{T, 1},
    init_currents::Array{T, 1},
) where {T, U}
    times_grown = _find_growth(currents, init_currents, ic.area.plasticity)
    # Update assembly synapses
    for (i, pow) in enumerate(times_grown)
        if pow > 0
            ic.currents[i] *= (1 + ic.area.plasticity) ^ pow
            neuron = ic.area.neurons[i]
            hebb_update!(neuron, attrib, pow=pow)
        end
    end
end

"""Computes the number of multiplicative updates applied to the ion
currents during the simulation.
"""
function _find_growth(
    currents::Array{T, 1},
    init_currents::Array{T, 1},
    plasticity::T
) where T
    # Find nonzero ion currents
    nonzero = findall(currents .!= T(0))
    nonzero == nothing && return zeros(T, size(currents))
    # Compute multiplicative change from the initial current
    change = currents[nonzero] ./ init_currents[nonzero]
    # Compute how many multiplicative updates occured given the plasticity
    nz_times_grown = log.(change) ./ log(1 + plasticity)
    times_grown = zeros(T, size(currents))
    times_grown[nonzero] .= nz_times_grown
    return times_grown
end


text:
# Syntactic sugar for arrow syntax
# we need version ≥ 1.3 in order to make use of multiple dispatch
# over abstract types

# This allows implicit: data |> machine
(mach::Machine{<:Unsupervised})(data) = transform(mach, data)
(mach::Machine{<:Unsupervised})(data::NTuple{<:Any,AbstractNode}) =
    transform(mach, data...)
(mach::Machine{<:Supervised})(data)   = predict(mach, data)
(mach::Machine)(data::AbstractMatrix) = data |> table |> mach

# This allows implicit: data |> Unsupervised
(m::Unsupervised)(data::AbstractNode) = data |> machine(m, data)
(m::Static)(data:: AbstractNode) = data |> machine(m)
(m::Static)(data::NTuple{<:Any, AbstractNode}) = data |> machine(m)
(m::Unsupervised)(data) = source(data) |> m
(m::Unsupervised)(data::AbstractMatrix) = data |> table |> m

# This allows implicit: data |> Supervised
(m::Supervised)(data::NTuple{2,AbstractNode}) = data[1] |> machine(m, data...)
(m::Supervised)(data::Tuple{AbstractNode,Any}) =
    (data[1], source(data[2])) |> m
(m::Supervised)(data::Tuple) = (source(data[1]), data[2]) |> m
(m::Supervised)(data::Tuple{AbstractMatrix,Any}) =
    (data[1] |> table, data[2]) |> m

# This allows implicit: data |> inverse_transform(node)
inverse_transform(node::Node{<:Machine{<:Unsupervised}}) =
    data -> inverse_transform(node.machine, data)


text:
get_all_plots_types() = Set([:fit, :residuals, :normal_checks, :cooksd, :leverage, :homoscedasticity])
get_needed_plots(s::String) = return get_needed_plots([s])
get_needed_plots(s::Symbol) = return get_needed_plots([s])
get_needed_plots(s::Vector{String}) = return get_needed_plots(Symbol.(lowercase.(s)))
get_needed_plots(::Vector{Any}) = return get_needed_plots([:none])
get_needed_plots(::Set{Any}) = return get_needed_plots([:none])
get_needed_plots(s::Set{Symbol}) = return get_needed_plots(collect(s))

function get_needed_plots(p::Vector{Symbol})

    needed_plots = Set{Symbol}()
    length(p) == 0 && return needed_plots
    :none in p && return needed_plots

    if :all in p
        return get_all_plots_types()
    end

    :fit in p && push!(needed_plots, :fit)
    :residuals in p && push!(needed_plots, :residuals)
    :normal_checks in p && push!(needed_plots, :normal_checks)
    :cooksd in p && push!(needed_plots, :cooksd)
    :leverage in p && push!(needed_plots, :leverage)
    :homoscedasticity in p && push!(needed_plots, :homoscedasticity)
    return needed_plots

end

"""
    function get_robust_cov_stats()

    Return all robust covariance estimators.
"""
get_all_robust_cov_stats() = Set([:white, :nw, :hc0, :hc1, :hc2, :hc3])
get_needed_robust_cov_stats(s::String) = return get_needed_robust_cov_stats([s])
get_needed_robust_cov_stats(s::Symbol) = return get_needed_robust_cov_stats([s])
get_needed_robust_cov_stats(s::Vector{String}) = return get_needed_robust_cov_stats(Symbol.(lowercase.(s)))
get_needed_robust_cov_stats(::Vector{Any}) = return get_needed_robust_cov_stats([:none])
get_needed_robust_cov_stats(::Set{Any}) = return get_needed_robust_cov_stats(Set([:none]))
get_needed_robust_cov_stats(s::Set{Symbol}) = return get_needed_robust_cov_stats(collect(s))

function get_needed_robust_cov_stats(s::Vector{Symbol})

    needed_white = Vector{Symbol}()
    needed_hac = Vector{Symbol}()

    length(s) == 0 && return (needed_white, needed_hac)
    :none in s && return (needed_white, needed_hac)
    if :all in s
        s = collect(get_all_robust_cov_stats())
    end

    :white in s && push!(needed_white, :white)
    :hc0 in s && push!(needed_white, :hc0)
    :hc1 in s && push!(needed_white, :hc1)
    :hc2 in s && push!(needed_white, :hc2)
    :hc3 in s && push!(needed_white, :hc3)

    :nw in s && push!(needed_hac, :nw)

    return (needed_white, needed_hac)

end

"""
    function get_all_model_stats()

    Returns all statistics availble for the fitted model.
"""
get_all_model_stats() = Set([:coefs, :sse, :mse, :sst, :rmse, :aic, :sigma, :t_statistic, :vif, :r2, :adjr2, :stderror, :t_values, :p_values, :ci,
                            :diag_normality, :diag_ks, :diag_ad, :diag_jb, :diag_heteroskedasticity, :diag_white, :diag_bp, :press,
                            :t1ss, :t2ss, :pcorr1, :pcorr2 , :scorr1, :scorr2 ])

get_needed_model_stats(req_stats::String) = return get_needed_model_stats([req_stats])
get_needed_model_stats(req_stats::Symbol) = return get_needed_model_stats(Set([req_stats]))
get_needed_model_stats(req_stats::Vector{String}) = return get_needed_model_stats(Symbol.(lowercase.(req_stats)))
get_needed_model_stats(::Vector{Any}) = return get_needed_model_stats([:none])
get_needed_model_stats(::Set{Any}) = return get_needed_model_stats(Set([:none]))
get_needed_model_stats(req_stats::Set{Symbol}) = get_needed_model_stats(collect(req_stats))

"""
    function get_needed_model_stats(req_stats::Vector{Symbol})

    return the list of needed statistics given the list of statistics about the model the caller wants.
"""
function get_needed_model_stats(req_stats::Vector{Symbol})
    needed = Set([:coefs, :sse, :mse])
    default = Set([:coefs, :sse, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    full = get_all_model_stats()
    unique!(req_stats)

    length(req_stats) == 0 && return needed
    :none in req_stats && return needed
    :all in req_stats && return full

    :default in req_stats && union!(needed, default)

    :sst in req_stats && push!(needed, :sst)
    :t1ss in req_stats && push!(needed, :t1ss)
    :t2ss in req_stats && push!(needed, :t2ss)
    :press in req_stats && push!(needed, :press)
    :rmse in req_stats && push!(needed, :rmse)
    :aic in req_stats && push!(needed, :aic)
    :sigma in req_stats && push!(needed, :sigma)
    :t_statistic in req_stats && push!(needed, :t_statistic)
    :vif in req_stats && push!(needed, :vif)
    :diag_ks in req_stats && push!(needed, :diag_ks)
    :diag_ad in req_stats && push!(needed, :diag_ad)
    :diag_jb in req_stats && push!(needed, :diag_jb)
    :diag_white in req_stats && push!(needed, :diag_white)
    :diag_bp in req_stats && push!(needed, :diag_bp)

    if :diag_normality in req_stats
        push!(needed, :diag_ks)
        push!(needed, :diag_ad)
        push!(needed, :diag_jb)
    end
    if :diag_heteroskedasticity in req_stats
        push!(needed, :diag_white)
        push!(needed, :diag_bp)
    end
    if :pcorr1 in req_stats
        push!(needed, :t1ss)
        push!(needed, :pcorr1)
    end
    if :pcorr2 in req_stats
        push!(needed, :t2ss)
        push!(needed, :pcorr2)
    end
    if :scorr1 in req_stats
        push!(needed, :sst)
        push!(needed, :t1ss)
        push!(needed, :scorr1)
    end
    if :scorr2 in req_stats
        push!(needed, :sst)
        push!(needed, :t2ss)
        push!(needed, :scorr2)
    end
    if :r2 in req_stats
        push!(needed, :sst)
        push!(needed, :r2)
    end
    if :adjr2 in req_stats
        push!(needed, :sst)
        push!(needed, :r2)
        push!(needed, :adjr2)
    end
    if :stderror in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
    end
    if :t_values in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_values)
    end
    if :p_values in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_values)
        push!(needed, :p_values)
    end
    if :ci in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_statistic)
        push!(needed, :ci)
    end

    return needed
end

"""
    function get_all_prediction_stats()

    get all the available statistics about the values predicted by a fitted model
"""
get_all_prediction_stats() = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])

get_prediction_stats(req_stats::String) = return get_prediction_stats([req_stats])
get_prediction_stats(req_stats::Vector{String}) = return get_prediction_stats(Symbol.(lowercase.(req_stats)))
get_prediction_stats(::Vector{Any}) = return get_prediction_stats([:none])
get_prediction_stats(::Set{Any}) = return get_prediction_stats(Set([:none]))
get_prediction_stats(req_stats::Set{Symbol}) = return get_prediction_stats(collect(req_stats))

"""
    function get_prediction_stats(req_stats::Vector{Symbol})

    return the list of needed statistics and the statistics that need to be presentd given the list of statistics about the predictions the caller wants.
"""
function get_prediction_stats(req_stats::Vector{Symbol})
    needed = Set([:predicted])
    full = get_all_prediction_stats()
    present = Set([:predicted])
    unique!(req_stats)

    length(req_stats) == 0 && return needed, present
    :none in req_stats && return needed, present
    :all in req_stats && return full, full

    :leverage in req_stats && push!(present, :leverage)
    :residuals in req_stats && push!(present, :residuals)
    if :stdp in req_stats
        push!(needed, :leverage)
        push!(present, :stdp)
    end
    if :stdi in req_stats
        push!(needed, :leverage)
        push!(present, :stdi)
    end
    if :stdr in req_stats
        push!(needed, :leverage)
        push!(present, :stdr)
    end
    if :student in req_stats
        push!(needed, :leverage)
        push!(needed, :residuals)
        push!(needed, :stdr)
        push!(present, :student)
    end
    if :rstudent in req_stats
        push!(needed, :leverage)
        push!(needed, :residuals)
        push!(needed, :stdr)
        push!(needed, :student)
        push!(present, :rstudent)
    end
    if :lcli in req_stats
        push!(needed, :leverage)
        push!(needed, :stdi)
        push!(present, :lcli)
    end
    if :ucli in req_stats
        push!(needed, :leverage)
        push!(needed, :stdi)
        push!(present, :ucli)
    end
    if :lclp in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(present, :lclp)
    end
    if :uclp in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(present, :uclp)
    end
    if :press in req_stats
        push!(needed, :residuals)
        push!(needed, :leverage)
        push!(present, :press)
    end
    if :cooksd in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(needed, :stdr)
        push!(needed, :residuals)
        push!(needed, :student)
        push!(present, :cooksd)
    end

    union!(needed, present)

    return needed, present

end

"""
    function encapsulate_string(s)

    (internal) Only used to encapsulate a string into an array.

    used exclusively to handle the function ```StatsBase.coefnames``` which sometime return an array or when there is only one element the element alone.
"""
function encapsulate_string(s::String)
    return [s]
end

"""
    function encapsulate_string(v)

    (internal) Only used to encapsulate a string into an array.

    used exclusively to handle the function ```StatsBase.coefnames``` which sometime return an array or when there is only one element the element alone.
"""
function encapsulate_string(v::Vector{String})
    return v
end

import Printf
"""
    macro gprintf(fmt::String)

    (internal) used to format with %g

    Taken from message published by user o314 at https://discourse.julialang.org/t/printf-with-variable-format-string/3805/6
"""
macro gprintf(fmt::String)
    :((io::IO, arg) -> Printf.@printf(io, $fmt, arg))
end

"""
    function fmt_pad(s::String, value, pad=0)
    (internal) helper to format and pad string for results display
"""
function fmt_pad(s::String, value, pad=0)
    fmt = @gprintf("%g")
    return rpad(s * sprint(fmt, value), pad)
end

using NamedArrays
"""
    function my_namedarray_print([io::IO = stdout], n::NamedArray)

    (internal) Print the NamedArray without the type annotation (on the first line).
"""
function my_namedarray_print(io::IO, n)
    tmpio = IOBuffer()
    show(tmpio, n)
    println(io, split(String(take!(tmpio)), "\n", limit=2)[2])
end
my_namedarray_print(n::NamedArray) = my_namedarray_print(stdout::IO, n)

"""
    function helper_print_table(io, title, stats::Vector, stats_name::Vector, updformula)

    (Internal) Convenience function to display a table of statistics to the user.
"""
function helper_print_table(io::IO, title, stats::Vector, stats_name::Vector, updformula)
    println(io, "\n$title")
    todelete = [i for (i, v) in enumerate(stats) if isnothing(v)]
    deleteat!(stats, todelete)
    deleteat!(stats_name, todelete)
    m_all_stats = reduce(hcat, stats)
    if m_all_stats isa Vector
        m_all_stats = reshape(m_all_stats, length(m_all_stats), 1)
    end
    na = NamedArray(m_all_stats)
    setnames!(na, encapsulate_string(string.(StatsBase.coefnames(updformula.rhs))), 1)
    setnames!(na, encapsulate_string(string.(stats_name)), 2)
    setdimnames!(na, ("Terms", "Stats"))
    my_namedarray_print(io, na)
end

function present_breusch_pagan_test(X, residuals, α)
    bpt = HypothesisTests.BreuschPaganTest(X, residuals)
    pval = pvalue(bpt)
    alpha_value= round((1 - α)*100, digits=3)
    topresent = string("Breush-Pagan Test (heteroskedasticity of residuals):\n  T*R² statistic: $(round(bpt.lm, sigdigits=6))    degrees of freedom: $(round(bpt.dof, digits=6))    p-value: $(round(pval, digits=6))\n")
    if pval > α
        topresent *= "  with $(alpha_value)% confidence: fail to reject null hyposthesis.\n"
    else
        topresent *= "  with $(alpha_value)% confidence: reject null hyposthesis.\n"
    end
    return topresent
end

function present_white_test(X, residuals, α)
    bpt = HypothesisTests.WhiteTest(X, residuals)
    pval = pvalue(bpt)
    alpha_value= round((1 - α)*100, digits=3)
    topresent = string("White Test (heteroskedasticity of residuals):\n  T*R² statistic: $(round(bpt.lm, sigdigits=6))    degrees of freedom: $(round(bpt.dof, digits=6))    p-value: $(round(pval, digits=6))\n")
    if pval > α
        topresent *= "  with $(alpha_value)% confidence: fail to reject null hyposthesis.\n"
    else
        topresent *= "  with $(alpha_value)% confidence: reject null hyposthesis.\n"
    end
    return topresent
end

function present_kolmogorov_smirnov_test(residuals, α)
    fitted_residuals = fit(Normal, residuals)
    kst = HypothesisTests.ApproximateOneSampleKSTest(residuals, fitted_residuals)
    pval = pvalue(kst)
    KS_stat = sqrt(kst.n)*kst.δ
    alpha_value= round((1 - α)*100, digits=3)
    topresent = string("Kolmogorov-Smirnov test (Normality of residuals):\n  KS statistic: $(round(KS_stat, sigdigits=6))    observations: $(kst.n)    p-value: $(round(pval, digits=6))\n")
    if pval > α
        topresent *= "  with $(alpha_value)% confidence: fail to reject null hyposthesis.\n"
    else
        topresent *= "  with $(alpha_value)% confidence: reject null hyposthesis.\n"
    end
end

function present_anderson_darling_test(residuals, α)
    fitted_residuals = fit(Normal, residuals)
    adt = HypothesisTests.OneSampleADTest(residuals, fitted_residuals)
    pval = pvalue(adt)
    alpha_value= round((1 - α)*100, digits=3)
    topresent = string("Anderson–Darling test (Normality of residuals):\n  A² statistic: $(round(adt.A², digits=6))    observations: $(adt.n)    p-value: $(round(pval, digits=6))\n")
    if pval > α
        topresent *= "  with $(alpha_value)% confidence: fail to reject null hyposthesis.\n"
    else
        topresent *= "  with $(alpha_value)% confidence: reject null hyposthesis.\n"
    end
end

function present_jarque_bera_test(residuals, α)
    jbt = HypothesisTests.JarqueBeraTest(residuals)
    pval = pvalue(jbt)
    alpha_value= round((1 - α)*100, digits=3)
    topresent = string("Jarque-Bera test (Normality of residuals):\n  JB statistic: $(round(jbt.JB, digits=6))    observations: $(jbt.n)    p-value: $(round(pval, digits=6))\n")
    if pval > α
        topresent *= "  with $(alpha_value)% confidence: fail to reject null hyposthesis.\n"
    else
        topresent *= "  with $(alpha_value)% confidence: reject null hyposthesis.\n"
    end
end

function warn_sigma(lm, stat)
    warn_sigma(lm.white_types, lm.hac_types , stat)
end

function warn_sigma(white_needed, hac_needed, stat)
    if length(white_needed) > 0 || length(hac_needed) > 0
        println(io, "The $(stat) statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
    end
end

function real_sqrt(x)
    return @. real(sqrt(complex(x, 0)))
end

isnotintercept(t::AbstractTerm) = t isa InterceptTerm ? false : true
iscontinuousterm(t::AbstractTerm) = t isa ContinuousTerm ? true : false
iscategorical(t::AbstractTerm) = t isa CategoricalTerm ? true : false

function check_cardinality(df::AbstractDataFrame, f, verbose=false)
    cate_terms = [a.sym for a in filter(iscategorical, terms(f.rhs))]
    if length(cate_terms) > 0
        freqt = freqtable(df, cate_terms...)
        if count(i -> (i == 0), freqt) > 0
            println("At least one group of categories have no observation. Use frequency tables to identify which one(s).")
            println(my_namedarray_print(freqt))
        elseif verbose == true
            println("No issue identified.")
        end
    end
end


text:
using Base.Test
using PGM
using Combinatorics

@testset "Utils" begin

@testset "range_rm" begin
    @test PGM.Utils.range_rm(1, 1) == Int[]
    @test PGM.Utils.range_rm(2, 1) == [2]
    @test PGM.Utils.range_rm(2, 2) == [1]
    @test PGM.Utils.range_rm(5, 1) == [2, 3, 4, 5]
    @test PGM.Utils.range_rm(5, 3) == [1, 2, 4, 5]
    @test PGM.Utils.range_rm(5, 5) == [1, 2, 3, 4]

    @test PGM.Utils.range_rm(1, [1]) == Int[]
    @test PGM.Utils.range_rm(2, [1]) == [2]
    @test PGM.Utils.range_rm(2, [2]) == [1]
    @test PGM.Utils.range_rm(2, []) == [1, 2]
    @test PGM.Utils.range_rm(2, [1,2]) == Int[]
    @test PGM.Utils.range_rm(5, [1]) == [2, 3, 4, 5]
    @test PGM.Utils.range_rm(5, [3]) == [1, 2, 4, 5]
    @test PGM.Utils.range_rm(5, [5]) == [1, 2, 3, 4]
    @test PGM.Utils.range_rm(5, [1, 3]) == [2, 4, 5]
    @test PGM.Utils.range_rm(5, [3, 5]) == [1, 2, 4]
    @test PGM.Utils.range_rm(5, [1, 5]) == [2, 3, 4]
    @test PGM.Utils.range_rm(5, 1:4) == [5]
    @test PGM.Utils.range_rm(5, 2:5) == [1]
    @test PGM.Utils.range_rm(5, 1:5) == Int[]
end

@testset "isortperm" begin
    for p in permutations(1:3)
        @test PGM.Utils.isortperm(p + 2) == p
    end

    for p in permutations(1:4)
        @test PGM.Utils.isortperm(p + 2) == p
    end

    for p in permutations(1:5)
        @test PGM.Utils.isortperm(p + 2) == p
    end
end

end  # Utils