text:
function principal_generator_eichler(I::AlgAssRelOrdIdl)
  O = left_order(I)
  @assert ismaximal(O)
  A = algebra(O)
  @assert iseichler(A)

  orders = representatives_of_maximal_orders(O)
  d = discriminant(O)

  # Compute r such that r orders[i] \subseteq orders[j] for all i, j
  r = fmpz(1)
  for i = 1:length(orders)
    for j = 1:length(orders)
      if i == j
        continue
      end

      # Consider orders[i] as an ideal.
      # The basis pseudo-matrix is probably not in HNF, but we don't need this.
      OO = ideal(A, basis_pmatrix(orders[i], copy = false), true)
      r = lcm(r, denominator(OO, orders[j]))
    end
  end
  rd = r*d

  y = integral_coprime_representative(O, I, rd)
  J = I*y

  N = normred(J, O)
  @assert denominator(N, copy = false) == 1 # J should be integral
  N = numerator(N, copy = false)
  OK = order(N)
  R, mR = ray_class_group(OK(1)*OK, ramified_infinite_places(A))
  @assert iszero(mR\N) "Ideal is not principal"

  primes = collect(keys(factor(N)))
  valuations = [ valuation(N, p) for p in primes ]
  w = _norm_equation_valuations_only(O, primes, valuations)
  w = evaluate(w)

  Ow = O*O(w)
  ww = _eichler_find_transforming_unit(J, Ow)
  return w*ww*inv(y)
end

# Requires nr(M) == nr(N) and that nr(N) is coprime to r*d where d is the
# discriminant of O and r in base_ring(O) such that r O_i \subseteq O_j for a
# system of representatives of the maximal orders.
# Returns t in O^\times such that M = Nt
function _eichler_find_transforming_unit_maximal(M::T, N::T) where { T <: Union{ AlgAssAbsOrdIdl, AlgAssRelOrdIdl } }
  A = algebra(M)
  O = left_order(M)
  @assert ismaximal(O)
  # We assume that left_order(N) == O, but testing this would be really expensive

  if M == N
    return one(O)
  end

  F = FieldOracle(A, [ O ])
  p = normred(M, O)
  @assert p == normred(N, O)
  @assert denominator(p) == 1
  p = numerator(p)
  OpO, toOpO = quo(O, p*O, p)
  B, toB = _as_matrix_algebra(OpO)

  I = ideal_from_gens(B, [ toB(toOpO(O(b))) for b in absolute_basis(M) ])
  J = ideal_from_gens(B, [ toB(toOpO(O(b))) for b in absolute_basis(N) ])

  # Compute the image of 1 under the canonical projections O -> O/M respectively O -> O/N
  Fq = base_ring(B)
  vv = mod(one(B), I)
  nonZeroCol = 0
  for i = 1:degree(B)
    for j = 1:degree(B)
      if !iszero(vv[i, j])
        nonZeroCol = j
        break
      end
    end
    if !iszero(nonZeroCol)
      break
    end
  end
  v = matrix(Fq, degree(B), 1, [ vv[i, nonZeroCol] for i = 1:degree(B) ])
  ww = mod(one(B), J)
  nonZeroCol = 0
  for i = 1:degree(B)
    for j = 1:degree(B)
      if !iszero(ww[i, j])
        nonZeroCol = j
        break
      end
    end
    if !iszero(nonZeroCol)
      break
    end
  end
  w = matrix(Fq, degree(B), 1, [ ww[i, nonZeroCol] for i = 1:degree(B) ])

  b = ceil(Int, degree(base_ring(B))*dim(B)*log2(BigInt(characteristic(base_ring(B)))))
  # A minimal set of generators of around b elements should generate B^\times
  units = _find_some_units(F, 1, b)
  push!(units, -one(O))
  generators = [ matrix(toB(toOpO(u)), copy = false) for u in units ]
  path_exists, path = find_path(generators, v, w)

  while !path_exists
    # Add 10 more generators
    new_units = _find_some_units(F, 1, 10)
    append!(units, new_units)
    append!(generators, [ matrix(toB(toOpO(u)), copy = false) for u in new_units ])
    path_exists, path = find_path(generators, v, w)
  end

  t = one(O)
  for i = 1:length(path)
    u = units[path[i]]
    t = u*t
  end
  @assert matrix(toB(toOpO(t)))*v == w

  return t
end

# Finds at least n units in the order F.maximal_orders[order_num]
function _find_some_units(F::FieldOracle{S, T, U, M}, order_num::Int, n::Int) where { S <: AbsAlgAss{nf_elem}, T, U, M }
  O = F.maximal_orders[order_num]
  units = Vector{elem_type(O)}()
  while length(units) < n
    LtoA = get_new_field(F, order_num)
    L = domain(LtoA)
    OL = maximal_order(L)
    K, KtoL, ktoK = simplified_absolute_field(L)
    OK = maximal_order_via_relative(K, KtoL)
    UK, mUK = unit_group(OK)
    for j = 1:ngens(UK)
      u = mUK(UK[j])
      if u == -one(OK)
        continue
      end
      push!(units, O(LtoA(KtoL(elem_in_nf(u, copy = false)))))
    end
  end
  return units
end

function _eichler_find_transforming_unit_recurse(I::S, J::S, primes::Vector{T}) where { S <: Union{ AlgAssAbsOrdIdl, AlgAssRelOrdIdl }, T <: Union{ Int, fmpz, NfAbsOrdIdl, NfRelOrdIdl } }
  if length(primes) == 1
    u = _eichler_find_transforming_unit_maximal(I, J)
    return elem_in_algebra(u, copy = false)
  end

  p = pop!(primes)
  IJ = I + J
  if valuation(normred(IJ, left_order(I)), p) != 0
    M = maximal_integral_ideal_containing(IJ, p, :left)
    I = divexact_left(I, M)
    J = divexact_left(J, M)
    return _eichler_find_transforming_unit_recurse(I, J, primes)
  end

  M = maximal_integral_ideal_containing(I, p, :left)
  N = maximal_integral_ideal_containing(J, p, :left)
  u = elem_in_algebra(_eichler_find_transforming_unit_maximal(M, N), copy = false)
  v = inv(u)
  I = divexact_left(I, M)
  OI = left_order(I) # We need it later anyway and it should be the same (as in ===) as the one of J

  J = divexact_left(J, N)
  J = v*J
  J = J*u
  J.left_order = OI

  t = _eichler_find_transforming_unit_recurse(I, J, primes)
  return u*t
end

# Finds x\in A^\times such that I == Jx.
# Requires nr(I) == nr(J) and that nr(I) is coprime to r*d where d is the
# discriminant of O and r in base_ring(O) such that r O_i \subseteq O_j for a
# system of representatives of the maximal orders.
function _eichler_find_transforming_unit(I::AlgAssRelOrdIdl, J::AlgAssRelOrdIdl)
  O = left_order(I)
  @assert ismaximal(O)
  # We assume that left_order(J) == O, but testing this would be really expensive

  n = normred(I, O)
  @assert n == normred(J, O)
  @assert denominator(n, copy = false) == 1
  n = numerator(n, copy = false)

  fac_n = factor(n)
  if isempty(fac_n)
    return one(algebra(order(O)))
  end

  primes = Vector{ideal_type(base_ring(O))}()
  for (p, e) in fac_n
    for i = 1:e
      push!(primes, p)
    end
  end
  t = _eichler_find_transforming_unit_recurse(I, J, primes)
  return t
end

degree(F::Union{ GaloisField, Generic.ResField{fmpz} }) = 1

function get_coeff_fmpz!(x::fq_nmod, n::Int, z::fmpz)
  ccall((:fmpz_set_ui, libflint), Nothing, (Ref{fmpz}, UInt), z, ccall((:nmod_poly_get_coeff_ui, libflint), UInt, (Ref{fq_nmod}, Int), x, n))
  return z
end

function get_coeff_fmpz!(x::fq, n::Int, z::fmpz)
  ccall((:fmpz_poly_get_coeff_fmpz, libflint), Nothing, (Ref{fmpz}, Ref{fq}, Int), z, x, n)
  return z
end

function lift!(x::gfp_elem, z::fmpz)
  ccall((:fmpz_set_ui, libflint), Nothing, (Ref{fmpz}, UInt), z, x.data)
  return z
end

function lift!(x::Generic.ResF{fmpz}, z::fmpz)
  ccall((:fmpz_set, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}), z, x.data)
  return z
end

