text:
function DFPermute(A::DF,v::Vector{Int})
    if length(A.var)!=length(v)
        error("Permute vector's dimension mismatch!")
    end

    valuespace = reshape(A.value, A.card...)
    newvaluespace = permutedims(valuespace, v)
    var = A.var[v]
    card = A.card[v]
    DF(var,card,newvaluespace[:])
end

function DFProduct(A::DF, B::DF)
    common_var = intersect(A.var, B.var)
    A_common_idx = indexin(common_var, A.var)
    B_common_idx = indexin(common_var, B.var)
    if A.card[A_common_idx] != B.card[B_common_idx]
        error("Inputs' common variables have mismatch cardinality!")
    end

    # Set the variables of C
    Cvar = union(A.var, B.var)

    Cvar_length = length(Cvar)

    mapA = indexin(A.var, Cvar)
    mapB = indexin(B.var, Cvar)

    # Set the cardinality of variables in C
    Ccard = zeros(Int, Cvar_length)
    Ccard[mapA] = A.card
    Ccard[mapB] = B.card

    Aarr0 = indexin(Cvar, A.var)
    Barr0 = indexin(Cvar, B.var)
    Aarr1 = Aarr0[Aarr0.>0]
    Barr1 = Barr0[Barr0.>0]

    A1 = DFPermute(A, Aarr1)
    B1 = DFPermute(B, Barr1)

    A1card = ones(Int, Cvar_length)
    B1card = ones(Int, Cvar_length)
    A1card[mapA] = A.card
    B1card[mapB] = B.card

    Avaluespace = reshape(A1.value, A1card...)
    Bvaluespace = reshape(B1.value, B1card...)

    Cvaluespace = Avaluespace .* Bvaluespace

    DF(Cvar, Ccard, Cvaluespace[:])
end

*(A::DF, B::DF) = DFProduct(A::DF, B::DF)

function DFMargin(A::DF, Remove_var::Vector{ASCIIString}, Remain_var::Vector{ASCIIString}, Remove_dims::Vector{Int}, Remain_dims::Vector{Int})
    Remain_card = A.card[Remain_dims]

    valuespace = reshape(A.value, A.card...)
    permute_dims = [Remain_dims,Remove_dims]
    permuted_valuespace = permutedims(valuespace, permute_dims)

    squeeze_dims = [length(Remain_dims)+1:length(permute_dims)]

    sumvaluespace = sum(permuted_valuespace, squeeze_dims)
    Remain_valuespace = squeeze(sumvaluespace, squeeze_dims)

    DF(Remain_var, Remain_card, Remain_valuespace[:])
end

function DFMarginDrop(A::DF, Remove_var::Vector{ASCIIString})
    Remove_dims = indexin(Remove_var, A.var)
    if any(Remove_dims==0)
        error("Wrong variable!")
    end

    Remain_var = symdiff(A.var, Remove_var)
    Remain_dims = indexin(Remain_var, A.var)

    DiscreteFactor.DFMargin(A, Remove_var, Remain_var, Remove_dims, Remain_dims)
end

function DFMarginKeep(A::DF, Remain_var::Vector{ASCIIString})
    Remain_dims = indexin(Remain_var, A.var)
    if any(Remain_dims==0)
        error("Wrong variable!")
    end

    Remove_var = symdiff(A.var, Remain_var)
    Remove_dims = indexin(Remove_var, A.var)

    DiscreteFactor.DFMargin(A, Remove_var, Remain_var, Remove_dims, Remain_dims)
end

function DFReduce(A::DF, Reduce_var::Vector{ASCIIString}, Reduce_idx::Vector{Int})
    Reduce_dims = indexin(Reduce_var, A.var)
    if any(Reduce_dims.==0)
        error("Wrong variable!")
    end
    if any(A.card[Reduce_dims].<Reduce_idx)
        error("Index larger than cardinality!")
    end

    Bvar = deepcopy(A.var)
    Bcard = deepcopy(A.card)
    Bvar[Reduce_dims] = Bvar[Reduce_dims].*["_{$i}" for i in Reduce_idx]
    Bcard[Reduce_dims] = 1

    valuespace = reshape(A.value, A.card...)

    for i = 1:length(Reduce_dims)
        valuespace = slicedim(valuespace, Reduce_dims[i], Reduce_idx[i])
    end
    DF(Bvar, Bcard, valuespace[:])
end

function DFNormalize(A::DF)
    Bvar = deepcopy(A.var)
    Bcard = deepcopy(A.card)
    tmpvalue = deepcopy(A.value)
    Z =  sum(tmpvalue)
    Bvalue = tmpvalue./Z
    DF(Bvar, Bcard, Bvalue)
end


text:
function Base.summary(io::IO, H::KPHamiltonian{D}) where D
    print(io, "KPHamiltonian{", D, "} for ", irdim(H.lgir), "D irrep (",
              Crystalline.formatirreplabel(label(H.lgir)),
              ") with ", length(H.cs), " basis elements:")
end

# ---------------------------------------------------------------------------------------- #

function format_matrix_element(v, i, j)
    if iszero(v)
        return "·"
    else
        if real(v) == v
            return format_scalar(real(v))
        elseif imag(v)*1im == v
            return format_scalar(imag(v))*"i"
        else
            return format_scalar(v)
        end
    end
end

function format_scalar(v)
    if isinteger(v)
        return string(round(Int, v))
    else
        return string(v)
    end
end

function print_matrix(io, A)
    pretty_table(io, A;
        tf = tf_matrix, noheader = true,
        formatters = format_matrix_element,
        highlighters = Highlighter((data,i,j) -> iszero(data[i,j]), foreground = :dark_gray)
    )
end

function Base.show(io::IO, ::MIME"text/plain", H::KPHamiltonian{D}) where D
    summary(io, H)
    println(io)

    N = length(H.hs)
    hs_rowstrs = split.(sprint.(print_matrix, H.hs; context=:color=>true), '\n')
    coefstrs = map(H.cs) do cₐ
        map(1:length(H.hs)) do n
            io′ = IOBuffer()
            write(io′, '(')
            count = 0
            for d in 1:D
                v = cₐ[d][n]
                if !iszero(v)
                    print(io′, Crystalline.signaschar(v))
                    if !(abs(v) ≈ one(v))
                        print(io′, round(abs(v), digits=3))
                    end
                    print(io′, "k", d == 1 ? "₁" : d == 2 ? "₂" : d == 3 ? "₃" : error("unexpected dimension $d"))
                    count += 1
                end
            end
            print(io′, ')')
            String(take!(io′))
        end
    end

    max_tw = displaysize(io)[2]
    cntr_row = div(irdim(H)+2, 2, RoundUp)
    for a in 1:length(H.cs)
        for row in 1:irdim(H)+2
            printstyled(io, row == 1 ? Crystalline.subscriptify(string(a))*"₎ " : "   ", color=:light_black)
            first_nonzero_term = true
            tw = 2
            for (n, h_rowstrs) in enumerate(hs_rowstrs)
                coefstrs[a][n] == "()" && continue
                if !first_nonzero_term
                    print(io, row == cntr_row ? " + " : "   ")
                end
                # calculate how much space the next matrix would take to display (need to
                # remove ANSI codes from `h_rowstrs` for this bit, since they are included
                # in `textwidth`'s count)
                tw += textwidth(replace(h_rowstrs[row], r"\e[^m]*m"=>"")) +
                      textwidth(coefstrs[a][n]) + 3*(n≠N)
                if tw > max_tw
                    row == cntr_row && print(io, "…")
                    break
                end
                print(io, h_rowstrs[row])
                first_nonzero_term = false
                if row == cntr_row
                    print(io, coefstrs[a][n])
                else
                    print(io, " "^textwidth(coefstrs[a][n]))
                end
            end
            println(io)
        end
    end
end

ext: jl


text:
map = open("local_geology.txt") do f
    readlines(f)
end

slopes = [(1,1), (3,1), (5,1), (7,1), (1,2)]
println(prod([count(tuple -> (((tuple[1] + 1) % slope[2] == 0) && (tuple[2][Int((tuple[1] - 1) * slope[1]/slope[2]) % length(tuple[2]) + 1] == '#')), enumerate(map)) for slope in slopes]))

# slope = (1,2)
# for (y, row) in enumerate(map)
#     if (y+1)  % slope[2] != 0
#         println(row)
#         continue
#     end
#     x = Int(((y-1)*slope[1]/slope[2]) % length(map) + 1)
#     if row[x] == '#'
#         rep = "X"
#     elseif row[x] == '.'
#         rep = "O"
#     end
#     println(row[begin:x-1], rep, row[x+1:end])
# end


text:
function show(io::IO, mC::MapSUnitModUnitGrpFacElem)
  @show_name(io, mC)
  io = IOContext(io, :compact => true)
  println(io, "SUnits (in factored form) mod Units map of ")
  show(io, codomain(mC))
  println(io, "for $(mC.idl)")
end

#Plan:
# find x_i s.th. I[i]*x[i] is FB-smooth
#  find T sth. T R = (I[i]*x[i])^d
#  saturate T|-d??

@doc Markdown.doc"""
    sunit_mod_units_group_fac_elem(I::Array{NfOrdIdl, 1}) -> GrpAb, Map
For an array $I$ of (coprime prime) ideals, find the $S$-unit group defined
by $I$, ie. the group of non-zero field elements which are only divisible
by ideals in $I$ modulo the units of the field.
The map will return elements in factored form.
"""
function sunit_mod_units_group_fac_elem(I::Array{NfOrdIdl, 1})
  #deal with trivial case somehow!!!
  @assert length(I) > 0
  O = order(I[1])
  I_in = I

  @vprint :ClassGroup 1 "calling sunit_mod_units_group_fac_elem with $(length(I)) ideals\n"

  c = _get_ClassGrpCtx_of_order(O, false)
  if c == nothing
    L = lll(maximal_order(nf(O)))
    class_group(L)
    c = _get_ClassGrpCtx_of_order(L)
    I = map(IdealSet(L), I)
  end
  module_trafo_assure(c.M)
  H = c.M.basis
  T = c.M.trafo

  U = Array{FacElem{nf_elem, Nemo.AnticNumberField}, 1}()

  X = Array{nf_elem, 1}()

  rr = sparse_matrix(FlintZZ)

  # To track the valuation of the S-units
  vals_of_rels = SRow{fmpz}[]

  @vprint :ClassGroup 1 "finding relations ...\n"
  @vtime :ClassGroup 1 for (i, A) = enumerate(I)
    @vprint :ClassGroup 2 "doin' $(i)/$(length(I)):\n$A\n"
    @vtime :ClassGroup 2 x, r = class_group_ideal_relation(A, c)
# TODO: write == for Idl and FracIdl
#    @assert prod([c.FB.ideals[p]^Int(v) for (p,v) = r]) == x*A
    push!(X, x)
    push!(rr, r)
    v = SRow(FlintZZ)
    # We only track the valuation of the prime ideals in S.
    # Even though S might intersect the class group factor base
    # non-trivially, this should still be correct.
    push!(vals_of_rels, sparse_row(FlintZZ, [(i, fmpz(-1))], sort = false))
  end

  @vprint :ClassGroup 1 "... done\n"

  @vprint :ClassGroup 1 "solving...\n"
  @vtime :ClassGroup 1 R, d = solve_ut(H, rr)
  Rd = hcat(d*identity_matrix(SMat, FlintZZ, nrows(R)), fmpz(-1)*R)
  @vprint :ClassGroup 1 ".. done, now saturating ...\n"
  @vtime :ClassGroup 1 S = hnf(saturate(Rd))
  @vprint :ClassGroup 1 " done\n"
  S1 = sub(S, 1:nrows(S), 1:nrows(S))
  S2 = sub(S, 1:nrows(S), (nrows(S) + 1):ncols(S))
  @assert nrows(S1) == nrows(S2) && nrows(S1) == nrows(S)

  g = vcat(c.R_gen, c.R_rel)

  valuations = SRow{fmpz}[]

  for s = 1:S.r
    rs = zeros(fmpz, c.M.bas_gens.r + c.M.rel_gens.r)
    for (p, v) = S2[s]
      rs[p] = v
    end

    for i in length(T):-1:1
      apply_right!(rs, T[i])
    end

    _val_vec = sparse_row(FlintZZ)

    e = FacElem(g, rs)
    for (p, v) = S1[s]
      _val_vec = _val_vec + v * vals_of_rels[p]
      if haskey(e.fac, X[p])
        e.fac[X[p]] += v
      else
        e.fac[X[p]] = v
      end
    end

    _val_vec = -_val_vec
    inv!(e)

    push!(valuations, _val_vec)
    push!(U, e)  # I don't understand the inv
  end
  @vprint :ClassGroup 1 "reducing mod units\n"
  @vtime :ClassGroup 1 U = reduce_mod_units(U, _get_UnitGrpCtx_of_order(order(c)))
  @vprint :ClassGroup 1 "Done!\n"

  #for j in 1:length(I)
  #  @assert (O(evaluate(U[j]))*O) == prod(I[i]^Int(valuations[j][i]) for i in 1:length(I))
  #end

  C = abelian_group(fmpz[0 for i=U])
  r = MapSUnitModUnitGrpFacElem()
  r.idl = I_in

  local exp
  let U = U
    function exp(a::GrpAbFinGenElem)
      b = U[1]^a.coeff[1, 1]
      for i = 2:length(U)
        if iszero(a.coeff[1, i])
          continue
        end
        mul!(b, b, U[i]^a.coeff[1, i])
      end
      return b
    end
  end

  local log
  let I = I, S1 = S1, C = C

    function log(a::FacElem{nf_elem, AnticNumberField})
      b = SRow{fmpz}()
      for i=1:length(I)
        v = valuation(a, I[i])
        if v != 0
          push!(b.pos, i)
          push!(b.values, v)
        end
      end
      s, d = solve_ut(S1, b)
      @assert d == 1  # this would indicate element is not in group...
      c = zeros(fmpz, length(I))
      for (p,v) = s
        c[p] = v
      end
      return C(c)
    end

    function log(a::nf_elem)
      return log(FacElem(a))
    end
  end

  r.header = MapHeader(C, FacElemMon(nf(O)), exp, log)
  r.valuations = valuations

  return C, r