# Let v and w be non-zero k \times 1-matrices over a finite field and let
# generators be a vector of invertible k \times k-matrices over this field.
# This functions finds a matrix g in the group generated by the elements in
# generators such that gv == w, if such a matrix exists.
# Returns false, rubbish if no such element exists or true and a vector I of
# indices such that g = prod([ generators[i] for i in reverse(I) ]) is the
# searched element.
# This really SLOW! (But it is a hard problem, I guess.)
function find_path(generators::Vector{T}, v::T, w::T) where { T <: MatElem }

  function _weight(v::MatElem{T}, n::fmpz, t::fmpz, jmax::Int) where { T <: Union{ fq, fq_nmod } }
    n = zero!(n)
    for i = 1:nrows(v)
      for j = 0:jmax
        n = add!(n, n, get_coeff_fmpz!(v[i, 1], j, t))
      end
    end
    return n
  end

  function _weight(v::MatElem{T}, n::fmpz, t::fmpz, jmax::Int) where { T <: Union{ gfp_elem, Generic.ResF{fmpz} } }
    n = zero!(n)
    for i = 1:nrows(v)
      n = add!(n, n, lift!(v[i, 1], t))
    end
    return n
  end

  # The main philosophy is to consider the graph G whose vertices are the elements
  # of F_q^k (k := nrows(v)) and in which there is an edge between two vertices
  # v_1 and v_2 if there exsists an element g in generators such that g v_1 = v_2.
  # We then do a graph traversal starting in v and stop as soon as we reach w
  # (or we visited all vertices and w is not in the graph).

  # Initialize the edges (we also add the products of generators)
  edges = Vector{T}()
  edges_in_generators = Vector{Vector{Int}}()
  set_of_edges = Set{T}()
  for i = 1:length(generators)
    e = generators[i]
    if e in set_of_edges
      continue
    end
    push!(edges, e)
    push!(set_of_edges, e)
    push!(edges_in_generators, [ i ])
  end
  for i = 1:length(generators)
    for j = 1:length(generators)
      e = generators[i]*generators[j]
      if e in set_of_edges
        continue
      end
      push!(edges, e)
      push!(set_of_edges, e)
      push!(edges_in_generators, [ j, i ]) # One has to reverse the order
    end
  end
  inv_edges = Dict{Int, T}()
  last_added_edge = Int[ 1, 0 ]

  # We first do a depth first search, where we only visit a vertex if it has
  # smaller weight (with respect to the above weight function).

  # Initialize weight function stuff
  Fq = base_ring(parent(v))
  p = characteristic(Fq)
  ngv = fmpz()
  t = fmpz()
  jmax = degree(Fq) - 1
  minusW = -w
  temp = deepcopy(v)

  res = Vector{Int}()
  visited = Dict{T, Int}()
  no_return = false
  visited[v] = 0
  temp = add!(temp, v, minusW)
  nv = deepcopy(_weight(temp, ngv, t, jmax))
  gv = deepcopy(v)

  # Start the depth first search
  while !iszero(nv)
    j = visited[v] + 1

    # Find the next vertex
    ngv = deepcopy(nv)
    while j <= length(edges)
      gv = mul!(gv, edges[j], v)
      if haskey(visited, gv)
        j += 1
        continue
      end

      # We have not been to gv before, check whether it has smaller weight
      temp = add!(temp, gv, minusW)
      ngv = _weight(temp, ngv, t, jmax)
      if ngv < nv
        break
      end
      j += 1
    end
    visited[v] = j

    if ngv < nv
      v = deepcopy(gv)
      visited[v] = 0
      nv = deepcopy(ngv)
      push!(res, j)
      # Go to gv!
      continue
    end

    # The search ended without a new vertex: we have to go back
    if isempty(res)
      # We are back at the start: add some edges
      if last_added_edge[2] == length(generators)
        if last_added_edge[1] == length(generators)
          # We already added all pairs and 3-tuples of generators, so we give up
          # and do breadth first search now.
          no_return = true
          break
        end
        last_added_edge[1] += 1
      else
        last_added_edge[2] += 1
      end
      e1 = generators[last_added_edge[1]]
      e2 = generators[last_added_edge[2]]
      e12 = e1*e2
      for i = 1:length(generators)
        e = e12*generators[i]
        if e in set_of_edges
          continue
        end
        push!(edges, e)
        push!(set_of_edges, e)
        push!(edges_in_generators, [ i, last_added_edge[2], last_added_edge[1] ])
      end
    else
      # We go one step back
      i = pop!(res)
      if !haskey(inv_edges, i)
        inv_edges[i] = inv(edges[i])
      end
      v = inv_edges[i]*v
      temp = add!(temp, v, minusW)
      nv = deepcopy(_weight(temp, ngv, t, jmax))
    end
  end

  if !no_return
    # We finished with a solution!
    res_in_generators = Vector{Int}()
    for i in res
      append!(res_in_generators, edges_in_generators[i])
    end

    return true, res_in_generators
  end

  # We could not find a path until now, so we do breadth first search (without
  # any weight function)
  Q = Vector{Tuple{T, Vector{Int}}}()
  # I would like to use a queue here, but they are in DataStructures.jl.
  #Q = Queue{Tuple{T, Vector{Int}}}()
  visited = Set{T}()
  push!(Q, (v, Int[]))
  #enqueue!(Q, (v, Int[]))
  push!(visited, v)
  while !isempty(Q)
    u, path_to_u = popfirst!(Q)
    #u, path_to_u = dequeue!(Q)
    for i = 1:length(edges)
      temp = mul!(temp, edges[i], u)
      if temp in visited
        continue
      end
      if temp == w
        path_in_generators = Vector{Int}()
        for j in path_to_u
          append!(path_in_generators, edges_in_generators[j])
        end
        append!(path_in_generators, edges_in_generators[i])
        return true, path_in_generators
      end
      gu = deepcopy(temp)
      push!(visited, gu)
      push!(Q, (gu, push!(copy(path_to_u), i)))
      #enqueue!(Q, (gu, push!(copy(path_to_u), i)))
    end
  end
  # There does not exist a path
  return false, Vector{Int}()
end


text:
function guess()
    number = dec(rand(1:10))
    print("Guess my number! ")
    while readline() != number
        print("Nope, try again... ")
    end
    println("Well guessed!")
end

guess()

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is part of Kpax3. License is MIT.

"""
  readfasta(ifile::String, protein::Bool, miss::Vector{UInt8}, l::Int,
verbose::Bool, verbosestep::Int)

Read data in FASTA format and convert it to an integer matrix. Sequences are
required to be aligned. Only polymorphic columns are stored.

## Arguments

- `ifile` Path to the input data file
- `protein` `true` if reading protein data or `false` if reading DNA data
- `miss` Characters (as `UInt8`) to be considered missing. Use
`miss = zeros(UInt8, 1)` if all characters are to be considered valid. Default
characters for `miss` are:

    - DNA data: _?, \\*, #, -, b, d, h, k, m, n, r, s, v, w, x, y, j, z_
    - Protein data: _?, \\*, #, -, b, j, x, z_

- `l` Sequence length. If unknown, it is better to choose a value which is
surely greater than the real sequence length. If `l` is found to be
insufficient, the array size is dynamically increased (not recommended from a
performance point of view). Default value should be sufficient for most
datasets
- `verbose` If `true`, print status reports
- `verbosestep` Print a status report every `verbosestep` read sequences

## Details

When computing evolutionary distances, don't put the gap symbol `-` among the
missing values. Indeed, indels are an important piece of information for
genetic distances.

FASTA data is encoded as standard 7-bit ASCII codes. The only exception is
Uracil which is given the same value 84 of Thymidine, i.e. each 'u' is silently
converted to 't' when reading DNA data. Conversion tables are the following:

+----------------------------------------+
|         Conversion table (DNA)         |
+----------------------------------------+
|          Nucleotide |  Code |  Integer |
+---------------------+-------+----------+
|           Adenosine |   A   |    97    |
|            Cytosine |   C   |    99    |
|             Guanine |   G   |   103    |
|           Thymidine |   T   |   116    |
|              Uracil |   U   |   116    |
|     Purine (A or G) |   R   |   114    |
| Pyrimidine (C or T) |   Y   |   121    |
|                Keto |   K   |   107    |
|               Amino |   M   |   109    |
|  Strong Interaction |   S   |   115    |
|    Weak Interaction |   W   |   119    |
|               Not A |   B   |    98    |
|               Not C |   D   |   100    |
|               Not G |   H   |   104    |
|          Not T or U |   V   |   118    |
|                 Any |   N   |   110    |
|                 Gap |   -   |    45    |
|              Masked |   X   |   120    |
+----------------------------------------+

+------------------------------------------------+
|           Conversion table (PROTEIN)           |
+------------------------------------------------+
|                  Amino Acid |  Code |  Integer |
+-----------------------------+-------+----------+
|                     Alanine |   A   |    97    |
|                    Arginine |   R   |   114    |
|                  Asparagine |   N   |   110    |
|               Aspartic acid |   D   |   100    |
|                    Cysteine |   C   |    99    |
|                   Glutamine |   Q   |   113    |
|               Glutamic acid |   E   |   101    |
|                     Glycine |   G   |   103    |
|                   Histidine |   H   |   104    |
|                  Isoleucine |   I   |   105    |
|                     Leucine |   L   |   108    |
|                      Lysine |   K   |   107    |
|                  Methionine |   M   |   109    |
|               Phenylalanine |   F   |   102    |
|                     Proline |   P   |   112    |
|                 Pyrrolysine |   O   |   111    |
|              Selenocysteine |   U   |   117    |
|                      Serine |   S   |   115    |
|                   Threonine |   T   |   116    |
|                  Tryptophan |   W   |   119    |
|                    Tyrosine |   Y   |   121    |
|                      Valine |   V   |   118    |
| Asparagine or Aspartic acid |   B   |    98    |
|  Glutamine or Glutamic acid |   Z   |   122    |
|       Leucine or Isoleucine |   J   |   106    |
|                         Gap |   -   |    45    |
|            Translation stop |   *   |    42    |
|                         Any |   X   |   120    |
+------------------------------------------------+

## Value

A tuple containing the following variables:

- `data` Multiple Sequence Alignment (MSA) encoded as a UInt8 matrix
- `id` Units' ids
- `ref` Reference sequence, i.e. a vector of the same length of the original
sequences storing the values of homogeneous sites. SNPs are instead represented
by a value of 46 ('.')
"""
function readfasta(ifile::String,
                   protein::Bool,
                   miss::Vector{UInt8},
                   l::Int,
                   verbose::Bool,
                   verbosestep::Int)
  # we read the file twice
  # the first time we check the total number of sequences, if each sequence
  # has the same length and at what location the SNPs are
  # the second time we store the data, keeping only the SNPs

  # note: this function has been written with a huge dataset in mind
  f = open(ifile, "r")

  s = strip(readuntil(f, '>', keep=true))
  if length(s) == 0
    close(f)
    throw(KFASTAError("No sequence was found."))
  elseif length(s) > 1
    close(f)
    throw(KFASTAError("First non empty row is not a sequence id."))
  end

  # we now know that there is a first sequence to read... but is the ID empty?
  sid = strip(readuntil(f, '\n'))

  if length(sid) == 0
    close(f)
    throw(KFASTAError("Missing sequence identifier. Sequence: 1."))
  end

  if !isascii(sid)
    close(f)
    throw(KFASTAError("FASTA file is not ASCII encoded."))
  end

  seqlen = 0
  tmpseqref = zeros(UInt8, l)
  tmpmissseqref = falses(length(tmpseqref))

  # support variables
  c = '\0'
  w = 0
  u = 0x00

  # start reading the first sequence
  keepgoing = true
  while keepgoing
    s = strip(readline(f))

    if !isascii(s)
      throw(KFASTAError("FASTA file is not ASCII encoded."))
    end

    if length(s) > 0
      if s[1] != '>'
        w = seqlen + length(s)

        # do we have enough space to store the first sequence?
        if w > length(tmpseqref)
          tmpseqref = copyto!(zeros(UInt8, w + l), tmpseqref)
          tmpmissseqref = copyto!(falses(w + l), tmpmissseqref)
        end

        for c in s
          u = UInt8(lowercase(c))

          if !((u == UInt8(32)) || (UInt8(9) <= u <= UInt8(13))) # skip blanks
            seqlen += 1

            if !in(u, miss)
              tmpseqref[seqlen] = u
              tmpmissseqref[seqlen] = false
            else
              tmpseqref[seqlen] = 0x00
              tmpmissseqref[seqlen] = true
            end
          end
        end
      else
        keepgoing = false
      end
    elseif eof(f)
      keepgoing = false
    end
  end

  if seqlen == 0
    close(f)
    throw(KFASTAError("Missing sequence. Sequence: 1."))
  end

  # at least a sequence has been found
  n = 1

  seqref = copyto!(zeros(UInt8, seqlen), 1, tmpseqref, 1, seqlen)
  missseqref = copyto!(falses(seqlen), 1, tmpmissseqref, 1, seqlen)

  if s[1] == '>'
    if length(s) > 1
      sid = lstrip(s, '>')
    else
      # there is only the '>' character
      close(f)
      throw(KFASTAError("Missing sequence identifier. Sequence: 2."))
    end
  end

  curlen = 0

  snp = falses(seqlen)
  seq = zeros(UInt8, seqlen)
  missseq = falses(seqlen)

  for line in eachline(f)
    s = strip(line)

    if !isascii(s)
      throw(KFASTAError("FASTA file is not ASCII encoded."))
    end

    if length(s) > 0
      if s[1] != '>'
        for c in s
          u = UInt8(lowercase(c))

          if !((u == UInt8(32)) || (UInt8(9) <= u <= UInt8(13))) # skip blanks
            curlen += 1

            if curlen > seqlen
              close(f)
              throw(KFASTAError(string("Different sequence length: sequence ",
                                       n + 1, " (", sid, ") ", "is longer ",
                                       "than expected.")))
            end

            if !in(u, miss)
              seq[curlen] = u
              missseq[curlen] = false
            else
              seq[curlen] = 0x00
              missseq[curlen] = true
            end
          end
        end
      else
        # we just finished scanning the previous sequence
        if curlen != seqlen
          close(f)
          throw(KFASTAError(string("Different sequence length: sequence ",
                                   n + 1, " (", sid, ") ", "is shorter than ",
                                   "expected.")))
        end

        n += 1

        for b in 1:seqlen
          if !missseq[b]
            if missseqref[b]
              # this sequence has a non-missing value where it is missing in
              # the reference sequence
              seqref[b] = seq[b]
              missseqref[b] = false
            else
              # to be compared, both values must be non-missing
              snp[b] = snp[b] || (seq[b] != seqref[b])
            end
          end
        end

        if verbose && (n % verbosestep == 0)
          println(n, " sequences have been pre-processed.")
        end

        if length(s) > 1
          sid = lstrip(s, '>')
          curlen = 0
        else
          # there is only the '>' character
          close(f)
          throw(KFASTAError(string("Missing identifier at sequence ", n + 1,
                                   ".")))
        end
      end
    end
  end

  # by construction, the last sequence has not been pre-processed
  if curlen != seqlen
    close(f)
    throw(KFASTAError(string("Different sequence length: sequence ", n + 1,
                             "(", sid, ") ", "is shorter than expected.")))
  end

  n += 1

  for b in 1:seqlen
    if !missseq[b]
      if missseqref[b]
        # this sequence has a non-missing value where it is missing in the
        # reference sequence
        seqref[b] = seq[b]
        missseqref[b] = false
      else
        # to be compared, both values must be non-missing
        snp[b] = snp[b] || (seq[b] != seqref[b])
      end
    end
  end

  m = 0
  for b in 1:seqlen
    if snp[b]
      m += 1
    end
  end

  if verbose
    println("Found ", n, " sequences: ", m, " SNPs out of ", seqlen,
            " total sites.\nProcessing data...")
  end

  data = zeros(UInt8, m, n)
  id = Array{String}(undef, n)
  enc = zeros(UInt8, 127)

  h1 = 0
  h2 = 0
  missing = false

  for h1 in 1:127
    missing = false
    h2 = 1

    while !missing && (h2 <= length(miss))
      missing = (UInt8(h1) == miss[h2])
      h2 += 1
    end

    enc[h1] = !missing ? UInt8(h1) : UInt8('?')
  end

  if !protein
    enc[Int('u')] = UInt8('t')
  end

  # go back at the beginning of the file and start again
  seekstart(f)

  # we already checked that the first non-empty element is the id
  i = 1

  s = strip(readuntil(f, '>'))
  id[i] = strip(readuntil(f, '\n'))

  i1 = 0
  i2 = 0
  for line in eachline(f)
    s = strip(line)

    if length(s) > 0
      if s[1] != '>'
        for c in s
          u = UInt8(lowercase(c))

          if !((u == UInt8(32)) || (UInt8(9) <= u <= UInt8(13))) # skip blanks
            i1 += 1

            if snp[i1]
              i2 += 1
              data[i2, i] = enc[u]
            end
          end
        end
      else
        i += 1

        if verbose && (i % verbosestep == 0)
          println(i, " sequences have been processed.")
        end

        # move on with the next sequence
        id[i] = lstrip(s, '>')

        i1 = 0
        i2 = 0
      end
    end
  end

  if verbose
    println("All ", i, " sequences have been processed.")
  end

  close(f)

  ref = fill(UInt8('.'), seqlen)
  for b in 1:seqlen
    if !snp[b]
      ref[b] = seqref[b]
    end
  end

  (data, id, ref)
end

function countfields(string::String,
                     delim::Char)
  ncol = 0

  idx = search(string, delim)
  while idx > 0
    ncol += 1
    idx = search(string, delim, idx + 1)
  end

  ncol
end

function readdata(ifile::String,
                  delim::Char,
                  miss::Vector{String},
                  verbose::Bool,
                  verbosestep::Int)
  # characters to strip from data values
  chars = ['"', ' ', '\t', '\n', '\r']

  f = open(ifile, "r")

  # find the first sequence (skip initial blank lines)
  s = ""
  keepgoing = true
  while keepgoing
    s = strip(readline(f))

    if length(s) > 0 || eof(f)
        keepgoing = false
    end
  end

  if length(s) == 0
    close(f)
    throw(KCSVError("No observation was found."))
  end

  obsref = map(x -> String(strip(x, chars)), split(s, ','))

  p = length(obsref)

  if p == 1
    close(f)
    throw(KCSVError(string("Delimiter '",  delim,
                           "' was not found in the first observation.")))
  end

  # at least a sequence has been found
  n = 1

  # polymorphic columns
  pol = falses(p)

  # missing values
  missobs = indexin(obsref, miss) .!= nothing

  # count the observations and check that they are all the same length
  for line in eachline(f)
    s = strip(line)

    if length(s) > 0
      obs = map(x -> String(strip(x, chars)), split(s, ','))

      if p != length(obs)
        close(f)
        throw(KCSVError(string("Different length. Failed at observation ",
                               n + 1, " (", obs[1], ").")))
      end

      for b in 2:p
        if !(obs[b] in miss)
          if missobs[b]
            # this observation has a non-missing value where it is missing in
            # the reference observation
            obsref[b] = obs[b]
            missobs[b] = false
          else
            # to be compared, both values must be non-missing
            pol[b] = pol[b] || (obs[b] != obsref[b])
          end
        end
      end

      n += 1

      if verbose && (n % verbosestep == 0)
        println(n, " observations have been pre-processed.")
      end
    end
  end

  m = 0
  for b in 2:p
    if pol[b]
      m += 1
    end
  end

  if verbose
    println("Found ", n, " observations: ", m, " polymorphisms out of ",
            p - 1, " total columns.\nProcessing data...")
  end

  data = Array{String}(undef, m, n)
  id = Array{String}(undef, n)

  # go back at the beginning of the file and start again
  seekstart(f)

  i = 0
  for line in eachline(f)
    s = strip(line)

    if length(s) > 0
      obs = map(x -> String(strip(x, chars)), split(s, ','))

      i += 1

      j = 0
      for b in 2:p
        if pol[b]
          data[j += 1, i] = !(obs[b] in miss) ? obs[b] : ""
        end
      end

      id[i] = obs[1]

      if verbose && (i % verbosestep == 0)
        println(i, " observations have been processed.")
      end
    end
  end

  if verbose
    println("All ", i, " sequences have been processed.")
  end

  close(f)

  ref = fill(".", p - 1)
  for a in 2:p
    if !pol[a]
      ref[a - 1] = obsref[a]
    end
  end

  (data, id, ref)