end

function show(io::IO, mC::MapSUnitGrpFacElem)
  @show_name(io, mC)
  print(io, "SUnits (in factored form) map of ")
  show(IOContext(io, :compact => true), codomain(mC))
  println(io, " for S of length ", length(mC.idl))
  if mC.isquotientmap != -1
    println(io, " This is the quotient modulo $(mC.isquotientmap)")
  end
end

@doc Markdown.doc"""
    sunit_group_fac_elem(I::Array{NfOrdIdl, 1}) -> GrpAb, Map
For an array $I$ of (coprime prime) ideals, find the $S$-unit group defined
by $I$, ie. the group of non-zero field elements which are only divisible
by ideals in $I$.
The map will return elements in factored form.
"""
function sunit_group_fac_elem(I::Array{NfOrdIdl, 1})
  O = order(I[1])
  S, mS = sunit_mod_units_group_fac_elem(I)
  U, mU = unit_group_fac_elem(O)

  G = abelian_group(vcat(U.snf, S.snf))

  r = MapSUnitGrpFacElem()
  r.valuations = Vector{SRow{fmpz}}(undef, ngens(G))
  for i = 1:ngens(U)
    r.valuations[i] = SRow{fmpz}()
  end
  for i = 1:ngens(S)
    r.valuations[i+ngens(U)] = mS.valuations[i]
  end
  r.idl = I

  local exp
  let mU = mU, mS = mS, U = U, G = G
    function exp(a::GrpAbFinGenElem)
      return image(mU, GrpAbFinGenElem(U, sub(a.coeff, 1:1, 1:length(U.snf))))*
             image(mS, GrpAbFinGenElem(S, sub(a.coeff, 1:1, length(U.snf)+1:length(G.snf))))
    end
  end

  local log
  let mS = mS, mU = mU, G = G
    function log(a::FacElem{nf_elem, AnticNumberField})
      a1 = preimage(mS, a)
      a2 = a*inv(image(mS, a1))
      #     @assert isunit(O(evaluate(a2)))
      a3 = preimage(mU, a2)
      return GrpAbFinGenElem(G, hcat(a3.coeff, a1.coeff))
    end

    function log(a::nf_elem)
      return log(FacElem(a))
    end
  end

  r.header = MapHeader(G, FacElemMon(nf(O)), exp, log)

  return G, r
end

function show(io::IO, mC::MapSUnitGrp)
  @show_name(io, mC)
  print(io, "SUnits  map of ")
  show(IOContext(io, :compact => true), codomain(mC))
  println(io, " for $(mC.idl)")
end

@doc Markdown.doc"""
    sunit_group(I::Array{NfOrdIdl, 1}) -> GrpAb, Map
For an array $I$ of (coprime prime) ideals, find the $S$-unit group defined
by $I$, ie. the group of non-zero field elements which are only divisible
by ideals in $I$.
"""
function sunit_group(I::Array{NfOrdIdl, 1})
  O = order(I[1])
  G, mG = sunit_group_fac_elem(I)

  r = MapSUnitGrp()
  r.idl = I

  local exp
  let mG = mG
    function exp(a::GrpAbFinGenElem)
      return evaluate(image(mG, a))
    end
  end

  local log
  let mG = mG
    function log(a::nf_elem)
      return preimage(mG, FacElem(a))
    end
  end

  r.header = MapHeader(G, nf(O), exp, log)

  return G, r
end


text:
export i_inv!, i_affine!

"""
    i_inv!(out!, A)

Get the inverse of `A`.

```note!!!
this function is implemented as a primitive.
```
"""
@i function i_inv!(out!::AbstractMatrix{T}, A::AbstractMatrix{T}) where T
    @invcheckoff invA ← inv(A)
    out! .+= invA
    @invcheckoff invA → inv(A)
end

@i function i_inv!(out!::AbstractMatrix{T}, A::AbstractMatrix{T}) where T<:GVar
    @routine @invcheckoff begin
        invA ← inv(value.(A))
        gA ← -transpose(invA) * grad(out!) * transpose(invA)
    end
    for i=1:length(out!)
        (out![i] |> value) -= invA[i]
    end
    for i=1:length(A)
        (A[i] |> grad) -= gA[i]
    end
    ~@routine
end

@i function :(-=)(det)(out!::T, A::AbstractMatrix{T}) where T<:GVar
    @routine @invcheckoff begin
        vA ← value.(A)
        detA ← det(vA)
        gA ← detA * grad(out!) * transpose(inv(vA))
    end
    (out! |> value) -= detA
    for i=1:length(A)
        (A[i] |> grad) += gA[i]
    end
    ~@routine
end

@i function :(-=)(logdet)(out!::T, A::AbstractMatrix{T}) where T<:GVar
    @routine @invcheckoff begin
        gA ← grad(out!) * transpose(inv(value.(A)))
    end
    (out! |> value) -= det(A |> grad)
    for i=1:length(A)
        (A[i] |> grad) += gA[i]
    end
    ~@routine
end