end


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bilevel_limitables" begin

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")
    function create_instance(load_1, load_2,
                            wind_1,
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "wind_2_1", PSCOPF.Networks.LIMITABLE,
        #                                         0., 100.,
        #                                         0., 1.,
        #                                         Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[1], "S1", wind_1)
        # PSCOPF.add_uncertainty!(uncertainties, ech, "wind_2_1", TS[1], "S1", wind_2)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", load_1)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", load_2)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    SortedDict{String,PSCOPF.GeneratorState}(), #gen_initial_state
                                    uncertainties, nothing, logs)

        return context
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|                      |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |        "1_2"         |
      S1: 40            |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 30

    Available limitable production : 40
    Demand : 40

    Market : use wind_1_1 at 40.
    This satisfies the EOD constraint
        e = p_capping = 0.
        lol = p_loss_of_load = 0.
    This does not induce any RSO constraint violation

    The TSO does not need to take any actions :
        e_min = p_global_capping = 0.
        lol_min = p_global_loss_of_load = 0.
        no limitation

    TSO cost : 0
    Market cost : 0
    =#
    @testset "no_problem" begin

        context = create_instance(10., 30.,
                                40.,
                                35.)

        tso = PSCOPF.TSOBilevel()
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #no limitation
        @test value(result.upper.limitable_model.p_limit["wind_1_1",TS[1],"S1"]) > 40. - 1e-09
        @test value(result.upper.limitable_model.b_is_limited["wind_1_1",TS[1],"S1"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #costs
        @test value(PSCOPF.get_upper_obj_expr(result)) < 1e-09
        @test value(PSCOPF.get_lower_obj_expr(result)) < 1e-09

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|                      |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |         "1_2"        |
      S1: 100           |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 30

    Available limitable production : 100
    Demand : 40

    Market :
        The market needs to satisfy EOD constraint => only option is to cap 60MW of limitable production.
            e = p_capping = 60.
            lol = p_loss_of_load = 0.

    The TSO locates the capping to avoid RSO constraints (even if there is no risk of violating RSO here)
        p_capping[wind_1_1]=60, plim[wind1]=40, p_injected[wind_1_1]=40

    The TSO takes a limitation action:
      e_min = p_global_capping = 0.
      lol_min = p_global_loss_of_load = 0.
      plim[wind1] = 40 (due to locating capping)

    TSO cost : 0.001 (one limitation)
    Market cost : 60 (capping)

    NOTE: here, there is only one limitable,
     normally the TSO should not have to limit the production of wind_1_1
     cause the market does not have other options
    FIXME? : TSO limits for EOD reasons.
    =#
    @testset "eod_overproduction_requires_market_capping_and_tso_limitation" begin
        context = create_instance(10., 30.,
                                100.,
                                35.,
                                "test_li")

        tso = PSCOPF.TSOBilevel()
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #TSO :
        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) <= 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) <= 1e-09

        # TSO needs to limit because it is him who locates capping
        @test 40. ≈ value(result.upper.limitable_model.p_limit["wind_1_1",TS[1],"S1"])
        @test 1. ≈ value(result.upper.limitable_model.b_is_limited["wind_1_1",TS[1],"S1"])
        @test 60. ≈ ( value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"]))

        #Market :
        #Market caps ENR for EOD reasons
        @test 60. ≈ value(result.lower.limitable_model.p_global_capping[TS[1],"S1"])
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) <= 1e-09

        #costs
        @test 0.001 == (1*tso.configs.TSO_LIMIT_PENALTY) #one limitation
        @test ≈(0.001, value(PSCOPF.get_upper_obj_expr(result)), atol=1e-04) #due to the 1e5 objective coeff on global LoL
        @test (60. *tso.configs.MARKET_CAPPING_COST) ≈ value(PSCOPF.get_lower_obj_expr(result))

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 40            |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 100           |                      | S1: 30

    Available limitable production : 40
    Demand : 130

    Market :
        The market needs to satisfy EOD => reduces the consumption to 40
            e = p_capping = 0.
            lol = p_loss_of_load = 90.

    The TSO locates the LoL to avoid RSO constraints
        Here, there is no risk of violating the RSO constraint => any option will do:
            p_loss_of_load[bus_1]+p_loss_of_load[bus_2] = 90

    The TSO takes no further actions
      e_min = p_global_capping = 0.
      lol_min = p_global_loss_of_load = 0.
      no limitation

    TSO cost : 0
    Market cost : 90*1e5 (capping)

    FIXME? : cost LoL locations differently
    =#
    @testset "eod_overload_requires_market_lol" begin
        context = create_instance(100., 30.,
                                40.,
                                35.)

        tso = PSCOPF.TSOBilevel()
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Market needed to cut conso
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO :
        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) <= 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) <= 1e-09

        #no limitation
        @test value(result.upper.limitable_model.p_limit["wind_1_1",TS[1],"S1"]) > 40. - 1e-09
        @test value(result.upper.limitable_model.b_is_limited["wind_1_1",TS[1],"S1"]) < 1e-09
        @test value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"]) < 1e-09

        #Market :
        #Market cuts conso for EOD reasons
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 90. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        #TSO distributes the cut conso (arbitrarily) while assuring RSO constraint
        @test 90. ≈ ( value(result.upper.lol_model.p_loss_of_load["bus_1", TS[1],"S1"])
                    + value(result.upper.lol_model.p_loss_of_load["bus_2", TS[1],"S1"]) )

        #costs
        @test (90*tso.configs.TSO_LOL_PENALTY) ≈ value(PSCOPF.get_upper_obj_expr(result))
        @test (90*tso.configs.MARKET_LOL_PENALTY) ≈ value(PSCOPF.get_lower_obj_expr(result))

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 50            |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 40

    Available limitable production : 50
    Demand : 50

    Market :
        The EOD is satisfied, the market would simply use the limitable production to satisfy the demand.
        => this will cause a flow[branch_1_2] = 40 > 35  ==> violates the RSO constraint

    The TSO needs to impose some additional constraints to oblige the market to respect the RSO constraint:
        option 1 : impose a capping of 5MW, limit wind_1_1's production to 45MW
        => cost : 5 + 0.001
            This will oblige the market to cut conso by 5.
            Since the TSO locates the cut conso, he can assure cutting conso on bus2 to assure the RSO constraint
        option 2 : impose reducing the consumption of bus2 by 5MW, limit wind_1_1's production to 45MW
            bus 2 will only require 35MW which assures RSO constraints
            TSO will need to limit the prod of wind_1_1
        => cost : 5e05 + 0.001

    TSO adopts the cheaper option : option 1
        e_min = p_global_capping = 5.
        lol_min = p_global_loss_of_load = 0.
        plimit[wind_1_1] = 45

    Market
    The market is now constrained by the TSO decisions:
        e = p_capping = 5. (due to TSO decisions)
        lol = p_loss_of_load = 5. (due to EOD)

    The TSO locates the capping and LoL:
        p_capping[wind_1_1]=5
        p_loss_of_load[bus_1]=0, p_loss_of_load[bus_2]=5

    TSO cost : 5.001
    Market cost : 5 + 5*1e5

    #FIXME? TSO obliges cuts for RSO reasons
    #normally both capping and lol are due to RSO constraints
    #however the TSO simply tells the market to cap some MW which only costs 5€
    =#
    @testset "rso_problem" begin
        context = create_instance(10., 40.,
                                50.,
                                35.)

        tso = PSCOPF.TSOBilevel()
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Market needed to cut conso
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO :
        #TSO imposes caping 5MW
        @test 5. ≈ value(result.upper.limitable_model.p_global_capping[TS[1],"S1"])
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #limit wind_1_1
        @test 45. ≈ value(result.upper.limitable_model.p_limit["wind_1_1",TS[1],"S1"])
        @test 1 ≈ value(result.upper.limitable_model.b_is_limited["wind_1_1",TS[1],"S1"])
        @test 5. ≈ value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"])

        #Market :
        #Market only cuts as required since EOD is OK
        @test 5. ≈ value(result.lower.limitable_model.p_global_capping[TS[1],"S1"])
        @test 5. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        #TSO decides the distribution :
        @test 5. ≈ ( value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"]))
        @test 5. ≈ ( value(result.upper.lol_model.p_loss_of_load["bus_1",TS[1],"S1"])
                    + value(result.upper.lol_model.p_loss_of_load["bus_2",TS[1],"S1"]) )
        @test value(result.upper.lol_model.p_loss_of_load["bus_1",TS[1],"S1"]) < 1e-09
        @test 5. ≈ value(result.upper.lol_model.p_loss_of_load["bus_2",TS[1],"S1"])

        #costs
        @test value(PSCOPF.get_upper_obj_expr(result)) ≈  (5*tso.configs.TSO_CAPPING_COST + 1*tso.configs.TSO_LIMIT_PENALTY
                                                            + 5*tso.configs.TSO_LOL_PENALTY )
        @test (5*tso.configs.MARKET_CAPPING_COST + 5*tso.configs.MARKET_LOL_PENALTY) ≈ value(PSCOPF.get_lower_obj_expr(result))
    end

    #=
    same as preceding : c.f. rso_problem

    option 1 : impose a capping of 5MW, limit wind_1_1's production to 45MW
    => default parameters cost : 5 + 0.001
    => cost : 5e03 + 0.001
    option 2 : impose reducing the consumption of bus2 by 5MW, limit wind_1_1's production to 45MW
    => default parameters cost : 5e05 + 0.001
    => cost : 5*3 + 0.001

    with new config, option 2 is more appealing
    =#
    @testset "rso_problem_with_modifed_costs" begin
        context = create_instance(10., 40.,
                                50.,
                                35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(
                                            TSO_LOL_PENALTY = 3.,
                                            TSO_CAPPING_COST = 1e3
                                ))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Market needed to cut conso
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO :
        #TSO imposes reducing the consumption
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 5. ≈ value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"])

        #limit wind_1_1
        @test 45. ≈ value(result.upper.limitable_model.p_limit["wind_1_1",TS[1],"S1"])
        @test 1 ≈ value(result.upper.limitable_model.b_is_limited["wind_1_1",TS[1],"S1"])
        @test 5. ≈ value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"])

        #Market :
        #Market only cuts as required since EOD is OK
        @test 5. ≈ value(result.lower.limitable_model.p_global_capping[TS[1],"S1"])
        @test 5. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        #TSO decides the distribution :
        @test 5. ≈ ( value(result.upper.limitable_model.p_capping["wind_1_1",TS[1],"S1"]))
        @test 5. ≈ ( value(result.upper.lol_model.p_loss_of_load["bus_1",TS[1],"S1"])
                    + value(result.upper.lol_model.p_loss_of_load["bus_2",TS[1],"S1"]) )
        @test value(result.upper.lol_model.p_loss_of_load["bus_1",TS[1],"S1"]) < 1e-09
        @test 5. ≈ value(result.upper.lol_model.p_loss_of_load["bus_2",TS[1],"S1"])

        #costs
        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ 30.001 ≈ (5*tso.configs.TSO_LOL_PENALTY + #due to TSO LoL
                                                                    1*tso.configs.TSO_LIMIT_PENALTY +
                                                                    5*tso.configs.TSO_LOL_PENALTY) #due to market LoL
        @test (5+5e05) == (5*tso.configs.MARKET_CAPPING_COST + 5*tso.configs.MARKET_LOL_PENALTY) ≈ value(PSCOPF.get_lower_obj_expr(result))
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using DifferentialEquations

# TODOS:
# macro generate dTdt_s! for MORE speed...
# clean up and doc ...

"""
EssentialModelComprehension:
Essential Properties of the different modell structs, needed to setup the solution
procedure
"""
struct EssentialModelComprehension{T1 <: AbstractFloat, T2 <: Integer}
    n_c::T2
    n_c_s_sub::T2
    T_in_f::T1
    T_in_s::T1
    dx_f::T1
    mf_f::T1
    mf_s::T1
    A_fs::T1
    A_f::T1
    htc_fs::HeatTransfer
    mp_f::IncompressibleFluidThermalMaterial
    mp_s::SolidThermalMaterial
    V_f::T1
    V_s::T1
    V_fs::T1
    scrap_cond::Bool
    dx_s_sub::T1
    A_s_sub::T1
    V_s_sub::T1
    cross_rad::Bool
    front_rad::Bool
    front_rad_model::Union{Nothing,HotSideRadiationModel}
    A_fs_rad::T1
    b::Vector{<:T1}
    k::T2
end

function solve_1d1d_diffeq(m::M1d1d)
    print_log(m1d1d_logfile, "Using diff eq solver!")

    # heat_flow_exchng_sf = zeros(Float64, m.mesh_settings.no_cells_bulk)
    emc = Vector{EssentialModelComprehension{Float64, Int64}}(undef, length(m.sim_phase))

    # Init Results according to used meshes
    results = Vector{Results1d1d}(undef,0)

    mesh_tmp = nothing

    phase_start_time = 0.0
    for pid=1:1:length(m.sim_phase)
        print_log(m1d1d_logfile, ["\nStarting Simulation of Phase ", pid])
        local T

        mesh = M1d1dMesh(m, pid)
        printInfo(mesh; tabs = 1)

        # just make sure that meshes do not change between the phases
        # because we than would have to interpolate between the meshes if they
        # change:
        if pid > 1
            if mesh_tmp != mesh
                error("Mesh change between different phases is not supported atm")
                print_log(m1d1d_logfile, "Error: Unsupported mesh change")
            else
                mesh_tmp = mesh
            end
        else
            mesh_tmp = mesh
        end

        emc[pid] = set_emc(m, pid, mesh)

        if pid > 1
            phase_start_time += m.sim_phase[pid-1].duration
            T = set_T(emc[pid], pid, results[end])
        else
            T = set_T(emc[pid], pid)
        end

        qalpha_ex = zeros(emc[pid].n_c)
        dTdt = zeros(size(T))

        tspan = (phase_start_time, phase_start_time + m.sim_phase[pid].duration)

        # p = [emc[pid], qalpha_ex]
        range = zeros(Int, 2*emc[pid].k)
        p = (emc[pid], qalpha_ex, emc[pid].mp_f.lambda, emc[pid].mp_f.cp_true,
                emc[pid].mp_f.rho, emc[pid].mp_s.lambda, emc[pid].mp_s.cp_true,
                emc[pid].mp_s.rho, emc[pid].mp_s.epsilon, range)
        @show T[1,:]
        @show tspan
        prob = ODEProblem(dTdt_model!,T , tspan, p)

        @time sol = solve(prob,Tsit5(), progress = true, saveat=m.sim_phase[pid].dt_save
                , reltol=1e-8, abstol=1e-8,dtmin=1e-6, maxiters=1e8)

        # test perf / allocs ...
        # fix many allocations of material functions pX_bounden and integrate_pX
        # @time dTdt_model!(dTdt,T,p,0.2)
        # @time dTdt_model2!(dTdt,T,p,0.2)
        push_sol_to_res!(results, m, mesh, emc[pid].mp_s.cp_true, emc[pid].mp_s.rho,
                emc[pid].V_s_sub, sol)
    end

    return results
 end

function push_sol_to_res!(res::Vector{Results1d1d},m::M1d1d, mesh::M1d1dMesh,
        cp_s::MaterialPropertyOfX, rho_s::MaterialPropertyOfX, V_s_sub_c::Real,
        sol)

    for t = 1:length(sol.t)
        Tres = M1d1dVals(mesh, init_val = 0.0)

        T_H = collect(0:0.05:2000)
        H_T_s = H_of_T.(cp_s, T_H)

        @inbounds for i = 1:1:size(sol.u[t],1)
            Tres.f[i] =  sol.u[t][i,1]

            if size(sol.u[t],2) > 2
                for j = 1:1:size(sol.u[t],2)
                    Tres.s_sub[i][j] = sol.u[t][i,j+1]
                end

                # determine Tmean of all T_s_subs if we use only T_s at a later
                # stage
                m_sum = 0.0
                H_T_sub_sum = 0.0
                @simd for j = 2:1:size(sol.u[t],2)
                    m_sum += pX_bounded(rho_s,sol.u[t][i,j])  * V_s_sub_c
                    H_T_sub_sum += H_of_T(cp_s, sol.u[t][i,j]) *
                            pX_bounded(rho_s, sol.u[t][i,j]) * V_s_sub_c
                end

                Tres.s[i] = T_of_H(H_T_s, T_H, H_T_sub_sum/(m_sum))
            else
                Tres.s[i] =  sol.u[t][i,2]

                if !isnothing(Tres.s_sub)
                # make sure that we can start with T_s_sub in later pid run
                    for j = 1:1:length(Tres.s_sub[1])
                        Tres.s_sub[i][j] = sol.u[t][i,2]
                    end
                end
            end
        end

        push!(res, Results1d1d(sol.t[t], Tres))
    end
end

function set_T(emc::EssentialModelComprehension{<:AbstractFloat, <:Integer},
        phaseid::Integer, Tprev::Union{Nothing, Results1d1d} = nothing)

    if emc.scrap_cond
        T = zeros(emc.n_c, emc.n_c_s_sub+1)
    else
        T = zeros(emc.n_c, 2)
    end

    if phaseid == 1
        T .= emc.T_in_s

    elseif !isnothing(Tprev)
        for  i = 1:emc.n_c
            T[i,1] = Tprev.T.f[i]
        end

        if emc.scrap_cond
            for i = 1:emc.n_c, j=2:1:emc.n_c_s_sub+1
                T[i,j] = Tprev.T.s_sub[i][j-1]
            end
        else
            for  i = 1:emc.n_c
                T[i,2] = Tprev.T.s[i]
            end
        end

    else
        error("No Tprev given!")
    end

    return T
 end

function dTdt_model!(dTdt::Matrix{<:Real}, T::Matrix{<:Real}, p::Tuple{EssentialModelComprehension{<:AbstractFloat, <:Integer},
    Vector{<:AbstractFloat}, MaterialPropertyOfX, MaterialPropertyOfX,MaterialPropertyOfX,
    MaterialPropertyOfX,MaterialPropertyOfX,MaterialPropertyOfX,MaterialPropertyOfX,Vector{<:Integer}}, t)
    emc = p[1]
    qalpha_ex = p[2]
    lam_f = p[3]
    cp_f = p[4]
    rho_f = p[5]
    lam_s = p[6]
    cp_s = p[7]
    rho_s = p[8]
    eps_s  = p[9]
    range = p[10]

    dTdt_f!(dTdt, T, qalpha_ex, emc, lam_f, cp_f, rho_f)
    dTdt_s!(dTdt, T, qalpha_ex, emc, lam_s, cp_s, rho_s, eps_s, range)
end

function dTdt_f!(dTdt::Matrix{<:AbstractFloat}, T::Matrix{<:AbstractFloat}, qalpha_ex::Vector{<:AbstractFloat}, emc::EssentialModelComprehension{<:AbstractFloat, <:Integer},
        lam_f::MaterialPropertyOfX, cp_f::MaterialPropertyOfX, rho_f::MaterialPropertyOfX)

    lambda_dT(TI, T2)::Float64 = harm_mean(pX_bounded(lam_f,T2),
            pX_bounded(lam_f, TI)) * (T2 - TI)
    f_dHdt_f_mvmnt(T1, T2)::Float64 = ( emc.mf_f * integrate_pX(cp_f, T1, T2) )

    @inbounds @simd for i = 1:1:emc.n_c
        local dHdT_alpha::Float64
        local dHdT_cond::Float64
        local dHdt_f_mvmnt::Float64
        local dHdt_f::Float64
        local dTdt_f::Float64

        dHdT_alpha = 0.0
        dHdT_alpha = (emc.A_fs * htc_fs(emc.htc_fs) *
                (T[i,2] - T[i,1]))

        qalpha_ex[i] = dHdT_alpha
        dHdT_cond = 0.0
        dHdt_f_mvmnt = 0.0

        if i != 1 && i != emc.n_c
            dHdT_cond = emc.A_f/emc.dx_f * ( lambda_dT(T[i,1], T[i-1,1]) +
                lambda_dT(T[i,1], T[i+1,1]) )

            dHdt_f_mvmnt = f_dHdt_f_mvmnt(T[i,1], T[i-1,1])
        elseif i == emc.n_c
            dHdT_cond = emc.A_f/emc.dx_f * lambda_dT(T[i,1], T[i-1,1])

            dHdt_f_mvmnt = f_dHdt_f_mvmnt(T[i,1], T[i-1,1])
        elseif i == 1
            dHdT_cond = emc.A_f/emc.dx_f * lambda_dT(T[i,1], T[i+1,1])

            dHdt_f_mvmnt = f_dHdt_f_mvmnt(T[i,1], emc.T_in_f)
            pX_bounded(lam_f, 100.0)
        end

        dHdt_f = (dHdt_f_mvmnt + dHdT_alpha + dHdT_cond)
        dTdt_f = 0.0
        dTdt_f = dHdt_f/(emc.V_f * pX_bounded(rho_f, T[i,1]) *
                    pX_bounded(cp_f, T[i,1]))

        dTdt[i,1] = dTdt_f
    end
end

"""
we could optimize that further if we use a macro to write a specific version
for this function ..
"""
function dTdt_s!(dTdt::Matrix{T1}, T::Matrix{T1}, qalpha_ex::Vector{T1},
    emc::EssentialModelComprehension{T1, T2}, lam_s::MaterialPropertyOfX,
    cp_s::MaterialPropertyOfX, rho_s::MaterialPropertyOfX,
    eps_s::MaterialPropertyOfX, range::Vector{T2}) where T1 <: AbstractFloat where T2 <: Integer

    lambda_dT(TI, T2)::T1 = ( harm_mean(pX_bounded(lam_s,T2),
    pX_bounded(lam_s, TI)) * (T2 - TI) )

    #  f(i,k) = vcat(collect(i-k:1:i-1), collect(i+1:1:i+k))
    function update_range!(range, i, k)
        @inbounds @simd for ii = 0:1:k
            range[ii] = (ii+i-1) - k
        end

        @inbounds @simd for ii = k+1:1:2*k
            range[ii] = (ii-k) + i
        end
    end

    @inbounds for i = emc.n_c:-1:1
        dHdT_alpha = -qalpha_ex[i]

        local dHdt_rad_front::T1
        local dHdt_rad_cross::T1
        local dHdt_cond::T1
        local dHdt_s_mvmnt::T1
        local dHdt_s::T1

        dHdt_rad_front = 0.0
        if emc.front_rad
            if i <= emc.k
                if i == 1
                    dHdt_rad_front= heat_flow_front_rm(
                            emc.front_rad_model,
                            pX_bounded(eps_s, T[i,2]), emc.b[i],
                            emc.A_fs_rad, T[i,2])
                else
                    dHdt_rad_front = heat_flow_front_rm(
                            emc.front_rad_model,
                            pX_bounded(eps_s, T[i,2]), emc.b[i+1],
                            emc.A_fs_rad, T[i,2])
                end
            end
        end

        # cross radiation model
        dHdt_rad_cross = 0.0
        if emc.cross_rad
            update_range!(range, i, emc.k)

            for ir = 1:1:2*emc.k
                local r
                r = range[ir]
                if r > 0 && r <= emc.n_c
                    T_r = T[r,2]
                    T_i = T[i,2]

                    r_dist = abs(i-r)

                    dHdt_rad_cross += emc.A_fs_rad *
                            hf_radiation_par_surf_b_weighted(
                            pX_bounded(eps_s, T_r), pX_bounded(eps_s, T_i),
                            emc.b[r_dist+1], T_r, T_i )
                end
            end
        end

        if emc.scrap_cond
            for j = 2:1:emc.n_c_s_sub+1
                dHdt_cond = 0.0
                scrap_cond_fac = emc.A_s_sub/emc.dx_s_sub
                if j != 2 && j != (emc.n_c_s_sub+1)
                    # CDS
                    dHdt_cond = scrap_cond_fac * ( lambda_dT(T[i,j], T[i,j-1]) +
                                lambda_dT(T[i,j], T[i,j+1]) )
                elseif j == (emc.n_c_s_sub+1) && j > 2
                    # symmetry only 1 dir (BDS + 0)
                        dHdt_cond = scrap_cond_fac * lambda_dT(T[i,j], T[i,j-1])
                elseif j == 2
                    # convection (FDS)
                    dHdt_cond = scrap_cond_fac * lambda_dT(T[i,j], T[i,j+1])
                end

                dHdt_s_mvmnt = 0.0
                if i != emc.n_c
                    dHdt_s_mvmnt = (emc.mf_s * integrate_pX(cp_s,
                            T[i,j], T[i+1,j]))
                else
                    dHdt_s_mvmnt = (emc.mf_s * integrate_pX(cp_s,
                            T[i,j], emc.T_in_s))
                end

                dHdt_s = 0.0
                if j == 2
                    dHdt_s = ( 0.5 * (dHdT_alpha +
                            dHdt_rad_front + dHdt_rad_cross) +
                            dHdt_cond + dHdt_s_mvmnt)
                else
                    dHdt_s = (dHdt_s_mvmnt + dHdt_cond)
                end

                dTdt_s = dHdt_s/(emc.V_s_sub * pX_bounded(rho_s, T[i,j]) *
                        pX_bounded(cp_s, T[i,j]))

                dTdt[i,j] = dTdt_s
            end
        else
            dHdt_s_mvmnt = 0.0
            if i != emc.n_c
                dHdt_s_mvmnt = (emc.mf_s * integrate_pX(cp_s,
                        T[i,2], T[i+1,2]))
            else
                dHdt_s_mvmnt = (emc.mf_s * integrate_pX(cp_s,
                        T[i,2], emc.T_in_s))
            end

            dHdt_s = (dHdT_alpha +  dHdt_rad_front + dHdt_rad_cross  + dHdt_s_mvmnt)

            dTdt_s = dHdt_s/(emc.V_s * pX_bounded(rho_s, T[i,2]) *
                    pX_bounded(cp_s, T[i,2]))

            dTdt[i,2] = dTdt_s
        end
    end
end

function set_emc(m::M1d1d, phaseid::Integer, mesh::M1d1dMesh)

    scrap_conduction = false
    scrap_conduction_temp = false
    A_fs_rad = 0.0

    no_cells = mesh.no_cells
    dx_f = mesh.gas_cells.len

    A_f = mesh.gas_cells.dA
    A_s = mesh.solid_cells.dA

    V_f = mesh.gas_cells.vol
    V_s = mesh.solid_cells.vol
    V_fs = V_f + V_s

    mf_f = m.offgas[phaseid].mass_flow
    mf_s = m.bulk[phaseid].mass_flow

    max_b_fac_layers = floor(Int64, no_cells/2)
    b_fac_over_layers = zeros(max_b_fac_layers)
    k = 0

    print_log(m1d1d_logfile, ["Mass flow fluid: ", mf_f])
    print_log(m1d1d_logfile, ["Mass flow solid: ", mf_s])

    A_fs =  mesh.A_fs
    htc_fs = m.heat_transfer[phaseid]

    print_log(m1d1d_logfile, ["Af_s: ", A_fs]);

    no_sub_cells = 0
    dx_s_sub = 0.0
    A_s_sub =  0.0
    V_s_sub = 0.0

    if m.sim_phase[phaseid].submodels.scrap_discretization_model_enabled
        print_log(m1d1d_logfile, "Scrap conduction enabled")

        if isnothing(mesh.solid_1d_submesh)
            error("No submesh generated for scrap part heating!")
        else
            no_sub_cells = mesh.solid_1d_submesh.no_cells
            dx_s_sub = mesh.solid_1d_submesh.sub_cells.len
            A_s_sub = mesh.solid_1d_submesh.sub_cells.dA
            V_s_sub = mesh.solid_1d_submesh.sub_cells.vol
            mf_s = mf_s / (no_sub_cells * 2.0)

            print_log(m1d1d_logfile, "Scrap conduction settings: ")
            print_log(m1d1d_logfile, ["\t no_sub_cells: ", no_sub_cells])
            print_log(m1d1d_logfile, ["\t dx_s_sub: ", dx_s_sub])
            print_log(m1d1d_logfile, ["\t A_s_sub: ", A_s_sub])
            print_log(m1d1d_logfile, ["\t A_s_sub/dx_s_sub: ", A_s_sub/dx_s_sub])
            print_log(m1d1d_logfile, ["\t V_s_sub: ", V_s_sub])
            print_log(m1d1d_logfile, ["\t mf_s: ", mf_s])

            if no_sub_cells < 3
                error("Not able to handle submesh with less than 3 cells")
            end
        end
    end

    print_log(m1d1d_logfile, "Radiation Settings: ")

    cross_radiation = m.sim_phase[phaseid].submodels.cross_radiation_model_enabled

    if cross_radiation
        print_log(m1d1d_logfile,"\t Cross radiation enabled")
    else
        print_log(m1d1d_logfile,"\t Cross radiation disabled")
    end

    front_radiation = !isnothing(
            m.sim_phase[phaseid].submodels.furnace_side_rad_model)

    if front_radiation
        print_log(m1d1d_logfile, "\t Front radiation enabled")
    else
        print_log(m1d1d_logfile, "\t Front radiation disabled")
    end

    if cross_radiation || front_radiation

        c = m.sim_phase[phaseid].submodels.cross_radiation_constant_c

        @assert c <= (1.0 + eps(typeof(c)))  && c >= (0.0 - eps(typeof(c)))

        A_fs_rad = c * 0.5 * A_fs

        A_cross = m.bulk[phaseid].cross_sect_surf_area
        eps_s = m.bulk[phaseid].scrap_props.thermal_props.epsilon

        if A_fs_rad > A_cross/pX_bounded(eps_s, Tmax)
            # TODO: we asume here pX_bounded(eps_s, Tmax) == epsilon_max
            # me may want to enhance the way this limit here is set ...
            warnmsg = ("A_fs_rad bigger than  A_cross/epsilon_max, this would lead " *
                    "to unphysical results, consider choosing a smaller " *
                    "cross radiation constant c. Limiting A_fs_rad to A_cross/epsilon_max")

            @warn   warnmsg

            print_log(m1d1d_logfile, "Warning: " * warnmsg, print_to_stdout=false)
            A_fs_rad =  A_cross/pX_bounded(eps_s, Tmax)

            # TODO: is that right because b[1] stands for self viewing? or not? ...
            k = 1
            b_fac_over_layers[1] = 1.0
            b_fac_over_layers[2] = 1.0

            b_fac_over_layers_print = b_fac_over_layers[b_fac_over_layers .> 1E-3]
            print_log(m1d1d_logfile, ["\t b_fac_over_layers: ", b_fac_over_layers_print])
            print_log(m1d1d_logfile, ["\t k: ",  k])

        elseif A_fs_rad > 0
            mean_overlap = calc_mean_overlap(A_cross, A_fs_rad, 1, 10_000)

            for j = 1:length(b_fac_over_layers)
                if j == 1
                    b_fac_over_layers[j] = 1 # with itself
                elseif j == 2
                    b_fac_over_layers[j] = mean_overlap  # with direct neighbor
                else
                    b_fac_over_layers[j] = (1-mean_overlap).^(j-2) *
                            mean_overlap # thrugh cells bfore
                end
            end

            b_fac_over_layers_print = b_fac_over_layers[b_fac_over_layers .> 1E-3]
            print_log(m1d1d_logfile, ["\t b_fac_over_layers: ", b_fac_over_layers_print])

            k = findfirst( x -> x < 1E-3, b_fac_over_layers)
            print_log(m1d1d_logfile, ["\t k: ",  k])
        else
            error("Wrong settings of cross_radiation_constant_c?")
        end

        if isnothing(k)

            warnmsg = ("b_fac_over_layers[" * string(max_b_fac_layers) *
                    "] is very high, maybe there is " *
                    "something wrong? k will be set to maximum of " *
                    string(max_b_fac_layers))

            @warn warnmsg
            print_log(m1d1d_logfile, "Warning: " * warnmsg, print_to_stdout=false)

            k = max_b_fac_layers
        end

        print_log(m1d1d_logfile, ["\t Af_s_rad: ", A_fs_rad]);
    end

    return EssentialModelComprehension(no_cells, no_sub_cells,
            m.offgas[phaseid].T_in, m.bulk[phaseid].T_charged, dx_f, mf_f, mf_s,
            A_f, A_fs, htc_fs, m.offgas[phaseid].thermal_props,
            m.bulk[phaseid].scrap_props.thermal_props, V_f, V_s, V_fs,
            m.sim_phase[phaseid].submodels.scrap_discretization_model_enabled,
            dx_s_sub, A_s_sub, V_s_sub,
            cross_radiation, front_radiation,
            m.sim_phase[phaseid].submodels.furnace_side_rad_model,
            A_fs_rad, b_fac_over_layers, k)
end


text:
using EHRAuthentication
using Documenter

DocMeta.setdocmeta!(EHRAuthentication, :DocTestSetup, :(using EHRAuthentication); recursive=true)

makedocs(;
    modules=[EHRAuthentication],
    authors="Dilum Aluthge, Pumas-AI Inc., and contributors",
    repo="https://github.com/JuliaHealth/EHRAuthentication.jl/blob/{commit}{path}#{line}",
    sitename="EHRAuthentication.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaHealth.github.io/EHRAuthentication.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaHealth/EHRAuthentication.jl",
    devbranch="main",
)


text:
include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")

# This works assuming the reference solution is on a larger discretization.
function relerr(UN::Matrix, UR::Matrix)
    m, n = size(UN)
    M, N = size(UR)
    ERR = deepcopy(UR)
    @inbounds for j = 1:n, i = 1:m
        ERR[i,j] -= UN[i,j]
    end
    return vecnorm(ERR)/vecnorm(UR)
end

using PyPlot

IMAX = 8
Ntimes = 2

n = 1024

θ = (0.5:n-0.5)/n
φ = (0:2n-2)*2/(2n-1)
F = [cos(10*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)) for θ in θ, φ in φ]
V = zero(F)
A_mul_B!(V, FastTransforms.plan_analysis(F), F)
U0 = fourier2sph(V; sketch = :none)

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.01, -0.5, 1.0), NonlinearOperator(u->u-u^3), U0)
UN = Matrix{Matrix{Float64}}(IMAX, Ntimes)
for j = 1:Ntimes
    T = 1.0*4^(j-1)
    for i = 1:IMAX
        UN[i,j] = ETDRK4(SLPDE, T, round(Int, 2.0^i*T))
        println("Done i = ", i," j = ", j)
    end