"""
    i_affine!(y!, W, b, x)

`affine!` transformation `y! += W*x + b`.
"""
@i function i_affine!(y!::AbstractVector{T}, W::AbstractMatrix{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T
    @safe @assert size(W) == (length(y!), length(x)) && length(b) == length(y!)
    @invcheckoff for j=1:size(W, 2)
        for i=1:size(W, 1)
            @inbounds y![i] += W[i,j]*x[j]
        end
    end
    @invcheckoff for i=1:size(W, 1)
        @inbounds y![i] += b[i]
    end
end


text:
"""
    legendre_p(ν::arb, μ::arb, z::arb)
> Returns the associated Legendre function of the first kind evaluated
  for degree ν, order μ, and argument z. When μ is zero, this reduces
  to the Legendre polynomial P_ν(z).
"""
function legendre_p(ν::arb, μ::arb, z::arb)
    res = parent(z)()
    ccall(
        (:arb_hypgeom_legendre_p, Nemo.libarb),
        Nothing,
        (Ref{arb}, Ref{arb}, Ref{arb}, Ref{arb}, Cint, Clong),
        res,
        ν,
        μ,
        z,
        0,
        parent(z).prec,
    )
    return res
end

function legendre_p_safe(ν::arb, μ::arb, z::arb)
    res = parent(z)(NaN)
    prec = parent(z).prec
    while !isfinite(res) && prec <= 16 * parent(z).prec
        ccall(
            (:arb_hypgeom_legendre_p, Nemo.libarb),
            Nothing,
            (Ref{arb}, Ref{arb}, Ref{arb}, Ref{arb}, Cint, Clong),
            res,
            ν,
            μ,
            z,
            0,
            prec,
        )
        prec *= 2
    end

    return res
end

function bessel_j(ν::arb, z::arb)
    res = parent(z)()
    ccall(
        (:arb_hypgeom_bessel_j, Nemo.libarb),
        Nothing,
        (Ref{arb}, Ref{arb}, Ref{arb}, Clong),
        res,
        ν,
        z,
        parent(z).prec,
    )
    return res
end

function bessel_y(ν::arb, z::arb)
    res = parent(z)()
    ccall(
        (:arb_hypgeom_bessel_y, Nemo.libarb),
        Nothing,
        (Ref{arb}, Ref{arb}, Ref{arb}, Clong),
        res,
        ν,
        z,
        parent(z).prec,
    )
    return res
end

"""
    getinterval(x::arb)
> Return an interval [a,b] containing the ball x.
"""
function getinterval(x::arb)
    getinterval(arb, x)
end

"""
    getinterval(::Type{arb}, x::arb)
> Return an interval [a,b] containing the ball x.
"""
function getinterval(::Type{arb}, x::arb)
    a, b = x.parent(), x.parent()
    a_mid = ccall((:arb_mid_ptr, Nemo.libarb), Ptr{Nemo.arf_struct}, (Ref{arb},), a)
    b_mid = ccall((:arb_mid_ptr, Nemo.libarb), Ptr{Nemo.arf_struct}, (Ref{arb},), b)
    ccall(
        (:arb_get_interval_arf, Nemo.libarb),
        Cvoid,
        (Ptr{Nemo.arf_struct}, Ptr{Nemo.arf_struct}, Ref{arb}, Clong),
        a_mid,
        b_mid,
        x,
        x.parent.prec,
    )

    (a, b)
end

"""
    getinterval(::Type{BigFloat}, x::arb)
> Return an interval [a,b] containing the ball x.
"""
function getinterval(::Type{BigFloat}, x::arb)
    a, b = BigFloat(), BigFloat()
    ccall(
        (:arb_get_interval_mpfr, Nemo.libarb),
        Cvoid,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{arb}),
        a,
        b,
        x,
    )

    (a, b)
end


text:
module Cells

# A cell -- a container holding data, no intelligence here

using Defs, Norms

import Base.show
import Norms.norm

export Cell
export axpy, norm, analytic, initial, boundary, error, rhs

immutable Cell
    x::Float64
    u::Float64
    rho::Float64
    v::Float64
end
Cell() = Cell(0.0, 0.0, 0.0, 0.0)

show(io::IO, c::Cell) =
    print(io, "cell: x=", c.x, " u=", c.u, " rho=", c.rho, " v=", c.v)

# Linear combination
axpy(a::Float64, x::Cell, y::Cell) =
    Cell(a*x.x + y.x, a*x.u + y.u, a*x.rho + y.rho, a*x.v + y.v)

# Norm
norm(c::Cell) = norm(c.u) + norm(c.rho) + norm(c.v)

# Analytic solution
analytic(t::Float64, x::Float64) =
    Cell(x,
         sin(2pi * t) * sin(2pi * x),
         2pi * cos(2pi * t) * sin(2pi * x),
         2pi * sin(2pi * t) * cos(2pi * x))

# Initial condition
initial(t::Float64, x::Float64) = analytic(t, x)

# Boundary condition
boundary(t::Float64, x::Float64) = analytic(t, x)

# Error
error(c::Cell, t::Float64) = axpy(-1.0, analytic(t, c.x), c)

# RHS
rhs(c::Cell, cm::Cell, cp::Cell) =
    Cell(0.0,                   # dx/dt
         c.rho,
         (cp.v - cm.v) / (2 * Defs.dx()),
         (cp.rho - cm.rho) / (2 * Defs.dx()))

end


text:
module MicroHH

## Export types and functions.
export Model
export prepare_model!, step_model!, save_model, load_model!

## Load packages.
using LoopVectorization
using Printf
using HDF5

## Include the necessary files.
include("StencilBuilder.jl")
include("Grid.jl")
include("Fields.jl")
include("Boundary.jl")
include("Timeloop.jl")
include("Dynamics.jl")
include("Pressure.jl")
include("Diagnostics.jl")
include("MultiDomain.jl")

## Global model data.
struct Model{TF <: Union{Float32, Float64}}
    name::String
    n_domains::Int
    last_measured_time::Ref{UInt64}

    grid::Vector{Grid}
    fields::Vector{Fields}
    boundary::Vector{Boundary}
    timeloop::Vector{Timeloop}
    pressure::Vector{Pressure}
    multidomain::Vector{MultiDomain}
end

function Model(name, n_domains, settings, TF)
    m = Model{TF}(name, n_domains, 0, [], [], [], [], [], [])
    for i in 1:n_domains
        push!(m.grid, Grid(settings[i]["grid"], TF))
        push!(m.fields, Fields(m.grid[i], settings[i]["fields"], TF))
        push!(m.boundary, Boundary(settings[i]["boundary"]))
        push!(m.timeloop, Timeloop(settings[i]["timeloop"]))
        push!(m.pressure, Pressure(m.grid[i], TF))
        push!(m.multidomain, MultiDomain(m.grid[i], settings[i]["multidomain"], TF))
    end

    return m
end

function calc_rhs!(m::Model, i)
    set_boundary!(m.fields[i], m.grid[i], m.boundary[i])
    calc_dynamics_tend!(m.fields[i], m.grid[i])
    calc_nudge_tend!(m.fields[i], m.grid[i], m.multidomain[i])
    calc_pressure_tend!(m.fields[i], m.grid[i], m.timeloop[i], m.pressure[i])
end

function prepare_model!(m::Model)
    # Prevent slow single precision performance due to subnormals.
    set_zero_subnormals(true)

    m.last_measured_time[] = time_ns()

    for i in 1:m.n_domains
        calc_rhs!(m, i)
    end

    check_model(m)

    return is_in_progress(m.timeloop[1])
end

function save_domain(m::Model, i)
    f = m.fields[i]
    g = m.grid[i]
    t = m.timeloop[i]
    b = m.boundary[i]

    # Update the boundary.
    set_boundary!(f, g, b)

    filename = @sprintf("%s.%02i.%08i.h5", m.name, i, round(t.time))
    h5open(filename, "w") do fid
        # Save the grid.
        write(fid, "x" , g.x[g.is:g.ie])
        write(fid, "xh", g.xh[g.is:g.ie])
        write(fid, "y" , g.y[g.js:g.je])
        write(fid, "yh", g.yh[g.js:g.je])
        write(fid, "z" , g.z[g.ks:g.ke])
        write(fid, "zh", g.zh[g.ks:g.keh])

        # Make the grid variables dimensions.
        HDF5.h5ds_set_scale(fid["x" ], "x" )
        HDF5.h5ds_set_scale(fid["xh"], "xh")
        HDF5.h5ds_set_scale(fid["y" ], "y" )
        HDF5.h5ds_set_scale(fid["yh"], "yh")
        HDF5.h5ds_set_scale(fid["z" ], "z" )
        HDF5.h5ds_set_scale(fid["zh"], "zh")

        # Save the fields.
        write(fid, "u", f.u[g.is:g.ie, g.js:g.je, g.ks:g.ke ])
        write(fid, "v", f.v[g.is:g.ie, g.js:g.je, g.ks:g.ke ])
        write(fid, "w", f.w[g.is:g.ie, g.js:g.je, g.ks:g.keh])
        write(fid, "s", f.s[g.is:g.ie, g.js:g.je, g.ks:g.ke ])
        write(fid, "s_bot", f.s_bot[g.is:g.ie, g.js:g.je])
        write(fid, "s_top", f.s_top[g.is:g.ie, g.js:g.je])
        write(fid, "s_gradbot", f.s_gradbot[g.is:g.ie, g.js:g.je])
        write(fid, "s_gradtop", f.s_gradtop[g.is:g.ie, g.js:g.je])

        # Attach the dimensions. Note the c-indexing.
        HDF5.h5ds_attach_scale(fid["u"], fid["xh"], 2)
        HDF5.h5ds_attach_scale(fid["u"], fid["y"], 1)
        HDF5.h5ds_attach_scale(fid["u"], fid["z"], 0)

        HDF5.h5ds_attach_scale(fid["v"], fid["x"], 2)
        HDF5.h5ds_attach_scale(fid["v"], fid["yh"], 1)
        HDF5.h5ds_attach_scale(fid["v"], fid["z"], 0)

        HDF5.h5ds_attach_scale(fid["w"], fid["x"], 2)
        HDF5.h5ds_attach_scale(fid["w"], fid["y"], 1)
        HDF5.h5ds_attach_scale(fid["w"], fid["zh"], 0)

        HDF5.h5ds_attach_scale(fid["s"], fid["x"], 2)
        HDF5.h5ds_attach_scale(fid["s"], fid["y"], 1)
        HDF5.h5ds_attach_scale(fid["s"], fid["z"], 0)

        HDF5.h5ds_attach_scale(fid["s_bot"], fid["x"], 1)
        HDF5.h5ds_attach_scale(fid["s_bot"], fid["y"], 0)

        HDF5.h5ds_attach_scale(fid["s_top"], fid["x"], 1)
        HDF5.h5ds_attach_scale(fid["s_top"], fid["y"], 0)

        HDF5.h5ds_attach_scale(fid["s_gradbot"], fid["x"], 1)
        HDF5.h5ds_attach_scale(fid["s_gradbot"], fid["y"], 0)

        HDF5.h5ds_attach_scale(fid["s_gradtop"], fid["x"], 1)
        HDF5.h5ds_attach_scale(fid["s_gradtop"], fid["y"], 0)
    end
end

function save_model(m::Model)
    for i in 1:m.n_domains
        save_domain(m, i)
    end
end

function load_domain!(m::Model, i)
    f = m.fields[i]
    g = m.grid[i]
    t = m.timeloop[i]
    b = m.boundary[i]

    filename = @sprintf("%s.%02i.%08i.h5", m.name, i, round(t.time))
    h5open(filename, "r") do fid
        f.u[g.is:g.ie, g.js:g.je, g.ks:g.ke ] = read(fid, "u")
        f.v[g.is:g.ie, g.js:g.je, g.ks:g.ke ] = read(fid, "v")
        f.w[g.is:g.ie, g.js:g.je, g.ks:g.keh] = read(fid, "w")
        f.s[g.is:g.ie, g.js:g.je, g.ks:g.ke ] = read(fid, "s")
        f.s_bot[g.is:g.ie, g.js:g.je] = read(fid, "s_bot")
        f.s_top[g.is:g.ie, g.js:g.je] = read(fid, "s_top")
        f.s_gradbot[g.is:g.ie, g.js:g.je] = read(fid, "s_gradbot")
        f.s_gradtop[g.is:g.ie, g.js:g.je] = read(fid, "s_gradtop")
    end

    set_boundary!(f, g, b)
end

function load_model!(m::Model)
    for i in 1:m.n_domains
        load_domain!(m, i)
    end
end

function step_model!(m::Model)
    # Prevent slow single precision performance due to subnormals.
    set_zero_subnormals(true)

    itime_next = m.timeloop[1].itime + m.timeloop[1].idt

    for i in 1:m.n_domains
        # Compute the nudging fields, which remain constant throughout the step.
        if i > 1
            calc_nudge_fields!(m.multidomain[i], m.fields[i], m.fields[i-1], m.grid[i], m.grid[i-1])
        end

        while (m.timeloop[i].itime < itime_next)
            integrate_time!(m.fields[i], m.grid[i], m.timeloop[i])
            step_time!(m.timeloop[i])

            if is_save_time(m.timeloop[i])
                save_domain(m, i)
            end

            calc_rhs!(m, i)
        end
    end

    if is_check_time(m.timeloop[1])
        check_model(m)
    end

    return is_in_progress(m.timeloop[1])
end

function check_model(m::Model)
    # Calculate the time since the last check.
    old_time = m.last_measured_time[]
    m.last_measured_time[] = time_ns()

    # First, print the model time and the wall clock since last sample.
    status_string = @sprintf("(%11.2f) Time = %8.3f",
        m.timeloop[1].time,
        (m.last_measured_time[] - old_time) * 1e-9)
    println(status_string)

    # Second, compute the divergence and CFL for each domain.
    for i in 1:m.n_domains
        status_string = @sprintf("  (%02i) Div = %6.3E, CFL = %6.3f, subnormals = %i",
            i,
            calc_divergence(m.fields[i], m.grid[i]),
            calc_cfl(m.fields[i], m.grid[i], m.timeloop[i]),
            sum(issubnormal.(m.fields[i].u)) + sum(issubnormal.(m.fields[i].v)) + sum(issubnormal.(m.fields[i].w)))
        println(status_string)
    end
end

end


text:
## interface to be implemented by specific problem types

abstract type AbstractMuseProblem end

struct Transformedθ end
struct UnTransformedθ end
is_transformed(::Transformedθ) = true
is_transformed(::UnTransformedθ) = false

@doc doc"""
    transform_θ(prob::AbstractMuseProblem, θ)

Map `θ` to a space where its domain is $(-\infty,\infty)$. Defaults to
identity function.
"""
transform_θ(prob::AbstractMuseProblem, θ) = θ

@doc doc"""
    inv_transform_θ(prob::AbstractMuseProblem, θ)

Map `θ` from the space where its domain is $(-\infty,\infty)$ back to
the original space. Defaults to identity function.
"""
inv_transform_θ(prob::AbstractMuseProblem, θ) = θ

@doc doc"""
Return the gradient of the joint log likelihood with respect to hyper
parameters `θ`, evaluated at data `x` and latent space `z`. The
signature of the function should be:

    ∇θ_logLike(prob::AbstractMuseProblem, x, z, θ)

If the problem needs a transformation of `θ` to map its domain to
$(-\infty,\infty)$, then it should instead implement:

    ∇θ_logLike(prob::AbstractMuseProblem, x, z, θ, θ_space)

where `θ_space` will be either `Transformedθ()` or `UnTransformedθ()`.
In this case, the `θ` argument will be passed in the space given by
`θ_space` and the gradient should be w.r.t. to `θ` in that space.

`z` must have domain $(-\infty,\infty)$. If a transformation is
required to make this the case, that should be handled internal to
this function and `z` will always refer to the transformed `z`.
"""
function ∇θ_logLike end
∇θ_logLike(prob::AbstractMuseProblem, x, z, θ, θ_space) =
    ∇θ_logLike(prob::AbstractMuseProblem, x, z, θ)

@doc doc"""
Return a tuple `(logLike, ∇z_logLike)` which give the log likelihood
and its gradient with respect to the latent space `z`, evaluated at
hyper parameters `θ` and data `x` . The signature of the function
should be:

    logLike_and_∇z_logLike(prob::AbstractMuseProblem, x, z, θ)

`z` must have domain $(-\infty,\infty)$. If a transformation is
required to make this the case, that should be handled internal to
this function and `z` will always refer to the transformed `z`.

The `θ` argument to this function will always be in the un-transfored
`θ` space.

!!! note

    Alternatively, custom problems can implement `ẑ_at_θ` directly and
    forego this method. The default `ẑ_at_θ` runs LBFGS with Optim.jl
    using `logLike_and_∇z_logLike`.
"""
function logLike_and_∇z_logLike end

@doc doc"""
Return a tuple `(x,z)` with data `x` and latent space `z` which are a
sample from the joint likelihood, given `θ`. The signature of the
function should be:

    sample_x_z(prob::AbstractMuseProblem, rng::AbstractRNG, θ)

Random numbers generated internally should use `rng`.

The `θ` argument to this function will always be in the un-transfored
`θ` space.
"""
function sample_x_z end

@doc doc"""
Return the log-prior at `θ`. The signature of the
function should be:

    logPriorθ(prob::AbstractMuseProblem, θ)

If the problem needs a transformation of `θ` to map its domain to
$(-\infty,\infty)$, then it should instead implement:

    logPriorθ(prob::AbstractMuseProblem, θ, θ_space)

where `θ_space` will be either `Transformedθ()` or `UnTransformedθ()`.
In this case, the `θ` argument will be passed in the space given by
`θ_space`.

Defaults to zero log-prior.
"""
logPriorθ(prob::AbstractMuseProblem, θ, θ_space) = logPriorθ(prob, θ)
logPriorθ(prob::AbstractMuseProblem, θ) = 0

@doc doc"""
Pre-process a user-provided `θ` into the data-structure used
internally in the computation. E.g. this allow the user to pass a
`NamedTuple` to functions like `muse` or `get_J!` while internally
converting it to a `ComponentVector`. The signature of the function
should be:

    standardizeθ(prob::AbstractMuseProblem, θ)

"""
standardizeθ(prob::AbstractMuseProblem, θ) = θ

@doc doc"""
Return the best-fit latent space `z` given data `x` and parameters
`θ`. The signature of the function should be:

    ẑ_at_θ(prob::AbstractMuseProblem, x, z₀, θ; ∇z_logLike_atol)

The return value should be `(ẑ, info)` where `info` can be any extra
diagonstic info which will be saved in the MUSE result.

The `θ` argument to this function will always be in the un-transfored
`θ` space.

The `z₀` should be used as a starting guess for the solution.

`z` must have domain $(-\infty,\infty)$. If a transformation is
required to make this the case, that should be handled internal to
this function, and the return value should refer to the transformed
`z`.

The default implementation of this method uses
[`logLike_and_∇z_logLike`](@ref) and Optim.jl's LBFGS to iteratively
maximize the log likelihood. Custom problems are free to override this
default if desired, in which case [`logLike_and_∇z_logLike`](@ref)
does not need to be implemented.
"""
function ẑ_at_θ(prob::AbstractMuseProblem, x, z₀, θ; ∇z_logLike_atol)
    soln = optimize(Optim.only_fg(z -> .-logLike_and_∇z_logLike(prob, x, z, θ)), z₀, Optim.LBFGS(), Optim.Options(g_tol=∇z_logLike_atol))
    _check_optim_soln(soln)
    soln.minimizer, soln
end

function _check_optim_soln(soln)
    Optim.converged(soln) || @warn("MAP solution did not converge within tolerance, result could be erroneous. Try tweaking `θ₀` or `∇z_logLike_atol` arguments to `muse` or fixing model.")
    isfinite(soln.minimum) || @error("MAP solution failed with logjoint(MAP)=$(soln.minimum).")
end

@doc doc"""
    check_self_consistency(
        prob,
        θ;
        fdm = central_fdm(3, 1),
        atol = 1e-3,
        rng = Random.default_rng(),
        has_volume_factor = true
    )

Checks the self-consistency of a defined problem at a given `θ`, e.g.
check that `inv_transform_θ(prob, transform_θ(prob, θ)) ≈ θ`, etc...
This is mostly useful as a diagonostic when implementing a new
`AbstractMuseProblem`.

A random `x` and `z` are sampled from `rng`. Finite differences are
computed using `fdm` and `atol` set the tolerance for `≈`.
`has_volume_factor` determines if the transformation includes the
logdet jacobian in the likelihood.
"""
function check_self_consistency(
    prob,
    θ;
    fdm = central_fdm(3, 1),
    atol = 1e-2,
    rng = Random.default_rng(),
    has_volume_factor = true
)

    θ = standardizeθ(prob, θ)
    x, z = sample_x_z(prob, rng, θ)
    # volume factor which is added by transformations. dont assume the
    # transformation is AD-able (eg it isnt for Turing)
    J(θ) = has_volume_factor ? FiniteDifferences.jacobian(fdm, θ -> transform_θ(prob, θ), θ)[1] : 1
    V(θ) = has_volume_factor ? logdet(J(θ)) : 0
    ∇θ_V(θ) = has_volume_factor ? FiniteDifferences.grad(fdm, V, θ)[1] : 0
    @testset "Self-consistency" begin
        @test inv_transform_θ(prob, transform_θ(prob, θ)) ≈ θ  atol=atol
        @test logPriorθ(prob, θ, UnTransformedθ()) ≈ logPriorθ(prob, transform_θ(prob, θ), Transformedθ()) .+ V(θ)  atol=atol
        @test ∇θ_logLike(prob, x, z, θ, UnTransformedθ()) ≈ J(θ)' * ∇θ_logLike(prob, x, z, transform_θ(prob, θ), Transformedθ()) .+ ∇θ_V(θ)  atol=atol
    end
end


text:
using CuArrays
using Distributions: pdf
using Random
using Flux
using BSON

RLBase.get_prob(p::AbstractPolicy, obs, ::RLBase.AbstractActionStyle, a) =
    pdf(get_prob(p, obs), a)

Random.rand(s::MultiContinuousSpace{<:CuArray}) = rand(CuArrays.CURAND.generator(), s)

# avoid fallback silently
Flux.testmode!(p::AbstractPolicy, mode = true) =
    @error "someone forgets to implement this method!!!"

function save(f::String, p::AbstractPolicy)
    policy = cpu(p)
    BSON.@save f policy
end

function load(f::String, ::Type{<:AbstractPolicy})
    BSON.@load f policy
    gpu(policy)
end


text:
function NRLMSISE00(arg0::NRLMSISE00InputParameters, arg1::PVCoordinatesProvider, arg2::BodyShape)
    return NRLMSISE00((NRLMSISE00InputParameters, PVCoordinatesProvider, BodyShape), arg0, arg1, arg2)
end

function NRLMSISE00(arg0::NRLMSISE00InputParameters, arg1::PVCoordinatesProvider, arg2::BodyShape, arg3::TimeScale)
    return NRLMSISE00((NRLMSISE00InputParameters, PVCoordinatesProvider, BodyShape, TimeScale), arg0, arg1, arg2, arg3)
end

function get_density(obj::NRLMSISE00, arg0::AbsoluteDate, arg1::Vector3D, arg2::Frame)
    return jcall(obj, "getDensity", jdouble, (AbsoluteDate, Vector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::NRLMSISE00, arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::Frame)
    return jcall(obj, "getDensity", RealFieldElement, (FieldAbsoluteDate, FieldVector3D, Frame), arg0, arg1, arg2)
end

function get_frame(obj::NRLMSISE00)
    return jcall(obj, "getFrame", Frame, ())
end

function with_switch(obj::NRLMSISE00, arg0::jint, arg1::jint)
    return jcall(obj, "withSwitch", NRLMSISE00, (jint, jint), arg0, arg1)
end


text:
using Documenter, AutomotiveDrivingModels

makedocs(
	 modules = [AutomotiveDrivingModels],
	 format = Documenter.HTML(),
	 sitename="AutomotiveDrivingModels.jl"
	 )

deploydocs(
    repo = "github.com/sisl/AutomotiveDrivingModels.jl.git",
)


text:
using Distributed, DelimitedFiles, JuMP, Gurobi, Statistics, StatsBase
include("../basics/init.jl")
include("../input_data/read_data.jl")
include("../input_data/data_manipulation.jl")
include("../basics/basic_functions.jl")
include("../basics/validations.jl")
include("../heuristics/indicators.jl")
include("../heuristics/construction_algorithms.jl")
include("../heuristics/local_search.jl")
include("./JuMP_functions.jl")
include("./JuMP_models.jl")
include("../solution_methods/GRASP_algorithm.jl")
include("../output_data/results.jl")

for file in reverse(readdir("../../data/Benchmark/"))
#for file in reverse(readdir("./data/pending/"))
    println("-------------")
    println(file)
    # if file != ".ipynb_checkpoints" && file != "00_DataDescription.txt"
    if file == "60C_10Type_LessDense_2QC.txt"
        doc = "../../data/Benchmark/"*file
        C,P,CP,Q,J,tt,d,ci,pj,bj,prejj,cpij = read_data(doc)

        possible_tasks = Array{LTask, 1}()
        tasks_by_position = Dict{Int, Array{LTask, 1}}()
        prec = Dict{Int, Array}()
        task_times = zeros(Int, P, C)

        for p = 1:P
            l=Array{LTask,1}()
            for t=1:CP
                target_bay = bj[cpij[t,1]+1]
                if p==P
                    #create array with all tasks
                    push!(possible_tasks, LTask(cpij[t,1]+1, target_bay, cpij[t,2]+1, cpij[t,3]))
                    #create matrix with all task times
                    task_times[cpij[t,1]+1, cpij[t,2]+1] = cpij[t,3]
                end
                #create dict with the tasks sorted by position
                if cpij[t,1]+1==p
                    push!(l,LTask(cpij[t,1]+1, target_bay, cpij[t,2]+1, cpij[t,3]))
                end
            end
            tasks_by_position[p]=sort!(l, by=x->x.t)

            #create dict with the position precedences
            s=Array{Int,1}()
            for j = 1:P
                if prejj[p,j]==1
                    push!(s,j)
                end
            end
            prec[p]=s

        end

        H = horizon_plan_new(C,P,J,Q,tt,d,bj,tasks_by_position)
        global CTS = Constants(C,P,J,Q,H,tt,d)
        global makespan, sol_x, sol_w
        #get minimum LoadingTime
        start = time()
        tasks_by_w = OperationalShippingProblem(0, task_times, bj, CTS)
        LT_min=0
        for (key,value) in tasks_by_w
            LT_min += value.t
        end
        LB = round((2*LT_min + CTS.tt*(maximum(bj)-minimum(bj)))/CTS.Q)
        # print("LT_min ="), print(LT_min), print(",  LB ="), println(LB)
        print("LB ="), println(LB)
        print("UB (H) ="), println(CTS.H)
        # start = time()
        # makespan, sol_x, sol_w = FlexibleShipLoadingProblem(1, 0, prec, task_times, bj, LB, CTS)
        exec_time = time() - start
        # println(makespan)
        println(exec_time)
    end
end

max_stars_repo_name: marcmartinezruiz/thesis

max_issues_repo_name: marcmartinezruiz/thesis

max_forks_repo_name: marcmartinezruiz/thesis


text:
using LinearAlgebra
using Rimu
using Rimu.RMPI
using Rimu.RMPI: sort_and_count!
using Test

@testset "RNG independence" begin
    m = n = 6
    add = near_uniform(BoseFS{n,m})
    svec = DVec(add => 2)
    dv = MPIData(svec)
    @test ConsistentRNG.check_crng_independence(dv) ==
        mpi_size()*Threads.nthreads()*fieldcount(ConsistentRNG.CRNG)
end

@testset "DistributeStrategies" begin
    # `DistributeStrategy`s
    ham = HubbardReal1D(BoseFS((1,2,3)))
    for setup in [RMPI.mpi_no_exchange, RMPI.mpi_all_to_all, RMPI.mpi_point_to_point]
        dv = DVec(starting_address(ham)=>10; style=IsDynamicSemistochastic())
        v = MPIData(dv; setup)
        df, state = lomc!(ham,v)
        @test size(df) == (100, 12)
    end
    # need to do mpi_one_sided separately
    dv = DVec(starting_address(ham)=>10; style=IsDynamicSemistochastic())
    v = RMPI.mpi_one_sided(dv; capacity = 1000)
    df, state = lomc!(ham,v)
    @test size(df) == (100, 12)
end

@testset "sort_and_count!" begin
    for l in (1, 2, 30, 1000)
        for k in (2, 10, 100)
            @testset "k=$k, l=$l" begin
                ordfun(x) = hash(x, hash(1)) % k
                vals = rand(1:10, l)
                counts = zeros(Int, k)
                displs = zeros(Int, k)

                RMPI.sort_and_count!(counts, displs, vals, ordfun.(vals), (0, k-1))
                @test issorted(vals, by=ordfun)
                @test sum(counts) == l

                for i in 0:(k - 1)
                    c = counts[i + 1]
                    d = displs[i + 1]
                    r = (1:c) .+ d
                    ords = ordfun.(vals)
                    @test all(ords[r] .== i)
                end
            end
        end
    end
end

@testset "MPIData" begin
    dv1 = MPIData(DVec(4 => 1, 3 => 2, 2 => 3, 1 => 4))
    dv2 = MPIData(empty(localpart(dv1)))

    @testset "Iteration and reductions" begin
        @test sort(collect(localpart(values(dv1)))) == 1:4

        @test sum(first, pairs(dv1)) == 10
        @test sum(last, pairs(dv1)) == 10
        @test prod(keys(dv1)) == 24
        @test sum(values(dv2)) == 0
    end
    @testset "Errors" begin
        @test_throws ErrorException [p for p in pairs(dv1)]
        @test_throws ErrorException [k for k in keys(dv1)]
        @test_throws ErrorException [v for v in values(dv1)]
    end
    @testset "Norms" begin
        @test walkernumber(dv1) ≡ 10.0
        @test norm(dv1, 1) ≡ 10.0
        @test norm(dv1, 2) ≈ 5.477225575
        @test norm(dv1, Inf) ≡ 4.0

        @test walkernumber(dv2) ≡ 0.0
        @test norm(dv2, 1) ≡ 0.0
        @test norm(dv2, 2) ≡ 0.0
        @test norm(dv2, Inf) ≡ 0.0

        @test_throws ErrorException norm(dv1, 3)
        @test_throws ErrorException norm(dv2, 15)
    end
    @testset "dot" begin
        @test dot(dv1, dv2) == 0
        @test dot(dv1, dv1) == dot(localpart(dv1), dv1)
        rand_ham = MatrixHamiltonian(rand(ComplexF64, 4,4))
        ldv1 = localpart(dv1)
        @test norm(dot(dv1, rand_ham, dv1)) ≈ norm(dot(ldv1, rand_ham, ldv1))
    end
end


text:
const MIN_SITE_FRACTION = 1.0e-15;
const R = 8.3145;
@variables P T


text:
export compute_R, build_R, build_R_acausal, parameter_search, get_data, fill_non_existing_entries

function compute_R(data::AbstractVector, windowsize::Int)
    denominator = data[1:end - windowsize]
    numerator = data[1 + windowsize:end]
    compute_R(numerator, denominator)
end

function compute_R(numerator::AbstractVector, denominator::AbstractVector)
    @assert length(numerator) == length(denominator) "inconsistent data lenghts"
    numerator ./ denominator
end

function build_R(df::DataFrame, past::Int, future::Int, k_gen::Int, pop_days::Int=0)
    @assert past >= 0 && future >= 0 "values for past = $past and future = $future must be non-negative"
    @assert k_gen > 0 "invalid value for k_gen ($(k_gen))"

    dates, data = df.days, df.cases

    window = 1 + past + future
    N = rollmean(data, window)
    days_N = dates[1 + past:end - future]

    R = compute_R(N, k_gen)
    days_R = dates[1 + k_gen + past:end - future]

    DataFrame(days = days_N, cases = N), DataFrame(days = days_R, R = R)
end

function build_R_acausal(df::DataFrame, past::Day, future::Day, k_gen::Day)
    df_N_temp, df_R = build_R(df, Dates.value(past), Dates.value(future), Dates.value(k_gen))
    if typeof(df_R.days) == Vector{Int}
       df_R.days = map(x -> Date(2020,2,3) + Day(x), df_R.days)
    end

    df_N_acausal = compute_cases_acausal(df, df_R, past, future, k_gen)
    relevant_days = df_N_acausal.days - k_gen

    num = df_N_acausal.cases
    den = df_N_temp[ first(relevant_days) .<= df_N_temp.days .<= last(relevant_days), :cases]
    R = compute_R(num, den)

    df_N_acausal, DataFrame(days = df_N_acausal.days, R = R)
end

build_R_acausal(df::DataFrame, past::Int, future::Int, k_gen::Int) = build_R_acausal(df, Day(past), Day(future), Day(k_gen))

function compute_cases_acausal(df_cases::DataFrame, df_reproduction::DataFrame, past::Day, future::Day, k_gen::Day)
    @assert typeof(df_reproduction.days) == Vector{Date} "Column `days` of df_cases  needs to have `Date` entries"
    start_date = first(df_reproduction).days - past
    end_date = last(df_reproduction).days + future

    @assert start_date ∈ df_cases.days "missing case numbers from the past (required start date: $(start_date))"
    @assert end_date ∈ df_cases.days "missing case numbers from the future (required end date: $(end_date))"

    N = Vector{Float64}()
    for row in eachrow(df_reproduction)
        day, R = row.days, row.R

        w = get_weights(R, past, future, k_gen)
        cases = get_values(df_cases, day, past, future, column = :cases)

        push!(N, mean(cases .* w))
    end

    DataFrame(days = df_reproduction.days .+ past, cases = N)
end

function get_data(df::DataFrame; days_col::String, data_col::String, kind::String)
    @assert kind ∈ ( "R", "cases" ) "kind $kind not supported"
    index = df[!, days_col]
    data = df[!, data_col]
    inds = .!ismissing.(data)

    DataFrame( Dict("days" => Vector(index[inds]), "$(kind)" => Vector(data[inds]) ) )
end

function parameter_search(df::DataFrame, past::Array{Int64,1}, future::Array{Int64,1}, k_gen:: Array{Int64,1}, data_col::String)
    result_r = DataFrame(NEUH = Float64[], NEUHA= Float64[], Past = Int[], Future = Int[], k_gen=Int[])
    result_n = DataFrame(Method = String[], Case= String[], Past = Int[], Future = Int[], k_gen=Int[])
    df_cases = get_data(df, days_col = "k", data_col = data_col, kind = "cases")
    df_cases.days = map(x -> Date(2020,2,3) + Day(x), df_cases.days)
    r_true = df[!, "true R"]
    for k in k_gen
        for i in past
            if i > k
                break
            end
            for j in future
                if j > k
                    break
                end
                neu_h_N, neu_h_R = build_R(df_cases,i, j,  k)
                error_h = mean(abs.(neu_h_R.R[max(1, 11 - i - j):length(neu_h_R.R)]- r_true[(max(1,11 - i -j)+ k + i):(length(r_true)- j)]))
                neu_ha_N, neu_ha_R = build_R_acausal(df_cases, i, j, k)
                error_ha = mean(abs.(neu_ha_R.R[max(1, 11 - i - j):length(neu_h_R.R)] -r_true[(max(1,11 - i -j)+ k + i):(length(r_true)- j)]))
                push!(result_r, [error_h, error_ha, i, j, k])
            end
        end
    end
    result_r
end

"""
function fill_non_existing_entries(x_ref::Vector, x::Vector, y::Vector; fill_with::Real=0)

    This function takes a pair `(x, y)` and compares it against the reference x_ref.
    In case there are missing values in `x` relative to `x_ref`, then `fill_with` is added at the corresponding index in `y` and returned (the value of `y` remains unchanged).

    ## Example

    ```julia
    x_ref = 1:10
    x = 1:2:10
    y = rand(length(x))

    y_ref = fill_non_existing_entries(x_ref, x, y)
    ```

    This yields.

    ```
    julia> [x y]
        5×2 Array{Float64,2}:
        1.0  0.245373
        3.0  0.682147
        5.0  0.919954
        7.0  0.153142
        9.0  0.235671

        julia> [x_ref y_ref]
        10×2 Array{Float64,2}:
        1.0  0.245373
        2.0  0.0
        3.0  0.682147
        4.0  0.0
        5.0  0.919954
        6.0  0.0
        7.0  0.153142
        8.0  0.0
        9.0  0.235671
        10.0  0.0
    ```
"""
function fill_non_existing_entries(x_ref, x, y::Vector; fill_with::Real=0)
    # sanity checks
    @assert length(x_ref) >= length(x)
    @assert length(x) == length(y)
    @assert eltype(x_ref) == eltype(x)

    # find all indicies of values of x in the reference x_ref
    inds_correct = findall(in(x), x_ref)
    # find all the other indicies, i.e. values that are not in x
    inds_to_correct = setdiff(1:length(x_ref), inds_correct)
    # generate target z
    z = zeros(eltype(y), length(x_ref))
    # ... and set correct values
    z[inds_correct] .= y
    # ... and set missing values
    z[inds_to_correct] .= fill_with
    # return
    z
end


text:
using RecipesBase, Distributions

@recipe function plot(X::SamplePath)
    X.tt, X.yy
end

@recipe function plot(X::VSamplePath)
    X.tt, X.yy'
end

@recipe function plot(X::Bridge.EstSamplePath{Float64}; quantile = 0.95)
    q = Distributions.quantile(Normal(), 1 - (1 - quantile)/2)
    ribbon --> q*sqrt.(X.vv)
    X.tt, X.yy
end

#@recipe function plot(X::Bridge.SamplePathBand)
#    ribbon --> (X.uu .- X.ll)
#    X.tt, (X.uu .+ X.ll)./2
#end


text:
"""

Computing GMD contributions for the parts of the solute

"""

function gmdget(topology :: String, gmd :: Array{}, sel1 :: String, sel2 :: String; vmd="vmd" )

  # gmd is the array resulting from reading the gmd contributions per atom
  # for the solute, produced by the gmd.f90 program of mdanalysis

  # The first selection corresponds the the total solute, the second
  # selection to the section of the solute that will be considered here

  sel1 = M3GTools.select(topology,sel1,vmd=vmd)
  sel2 = M3GTools.select(topology,sel2,vmd=vmd)

  # To do so, define a vector that contains the index in "solute" of the atoms
  # of the "backbone" selection. We sum 2 because this will correspond directly
  # to the columns of the gmd_solute array.

  iatoms = similar(sel2);
  for i in 1:length(sel2)
    iatoms[i] = findfirst(isequal(sel2[i]),sel1) + 2
  end

  # Sum all contributions of the gmd_solute array for iatoms:
  gmd_get = zeros(length(gmd[:,1]));
  for i in 1:length(gmd[:,1])
    for j in 1:length(sel2)
      gmd_get[i] = gmd_get[i] + gmd[i,iatoms[j]]
    end
  end

  return gmd_get

end

gmdget(topology :: String, gmd :: Array{}; data = "none", get = "none", vmd="vmd" ) = gmdget(topology, gmd, data, get, vmd=vmd)

gmdget(simulation :: Simulation, gmd :: Array{}, data, get) = gmdget(simulation, gmd, data, get, vmd=simulation.vmd)
gmdget(simulation :: Simulation, gmd :: Array{}; data = "none", get = "none" ) = gmdget(simulation.psf, gmd, data, get, vmd=simulation.vmd)


text:
function hold(cur_player::AbstractPlayer)

  # ------------
  #  start hold
  # ------------

  can_move(cur_player) || return

  cur_player.state.can_hold || return
  cur_player.state.can_hold = false

  score!(cur_player, "hold", -1)

  # ------------
  #  clean hold
  # ------------

  cur_js = """
    \$(".js-hold-piece td").removeClass();
    var tmp_cell;
  """

  if is_repl
    cur_coords = map(
      cur_block -> calc_block_coords(cur_block),
      cur_player.piece.blocks
    )

    cur_string = []

    row_info = map(first, cur_coords)
    col_info = map(last, cur_coords)

    max_row = min(cur_player.grid.rows, maximum(row_info) + cur_player.piece.width )

    min_col = max(1, minimum(col_info) - cur_player.piece.width )
    max_col = min(cur_player.grid.cols, maximum(col_info) + cur_player.piece.width )

    for cur_row in 1:max_row
      for cur_col in min_col:max_col
        has_piece = cur_player.grid.table[cur_row, cur_col] != ""
        has_piece && continue

        push!(cur_string, "\x1b[$(cur_player.grid.rows-cur_row+4);$(29+2*cur_col)H")
        push!(cur_string, crayon_dict["invisible"])
        push!(cur_string, "  ")
        push!(cur_string, inv(crayon_dict["invisible"]))
      end
    end
  end

  # -------------------
  #  show hold preview
  # -------------------

  work_piece, cur_player.piece =
    cur_player.piece, cur_player.hold

  work_piece.owner = Nullable{AbstractContainer}()

  make_piece!(
    cur_player,
    work_piece.name,
    :hold
  )

  cur_game = cur_player.game
  cur_piece = cur_player.hold

  if is_repl
    for cur_row in 1:2
      for cur_col in 1:4
        push!(cur_string, "\x1b[$(21-cur_row);$(9+2*cur_col)H")
        push!(cur_string, crayon_dict["invisible"])
        push!(cur_string, "  ")
        push!(cur_string, inv(crayon_dict["invisible"]))
      end
    end
  end

  for cur_block in cur_piece.blocks
    (cur_row, cur_col) = calc_block_coords(cur_block)

    cur_row += 2
    cur_col += 3

    ( cur_piece.name == 'i' ) && ( cur_row -= 1 )

    cur_js *= """
      tmp_cell = \$(".js-hold-piece .cs-row-$(cur_row) td:nth-child($(cur_col))");
      tmp_cell.addClass("cs-color cs-$(cur_piece.color)");
    """

    if is_repl
      push!(cur_string, "\x1b[$(21-cur_row);$(9+2*cur_col)H")
      push!(cur_string, crayon_dict["$(cur_piece.color)_shadow"])
      push!(cur_string, "  ")
      push!(cur_string, inv(crayon_dict["$(cur_piece.color)_shadow"]))
    end
  end

  # --------------
  #  wrap-up hold
  # --------------

  if is_repl
    push!(cur_string, "\x1b[u")
    print(cur_string...)
  end

  tetris_js(
    cur_game.scope,
    JSString(cur_js)
  )

  isnull(cur_player.piece) && return
  pin_piece!(cur_player.piece)

  reset_shadow!(cur_player)
  move!(cur_player)

  step(cur_player)

  return

end


text:
# Unfortunately the default wandb agent API causes a segfault.

# function sweep(sweep_config::Dict; entity = nothing, project = nothing)
#     return wandb.sweep(sweep_config, entity, project)
# end

# function agent(
#     sweep_id;
#     func = nothing,
#     entity = nothing,
#     project = nothing,
#     count = nothing,
# )
#     return wandb.agent(sweep_id, func, entity, project, count)
# end

mutable struct WandbHyperParameterSweep
    sweep_tag::String
end

WandbHyperParameterSweep() = WandbHyperParameterSweep(randstring(12))

function (hpsweep::WandbHyperParameterSweep)(func, config,
                                             # For Compat with FluxTraining and other Integrations
                                             logger=WandbLogger, args...; func_args=(), func_kwargs=(;), kwargs...)
    lg = logger(args...; tags=[hpsweep.sweep_tag], kwargs...)

    try
        @info "Logging to Wandb" logger = lg

        # The user can pass other global configuration options to the sweep
        update_config!(lg, config)

        res = func(lg, Dict(get_config(lg).__dict__["_items"]), func_args...; func_kwargs...)

        close(lg)

        return res
    catch
        # Close the logger else we have to restart the julia session to get things
        # back into sync
        close(lg)
        rethrow()
    end
end


text:
using testDoc
using Test

@testset "testDoc.jl" begin
    # Write your tests here.
end


text:
import MPI
include("04-sendrecv-impl.jl")

function main()
    MPI.Init()

    do_sendrecv()

    MPI.Finalize()
end

main()


text:
module MARY_fRG

include("basics/basics.jl")

include("helpers/drawers.jl")

include("helpers/triangulated.jl")

include("helpers/refine.jl")

include("helpers/interactions.jl")

include("fermi/fermi.jl")

include("flowequ/flow_equation.jl")

end # module


text:
# TODO: make all insertions to be carried via one method (?)
#       Is that possible? Yes

# TODO: make immutable
#=
    stores hash of a monomial,
    corresponding divmask to speed up divisibility checks,
    index for position matrix (defaults to zero),
    and the todal degree
=#
mutable struct Hashvalue
    hash::UInt32

    #=

    =#
    divmask::UInt32

    idx::Int
    deg::Int
end

#=
    Hashtable implementing linear probing
    and designed to store and operate with monomials
=#
mutable struct MonomialHashtable
    # stores monomial exponent vectors,
    # assumes degrees are < 2^16
    exponents::Vector{Vector{UInt16}}

    # maps exponent hash to its position in exponents array
    hashtable::Vector{Int}

    # stores hashes, division masks,
    # and other valuable info
    # for each hashtable enrty
    hashdata::Vector{Hashvalue}

    # values to hash exponents with, i.e
    # hash(e) = sum(hasher .* e)
    hasher::Vector{UInt32}

    #= Ring information =#
    # number of variables
    nvars::Int
    # raw length of exponent vector
    explen::Int
    # ring monomial ordering
    ord::Symbol

    #= Divisibility =#
    # divisor map to check divisibility faster
    divmap::Vector{UInt32}
    # variables selected for divmap
    divvars::Vector{Int}
    # count of divmap variables
    ndivvars::Int
    # bits per div variable
    ndivbits::Int

    size::Int
    # elements added
    load::Int
    #
    offset::Int
end

#------------------------------------------------------------------------------

# initialize and set fields for basis hashtable
# TODO: Initial size is a nice parameter to play with
function initialize_basis_hash_table(
        ring::PolyRing,
        rng::Random.AbstractRNG;
        initial_size::Int=2^16)

    # TODO: adapive choice of `initial_size`
    exponents = Vector{Vector{UInt16}}(undef, initial_size)
    hashdata  = Vector{Hashvalue}(undef, initial_size)
    hashtable = zeros(Int, initial_size)

    nvars = ring.nvars
    explen = ring.explen
    ord = ring.ord

    # initialize hashing vector
    hasher = zeros(UInt32, explen)
    for i in 1:explen
        # we don't want hash vector components to be zero
        while hasher[i] == 0
            # TODO: pass random generator
            hasher[i] = rand(rng, UInt32)
        end
    end

    # exponents[1:load] cover all stored exponents
    # , also exponents[1] is zeroed by default
    load = 1
    size = initial_size

    # exponents array starts from index offset,
    # We store buffer array at index 1
    offset = 2

    # initialize fast divisibility params
    charbit    = 8 # TODO ??
    int32bits  = charbit * sizeof(Int32)
    int32bits != 32 && error("Strange story with ints")
    ndivbits   = div(int32bits, nvars)
    # division mask stores at least 1 bit
    # per each of first charbit*sizeof(Int32) variables
    ndivbits == 0 && (ndivbits += 1)
    ndivvars = nvars < int32bits ? nvars : int32bits
    divvars  = Vector{Int}(undef, ndivvars)
    divmap   = Vector{UInt32}(undef, ndivvars * ndivbits)
    # count only first ndivvars variables for divisibility checks
    for i in 1:ndivvars
        divvars[i] = i
    end

    # first stored exponent used as buffer lately
    exponents[1] = zeros(UInt16, explen)

    MonomialHashtable(
        exponents, hashtable, hashdata, hasher,
        nvars, explen, ord,
        divmap, divvars, ndivvars, ndivbits,
        size, load, offset)
end

# initialize hashtable either for `symbolic_preprocessing` or for `update` functions
# These are of the same purpose and structure as basis hashtable,
# but are more local oriented
function initialize_secondary_hash_table(basis_ht::MonomialHashtable)

    # hmm TODO
    initial_size = 2^6

    exponents = Vector{Vector{UInt16}}(undef, initial_size)
    hashdata  = Vector{Hashvalue}(undef, initial_size)
    hashtable = zeros(Int, initial_size)

    # preserve ring info
    explen = basis_ht.explen
    nvars  = basis_ht.nvars
    ord    = basis_ht.ord

    # preserve division info
    divmap   = basis_ht.divmap
    divvars  = basis_ht.divvars
    ndivvars = basis_ht.ndivvars
    ndivbits = basis_ht.ndivbits

    # TODO: preserve hasher?
    # yes
    hasher = basis_ht.hasher

    load = 1
    size = initial_size
    offset = 2

    exponents[1] = zeros(UInt16, explen)

    MonomialHashtable(
        exponents, hashtable, hashdata, hasher,
        nvars, explen, ord,
        divmap, divvars, ndivvars, ndivbits,
        size, load, offset)
end

#------------------------------------------------------------------------------

# resizes (if needed) ht so that it can store `size` elements,
# and clears all previoud data
function reinitialize_hash_table!(ht::MonomialHashtable, size::Int)
    if size > ht.size
        while size > ht.size
            ht.size *= 2
        end
        # TODO
        resize!(ht.hashdata, ht.size)
        resize!(ht.exponents, ht.size)
    end
    ht.hashtable = zeros(Int, ht.size)
    ht.load = 1
end

# doubles the size of storage in `ht`,
# and rehashes all elements
function enlarge_hash_table!(ht::MonomialHashtable)
    ht.size *= 2
    resize!(ht.hashdata,  ht.size)
    resize!(ht.exponents, ht.size)

    ht.hashtable = zeros(Int, ht.size)

    mod = UInt32(ht.size - 1)
    for i in ht.offset:ht.load  # TODO: hardcoding 2 is bad
        # hash for this elem is already computed
        he = ht.hashdata[i].hash
        hidx = he
        for j in UInt32(1):UInt32(ht.size)
            hidx = (he + j) & mod + UInt32(1)
            ht.hashtable[hidx] != 0 && continue
            ht.hashtable[hidx] = i
            break
        end
    end
end

#------------------------------------------------------------------------------

function insert_in_hash_table!(ht::MonomialHashtable, e::Vector{UInt16})
    # generate hash
    he = UInt32(0)
    for i in 1:ht.explen
        he += ht.hasher[i] * e[i]
    end
    # find new elem position in the table
    hidx = Int(he)  # Int for type stability
    # power of twoooo
    mod = UInt32(ht.size - 1)
    i = UInt32(1)

    @label Restart
    while i < ht.size
        # TODO: & instead of %
        hidx = (he + i) & mod + UInt32(1)
        @inbounds vidx  = ht.hashtable[hidx]

        # if free
        vidx == 0 && break

        # if not free and not same hash
        if ht.hashdata[vidx].hash != he
            i += UInt32(1)
            continue
        end

        present = ht.exponents[vidx]
        for j in 1:ht.explen
            # if hash collision
            if present[j] != e[j]
                i += UInt32(1)
                @goto Restart
            end
        end

        # already present in hashtable
        return vidx
    end

    # add its position to hashtable, and insert exponent to that position
    vidx = ht.load + 1
    ht.hashtable[hidx] = vidx
    ht.exponents[vidx] = copy(e)

    divmask = generate_monomial_divmask(e, ht)
    ht.hashdata[vidx] = Hashvalue(he, divmask, 0, e[end])

    ht.load += 1

    return vidx
end

#------------------------------------------------------------------------------

#=
    Having `ht` filled with monomials from input polys,
    computes ht.divmap and divmask for each of already stored monomials
=#
function fill_divmask!(ht::MonomialHashtable)
    ndivvars = ht.ndivvars
    divvars  = ht.divvars

    min_exp  = Vector{UInt16}(undef, ndivvars)
    max_exp  = Vector{UInt16}(undef, ndivvars)

    e = ht.exponents[ht.offset]
    for i in 1:ndivvars
        min_exp[i] = e[divvars[i]]
        max_exp[i] = e[divvars[i]]
    end

    for i in ht.offset:ht.load # TODO: offset
        e = ht.exponents[i]
        for j in 1:ndivvars
            if e[divvars[j]] > max_exp[j]
                max_exp[j] = e[divvars[j]]
                continue
            end
            if e[divvars[j]] < min_exp[j]
                min_exp[j] = e[divvars[j]]
            end
        end
    end

    ctr   = 1
    steps = UInt32(0)
    for i in 1:ndivvars
        steps = div(max_exp[i] - min_exp[i], ht.ndivbits)
        (steps == 0) && (steps += 1)
        for j in 1:ht.ndivbits
            ht.divmap[ctr] = steps
            steps += 1
            ctr += 1
        end
    end

    for vidx in ht.offset:ht.load
        unmasked = ht.hashdata[vidx]
        e = ht.exponents[vidx]
        divmask  = generate_monomial_divmask(e, ht)
        ht.hashdata[vidx] = Hashvalue(unmasked.hash, divmask, 0, e[end])
    end
end

#=
    TODO

=#
function generate_monomial_divmask(
            e::Vector{UInt16},
            ht::MonomialHashtable)

    divvars = ht.divvars
    divmap  = ht.divmap

    ctr = UInt32(1)
    res = UInt32(0)
    for i in 1:ht.ndivvars
        for j in 1:ht.ndivbits
            @inbounds if e[divvars[i]] >= divmap[ctr]
                res |= UInt32(1) << (ctr - 1)
            end
            ctr += UInt32(1)  # for typ stability
        end
    end

    res
end

#------------------------------------------------------------------------------

# h1 divisible by h2
function is_monom_divisible(h1::Int, h2::Int, ht::MonomialHashtable)

    if (ht.hashdata[h2].divmask & ~ht.hashdata[h1].divmask) != 0
        return false
    end

    e1 = ht.exponents[h1]
    e2 = ht.exponents[h2]
    # TODO: one less iteration is possible
    for i in 1:ht.explen
        if e1[i] < e2[i]
            return false
        end
    end

    return true
end

function is_gcd_const(h1::Int, h2::Int, ht::MonomialHashtable)
    e1 = ht.exponents[h1]
    e2 = ht.exponents[h2]

    for i in 1:ht.explen - 1
        if e1[i] != 0 && e2[i] != 0
            return false
        end
    end

    return true
end

#------------------------------------------------------------------------------

# compare pairwise divisibility of lcms from a[first:last] with lcm
function check_monomial_division_in_update(
            a::Vector{Int}, first::Int, last::Int,
            lcm::Int, ht::MonomialHashtable)

    # pairs are sorted, we only need to check entries above starting point

    divmask = ht.hashdata[lcm].divmask
    lcmexp  = ht.exponents[lcm]

    j = first
    @label Restart
    while j <= last
        # bad lcm
        if a[j] == 0
            j += 1
            continue
        end
        # fast division check
        if (~ht.hashdata[a[j]].divmask & divmask) != 0
            j += 1
            continue
        end
        ea = ht.exponents[a[j]]
        for i in 1:ht.explen
            if ea[i] < lcmexp[i]
                j += 1
                @goto Restart
            end
        end
        # mark as redundant
        a[j] = 0
    end

end

#------------------------------------------------------------------------------

# add monomials from `poly` multiplied by exponent vector `etmp`
# with hash `htmp` to hashtable `symbol_ht`,
# and substitute hashes in row
function insert_multiplied_poly_in_hash_table!(
        row::Vector{Int},
        htmp::UInt32,
        etmp::Vector{UInt16},
        poly::Vector{Int},
        ht::MonomialHashtable,
        symbol_ht::MonomialHashtable)

    # oof

    # length of poly to add
    len    = length(poly)
    explen = ht.explen

    mod = UInt32(symbol_ht.size - 1)

    bexps = ht.exponents
    bdata = ht.hashdata

    sexps = symbol_ht.exponents
    sdata = symbol_ht.hashdata

    l = 1 # hardcoding 1 does not seem nice =(
    @label Letsgo
    while l <= len
        # we iterate over all monoms of the given poly,
        # multiplying them by htmp/etmp,
        # and inserting into symbolic hashtable

        # hash is linear, so that
        # hash(e1 + e2) = hash(e1) + hash(e2)
        # We also assume that the hashing vector is shared same
        # between all created hashtables
        h = htmp + bdata[poly[l]].hash
        e = bexps[poly[l]]
        # println("monom of index $(poly[l]) : $e")

        lastidx = symbol_ht.load + 1
        #=
        if !isassigned(sexps, lastidx)
            sexps[lastidx] = Vector{UInt16}(undef, explen)
        end
        enew = sexps[lastidx]
        =#
        @inbounds enew = sexps[1]
        @inbounds for j in 1:explen
            # multiplied monom exponent
            enew[j] = etmp[j] + e[j]
        end

        # now insert into hashtable
        k = h
        i = UInt32(1)
        @label Restart
        while i <= symbol_ht.size  # TODO: < or <= ?
            # @assert ispow2(mod + 1)
            k = (h + i) & mod + UInt32(1)
            @inbounds vidx = symbol_ht.hashtable[k]
            # if index is free
            vidx == 0 && break
            # if different exponent is stored here
            if sdata[vidx].hash != h
                i += UInt32(1)
                continue
            end

            estored = sexps[vidx]
            @inbounds for j in 1:explen
                # hash collision, restarting search
                if estored[j] != enew[j]
                    i += UInt32(1)
                    @goto Restart
                end
            end

            row[l] = vidx
            l += 1

            @goto Letsgo
        end

        # add multiplied exponent to hash table
        if !isassigned(sexps, lastidx)
            sexps[lastidx] = Vector{UInt16}(undef, explen)
        end
        sexpsnew = sexps[lastidx]
        for j in 1:explen
            # multiplied monom exponent
            @inbounds sexpsnew[j] = enew[j]
        end
        symbol_ht.hashtable[k] = lastidx

        divmask = generate_monomial_divmask(enew, symbol_ht)
        sdata[lastidx] = Hashvalue(h, divmask, 0, enew[end])

        row[l] = lastidx
        l += 1
        symbol_ht.load += 1
    end

    row
end

function multiplied_poly_to_matrix_row!(
        symbolic_ht::MonomialHashtable, basis_ht::MonomialHashtable,
        htmp::UInt32, etmp::Vector{UInt16}, poly::Vector{Int})

    row = copy(poly)
    while symbolic_ht.load + length(poly) >= symbolic_ht.size
        enlarge_hash_table!(symbolic_ht)
    end
    insert_multiplied_poly_in_hash_table!(row, htmp, etmp, poly, basis_ht, symbolic_ht)
end

#------------------------------------------------------------------------------

function insert_in_basis_hash_table_pivots(
        row::Vector{Int},
        ht::MonomialHashtable,
        symbol_ht::MonomialHashtable,
        col2hash::Vector{Int})

    while ht.size - ht.load <= length(row)
        enlarge_hash_table!(ht)
    end

    sdata = symbol_ht.hashdata
    sexps = symbol_ht.exponents

    mod    = UInt32(ht.size - 1)
    explen = ht.explen
    bdata  = ht.hashdata
    bexps  = ht.exponents
    bhash  = ht.hashtable

    l = 1
    @label Letsgo
    while l <= length(row)
        hidx = col2hash[row[l]]

        # symbolic hash
        h = sdata[hidx].hash

        lastidx = ht.load + 1
        # TODO: speed this up
        bexps[lastidx] = copy(sexps[hidx])
        e = bexps[lastidx]

        k = h
        i = UInt32(1)
        @label Restart
        while i <= ht.size
            k = (h + i) & mod + UInt32(1)
            hm = bhash[k]

            hm == 0 && break
            if bdata[hm].hash != h
                i += UInt32(1)
                continue
            end

            ehm = bexps[hm]
            for j in 1:explen
                if e[j] != ehm[j]
                    i += UInt32(1)
                    @goto Restart
                end
            end

            row[l] = hm
            l += 1
            @goto Letsgo
        end

        bhash[k] = pos = lastidx
        row[l] = pos
        l += 1

        bdata[pos] = Hashvalue(h, sdata[hidx].divmask,
                             sdata[hidx].idx, sdata[hidx].deg)

        ht.load += 1
    end
end

function insert_plcms_in_basis_hash_table!(
        pairset::Pairset,
        off::Int,
        ht::MonomialHashtable,
        update_ht::MonomialHashtable,
        basis::Basis,
        plcm::Vector{Int},
        ifirst::Int, ilast::Int)

    # including ifirst and not including ilast

    gens  = basis.gens
    mod   = UInt32(ht.size - 1)
    ps    = pairset.pairs

    m = ifirst
    l = 1
    @label Letsgo
    while l < ilast
        if plcm[l] == 0
            l += 1
            continue
        end

        if is_gcd_const(gens[ps[off + l].poly1][1], gens[ps[off + 1].poly2][1], ht)
            l += 1
            continue
        end

        ps[m] = ps[off + l]

        # TODO: IT IS NOT CORRECT
        # upd: it is, but it can be done better
        h = update_ht.hashdata[plcm[l]].hash
        ht.exponents[ht.load+1] = copy(update_ht.exponents[plcm[l]])
        n = ht.exponents[ht.load+1]

        k = h
        i = UInt32(1)
        @label Restart
        while i <= ht.size
            k = (h + i) & mod + UInt32(1)
            hm = ht.hashtable[k]

            hm == 0 && break
            if ht.hashdata[hm].hash != h
                i += UInt32(1)
                continue
            end

            ehm = ht.exponents[hm]

            # @info "SO, we have " n ehm m
            for j in 1:ht.explen
                if ehm[j] != n[j]
                    i += 1
                    @goto Restart
                end
            end

            ps[m] = SPair(ps[m].poly1, ps[m].poly2, hm, ps[m].deg)
            m += 1
            l += 1
            @goto Letsgo
        end

        ht.hashtable[k] = pos = ht.load + 1

        uhd = update_ht.hashdata

        ll = plcm[l]
        ht.hashdata[ht.load + 1] = Hashvalue(h, uhd[ll].divmask, 0, uhd[ll].deg)

        ht.load += 1
        ps[m] = SPair(ps[m].poly1, ps[m].poly2, pos, ps[m].deg)
        m += 1
        l += 1
    end

    pairset.load = m - 1
end

#------------------------------------------------------------------------------

# computes lcm of he1 and he2 as exponent vectors from ht1
# and inserts it in ht2
function get_lcm(he1::Int, he2::Int,
                    ht1::MonomialHashtable, ht2::MonomialHashtable)

    e1   = ht1.exponents[he1]
    e2   = ht1.exponents[he2]
    etmp = ht1.exponents[1]

    # TODO: degrevlex only
    etmp[end] = 0
    for i in 1:ht1.explen-1
        etmp[i]    = max(e1[i], e2[i])
        etmp[end] += etmp[i]
    end

    insert_in_hash_table!(ht2, etmp)
end

#------------------------------------------------------------------------------


text:
"Maps some argument types from hsa_*_t to special wrapper types in the julia API"
function map_argtypes(obuf)
    for exu in obuf
        if isa(exu, Expr) && exu.head == :function
            # Go through all generated functions and extract
            # the arguments from the signature
            sig = exu.args[1]
            args = sig.args[2:end]

            for a in args
                a_type = a.args[2]

                if haskey(argtype_map, a_type)
                    a.args[2] = argtype_map[a_type]
                end
            end
        end
    end
end


text:
using Plots
using LinearAlgebra

σ11 = σ22 = 1.0
σ12 = σ21 = -0.7

Σ = [
    σ11 σ12
    σ21 σ22
]

L = cholesky(Σ).L
cL = size(L, 2)

boxmuller(r1, r2) = √(-2 * log(r1)) * sin(2π * r2)
univariate(n...) = boxmuller.(rand(n...), rand(n...))

function multivariate(n)
    x = univariate(cL, n)
    return L * x
end

xy = multivariate(100)
scatter(xy[1, :], xy[2, :])

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module Calculators

export AbstractCalculator,
       AbstractCalculationProgram,
       bsse_corrected_energy,
       calculate_energy,
       Calculator,
       clean_calculation_files,
       getBSSEsteps,
       Orca,
       read_energy,
       write_input

using ..Clusters
using Distributed

"""
print_orca_xyz(io::IO, c::AbstractClusterWithSymbols; ghost=undef)

Prints ORCA input style xyz input. Idexes in `ghost` are marked as ghost atoms.
"""
function print_orca_xyz(io::IO, c::AbstractClusterWithSymbols; ghost=undef)
    #TODO remove expicit referenzes to Cluster.xyz
    if ghost == undef
        for i in 1:length(c)
            println(io, c.atoms[i].id, "   ", c.xyz[1,i], "  ", c.xyz[2,i], "  ", c.xyz[3,i])
        end
    else
        for i in 1:length(c)
            if i in ghost
                println(io, c.atoms[i].id, ":   ", c.xyz[1,i], "  ", c.xyz[2,i], "  ", c.xyz[3,i])
            else
                println(io, c.atoms[i].id, "   ", c.xyz[1,i], "  ", c.xyz[2,i], "  ", c.xyz[3,i])
            end
        end
    end
end

abstract type AbstractCalculationProgram end

"""
    mutable struct Orca <: AbstractCalculationProgram

Holds information of how to use ORCA program.

Note that if you are using more than 1 core you need to specify exact
path to orca executable.

# Fields
- `executable` : path to "orca" executable - default = `"orca"`
- `ncore::Uint` : number of cores used by orca - default = `1`
- `maxmem:.Uint` : maximum memory per core for orca in mega bytes - default = `1000`
- `tmp_dir` : directory where calculations are done - default = `mktempdir()`

# Examples
```julia
julia> Orca()
Orca("orca", 0x0000000000000001, 0x00000000000003e8, "/tmp/jl_k4hnkY")

julia> Orca(executable="/opt/Orca/bin/orca", maxmem=4000, tmp_dir="/tmp/calculations/")
Orca("/opt/Orca/bin/orca", 0x0000000000000001, 0x0000000000000fa0, "/tmp/calculations/")
```
"""
mutable struct Orca <: AbstractCalculationProgram
    "path for orca excecutable"
    executable
    "number of cores in calculation"
    ncore::UInt
    "maximum memory per core"
    memcore::UInt
    "directory where calculation is done"
    tmp_dir
    function Orca(;executable="orca",
                   ncore=1, maxmem=1000, tmp_dir=mktempdir())
        new(executable, ncore, maxmem, tmp_dir)
    end
end

"""
    getBSSEsteps(cal::AbstractCalculationProgram)

Used to set up progressbar. Returns 1, if calculator program has a method to calculate
counter poise correction with a single command and 3, if program needs to be called
multiple times to make counter poise correction.
"""
getBSSEsteps(cal::AbstractCalculationProgram) = 1
getBSSEsteps(cal::Orca) = 3

abstract type AbtractCalculator{T<:AbstractCalculationProgram} end

"""
    mutable struct Calculator

Struct used to hold together different bits of information needed in calculations.

# Fields
- `method` : calculation method informaton
- `basis` : basis information
- `calculator` : calculation program information

# Examples
```julia
julia> Calculator("method", "basis", Orca())
Calculator{Orca}("method", "basis", Orca("orca", 0x0000000000000001, 0x00000000000003e8, "/tmp/tmpVg959k"))
```
"""
mutable struct Calculator{T} <: AbtractCalculator{T}
    method
    basis
    calculator::T
    function Calculator(method::AbstractString, basis::AbstractString, calculator::T) where T<:AbstractCalculationProgram
        new{T}(method, basis, calculator)
    end
end

"""
    write_input(io::IO, cal::Calculator,
                     c::AbstractClusterWithSymbols; ghost=undef)

Writes input files for calculation program. Only ORCA input is supported now.

# Arguments
- `io::IO` : stream where writing is done
- `cal::Calculator` : calculation information
- `c::AbstractClusterWithSymbols` : molecule cluster that is calculated
- `ghost=undef` : collection for atom idexes that are considered ghost atoms

"""
function write_input(io::IO, cal::Calculator{Orca},
                     c::AbstractClusterWithSymbols; ghost=undef)
    if typeof(cal.calculator) == Orca
        println(io, "! ", cal.method)
        println(io, "! ", cal.basis)
        println(io, "! MINIPRINT")
        if cal.calculator.ncore > 1
            println(io, "%pal nprocs $(cal.calculator.ncore) end")
        end
        println(io, "%maxcore $(cal.calculator.memcore)")
        println(io, "* xyz 0 1")
        print_orca_xyz(io, c, ghost=ghost)
        println(io,"*")
    else
        error("Calculator type not recogniced")
    end
end

"""
    read_energy(fname)

Read ORCA output for final energy. Throws an error if calculation failed.
"""
function read_energy(fname)
    lines = ""
    open(fname, "r") do f
        lines = readlines(f)
    end
    if ! occursin("****ORCA TERMINATED NORMALLY****", lines[end-1])
        error("Orca failed : $(lines[end-1])")
    end
    for l in reverse(lines)
        if occursin("FINAL SINGLE POINT ENERGY", l)
            return parse(Float64, split(l)[end])
        end
    end
end

"""
    calculate_energy(cal::Calculator, points; basename="base", ghost=undef, id="", pchannel=pchannel)

Calculates energy for given clusters.

# Arguments
- `cal::Calculator` : calculation information
- `points` : collection of [`Cluster`](@ref) which energy is calculated

# Keywords
- `basename="base"` : base name for input/ouput files
- `ghost=undef` : indices for atoms considered as ghost
- `id=""` : additional information for calculator - needed for multiprocessing in same folder
- `pchannel=undef`  : (Remote)Channel where progess information is added
"""
function calculate_energy(cal::Calculator{Orca}, points; basename="base", ghost=undef, id="", pchannel=undef)
    clean_calculation_files(basename=basename)
    inname = "$(basename).inp"
    outname= "$(basename).out"
    cmd = pipeline(`$(cal.calculator.executable) $(inname)`, outname)
    function do_calculation(point)
        ts = time()
        open(inname,"w") do io
            write_input(io, cal, point, ghost=ghost)
        end
        run(cmd)
        out = read_energy(outname)
        te = time()
        @debug "Calculation done in $(round(te-ts, digits=1)) seconds"
        pchannel != undef && put!(pchannel,true)
        return out
    end
    return map(x -> do_calculation(x), points)
end

"""
    calculate_energy(cal::Calculator, point::Cluster; basename="base", ghost=undef, id="", pchannel=pchannel)

Calculates energy for given cluster.

# Arguments
- `cal::Calculator` : calculation information
- `point::Cluster` : cluster which energy is calculated

# Keywords
- `basename="base"` : base name for input/ouput files
- `ghost=undef` : indices for atoms considered as ghost
- `id=""` : additional information for calculator - needed for multiprocessing in same folder
- `pchannel=undef`  : (Remote)Channel where progess information is added
"""
function calculate_energy(cal::Calculator{Orca}, point::Cluster; basename="base", ghost=undef, id="", pchannel=undef)
    clean_calculation_files(basename=basename)
    inname = "$(basename).inp"
    outname= "$(basename).out"
    cmd = pipeline(`$(cal.calculator.executable) $(inname)`, outname)

    ts = time()
    open(inname,"w") do io
        write_input(io, cal, point, ghost=ghost)
    end
    run(cmd)
    out = read_energy(outname)
    te = time()
    @debug "Calculation done in $(round(te-ts, digits=1)) seconds"
    pchannel != undef && put!(pchannel,true)
    return out
end

"""
    bsse_corrected_energy(cal::Calculator, c1, c2; basename="base", id="", pchannel=pchannel)

Calculates energy of combined clusters taking into account basis set superposition error

# Arguments
- `cal::Calculator` : calcualation information
- `c1` : collection of clusters
- `c2` : collection of clusters

# Keywords
- `basename="base"` : base name for input/ouput files
- `id=""` : additional information for calculator - needed for multiprocessing in same folder
- `pchannel=undef`  : (Remote)Channel where progess information is added
"""
function bsse_corrected_energy(cal::Calculator{Orca}, c1, c2; basename="base", id="",
                                pchannel=undef)
    # expects ORCA calculator
    points = c1 .+ c2
    e = calculate_energy(cal, points, basename=basename, id=id, pchannel=pchannel)
    l1 = length(c1[1])
    l2 = length(c2[1]) + l1
    bsse1 = calculate_energy(cal, points, basename=basename, ghost=1:l1, id=id, pchannel=pchannel)
    bsse2 = calculate_energy(cal, points, basename=basename, ghost=(l1+1):l2, id=id, pchannel=pchannel)
    return e .- bsse1 .- bsse2
end

"""
    bsse_corrected_energy(cal::Calculator, c1::Cluster, c2::Cluster; basename="base", id="", pchannel=pchannel)

Calculates energy of combined cluster taking into account basis set superposition error

# Arguments
- `cal::Calculator` : calcualation information
- `c1::Cluster` : cluster
- `c2::Cluster` : cluster

# Keywords
- `basename="base"` : base name for input/ouput files
- `id=""` : additional information for calculator - needed for multiprocessing in same folder
- `pchannel=undef`  : (Remote)Channel where progess information is added
"""
function bsse_corrected_energy(cal::Calculator{Orca}, c1::Cluster, c2::Cluster;
                               basename="base", id="", pchannel=undef)
    # expects ORCA calculator
    points = c1 + c2
    e = calculate_energy(cal, points, basename=basename, id=id, pchannel=pchannel)
    l1 = length(c1)
    l2 = length(c2) + l1
    bsse1 = calculate_energy(cal, points, basename=basename, ghost=1:l1, id=id, pchannel=pchannel)
    bsse2 = calculate_energy(cal, points, basename=basename, ghost=(l1+1):l2, id=id, pchannel=pchannel)
    return e .- bsse1 .- bsse2
end

"""
    clean_calculation_files(;dir=".", basename="base")

Deletes files that calculations produced

# Arguments
- `dir="."` : directory to be cleaned
- `basename="base"` : base name used in calculations
"""
function clean_calculation_files(;dir=".", basename="base")
    filenames=readdir(dir)
    # There is a  possible exploit here related to basename
    i = map( x -> occursin(basename*".", x) || occursin(basename*"_", x), filenames)
    rm.(filenames[i])
end

end #module


text:
import Base: hash, ==

mutable struct Partial{T}
  value
end

struct Const{T}
  value::T
  Const(x) = new{Core.Typeof(x)}(x)
end

mutable struct Shape{T}
  size::NTuple{N,Int} where N
end

struct Node{T}
  value::Variable
end

for T in :[Partial, Const, Shape, Node].args
  @eval hash(x::$T, h::UInt64) = hash($T, hash(getfield(x, 1), h))
  @eval x::$T == y::$T = getfield(x, 1) == getfield(y, 1)
end

const AType{T} = Union{Type{T},Const{T},Shape{T},Partial{T},Node{T}}

Base.show(io::IO, c::Const) = print(io, "const(", c.value, ")")

widen(::AType{T}) where T = T
Base.eltype(x::Shape) = eltype(widen(x))
Base.size(x::Shape) = x.size
Base.ndims(x::Shape) = length(x.size)

Base.size(x::Const{<:Array}) = size(x.value)
Base.eltype(x::Const{<:Array}) = eltype(x.value)
Base.ndims(x::Const{<:Array}) = ndims(x.value)
Base.size(::AType{<:Number}) = ()
Base.eltype(T::AType{<:Number}) = widen(T)
Base.ndims(T::AType{<:Number}) = 0

ptuple() = Const(())
ptuple(x::Const...) = Const(map(x -> x.value, x))
ptuple(x::Type...) = Tuple{x...}
ptuple(x...) = Partial{Tuple{widen.(x)...}}((x...,))

_union(::Type{Union{}}, T) = T
_union(S::Const, T::Const) = S == T ? S : Union{widen(S),widen(T)}
_union(S::Partial, T::Partial) = S == T ? S : Union{widen(S),widen(T)}
_union(S, T) = Union{widen(S), widen(T)}
_issubtype(S, T::Type) = widen(S) <: T
_issubtype(S, T) = S == T

function instrument(ex)
  isexpr(ex, :new) ? xcall(Mjolnir, :__new__, ex.args...) :
  isexpr(ex, :splatnew) ? xcall(Mjolnir, :__splatnew__, ex.args...) :
  ex
end

function prepare_ir!(ir)
  ir |> IRTools.expand! |> IRTools.explicitbranch!
  for b in ir.blocks
    b.argtypes .= Union{}
    for i in 1:length(b.stmts)
      st = b.stmts[i]
      b.stmts[i] = stmt(st, expr = instrument(st.expr), type = Union{})
    end
  end
  return ir
end

function blockargs!(b, args)
  isempty(args) && return true # TODO be more fine grained
  changed = false
  for i = 1:length(argtypes(b))
    _issubtype(args[i], argtypes(b)[i]) && continue
    argtypes(b)[i] = _union(argtypes(b)[i], args[i])
    changed = true
  end
  return changed
end

# TODO clean up edges when blocks are reset
mutable struct Frame
  ir::IR
  inlined::Vector{Vector{IR}}
  edges::Vector{Any}
  stmts::Vector{Vector{Variable}}
  rettype::AType
end

getblock(fr::Frame, b, follow) =
  follow == 0 ? (block(fr.ir, b), fr.stmts[b]) :
    (block(fr.inlined[b][follow], 1), keys(fr.inlined[b][follow]))

# TODO clear up inlined edges
uninline!(fr::Frame, b, follow) = deleteat!(fr.inlined[b], follow+1:length(fr.inlined[b]))

Frame(ir::IR) = Frame(ir, [IR[] for _ = 1:length(blocks(ir))], [], keys.(blocks(ir)), Union{})

function frame(ir::IR, args...)
  prepare_ir!(ir)
  argtypes(ir) .= args
  return Frame(ir)
end

function inline_bbs!(fr)
  for b in blocks(fr.ir)
    while !isempty(fr.inlined[b.id])
      follow = popfirst!(fr.inlined[b.id])
      @assert length(openbranches(b)) == 1
      args = arguments(openbranches(b)[1])
      env = Dict(zip(arguments(follow), args))
      for (v, st) in follow
        env[v] = push!(b, rename(env, st))
      end
      empty!(branches(b))
      append!(branches(b), rename.((env,), branches(block(follow, 1))))
    end
  end
end

struct Inference
  frames::Dict{Vector{AType},Frame}
  queue::WorkQueue{Any}
  edges::IdDict{Partial,Any}
  primitives
end

struct MutCtx
  inf::Inference
  ip
end

function edge!(cx::MutCtx, x::Partial) # TODO: clear edges?
  push!(get!(cx.inf.edges, x, Set()), cx.ip)
  return
end

function visit!(cx::MutCtx, x::Partial) # TODO: can do this by key
  if haskey(cx.inf.edges, x)
    push!(cx.inf.queue, cx.inf.edges[x]...)
  end
end

exprtype(ir, x) = IRTools.exprtype(ir, x, typeof = Const)
exprtype(ir, x::GlobalRef) = Const(getproperty(x.mod, x.name))

function infercall!(inf, loc, block, ex)
  Ts = exprtype.((block.ir,), ex.args)
  Ts[1] === Const(Base.not_int) && (Ts[1] = Const(!))
  T = mutate(inf.primitives, MutCtx(inf, loc), Ts...)
  T == nothing || return T
  ir = IR(widen.(Ts)...)
  ir == nothing && error("No IR for $(Tuple{widen.(Ts)...})")
  if !haskey(inf.frames, Ts)
    fr = inf.frames[Ts] = frame(ir, Ts...)
    push!(inf.queue, (fr, 1, 0, 1))
  else
    fr = inf.frames[Ts]
  end
  push!(fr.edges, loc)
  return fr.rettype
end

function inferbranch!(inf, fr, block, follow, br)
  ir = IR(IRTools.block(fr.ir, br.block))
  argtypes(ir) .= exprtype.((getblock(fr, block, follow)[1].ir,), arguments(br))
  push!(fr.inlined[block], ir)
  push!(inf.queue, (fr, block, follow+1, 1))
end

function openbranches(bl)
  brs = []
  for br in branches(bl)
    br.condition == nothing && (push!(brs, br); break)
    cond = exprtype(bl.ir, br.condition)
    cond == Const(true) && continue
    cond == Const(false) && (push!(brs, br); break)
    push!(brs, br)
  end
  return brs
end

function step!(inf::Inference)
    frame, b, f, ip = pop!(inf.queue)
    block, stmts = getblock(frame, b, f)
    uninline!(frame, b, f)
    if ip <= length(stmts)
        var = stmts[ip]
        st = block[var]
        st.expr isa QuoteNode && return
        if isexpr(st.expr, :call)
            T = infercall!(inf, (frame, b, f, ip), block, st.expr)
            if T != Union{}
                block.ir[var] = stmt(block[var], type = _union(st.type, T))
                push!(inf.queue, (frame, b, f, ip+1))
            end
        elseif isexpr(st.expr, :inbounds)
            push!(inf.queue, (frame, b, f, ip+1))
        else
            error("Unrecognised expression $(st.expr)")
        end
    elseif (brs = openbranches(block); length(brs) == 1 && !isreturn(brs[1])
            && !(brs[1].block == length(frame.ir.blocks)))
        inferbranch!(inf, frame, b, f, brs[1])
    else
        for br in brs
            if isreturn(br)
                T = exprtype(block.ir, IRTools.returnvalue(block))
                _issubtype(T, frame.rettype) && return
                frame.rettype = _union(frame.rettype, T)
                foreach(loc -> push!(inf.queue, loc), frame.edges)
            else
                args = exprtype.((block.ir,), arguments(br))
                if blockargs!(IRTools.block(frame.ir, br.block), args)
                    push!(inf.queue, (frame, br.block, 0, 1))
                end
            end
        end
    end
    return
end

function infer!(inf::Inference)
    while !isempty(inf.queue)
        step!(inf)
    end
    for (_, fr) in inf.frames
        inline_bbs!(fr)
        fr.ir |> IRTools.Inner.trimblocks!
    end
    return inf
end

function Inference(fr::Frame, P)
    q = WorkQueue{Any}()
    push!(q, (fr, 1, 0, 1))
    Inference(Dict(argtypes(fr.ir)=>fr), q, IdDict(), P)
end

function infer!(P, ir::IR, args...)
    fr = frame(ir, args...)
    inf = Inference(fr, P)
    infer!(inf)
end

function return_type(ir::IR, args...)
    fr = frame(copy(ir), args...)
    inf = Inference(fr, Defaults())
    infer!(inf)
    return fr.rettype
end


text:
"""
    Native binary RData format IO stream wrapper.

    TODO write readers
"""
type NativeIO{T<:IO} <: RDAIO
    sub::T               # underlying IO stream

    NativeIO(io::T) = new(io)
    @compat (::Type{NativeIO}){T<:IO}(io::T) = new{T}(io)
end


text:
import JSON
using  Printf

elements = JSON.Parser.parsefile("e:/dev-jl/pano/data-osm/osm-cz-sk.json ")["elements"]
@printf("Found %d elements.", size(elements)[1])

io = open("e:/dev-jl/pano/data-osm/osm-cz-sk.tsv", "w")
@printf(io,"\"Summit\"\t\"Elevation\"\t\"Latitude\"\t\"Longitude\"\n")

for i in 1:size(elements)[1]
    ele::String = get(elements[i]["tags"], "ele", "")
    if length(ele) == 0
        continue
    end
    name::String = get(elements[i]["tags"], "name", "")
    if length(name) == 0
        continue
    end
    lat = elements[i]["lat"]
    lon = elements[i]["lon"]
    @printf(io, "\"%s\"\t%0.1f\t%0.5f\t%0.5f\n", name, parse(Float64, ele), lat, lon)
end
close(io)

max_stars_repo_name: pavel-perina/panorama-jl

max_issues_repo_name: pavel-perina/panorama-jl

max_forks_repo_name: pavel-perina/panorama-jl


text:
devcount() = (ret = Cint[0]; rt.cudaGetDeviceCount(ret); int(ret[1]))

device() = (ret = Cint[0]; rt.cudaGetDevice(ret); int(ret[1]))
device(dev::Integer) = (rt.cudaSetDevice(dev); dev)

device_reset() = device_reset(device())
if VERSION < v"0.4.0-dev"
    finalize(x) = nothing
end
function device_reset(dev::Integer)
    # Clear all items on this device from cuda_ptrs, so they don't get
    # freed later
    todelete = Any[]
    for (p,pdev) in cuda_ptrs
        if pdev == dev
            finalize(p)
            push!(todelete, p)
        end
    end
    for p in todelete
        delete!(cuda_ptrs, p)
    end
    # Reset the device
    device(dev)
    rt.cudaDeviceReset()
end

device_synchronize() = rt.cudaDeviceSynchronize()

device_properties(dev::Integer) = (aprop = Array(rt.cudaDeviceProp, 1); rt.cudaGetDeviceProperties(aprop, dev); aprop[1])

attribute(dev::Integer, code::Integer) = (ret = Cint[0]; rt.cudaDeviceGetAttribute(ret, code, dev); int(ret[1]))

capability(dev::Integer) = (attribute(dev,rt.cudaDevAttrComputeCapabilityMajor),
                            attribute(dev,rt.cudaDevAttrComputeCapabilityMinor))

name(p::rt.cudaDeviceProp) = bytestring(convert(Ptr{Uint8}, pointer([p.name])))

# criteria = dev -> Bool
function devices(criteria::Function; nmax::Integer = typemax(Int))
    devlist = find(map(criteria, 0:devcount().-1)).-1
    isempty(devlist) && error("No suitable devices found")
    devlist[1:min(nmax, length(devlist))]
end

# do syntax, f(devlist)
function devices(f::Function, criteria::Function; nmax::Integer = typemax(Int))
    devlist = devices(criteria, nmax=nmax)
    devices(f, devlist)
end

function devices(f::Function, devlist::Union(Integer,AbstractVector))
    local ret
    mdutils = [CuModule() for i = 1:length(devlist)]
    try
        init!(mdutils, devlist)
        ret = f(devlist)
    finally
        close!(mdutils, devlist)
    end
    ret
end

# A cache of useful CUDA kernels that gets loaded and closed
# by devices(f, devlist)
const global ptxdict = Dict()

function init!(mdutils::Array{CuModule}, devlist)
    funcnames = ["fill_contiguous", "fill_pitched"]
    funcexts  = ["double","float","int64","uint64","int32","uint32","int16","uint16","int8","uint8"]
    datatypes = [Float64,Float32,Int64,Uint64,Int32,Uint32,Int16,Uint16,Int8,Uint8]
    utilfile  = joinpath(Pkg.dir(), "CUDArt/deps/utils.ptx")
    # initialize all devices
    for idev = 1:length(devlist)
        dev = devlist[idev]
        device(dev)
        # allocate and destroy memory to force initialization
        free(malloc(Uint8, 1))
        # Load the utility functions
        mdutils[idev] = CuModule(utilfile, false)
        for func in funcnames
            for i = 1:length(funcexts)
                ptxdict[(dev, func, datatypes[i])] = CuFunction(mdutils[idev], func*"_"*funcexts[i])
            end
        end
        ptxdict[(dev, "clock_block")] = CuFunction(mdutils[idev], "clock_block")
    end
end

function close!(mdutils::Array{CuModule}, devlist)
    for idev = 1:length(devlist)
        if mdutils[idev].handle != C_NULL
            unload(mdutils[idev])
        end
        device_reset(devlist[idev])
    end
end


text:
function t_data_cells()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

# Simple cell checks
Test.@test A01 isa j4pr.CellDataVec
Test.@test A02 isa j4pr.CellDataVecVec
Test.@test A03 isa j4pr.CellDataVecMat

Test.@test A04 isa j4pr.CellDataMat
Test.@test A05 isa j4pr.CellDataMatVec
Test.@test A06 isa j4pr.CellDataMatMat

Test.@test (A01 isa j4pr.CellDataU) && !(A01 isa j4pr.CellDataL) && !(A01 isa j4pr.CellDataLL)
Test.@test (A04 isa j4pr.CellDataU) && !(A04 isa j4pr.CellDataL) && !(A04 isa j4pr.CellDataLL)

Test.@test !(A02 isa j4pr.CellDataU) && (A02 isa j4pr.CellDataL) && !(A02 isa j4pr.CellDataLL)
Test.@test !(A05 isa j4pr.CellDataU) && (A05 isa j4pr.CellDataL) && !(A05 isa j4pr.CellDataLL)

Test.@test !(A03 isa j4pr.CellDataU) && !(A03 isa j4pr.CellDataL) && (A03 isa j4pr.CellDataLL)
Test.@test !(A06 isa j4pr.CellDataU) && !(A06 isa j4pr.CellDataL) && (A06 isa j4pr.CellDataLL)

# Test data/observation concatenation
v=[1,2,3]; m=rand(3,3)
Test.@test j4pr.dcat((nothing,nothing)) == nothing
Test.@test j4pr.dcat(nothing) == nothing
Test.@test j4pr.dcat(nothing,nothing) == nothing
Test.@test j4pr.dcat(nothing,v) == v
Test.@test j4pr.dcat(v,nothing) == v
Test.@test j4pr.dcat(v) == v
Test.@test j4pr.dcat(view(v,:)) == v
Test.@test j4pr.dcat(m) == m
Test.@test j4pr.dcat(view(m,:,:)) == m
Test.@test j4pr.dcat(v,v) == [v';v']
Test.@test j4pr.dcat(v,m) == [v';m]
Test.@test j4pr.dcat(m,v) == [m;v']
Test.@test j4pr.dcat(m,m) == [m;m]
Test.@test j4pr.dcat(m,m,v) == [m;m;v']

Test.@test j4pr.ocat((nothing,nothing)) == nothing
Test.@test j4pr.ocat(nothing) == nothing
Test.@test j4pr.ocat(nothing,nothing) == nothing
Test.@test j4pr.ocat(nothing,v) == v
Test.@test j4pr.ocat(v,nothing) == v
Test.@test j4pr.ocat(v) == v
Test.@test j4pr.ocat(view(v,:)) == v
Test.@test j4pr.ocat(m) == m
Test.@test j4pr.ocat(view(m,:,:)) == m
Test.@test j4pr.ocat(v,v) == [v;v]
Test.@test j4pr.ocat(v,m) == [v m]
Test.@test j4pr.ocat(m,v) == [m v]
Test.@test j4pr.ocat(m,m) == [m m]
Test.@test j4pr.ocat(m,m,v) == [m m v]

# Test iteration
Test.@test try
	for i in j4pr.datacell(m,v) i end
	true
catch
	false
end

end