end

n = 2048

θ = (0.5:n-0.5)/n
φ = (0:2n-2)*2/(2n-1)
F = [cos(10*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)) for θ in θ, φ in φ]
V = zero(F)
A_mul_B!(V, FastTransforms.plan_analysis(F), F)
U0 = fourier2sph(V; sketch = :none)

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.01, -0.5, 1.0), NonlinearOperator(u->u-u^3), U0)
UREF = Matrix{Matrix{Float64}}(1, Ntimes)
for j = 1:Ntimes
    T = 1.0*4^(j-1)
    UREF[1,j] = ETDRK4(SLPDE, T, round(Int, 2.0^(IMAX+1)*T))
    println("Done j = ", j)
end

errt = zeros(IMAX, Ntimes)
for j = 1:Ntimes
    for i = 1:IMAX
        errt[i,j] = relerr(UN[i,j], UREF[1,j])
        println("i = ",i,", j = ", j," and the vector 2-norm relative error: ", errt[i,j])
    end
end

h = 1./(2.0.^(1:IMAX))

clf()
loglog(h, errt[:,1], ".k", h, errt[:,2], "xk")
loglog(h, h.^4./60, "-k")
xlabel("\$h\$"); ylabel("Relative Error"); grid()
ylim((1e-13,1e-0))
legend(["\$t=1\$","\$t=4\$","\$\\mathcal{O}(h^4)\$"], loc = "lower right")
savefig("NAC_temporal_error.pdf")

IMAX = 8
Ntimes = 2

UN = Matrix{Matrix{Float64}}(IMAX, Ntimes)
for j = 1:Ntimes
    T = 1.0*4^(j-1)
    for i = 1:IMAX
        n = 64*i
        θ = (0.5:n-0.5)/n
        φ = (0:2n-2)*2/(2n-1)
        F = [cos(10*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)) for θ in θ, φ in φ]
        V = zero(F)
        A_mul_B!(V, FastTransforms.plan_analysis(F), F)
        U0 = fourier2sph(V; sketch = :none)

        SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.01, -0.5, 1.0), NonlinearOperator(u->u-u^3), U0)

        UN[i,j] = ETDRK4(SLPDE, T, round(Int, 2.0^8*T))
        println("Done i = ", i," j = ", j)
    end
end

n = 64*(2*IMAX)

θ = (0.5:n-0.5)/n
φ = (0:2n-2)*2/(2n-1)
F = [cos(10*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)) for θ in θ, φ in φ]
V = zero(F)
A_mul_B!(V, FastTransforms.plan_analysis(F), F)
U0 = fourier2sph(V; sketch = :none)

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.01, -0.5, 1.0), NonlinearOperator(u->u-u^3), U0)
UREF = Matrix{Matrix{Float64}}(1, Ntimes)
for j = 1:Ntimes
    T = 1.0*4^(j-1)
    UREF[1,j] = ETDRK4(SLPDE, T, round(Int, 2.0^9*T))
    println("Done j = ", j)
end

errs = zeros(IMAX, Ntimes)
for j = 1:Ntimes
    for i = 1:IMAX
        errs[i,j] = relerr(UN[i,j], UREF[1,j])
        println("i = ",i,", j = ", j," and the vector 2-norm relative error: ", errs[i,j])
    end
end

n = 64*(1:IMAX)-1

clf()
semilogy(n, errs[:,1], "--.k", n, errs[:,2], "--xk")
xlabel("\$n\$"); ylabel("Relative Error"); grid()
ylim((1e-13,1e-0))
legend(["\$t=1\$","\$t=4\$"], loc = "lower left")
savefig("NAC_spatial_error.pdf")


text:
@inline function abs_fpdb(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(abs_fp_dd(x))
end

@inline function negabs_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(negabs_fp_dd(x))
end

@inline function neg_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(neg_fp_dd(x))
end

@inline function inv_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(inv_fp_dd(x))
end

@inline function powr2_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(powr2_fp_dd(x))
end

@inline function powr3_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(powr3_fp_dd(x))
end

@inline function powr4_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(powr4_fp_dd(x))
end

@inline function powr5_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(powr5_fp_dd(x))
end

@inline function powr6_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(powr6_fp_dd(x))
end

@inline function root2_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(root2_fp_dd(x))
end

@inline function root3_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(root3_fp_dd(x))
end

@inline function root4_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(root4_fp_dd(x))
end

@inline function root5_fp_db(x::T) where {T<:IEEEFloat}
    throw(ErrorException("not implemented"))
end

@inline function root6_fp_db(x::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(root6_fp_dd(x))
end


text:
## script to test function checkRootPlace
## Claudia May 2016

include("../src/types.jl")
include("../src/functions.jl")

net = readTopologyLevel1("(6,((5,#H7:0.0::0.29999):10.0,(((2,#H9:0.0::0.15902):0.02332,(1,(3)#H9:0.9417::0.84098):0.12997):1.2018,(4)#H7:0.01722::0.70001):9.99428):0.24593);")
plot(net)
checkRootPlace!(net,outgroup="6")
plot(net)
checkRootPlace!(net,outgroup="2")
plot(net)
checkRootPlace!(net,outgroup="1")
plot(net)
checkRootPlace!(net,outgroup="3") ## can't do it, does not change net
plot(net)


text:
# Solution of a Class of Multistage Dynamic Optimization Problems.
# 2.Problems with Path Constraints

function van_der_pol(; store_results=false::Bool)
    datadir = nothing
    if store_results
        datadir = generate_data_subdir(@__FILE__)
    end

    # initial conditions and timepoints
    t0 = 0.0f0
    tf = 5.0f0
    u0 = [0.0f0, 1.0f0, 0.0f0]
    tspan = (t0, tf)
    Δt = 0.1f0

    function system!(du, u, p, t, controller; input=:state)
        @argcheck input in (:state, :time)

        # neural network outputs the controls taken by the system
        x1, x2, x3 = u

        if input == :state
            c1 = controller(u, p)[1]  # control based on state and parameters
        elseif input == :time
            c1 = controller(t, p)[1]  # control based on time and parameters
        end

        # dynamics of the controlled system
        x1_prime = (1 - x2^2) * x1 - x2 + c1
        x2_prime = x1
        x3_prime = x1^2 + x2^2 + c1^2

        # update in-place
        @inbounds begin
            du[1] = x1_prime
            du[2] = x2_prime
            du[3] = x3_prime
        end
    end

    # set arquitecture of neural network controller
    controller = FastChain(
        FastDense(3, 16, tanh_fast),
        FastDense(16, 16, tanh_fast),
        FastDense(16, 1),
        (x, p) -> (1.3f0 .* sigmoid_fast.(x)) .- 0.3f0,
    )

    controlODE = ControlODE(controller, system!, u0, tspan; Δt)

    θ = initial_params(controlODE.controller)

    function collocation(
        u0;
        num_supports::Integer=length(controlODE.tsteps),
        nodes_per_element::Integer=4,
        constrain_states::Bool=false,
    )
        optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
        model = InfiniteModel(optimizer)
        method = OrthogonalCollocation(nodes_per_element)
        @infinite_parameter(
            model, t in [t0, tf], num_supports = num_supports, derivative_method = method
        )

        @variables(
            model,
            begin  # "start" sets the initial guess values
                # state variables
                x[1:3], Infinite(t)
                # control variables
                c[1], Infinite(t)
            end
        )

        # initial conditions
        @constraint(model, [i = 1:3], x[i](0) == u0[i])

        # control range
        @constraint(model, -0.3 <= c[1] <= 1.0)

        if constrain_states
            @constraint(model, -0.4 <= x[1])
        end

        # dynamic equations
        @constraints(
            model,
            begin
                ∂(x[1], t) == (1 - x[2]^2) * x[1] - x[2] + c[1]
                ∂(x[2], t) == x[1]
                ∂(x[3], t) == x[1]^2 + x[2]^2 + c[1]^2
            end
        )

        @objective(model, Min, x[3](tf))

        optimize!(model)

        model |> optimizer_model |> solution_summary
        states = hcat(value.(x)...) |> permutedims
        controls = hcat(value.(c)...) |> permutedims

        return model, supports(t), states, controls
    end

    phase_time = 0.0f0
    function square_bounds(u0, arista)
        low_bounds = u0 .- repeat([arista / 2], length(u0))
        high_bounds = u0 .+ repeat([arista / 2], length(u0))
        bounds = [(l, h) for (l, h) in zip(low_bounds, high_bounds)]
        return bounds
    end

    _, states_raw, _ = run_simulation(controlODE, θ)

    start_mark = InitialState(; points=states_raw[:, 1])
    marker_path = IntegrationPath(; points=states_raw)
    final_mark = FinalState(; points=states_raw[:, end])
    phase_portrait(
        controlODE,
        θ,
        phase_time,
        square_bounds(u0, 7);
        dimension=3,
        projection=[1, 2],
        markers=[marker_path, start_mark, final_mark],
        # start_points_x, start_points_y,
        # start_points=reshape(u0 .+ repeat([-1e-4], 3), 1, 3),
        title="Initial policy",
    )

    control_profile, infopt_model, times_collocation, states_collocation, controls_collocation = collocation_preconditioner(
        u0, collocation; plot=true
    )

    θ = preconditioner(
        controlODE,
        control_profile;
        #control_range_scaling=[maximum(controls_collocation) - minimum(controls_collocation)],
    )

    plot_simulation(controlODE, θ; only=:controls)
    store_simulation("precondition", controlODE, θ; datadir)

    _, states_raw, _ = run_simulation(controlODE, θ)
    start_mark = InitialState(; points=states_raw[:, 1])
    marker_path = IntegrationPath(; points=states_raw)
    final_mark = FinalState(; points=states_raw[:, end])
    phase_portrait(
        controlODE,
        θ,
        phase_time,
        square_bounds(u0, 7);
        dimension=3,
        projection=[1, 2],
        markers=[marker_path, start_mark, final_mark],
        # start_points_x, start_points_y,
        # start_points=reshape(u0 .+ repeat([-1e-4], 3), 1, 3),
        title="Preconditioned policy",
    )

    # closures to comply with required interface
    function plotting_callback(params, loss)
        return plot_simulation(
            controlODE, params; only=:states, vars=[1], show=loss
        )
    end

    ### define objective function to optimize
    function loss(controlODE, params; kwargs...)
        sol = solve(controlODE, params; kwargs...)
        return Array(sol)[3, end]  # return last value of third variable ...to be minimized
    end
    loss(params) = loss(controlODE, params)

    @info "Training..."
    @infiltrate
    result = sciml_train(
        loss,
        θ,
        LBFGS(; linesearch=BackTracking());
        # cb=plotting_callback,
        allow_f_increases=true,
    )
    θ = result.minimizer

    store_simulation(
        "unconstrained",
        controlODE,
        θ;
        datadir,
        metadata=Dict(:loss => loss(θ), :constraint => "none"),
    )

    ### now add state constraint x1(t) > -0.4 with
    function penalty_loss(controlODE, params; α=10.0f0, kwargs...)
        # integrate ODE system
        sol = Array(solve(controlODE, params; kwargs...))
        fault = min.(sol[1, 1:end] .+ 0.4f0, 0.0f0)
        penalty = α * Δt * sum(fault .^ 2)  # quadratic penalty
        return sol[3, end] + penalty
    end
    @info "Enforcing constraints..."
    penalty_coefficients = [10.0f0, 10.0f1, 10.0f2, 10.0f3]
    prog = Progress(
        length(penalty_coefficients);
        desc="Fiacco-McCormick iterations",
        dt=0.5,
        showspeed=true,
        enabled=true,
    )
    for α in penalty_coefficients

        # function plotting_callback(params, loss)
        #     return plot_simulation(
        #         controlODE, params; only=:states, vars=[1], show=loss
        #     )
        # end

        # closures to comply with interface
        penalty_loss_(params) = penalty_loss(controlODE, params; α)

        # plot_simulation(
        #     controlODE,
        #     θ;
        #     only=:controls,
        #     show=penalty_loss_(θ),
        # )

        result = sciml_train(
            penalty_loss_,
            θ,
            LBFGS(; linesearch=BackTracking(; iterations=20));
            iterations=50,
            allow_f_increases=true,
            # cb=plotting_callback,
        )
        θ = result.minimizer
        next!(prog; showvalues=[(:α, α), (:loss, penalty_loss_(θ))])
    end
    θ_opt = result.minimizer
    optimal = result.minimum

    # penalty_loss(result.minimizer, constrained_prob, tsteps; α=penalty_coefficients[end])
    plot_simulation(controlODE, θ_opt; only=:controls)

    store_simulation(
        "constrained",
        controlODE,
        θ_opt;
        datadir,
        metadata=Dict(
            :loss =>
                penalty_loss(controlODE, θ_opt; α=penalty_coefficients[end]),
            :constraint => "quadratic x2(t) > -0.4",
        ),
    )

    _, states_opt, _ = run_simulation(controlODE, θ_opt)
    start_mark = InitialState(; points=states_opt[:, 1])
    marker_path = IntegrationPath(; points=states_opt)
    final_mark = FinalState(; points=states_opt[:, end])
    shader = ShadeConf(; indicator=function (x, y)
        if x > -0.4
            return true
        end
        return false
    end)
    phase_portrait(
        controlODE,
        θ_opt,
        phase_time,
        square_bounds(u0, 7);
        shader,
        dimension=3,
        projection=[1, 2],
        markers=[marker_path, start_mark, final_mark],
        # start_points_x, start_points_y,
        # start_points=reshape(u0 .+ repeat([-1e-4], 3), 1, 3),
        title="Optimized policy",
    )

    # u0 = [0f0, 1f0, 0f0]
    perturbation_specs = [
        (variable=1, type=:positive, scale=1.0f0, samples=8, percentage=2.0f-2)
        (variable=2, type=:negative, scale=1.0f0, samples=8, percentage=2.0f-2)
        (variable=3, type=:positive, scale=20.0f0, samples=8, percentage=2.0f-2)
    ]
    constraint_spec = ConstRef(; val=-0.4, direction=:horizontal, class=:state, var=1)

    plot_initial_perturbations_collocation(
        controlODE,
        θ_opt,
        perturbation_specs,
        collocation;
        refs=[constraint_spec],
        storedir=generate_data_subdir(@__FILE__),
    )
    return optimal
end  # wrap


text:
@testset "\e[34mProgress" begin
    import Term.progress: track, ProgressBar, update, track, start, stop

    @test_nothrow begin
        for i in track(1:10)
            sleep(0.001)
        end
    end

    @test_nothrow begin
        vec = collect(1:10)
        for i in track(vec; description="[red]")
            sleep(0.001)
        end
    end

    for level in (:minimal, :default, :extensive)
        @test_nothrow begin
            for i in track(1:10; columns=level)
                sleep(0.001)
            end
        end
    end

    for redirect in (true, false)
        @test_nothrow begin
            for i in track(1:10; redirectstdout=redirect)
                sleep(0.001)
                print("test")
            end
        end
    end

    import Term.progress: DescriptionColumn, BarColumn, PercentageColumn
    cols = [DescriptionColumn("test"), BarColumn(), PercentageColumn()]
    @test_nothrow begin
        for i in track(1:10; columns=cols)
            sleep(0.001)
            print("test")
        end
    end

end


text:
struct Noop <: Codec end

const NoopStream{S} = TranscodingStream{Noop,S} where S<:IO

NoopStream(stream::IO; kwargs...) = TranscodingStream(Noop(), stream; kwargs...)

function TranscodingStream(codec::Noop, stream::IO;
                           bufsize::Integer=DEFAULT_BUFFER_SIZE,
                           sharedbuf::Bool=(stream isa TranscodingStream))
    checkbufsize(bufsize)
    checksharedbuf(sharedbuf, stream)
    if sharedbuf; buffer = stream.state.buffer1
    else; buffer = Buffer(bufsize)
    end
    TranscodingStream(codec, stream, State(buffer, buffer))
end

"""
    position(stream::NoopStream)

Get the current position of `stream`.

Note that this method may return a wrong position when
- some data have been inserted by `TranscodingStreams.unread`, or
- the position of the wrapped stream has been changed outside of this package.
"""
Base.position(stream::NoopStream) = position(stream.stream) - buffersize(stream.state.buffer1)

Base.seek(stream::NoopStream, pos::Integer) = (seek(stream.stream, pos); initbuffer!(stream.state.buffer1))
Base.seekstart(stream::NoopStream) = (seekstart(stream.stream); initbuffer!(stream.state.buffer1))
Base.seekend(stream::NoopStream) = (seekend(stream.stream); initbuffer!(stream.state.buffer1))

function Base.unsafe_read(stream::NoopStream, output::Ptr{UInt8}, nbytes::UInt)
    changemode!(stream, :read)
    buffer = stream.state.buffer1
    p = output
    p_end = output + nbytes
    while p < p_end && !eof(stream)
        if buffersize(buffer) > 0
            m = min(buffersize(buffer), p_end - p)
            copydata!(p, buffer, m)
        else
            # directly read data from the underlying stream
            m = p_end - p
            Base.unsafe_read(stream.stream, p, m)
        end
        p += m
    end
    if p < p_end && eof(stream); throw(EOFError()) end
end

function Base.unsafe_write(stream::NoopStream, input::Ptr{UInt8}, nbytes::UInt)
    changemode!(stream, :write)
    buffer = stream.state.buffer1
    if marginsize(buffer) ≥ nbytes;     copydata!(buffer, input, nbytes); Int(nbytes)
    else;                               flushbuffer(stream); unsafe_write(stream.stream, input, nbytes)
    end
end

Base.transcode(::Type{Noop}, data::ByteData) = transcode(Noop(), data)
# Copy data because the caller may expect the return object is not the same as from the input.
Base.transcode(::Noop, data::ByteData) = Vector{UInt8}(data)

# Stats
# -----

function stats(s::NoopStream)
    state,mode = s.state, state.mode
    @checkmode (:idle, :read, :write)
    buffer = state.buffer1
    @assert buffer == s.state.buffer2
    if mode == :idle;       consumed = supplied = 0
    elseif mode == :read;   supplied = buffer.total; consumed = supplied-buffersize(buffer)
    elseif mode == :write;  supplied = buffer.total+buffersize(buffer); consumed = buffer.total
    else;                   assert(false)
    end
    return Stats(consumed, supplied, supplied, supplied)
end

# Buffering
# ---------
#
# These methods are overloaded for the `Noop` codec because it has only one
# buffer for efficiency.

function fillbuffer(s::NoopStream)
    changemode!(s, :read)
    buffer = s.state.buffer1
    @assert buffer === s.state.buffer2
    if s.stream isa TranscodingStream && buffer === s.stream.state.buffer1
        # Delegate the operation when buffers are shared.
        return fillbuffer(s.stream)
    end
    nfilled::Int = 0
    while buffersize(buffer) == 0 && !eof(s.stream)
        makemargin!(buffer, 1)
        nfilled += readdata!(s.stream, buffer)
    end
    buffer.total += nfilled
    nfilled
end

function flushbuffer(s::NoopStream, all::Bool=false)
    changemode!(s, :write)
    buffer = s.state.buffer1
    @assert buffer === s.state.buffer2
    nflushed::Int = 0
    if all
        while buffersize(buffer) > 0; nflushed += writedata!(s.stream, buffer) end
    else
        nflushed += writedata!(s.stream, buffer)
        makemargin!(buffer, 0)
    end
    buffer.total += nflushed
    return nflushed
end

flushuntilend(s::NoopStream) = s.state.buffer1.total += writedata!(s.stream, s.state.buffer1)


text:
# For recording error information
#

# This type is for recording the minimum and maximum of some online sequence
# alongside the average magnitude of the input. This is specifically for
# recording propagation errors and negative or imagine probabilities. In that
# case it's much more interesting to know the average magnitude instead of true
# average.
mutable struct MagnitudeStats
    max::Float64
    min::Float64
    sum::Float64
    count::Int64
end

MagnitudeStats() = MagnitudeStats(-Inf, +Inf, 0.0, 0)

function Base.push!(stat::MagnitudeStats, value)
    v = log10(abs(value))
    stat.max = max(stat.max, v)
    stat.min = min(stat.min, v)
    stat.sum += v
    stat.count += 1
end

Base.min(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.min) : 0.0
Base.max(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.max) : 0.0
Statistics.mean(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.sum / s.count) : 0.0
Base.length(s::MagnitudeStats) = s.count

function Base.show(io::IO, s::MagnitudeStats)
    println(io, "MagnitudeStats: ($(s.count) Values)")
    println(io, "\tmin = $(min(s))")
    println(io, "\tmean = $(mean(s))")
    println(io, "\tmax = $(max(s))")
end

"""
    DQMCAnalysis

Analysis data of determinant quantum Monte Carlo (DQMC) simulations.
"""
@with_kw mutable struct DQMCAnalysis
    th_runtime::Float64 = 0.0
    me_runtime::Float64 = 0.0
    imaginary_probability::MagnitudeStats = MagnitudeStats()
    negative_probability::MagnitudeStats = MagnitudeStats()
    propagation_error::MagnitudeStats = MagnitudeStats()
end


text:
"""
    is_segment(node)

Checks if a node (n) has only one child (n+1). This is usefull to simplify a complex mtg to
become an mtg with nodes only at the branching points, has it is often measured on the field.

The function also takes care of passing the link of the node (n) to its child (n+1) if the
node (n) branches or decompose its parent (n-1). This allows a conservation of the relationships
as they previously were in the mtg.

See [`delete_nodes!`](@ref) for an example of application.
"""
function is_segment!(node)
    if !isleaf(node) && !isroot(node) && !isroot(node.parent) && length(node.children) == 1
        # We keep the root and the leaves, but want to delete the nodes with no branching.
        # We recognise them because they have only one child. Also we want to keep the very
        # first node even if it has only one child.

        # If there is only one child but it is branching:
        if node[1].MTG.link == "+"
            return false
        end

        # If it's a node that branches, set its unique child as the branching node instead:
        if node.MTG.link == "+"
            node_MTG = node[1].MTG
            node[1].MTG = typeof(node_MTG)("+", node_MTG.symbol, node_MTG.index, node_MTG.scale)
        end

        # If it's a node that decompose ("/"), set its unique child as the decomposing node:
        if node.MTG.link == "/"
            node_MTG = node[1].MTG
            node[1].MTG = typeof(node_MTG)("/", node_MTG.symbol, node_MTG.index, node_MTG.scale)
        end

        # And return true to delete it:
        return true
    end
    return false
end

"""
    filter_fun_nothing(filter_fun, ignore_nothing, attr_keys)

Returns a new filtering function that adds a filter on the keys value for `nothing` if
`ignore_nothing` is `true`
"""
function filter_fun_nothing(filter_fun, ignore_nothing, attr_keys)
    # Change the filtering function if we also want to remove nodes with nothing values.
    if ignore_nothing
        if filter_fun !== nothing
            filter_fun_ =
                function (node)
                    all([unsafe_getindex(node, i) !== nothing for i in attr_keys]) && filter_fun(node)
                end
        else
            filter_fun_ =
                function (node)
                    all([unsafe_getindex(node, i) !== nothing for i in attr_keys])
                end
        end
    else
        filter_fun_ = filter_fun
    end

    filter_fun_
end

function filter_fun_nothing(filter_fun, ignore_nothing, attr_keys::T) where {T<:Union{Symbol,String}}
    filter_fun_nothing(filter_fun, ignore_nothing, [attr_keys])
end


text:
using Documenter
include("../src/LightGraphs.jl")
using LightGraphs

# index is equal to the README for the time being
cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/index.md"); remove_destination=true)

# same for contributing and license
cp(normpath(@__FILE__, "../../CONTRIBUTING.md"), normpath(@__FILE__, "../src/contributing.md"); remove_destination=true)
cp(normpath(@__FILE__, "../../LICENSE.md"), normpath(@__FILE__, "../src/license.md"); remove_destination=true)

makedocs(
    modules     = [LightGraphs],
    format      = :html,
    sitename    = "LightGraphs",
    doctest     = false,
    pages       = Any[
        "Getting Started"                   => "index.md",
        "Basic Functions"                   => "basicmeasures.md",
        "Operators"                         => "operators.md",
        "Path and Traversal"                => "pathing.md",
        "Distance"                          => "distance.md",
        "Centrality Measures"               => "centrality.md",
        "Linear Algebra"                    => "linalg.md",
        "Matching"                          => "matching.md",
        "Community Structures"              => "community.md",
        "Degeneracy"                        => "degeneracy.md",
        "Flow and Cut"                      => "flowcut.md",
        "Graph Generators"                  => "generators.md",
        "Reading / Writing Graphs"          => "persistence.md",
        "Integration with other packages"   => "integration.md",
        "Contributing"                      => "contributing.md",
        "License Information"               => "license.md"
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/JuliaGraphs/LightGraphs.jl.git",
    target      = "build",
    julia       = "0.6",
    osname      = "linux"
)

rm(normpath(@__FILE__, "../src/index.md"))
rm(normpath(@__FILE__, "../src/contributing.md"))
rm(normpath(@__FILE__, "../src/license.md"))