text:
function compute_error(s, cpt_err)
    if cpt_err[1] == :across
        GroupedErrors._across(s, cpt_err[2])
    elseif cpt_err[1] == :bootstrap
        GroupedErrors._bootstrap(s, cpt_err[2])
    else
        error("Compute error $cpt_err not supported")
    end
end

compute_error(s, ::Void) = s

function smoothing_kwargs(a::Analysis)
    a.axis_type != :continuous && return []
    a.y in colnames(a.data.table) && return [(:span, (a.smoother+1.0)/100)]
    a.y in [:density, :hazard] && return [(:bandwidth, (a.smoother+1.0)*std(column(a.data.table, a.x))/200)]
    return []
end

function process(::Type{GroupedError}, a::Analysis)
    s = GroupedErrors.ColumnSelector(a.data.table)
    s = GroupedErrors._splitby(s, Symbol[a.data.splitby...])
    s = compute_error(s, a.compute_error)
    if a.axis_type != :pointbypoint
        maybe_nbins = a.axis_type == :binned ? (round(Int64, 101-a.smoother),) : ()
        s = GroupedErrors._x(s, a.x, a.axis_type, maybe_nbins...)
        y = a.y in colnames(a.data.table) ? (:locreg, a.y) : (a.y,)
        s = GroupedErrors._y(s, y...; smoothing_kwargs(a)...)
    else
        s = GroupedErrors._x(s, a.x, a.xfunc)
        s = GroupedErrors._y(s, a.y, a.yfunc)
    end
    @plot s a.plot()
end


text:
function loss()
    mse(ŷ, y, ∂) = loss_mse(ŷ, y, ∂)
    quadratic(ŷ, y, ∂ = false) = loss_quadratic(ŷ, y, ∂)
    ce(ŷ, y, ∂ = false) = loss_cross_entropy(ŷ, y, ∂)
    exp(ŷ, y, ∂ = false) = loss_exponential(ŷ, y, ∂)
    () -> (mse, quadratic, ce, exp)
end

function loss_mse(Ŷ, Y, ∂ = false)
    if ∂ == true
        if length(size(Y)) > 1
            return (2/size(Y)[2] * (Ŷ - Y))
        else
            return (2/size(Y)[1] * (Ŷ - Y))
        end
    else
        return 1/length(Y) * sum(abs2, Y .- Ŷ)
    end
end

function loss_quadratic(Ŷ, Y, ∂ = false)
    return 1/2 * sum(abs2, Y .- Ŷ)
end

function loss_cross_entropy(Ŷ, Y, ∂ = false)
    return 0
end

function loss_exponential(Ŷ, Y, ∂ = false)
    return 0
end

ext: jl


text:
using Logging, WeightsAndBiasLogger

args = (n_epochs=1_000, lr=1e-3)

logger = WBLogger(project="sample-project")

config!(logger, args)

with(logger) do
    loss = 0
    for i in 1:args.n_epochs
        loss += randn() * args.lr
        @info "train" i=i loss=loss
    end
end


text:
export DFT, IDFT

abstract type FourierTransform{N,C,D,T1,T2} <: LinearOperator end

"""
`DFT([domainType=Float64::Type,] dim_in::Tuple [,dims])`

`DFT(dim_in...)`

`DFT(x::AbstractArray [,dims])`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Discrete Fourier Transform over dimensions `dims` of `x`.

```julia
julia> DFT(Complex{Float64},(10,10))
ℱ  ℂ^(10, 10) -> ℂ^(10, 10)

julia> DFT(10,10)
ℱ  ℝ^(10, 10) -> ℂ^(10, 10)

julia> A = DFT(ones(3))
ℱ  ℝ^3 -> ℂ^3

julia> A*ones(3)
3-element Array{Complex{Float64},1}:
 3.0+0.0im
 0.0+0.0im
 0.0+0.0im
```

"""
struct DFT{N,
	      C<:RealOrComplex,
	      D<:RealOrComplex,
	      T1<:AbstractFFTs.Plan,
	      T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

"""
`IDFT([domainType=Float64::Type,] dim_in::Tuple [,dims])`

`IDFT(dim_in...)`

`IDFT(x::AbstractArray [,dims])`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Inverse Discrete Fourier Transform over dimensions `dims` of `x`.

```julia
julia> IDFT(Complex{Float64},(10,10))
ℱ⁻¹  ℂ^(10, 10) -> ℂ^(10, 10)

julia> IDFT(10,10)
ℱ⁻¹ ℝ^(10, 10) -> ℂ^(10, 10)

julia> A = IDFT(ones(3))
ℱ⁻¹  ℝ^3 -> ℂ^3

julia> A*ones(3)
3-element Array{Complex{Float64},1}:
 1.0+0.0im
 0.0+0.0im
 0.0+0.0im
```

"""
struct IDFT{N,
	       C<:RealOrComplex,
	       D<:RealOrComplex,
	       T1<:AbstractFFTs.Plan,
	       T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

# Constructors
#standard constructor
DFT(dim_in::NTuple{N,Int},dims=1:N) where {N} = DFT(zeros(dim_in),dims)

function DFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Real}
	A,At = plan_fft(x,dims), plan_bfft(fft(x,dims),dims)
	DFT{N,Complex{D},D,typeof(A),typeof(At)}(size(x),A,At)
end

function DFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Complex}
	A,At = plan_fft(x,dims), plan_bfft(fft(x,dims),dims)
	DFT{N,D,D,typeof(A),typeof(At)}(size(x),A,At)
end

DFT(T::Type,dim_in::NTuple{N,Int},dims=1:N) where {N} = DFT(zeros(T,dim_in),dims)
DFT(dim_in::Vararg{Int}) = DFT(dim_in)
DFT(T::Type,dim_in::Vararg{Int}) = DFT(T,dim_in)

function IDFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Real}
	A,At = plan_ifft(x,dims), plan_fft(ifft(x,dims),dims)
	IDFT{N,Complex{D},D,typeof(A),typeof(At)}(size(x),A,At)
end

#standard constructor
IDFT(T::Type,dim_in::NTuple{N,Int},dims=1:N) where {N} = IDFT(zeros(T,dim_in),dims)

function IDFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Complex}
	A,At = plan_ifft(x,dims), plan_fft(ifft(x,dims),dims)
	IDFT{N,D,D,typeof(A),typeof(At)}(size(x),A,At)
end

IDFT(dim_in::NTuple{N,Int},dims=1:N) where {N} = IDFT(zeros(dim_in),dims)
IDFT(dim_in::Vararg{Int}) = IDFT(dim_in)
IDFT(T::Type,dim_in::Vararg{Int}) = IDFT(T,dim_in)

# Mappings

function mul!(y::AbstractArray{C,N},
              L::DFT{N,C,C,T1,T2},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A,b)
end

function mul!(y::AbstractArray{C,N},
              L::DFT{N,C,D,T1,T2},
              b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2}
	mul!(y,L.A,complex(b))
end

function mul!(y::AbstractArray{C,N},
              L::AdjointOperator{DFT{N,C,C,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A.At,b)
end

function mul!(y::AbstractArray{D,N},
              L::AdjointOperator{DFT{N,C,D,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2}
	y2 = complex(y)
	mul!(y2,L.A.At,b)
	y .= real.(y2)
end

function mul!(y::AbstractArray{C,N},
              L::IDFT{N,C,C,T1,T2},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A,b)
end

function mul!(y::AbstractArray{C,N},
              L::IDFT{N,C,D,T1,T2},
              b::AbstractArray{D,N}) where {N,C<:Complex,D<:Real,T1,T2}
	mul!(y,L.A,complex(b))
end

function mul!(y::AbstractArray{C,N},
              L::AdjointOperator{IDFT{N,C,C,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A.At,b)
	y .*= L.A.A.scale
end

function mul!(y::AbstractArray{D,N},
              L::AdjointOperator{IDFT{N,C,D,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,D<:Real,T1,T2}

	y2 = complex(y)
	mul!(y2,L.A.At,b)
	y .= real.(y2).*L.A.A.scale
end

# Properties

size(L::FourierTransform) = (L.dim_in,L.dim_in)

fun_name(A::DFT) = "ℱ"
fun_name(A::IDFT) = "ℱ⁻¹"

domainType(L::FourierTransform{N,C,D,T1,T2}) where {N,C,D,T1,T2} = D
codomainType(L::FourierTransform{N,C,D,T1,T2}) where {N,C,D,T1,T2} = C

is_AcA_diagonal(L::FourierTransform)     = true
is_AAc_diagonal(L::FourierTransform)     = true
is_invertible(L::FourierTransform)       = true
is_full_row_rank(L::FourierTransform)    = true
is_full_column_rank(L::FourierTransform) = true

diag_AcA(L::DFT) = float(prod(size(L,1)))
diag_AAc(L::DFT) = float(prod(size(L,1)))

diag_AcA(L::IDFT) = 1/prod(size(L,1))
diag_AAc(L::IDFT) = 1/prod(size(L,1))


text:
################################################################################
#
# convenience ...
#
################################################################################

#CF: we have to "extend" AnticNumberField as NumberField is just defined by
#    NumberField = AnticNumberField in Nemo.
#    Possibly should be replaced by having optional 2nd arguments?
doc"""
***
  NumberField(f::fmpq_poly)

>  The number field Q[x]/f generated by f.
"""
function AnticNumberField(f::fmpq_poly)
  return NumberField(f, "_a")
end

function AnticNumberField(f::fmpz_poly, s::Symbol)
  Qx, x = PolynomialRing(QQ, string(parent(f).S))
  return NumberField(Qx(f), s)
end

function AnticNumberField(f::fmpz_poly, s::AbstractString)
  Qx, x = PolynomialRing(QQ, string(parent(f).S))
  return NumberField(Qx(f), s)
end

function AnticNumberField(f::fmpz_poly)
  Qx, x = PolynomialRing(QQ, string(parent(f).S))
  return NumberField(Qx(f))
end

################################################################################
#
#  Base case for dot products
#
################################################################################

dot(x::fmpz, y::nf_elem) = x*y

dot(x::nf_elem, y::fmpz) = x*y

################################################################################
# given a basis (an array of elements), get a linear combination with
# random integral coefficients
################################################################################
doc"""
***
  rand(b::Array{nf_elem,1}, r::UnitRange)

> A random linear combination of elements in b with coefficients in r
"""
function rand(b::Array{nf_elem,1}, r::UnitRange)
  length(b) == 0 && error("Array must not be empty")

  s = zero(b[1].parent)
  t = zero(b[1].parent)

  for i = 1:length(b)
    mul!(t, b[i], rand(r))
    add!(s, t, s)
  end
  return s
end

function rand!(c::nf_elem, b::Array{nf_elem,1}, r::UnitRange)
  length(b) == 0 && error("Array must not be empty")

  mul!(c, b[1], rand(r))
  t = zero(b[1].parent)

  for i = 2:length(b)
    mul!(t, b[i], rand(r))
    add!(c, t, c)
  end
  nothing
end

################################################################################
#
#  fmpq_poly with denominator 1 to fmpz_poly
#
################################################################################

#function Base.call(a::Hecke.FmpqPolyRing, b::String)
#  return a(1)
#end

function Base.call(a::FmpzPolyRing, b::fmpq_poly)
  (den(b) != 1) && error("denominator has to be 1")
  z = a()
  ccall((:fmpq_poly_get_numerator, :libflint), Void,
              (Ptr{fmpz_poly}, Ptr{fmpq_poly}), &z, &b)
  return z
end

doc"""
  basis(K::AnticNumberField)

> A Q-basis for K, ie. 1, x, x^2, ... as elements of K
"""
function basis(K::AnticNumberField)
  n = degree(K)
  g = gen(K);
  d = Array(typeof(g), n)
  b = K(1)
  for i = 1:n-1
    d[i] = b
    b *= g
  end
  d[n] = b
  return d
end

doc"""
  representation_mat(a::nf_elem) -> fmpz_mat

> The right regular representation of a, i.e. the matrix representing
> the multiplication by a map on the number field.
> den(a) must be one
"""
function representation_mat(a::nf_elem)
  @assert den(a) == 1
  dummy = fmpz(0)
  n = degree(a.parent)
  M = MatrixSpace(FlintZZ, n,n)()::fmpz_mat
  t = gen(a.parent)
  b = deepcopy(a)
  for i = 1:n-1
    elem_to_mat_row!(M, i, dummy, b)
    mul!(b, b, t) ## CF: should write and use mul_gen which is faster
  end
  elem_to_mat_row!(M, n, dummy, b)
  return M
end

function basis_mat(A::Array{nf_elem, 1})
  @assert length(A) > 0
  n = length(A)
  d = degree(parent(A[1]))

  M = MatrixSpace(FlintZZ, n, d)()

  deno = one(FlintZZ)
  dummy = one(FlintZZ)

  for i in 1:n
    deno = lcm(deno, den(A[i]))
  end

  for i in 1:n
    elem_to_mat_row!(M, i, dummy, A[i])
    for j in 1:d
      M[i, j] = divexact(deno, dummy) * M[i, j]
    end
  end
  return FakeFmpqMat(M, deno)
end

function set_den!(a::nf_elem, d::fmpz)
  ccall((:nf_elem_set_den, :libflint),
        Void,
       (Ptr{Nemo.nf_elem}, Ptr{Nemo.fmpz}, Ptr{Nemo.AnticNumberField}),
       &a, &d, &parent(a))
end

doc"""
***
  charpoly(a::nf_elem) -> fmpz_poly

> The characteristic polynomial of a.
"""
function charpoly(a::nf_elem)
  d = den(a)
  Zx = PolynomialRing(ZZ, string(parent(parent(a).pol).S))[1]
  f = charpoly(Zx, representation_mat(d*a))
  return f(gen(parent(f))*d)
end

doc"""
***
  minpoy(a::nf_elem) -> fmpz_poly

> The minimal polynomial of a.
"""
function minpoly(a::nf_elem)
  d = den(a)
  Zx = PolynomialRing(ZZ, string(parent(parent(a).pol).S))[1]
  f = minpoly(Zx, representation_mat(d*a))
  return f(gen(parent(f))*d)
end
###########################################################################
###########################################################################

function nf_poly_to_xy(f::PolyElem{Nemo.nf_elem}, x::PolyElem, y::PolyElem)
   K = base_ring(f)
   Qy = parent(K.pol)

   res = zero(parent(y))
   for i=degree(f):-1:0
     res *= y
     res += evaluate(Qy(coeff(f, i)), x)
   end
   return res
 end

doc"""
  norm(f::PolyElem{nf_elem}) -> fmpq_poly

> The norm of f, i.e. the product of all conjugates of f taken coefficientwise.
"""
function norm(f::PolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(f)
  Qy = parent(K.pol)
  y = gen(Qy)
  Qyx, x = PolynomialRing(Qy, "x")

  Qx = PolynomialRing(QQ, "x")[1]
  Qxy = PolynomialRing(Qx, "y")[1]

  T = evaluate(K.pol, gen(Qxy))
  h = nf_poly_to_xy(f, gen(Qxy), gen(Qx))
  return resultant(T, h)
end

doc"""
  factor(f::PolyElem{nf_elem}) -> Dict{PolyElem{nf_elem}, Int}

> The factorisation of f (using Trager's method).
"""
function factor(f::PolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(f)

  k = 0
  g = f
  N = 0

  while true
    N = norm(f)

    if !is_constant(N) && is_squarefree(N)
      break
    end

    k = k + 1
    g = compose(f, gen(Kx) - k*gen(K))
  end

  fac = factor(N)

  res = Dict{PolyElem{nf_elem}, Int64}()

  for i in keys(fac)
    t = zero(Kx)
    for j in 0:degree(i)
      t = t + K(coeff(i, j))*gen(Kx)^j
    end
    t = compose(t, gen(Kx) + k*gen(K))
    res[gcd(f, t)] = 1
  end

  return res
end

################################################################################
#
# Operations for nf_elem
#
################################################################################

function gen!(r::nf_elem)
   a = parent(r)
   ccall((:nf_elem_gen, :libflint), Void,
         (Ptr{nf_elem}, Ptr{AnticNumberField}), &r, &a)
   return r
end

function one!(r::nf_elem)
   a = parent(r)
   ccall((:nf_elem_one, :libflint), Void,
         (Ptr{nf_elem}, Ptr{AnticNumberField}), &r, &a)
   return r
end

function zero!(r::nf_elem)
   a = parent(r)
   ccall((:nf_elem_zero, :libflint), Void,
         (Ptr{nf_elem}, Ptr{AnticNumberField}), &r, &a)
   return r
end

function one(r::nf_elem)
   a = parent(r)
   return one(a)
end

function zero(r::nf_elem)
   return zero(parent(r))
end

*(a::nf_elem, b::Integer) = a * fmpz(b)

doc"""
***
   norm_div(a::nf_elem, d::fmpz, nb::Int) -> fmpz

> Computes divexact(norm(a), d) provided the result has at most nb bits.
> Typically, a is in some ideal and d is the norm of the ideal.
"""
function norm_div(a::nf_elem, d::fmpz, nb::Int)
   z = fmpq()
   #CF the resultant code has trouble with denominators,
   #   this "solves" the problem, but it should probably be
   #   adressed in c
   de = den(a)
   n = degree(parent(a))
   ccall((:nf_elem_norm_div, :libflint), Void,
         (Ptr{fmpq}, Ptr{nf_elem}, Ptr{AnticNumberField}, Ptr{fmpz}, UInt),
         &z, &(a*de), &a.parent, &(d*de^n), UInt(nb))
   return z
end

function sub!(a::nf_elem, b::nf_elem, c::nf_elem)
   ccall((:nf_elem_sub, :libflint), Void,
         (Ptr{nf_elem}, Ptr{nf_elem}, Ptr{nf_elem}, Ptr{AnticNumberField}),

         &a, &b, &c, &a.parent)
end

function ^(x::nf_elem, y::fmpz)
  if y < 0
    return inv(x)^(-y)
  elseif y == 0
    return parent(x)(1)
  elseif y == 1
    return deepcopy(x)
  elseif mod(y, 2) == 0
    z = x^(div(y, 2))
    return z*z
  elseif mod(y, 2) == 1
    return x^(y-1) * x
  end
end

doc"""
***
    root(a::nf_elem, n::Int) -> Bool, nf_elem

> Determines whether $a$ has an $n$-th root. If this is the case,
> the root is returned.
"""
function root(a::nf_elem, n::Int)
  #println("Compute $(n)th root of $a")
  Kx, x = PolynomialRing(parent(a), "x")

  f = x^n - a

  fac = factor(f)
  #println("factorization is $fac")

  for i in keys(fac)
    if degree(i) == 1
      return (true, -coeff(i, 0)//coeff(i, 1))
    end
  end

  return (false, zero(parent(a)))
end

function num(a::nf_elem)
   const _one = fmpz(1)
   z = copy(a)
   ccall((:nf_elem_set_den, :libflint), Void,
         (Ptr{nf_elem}, Ptr{fmpz}, Ptr{AnticNumberField}),
         &z, &_one, &a.parent)
   return z
end

copy(d::nf_elem) = deepcopy(d)

################################################################################
#
#  Minkowski map
#
################################################################################

doc"""
***
    minkowski_map(a::nf_elem, abs_tol::Int) -> Array{arb, 1}

> Returns the image of $a$ under the Minkowski embedding.
> Every entry of the array returned is of type `arb` with radius less then
> `2^(-abs_tol)`.
"""
function minkowski_map(a::nf_elem, abs_tol::Int = 32)
  K = parent(a)
  A = Array(arb, degree(K))
  r, s = signature(K)
  c = conjugate_data_arb(K)
  R = PolynomialRing(AcbField(c.prec), "x")[1]
  f = R(parent(K.pol)(a))
  CC = AcbField(c.prec)
  T = PolynomialRing(CC, "x")[1]
  g = T(f)

  for i in 1:r
    t = evaluate(g, c.real_roots[i])
    @assert isreal(t)
    A[i] = real(t)
    if !radiuslttwopower(A[i], -abs_tol)
      refine(c)
      return minkowski_map(a, abs_tol)
    end
  end

  t = base_ring(g)()

  for i in 1:s
    t = evaluate(g, c.complex_roots[i])
    t = sqrt(CC(2))*t
    if !radiuslttwopower(t, -abs_tol)
      refine(c)
      return minkowski_map(a, abs_tol)
    end
    A[r + 2*i - 1] = real(t)
    A[r + 2*i] = imag(t)
  end

  return A
end

function t2{T}(x::nf_elem, abs_tol::Int = 32, ::Type{T} = arb)
  p = 2*abs_tol
  z = mapreduce(y -> y^2, +, minkowski_map(x, p))
  while !radiuslttwopower(z, -abs_tol)
    p = 2 * p
    z = mapreduce(y -> y^2, +, minkowski_map(x, p))
  end
  return z
end

################################################################################
#
#  Conjugates and real embeddings
#
################################################################################

doc"""
***
    conjugates_arb(x::nf_elem, abs_tol::Int) -> Array{acb, 1}

> Compute the the conjugates of `x` as elements of type `acb`.
> Recall that we order the complex conjugates
> $\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
> $\sigma_{i}(x) = \overline{sigma_{i + s}(x)}$ for $r + 1 \leq i \leq r + s$.
>
> Every entry `y` of the array returned satisfies
> `radius(real(y)) < 2^abs_tol` and `radius(imag(y)) < 2^abs_tol` respectively.
"""
function conjugates_arb(x::nf_elem, abs_tol::Int = 32)
  K = parent(x)
  d = degree(K)
  c = conjugate_data_arb(K)
  r, s = signature(K)
  conjugates = Array(acb, r + 2*s)
  CC = AcbField(c.prec)

  for i in 1:r
    conjugates[i] = CC(evaluate(parent(K.pol)(x), c.real_roots[i]))
    if !isfinite(conjugates[i]) || !radiuslttwopower(conjugates[i], -abs_tol)
      refine(c)
      return conjugates_arb(x, abs_tol)
    end
  end

  for i in 1:s
    conjugates[r + i] = evaluate(parent(K.pol)(x), c.complex_roots[i])
    if !isfinite(conjugates[i]) || !radiuslttwopower(conjugates[i], -abs_tol)
      refine(c)
      return conjugates_arb(x, abs_tol)
    end
    conjugates[r + i + s] = Nemo.conj(conjugates[r + i])
  end

  return conjugates
end

doc"""
***
    conjugates_arb_real(x::nf_elem, abs_tol::Int) -> Array{arb, 1}

> Compute the the real conjugates of `x` as elements of type `arb`.
>
> Every entry `y` of the array returned satisfies
> `radius(y) < 2^abs_tol`.
"""
function conjugates_arb_real(x::nf_elem, abs_tol::Int = 32)
  r1, r2 = signature(parent(x))
  c = conjugates_arb(x, abs_tol)
  z = Array(arb, r1)

  for i in 1:r
    z[i] = real(c[i])
  end

  return z
end

doc"""
***
    conjugates_arb_complex(x::nf_elem, abs_tol::Int) -> Array{acb, 1}

> Compute the the complex conjugates of `x` as elements of type `acb`.
> Recall that we order the complex conjugates
> $\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
> $\sigma_{i}(x) = \overline{sigma_{i + s}(x)}$ for $r + 1 \leq i \leq r + s$.
>
> Every entry `y` of the array returned satisfies
> `radius(real(y)) < 2^abs_tol` and `radius(imag(y)) < 2^abs_tol`.
"""
function conjugates_arb_complex(x::nf_elem, abs_tol::Int)
end

doc"""
***
    conjugates_arb_log(x::nf_elem, abs_tol::Int) -> Array{arb, 1}

> Returns the elements
> $(\log(\lvert \sigma_1(x) \rvert),\dotsc,\log(\lvert\sigma_r(x) \rvert),
> \dotsc,2\log(\lvert \sigma_{r+1}(x) \rvert),\dotsc,
> 2\log(\lvert \sigma_{r+s}(x)\rvert))$ as elements of type `arb` radius
> less then `2^abs_tol`.
"""
function conjugates_arb_log(x::nf_elem, abs_tol::Int)
  K = parent(x)
  d = degree(K)
  r1, r2 = signature(K)
  c = conjugate_data_arb(K)

  # We should replace this using multipoint evaluation of libarb
  z = Array(arb, r1 + r2)
  for i in 1:r1
    z[i] = log(abs(evaluate(parent(K.pol)(x),c.real_roots[i])))
    if !isfinite(z[i]) || !radiuslttwopower(z[i], -abs_tol)
      refine(c)
      return conjugates_arb_log(x, abs_tol)
    end
  end
  for i in 1:r2
    z[r1 + i] = 2*log(abs(evaluate(parent(K.pol)(x), c.complex_roots[i])))
    if !isfinite(z[r1 + i]) || !radiuslttwopower(z[r1 + i], -abs_tol)
      refine(c)
      return conjugates_arb_log(x, abs_tol)
    end
  end
  return z
end

function conjugates_arb_log(x::nf_elem, R::ArbField)
  z = conjugates_arb_log(x, R.prec)
  return map(R, z)
end

################################################################################
#
#  Torsion units and related functions
#
################################################################################

doc"""
***
    is_torsion_unit(x::nf_elem, checkisunit::Bool = false) -> Bool

> Returns whether $x$ is a torsion unit, that is, whether there exists $n$ such
> that $x^n = 1$.
>
> If `checkisunit` is `true`, it is first checked whether $x$ is a unit of the
> maximal order of the number field $x$ is lying in.
"""
function is_torsion_unit(x::nf_elem, checkisunit::Bool = false)
  if checkisunit
    _is_unit(x) ? nothing : return false
  end

  K = parent(x)
  d = degree(K)
  c = conjugate_data_arb(K)
  r, s = signature(K)

  while true
    @vprint :UnitGroup 2 "Precision is now $(c.prec) \n"
    l = 0
    @vprint :UnitGroup 2 "Computing conjugates ... \n"
    cx = conjugates_arb(x, c.prec)
    A = ArbField(c.prec)
    for i in 1:r
      k = abs(cx[i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end
    for i in 1:s
      k = abs(cx[r + i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end

    if l == r + s
      return true
    end
    refine(c)
  end
end

doc"""
***
    torsion_unit_order(x::nf_elem, n::Int)

> Given a torsion unit $x$ together with a multiple $n$ of its order, compute
> the order of $x$, that is, the smallest $k \in \mathbb Z_{\geq 1}$ such
> that $x^`k` = 1$.
>
> It is not checked whether $x$ is a torsion unit.
"""
function torsion_unit_order(x::nf_elem, n::Int)
  # This is lazy
  # Someone please change this
  y = deepcopy(x)
  for i in 1:n
    if y == 1
      return i
    end
    mul!(y, y, x)
  end
  error("Something odd in the torsion unit order computation")
end

################################################################################
#
#  Serialization
#
################################################################################

# This function can be improved by directly accessing the numerator
# of the fmpq_poly representing the nf_elem
doc"""
***
    write(io::IO, A::Array{nf_elem, 1}) -> Void

> Writes the elements of `A` to `io`. The first line are the coefficients of
> the defining polynomial of the ambient number field. The following lines
> contain the coefficients of the elements of `A` with respect to the power
> basis of the ambient number field.
"""
function write(io::IO, A::Array{nf_elem, 1})
  if length(A) == 0
    return
  else
    # print some useful(?) information
    print(io, "# File created by Hecke $VERSION_NUMBER, $(Base.Dates.now()), by function 'write'\n")
    K = parent(A[1])
    polring = parent(K.pol)

    # print the defining polynomial
    g = K.pol
    d = den(g)

    for j in 0:degree(g)
      print(io, coeff(g, j)*d)
      print(io, " ")
    end
    print(io, d)
    print(io, "\n")

    # print the elements
    for i in 1:length(A)

      f = polring(A[i])
      d = den(f)

      for j in 0:degree(K)-1
        print(io, coeff(f, j)*d)
        print(io, " ")
      end

      print(io, d)

      print(io, "\n")
    end
  end
end

doc"""
***
    write(file::ASCIIString, A::Array{nf_elem, 1}, flag::ASCIString = "w") -> Void

> Writes the elements of `A` to the file `file`. The first line are the coefficients of
> the defining polynomial of the ambient number field. The following lines
> contain the coefficients of the elements of `A` with respect to the power
> basis of the ambient number field.
>
> Unless otherwise specified by the parameter `flag`, the content of `file` will be
> overwritten.
"""
function write(file::ASCIIString, A::Array{nf_elem, 1}, flag::ASCIIString = "w")
  f = open(file, flag)
  write(f, A)
  close(f)
end

# This function has a bad memory footprint
doc"""
***
    read(io::IO, K::AnticNumberField, ::Type{nf_elem}) -> Array{nf_elem, 1}

> Given a file with content adhering the format of the `write` procedure,
> this functions returns the corresponding object of type `Array{nf_elem, 1}` such that
> all elements have parent $K$.

**Example**

    julia> Qx, x = QQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(io::IO, K::AnticNumberField, ::Type{Hecke.nf_elem})
  Qx = parent(K.pol)

  A = Array{nf_elem, 1}()

  i = 1

  for ln in eachline(io)
    if ln[1] == '#'
      continue
    elseif i == 1
      # the first line read should contain the number field and will be ignored
      i = i + 1
    else
      coe = map(Hecke.fmpz, split(ln, " "))
      t = fmpz_poly(Array(slice(coe, 1:(length(coe) - 1))))
      t = Qx(t)
      t = divexact(t, coe[end])
      push!(A, K(t))
      i = i + 1
    end
  end

  return A
end

doc"""
***
    read(file::ASCIIString, K::AnticNumberField, ::Type{nf_elem}) -> Array{nf_elem, 1}

> Given a file with content adhering the format of the `write` procedure,
> this functions returns the corresponding object of type `Array{nf_elem, 1}` such that
> all elements have parent $K$.

**Example**

    julia> Qx, x = QQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(file::ASCIIString, K::AnticNumberField, ::Type{Hecke.nf_elem})
  f = open(file, "r")
  A = read(f, K, Hecke.nf_elem)
  close(f)
  return A
end

function dot(a::Array{nf_elem, 1}, b::Array{fmpz, 1})
  d = zero(parent(a[1]))
  t = zero(d)
  for i=1:length(a)
    Nemo.mul!(t, a[i], b[i])
    Nemo.add!(d, d, t)
  end
  return d
end

type nf_elem_deg_1_raw
  num::Int  ## fmpz!
  den::Int
end

type nf_elem_deg_2_raw
  nu0::Int  ## fmpz - actually an fmpz[3]
  nu1::Int
  nu2::Int
  den::Int
end

type nf_elem_deg_n_raw  #actually an fmpq_poly_raw
  A::Ptr{Int} # fmpz
  den::Int # fmpz
  alloc::Int
  len::Int
end

type nmod_t
  n::Int
  ni::Int
  norm::Int
end

#nf_elem is a union of the three types above
#ignores the denominator completely

function nf_elem_to_nmod_poly_no_den!(r::nmod_poly, a::nf_elem)
  d = degree(a.parent)
  zero!(r)
  p = r.mod_n
  if d == 1
    ra = pointer_from_objref(a)
    s = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Void}, UInt), &, p)
    ccall((:nmod_poly_set_coeff_ui, :libflint), Void, (Ptr{nmod_poly}, Int, UInt), &r, 0, s)
  elseif d == 2
    ra = pointer_from_objref(a)
    s = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Void}, UInt), ra, p)
    ccall((:nmod_poly_set_coeff_ui, :libflint), Void, (Ptr{nmod_poly}, Int, UInt), &r, 0, s)
    s = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Void}, UInt), ra + sizeof(Int), p)
    ccall((:nmod_poly_set_coeff_ui, :libflint), Void, (Ptr{nmod_poly}, Int, UInt), &r, 1, s)
  else
    ccall((:_fmpz_vec_get_nmod_poly, :libhecke), Void, (Ptr{nmod_poly}, Ptr{Int}, Int), &r, a.elem_coeffs, a.elem_length)
# this works without libhecke:
#    ccall((:nmod_poly_fit_length, :libflint), Void, (Ptr{nmod_poly}, Int), &r, a.elem_length)
#    ccall((:_fmpz_vec_get_nmod_vec, :libflint), Void, (Ptr{Void}, Ptr{Void}, Int, nmod_t), r._coeffs, a.elem_coeffs, a.elem_length, nmod_t(p, 0, 0))
#    r._length = a.elem_length
#    ccall((:_nmod_poly_normalise, :libflint), Void, (Ptr{nmod_poly}, ), &r)
  end
end

function nf_elem_to_nmod_poly_den!(r::nmod_poly, a::nf_elem)
  d = degree(a.parent)
  p = r.mod_n
  if d == 1
    ra = pointer_from_objref(a)
    den = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Void}, UInt), ra + sizeof(Int), p)
  elseif d == 2
    ra = pointer_from_objref(a)
    den = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Void}, UInt), ra + 3*sizeof(Int), p)
  else
    den = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Int}, UInt), &a.elem_den, p)
  end
  den = ccall((:n_invmod, :libflint), UInt, (UInt, UInt), den, p)
  nf_elem_to_nmod_poly_no_den!(r, a)
  mul!(r, r, den)
end

function nf_elem_to_nmod_poly(Rx::Nemo.NmodPolyRing, a::nf_elem)
  r = Rx()
  nf_elem_to_nmod_poly_den!(r, a)
  return r
end

Base.call(R::Nemo.NmodPolyRing, a::nf_elem) = nf_elem_to_nmod_poly(R, a)


text:
"""
    struct UnstructuredGridTopology{Dc,Dp,T,O} <: GridTopology{Dc,Dp}
      # private fields
    end
"""
struct UnstructuredGridTopology{Dc,Dp,T,O} <: GridTopology{Dc,Dp}
  vertex_coordinates::Vector{Point{Dp,T}}
  n_m_to_nface_to_mfaces::Matrix{Table{Int,Vector{Int},Vector{Int32}}}
  cell_type::Vector{Int8}
  polytopes::Vector{Polytope{Dc}}
end

# Constructors

"""
    UnstructuredGridTopology(
      vertex_coordinates::Vector{<:Point},
      cell_vertices::Table,
      cell_type::Vector{<:Integer},
      polytopes::Vector{<:Polytope},
      orientation::Val{O}=Val{false}()) where O
"""
function UnstructuredGridTopology(
  vertex_coordinates::Vector{<:Point},
  cell_vertices::Table,
  cell_type::Vector{<:Integer},
  polytopes::Vector{<:Polytope},
  orientation::Val{O}=Val{false}()) where O

  D = num_dims(first(polytopes))
  n = D+1
  n_m_to_nface_to_mfaces = Matrix{Table{Int,Vector{Int},Vector{Int32}}}(undef,n,n)
  n_m_to_nface_to_mfaces[D+1,0+1] = cell_vertices
  vertex_cells = generate_cells_around(cell_vertices,length(vertex_coordinates))
  n_m_to_nface_to_mfaces[0+1,D+1] = vertex_cells

  P = eltype(vertex_coordinates)
  Dp = length(P)
  T = eltype(P)

  UnstructuredGridTopology{D,Dp,T,O}(
    vertex_coordinates,
    n_m_to_nface_to_mfaces,
    cell_type,
    polytopes)

end

"""
    UnstructuredGridTopology(
      vertex_coordinates::Vector{<:Point},
      d_to_dface_vertices::Vector{<:Table},
      cell_type::Vector{<:Integer},
      polytopes::Vector{<:Polytope},
      orientation::Val{O}=Val{false}()) where O
"""
function UnstructuredGridTopology(
  vertex_coordinates::Vector{<:Point},
  d_to_dface_vertices::Vector{<:Table},
  cell_type::Vector{<:Integer},
  polytopes::Vector{<:Polytope},
  orientation::Val{O}=Val{false}()) where O

  D = num_dims(first(polytopes))
  n = D+1
  n_m_to_nface_to_mfaces = Matrix{Table{Int,Vector{Int},Vector{Int32}}}(undef,n,n)
  nvertices = length(vertex_coordinates)
  for d in 0:D
    dface_to_vertices = d_to_dface_vertices[d+1]
    n_m_to_nface_to_mfaces[d+1,0+1] = dface_to_vertices
    vertex_to_dfaces = generate_cells_around(dface_to_vertices,nvertices)
    n_m_to_nface_to_mfaces[0+1,d+1] = vertex_to_dfaces
  end

  P = eltype(vertex_coordinates)
  Dp = length(P)
  T = eltype(P)

  UnstructuredGridTopology{D,Dp,T,O}(
    vertex_coordinates,
    n_m_to_nface_to_mfaces,
    cell_type,
    polytopes)

end

"""
    UnstructuredGridTopology(topo::GridTopology)
"""
function UnstructuredGridTopology(topo::GridTopology)

  vertex_coordinates = collect1d(get_vertex_coordinates(topo))
  cell_type = collect1d(get_cell_type(topo))
  polytopes = get_polytopes(topo)
  D = num_cell_dims(topo)
  d_to_dface_vertices = [ Table(get_faces(topo,d,0)) for d in 0:D ]
  orientation = OrientationStyle(topo)

  UnstructuredGridTopology(
    vertex_coordinates,
    d_to_dface_vertices,
    cell_type,
    polytopes,
    orientation)

end

function UnstructuredGridTopology(topo::UnstructuredGridTopology)
  topo
end

function GridTopology(::Type{<:Polytope{D}},topo::GridTopology{D}) where D
  topo
end

function GridTopology(::Type{<:Polytope{D}},topo::GridTopology) where D
  _topo = UnstructuredGridTopology(topo)
  GridTopology(Polytope{D},_topo)
end

function GridTopology(::Type{<:Polytope{D}},topo::UnstructuredGridTopology) where D

  vertex_coordinates = get_vertex_coordinates(topo)
  cell_type = get_face_type(topo,D)
  polytopes = get_reffaces(Polytope{D},topo)

  d_to_dface_vertices = [ get_faces(topo,d,0) for d in 0:D ]
  orientation = OrientationStyle(topo)

  topo_d = UnstructuredGridTopology(
    vertex_coordinates,
    d_to_dface_vertices,
    cell_type,
    polytopes,
    orientation)

end

# Needed, do not remove
function num_faces(g::UnstructuredGridTopology,d::Integer)
  if d == 0
    length(g.vertex_coordinates)
  elseif d == num_cell_dims(g)
    length(g.cell_type)
  else
    D = num_cell_dims(g)
    face_to_cells = get_faces(g,d,D)
    length(face_to_cells)
  end
end

# Implementation of abstract API

OrientationStyle(::Type{UnstructuredGridTopology{Dc,Dp,T,O}}) where {Dc,Dp,T,O} = Val{O}()

get_vertex_coordinates(g::UnstructuredGridTopology) = g.vertex_coordinates

get_cell_type(g::UnstructuredGridTopology) = g.cell_type

get_polytopes(g::UnstructuredGridTopology) = collect1d(g.polytopes)

function get_faces(g::UnstructuredGridTopology, dimfrom::Integer, dimto::Integer)
  _setup_faces!(g,dimfrom,dimto)
  g.n_m_to_nface_to_mfaces[dimfrom+1,dimto+1]
end

function _setup_faces!(g,dimfrom,dimto)
  if isassigned(g.n_m_to_nface_to_mfaces,dimfrom+1,dimto+1)
    return nothing
  end
  D = num_cell_dims(g)
  if dimfrom==0
    if dimto==0
      _setup_face_to_face!(g,dimfrom)
    elseif dimto==D
      nothing
    else
      _setup_face_to_vertices!(g,dimto)
    end
  elseif dimfrom == D
    if dimto==0
      nothing
    elseif dimto==D
      _setup_face_to_face!(g,dimfrom)
    else
      _setup_cell_to_faces!(g,dimto)
    end
  else
    if dimto==0
      _setup_face_to_vertices!(g,dimfrom)
    elseif dimto==D
      _setup_cell_to_faces!(g,dimfrom)
    elseif dimto==dimfrom
      _setup_face_to_face!(g,dimfrom)
    elseif dimfrom > dimto
      _setup_nface_to_mface!(g,dimfrom,dimto)
    else
      _setup_nface_to_mface!(g,dimto,dimfrom)
    end
  end
  nothing
end

function _setup_cell_to_faces!(model,dimto)

  D = num_cell_dims(model)

  if isassigned(model.n_m_to_nface_to_mfaces,D+1,dimto+1)
    return
  end

  if isassigned(model.n_m_to_nface_to_mfaces,dimto+1,0+1)

    cell_to_vertices = get_faces(model,D,0)
    vertex_to_faces = get_faces(model,0,dimto)
    cell_to_cell_type = get_cell_type(model)
    polytopes = get_polytopes(model)
    cell_type_to_lface_to_lvertices = map( (p)->get_faces(p,dimto,0), polytopes )

    cell_to_faces = find_cell_to_faces(
      cell_to_vertices,
      cell_type_to_lface_to_lvertices,
      cell_to_cell_type,
      vertex_to_faces)

  else

    cell_to_vertices = get_faces(model,D,0)
    vertex_to_cells = get_faces(model,0,D)
    cell_to_cell_type = get_cell_type(model)
    polytopes = get_polytopes(model)
    cell_type_to_lface_to_lvertices = map( (p)->get_faces(p,dimto,0), polytopes )

    cell_to_faces = generate_cell_to_faces(
        cell_to_vertices,
        cell_type_to_lface_to_lvertices,
        cell_to_cell_type,
        vertex_to_cells)

  end

  faces_to_cells = generate_cells_around(cell_to_faces)

  model.n_m_to_nface_to_mfaces[D+1,dimto+1] = cell_to_faces
  model.n_m_to_nface_to_mfaces[dimto+1,D+1] = faces_to_cells

  nothing

end

function _setup_face_to_vertices!(model,dimfrom)

  if isassigned(model.n_m_to_nface_to_mfaces,dimfrom+1,0+1)
    return
  end

  D = num_cell_dims(model)

  cell_to_vertices = get_faces(model,D,0)
  cell_to_faces = get_faces(model,D,dimfrom)
  cell_to_ctype = get_cell_type(model)
  polytopes = get_polytopes(model)
  ctype_to_lface_to_lvertices = map( (p)->get_faces(p,dimfrom,0), polytopes )

  nfaces = num_faces(model,dimfrom)

  face_to_vertices = generate_face_to_vertices(
    cell_to_vertices,
    cell_to_faces,
    cell_to_ctype,
    ctype_to_lface_to_lvertices,
    nfaces)

  nvertices = num_faces(model,0)

  vertex_to_faces = generate_cells_around(face_to_vertices,nvertices)

  model.n_m_to_nface_to_mfaces[dimfrom+1,0+1] = face_to_vertices
  model.n_m_to_nface_to_mfaces[0+1,dimfrom+1] = vertex_to_faces

  nothing

end

function _setup_face_to_face!(model,d)

  if isassigned(model.n_m_to_nface_to_mfaces,d+1,d+1)
    return
  end

  nfaces = num_faces(model,d)
  id = identity_table(Int,Int32,nfaces)
  model.n_m_to_nface_to_mfaces[d+1,d+1] = id

  return

end

function _setup_nface_to_mface!(model,n,m)
  @assert n > m

  if isassigned(model.n_m_to_nface_to_mfaces,n+1,m+1)
    return
  end

  nface_to_vertices = get_faces(model,n,0)
  vertex_to_mfaces = get_faces(model,0,m)
  nface_to_nftype = get_face_type(model,n)

  polytopes = get_reffaces(Polytope{n},model)
  nftype_to_lmface_to_lvertices = map( (p)->get_faces(p,m,0), polytopes)

  nface_to_mfaces = find_cell_to_faces(
    nface_to_vertices,
    nftype_to_lmface_to_lvertices,
    nface_to_nftype,
    vertex_to_mfaces)

  num_mfaces = num_faces(model,m)
  mface_to_nfaces = generate_cells_around(nface_to_mfaces,num_mfaces)

  model.n_m_to_nface_to_mfaces[n+1,m+1] = nface_to_mfaces
  model.n_m_to_nface_to_mfaces[m+1,n+1] = mface_to_nfaces

  nothing

end


text:
#set OMP_NUM_THREADS=<number of threads to use>.
using MAT
using jInv.Mesh
using DivSigGrad
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.Utils
using ArgParse
using KrylovMethods
using Multigrid

include("getDCResistivitySourcesAMG.jl")

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin

	"--n1"
		help = "number of cells in x1 direction "
		arg_type = Int
		default = 48
	"--n2"
		help = "number of cells in x2 direction "
		arg_type = Int
		default = 48
	"--n3"
		help = "number of cells in x3 direction "
		arg_type = Int
		default = 24
	"--srcSpacing"
		help = "spacing between sources"
		arg_type = Int
		default = 2
	"--nthreads"
		help = "number of threads"
		arg_type = Int
		default = 1
	"--out"
		help ="file to write into"
		default = "times.csv"
	"--solver"
		help = "PDE solver: 1->PCG, 2->BlockPCG, 3->MUMPS"
		arg_type = Int
		default = 1
	end
	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()

	matfile   = matread("../3Dseg12812864.mat")
	m         = matfile["VELc"]
	nInv      = [128;128;64]

	n1 = parsed_args["n1"]
	n2 = parsed_args["n2"]
	n3 = parsed_args["n3"]
	ns = parsed_args["srcSpacing"]
	nthreads = parsed_args["nthreads"]

	out = parsed_args["out"]

	# set number of threads for openblac
	set_num_threads(nthreads)

	solver = parsed_args["solver"]

	domain  = [0;13.5;0;13.5;0;4.2]
	n       = [n1;n2;n3]

	M    = getRegularMesh(domain,n)
	MInv = getRegularMesh(domain,nInv)
	Q,P  = getDCResistivitySourcesAndRecAll(M,srcSpacing=[ns,ns]);
	nsrc = size(Q,2)

	@printf "problem size: n1=%d\tn2=%d\tn3=%d\tns=%d\n" n1 n2 n3 nsrc
	@printf "nworkers()=%d\n" nworkers()
	@printf "nthreads=%d \n" nthreads

	if solver==3
		@printf "solver=MUMPS\n\n"
		@printf "warmup\ttotal\tfactorization\tsolve\n"
	elseif solver==2
		@printf "solver=BlockPCG\n\n"
		@printf "warmup\ttotal\titer\tprecond\tcgTime\tmatvec\n"
	elseif solver==1
		@printf "solver=PCG\n\n"
		@printf "warmup\ttotal\titer\tprecond\tcgTime\tmatvec\n"
	end
	Ainv = []
	if solver==1
		Ainv   = getIterativeSolver(KrylovMethods.cg,out=-1,sym=1,PC=:jac);
		Ainv.tol=1e-8
		Ainv.nthreads=nthreads
	elseif solver==2
 		Ainv         = getBlockIterativeSolver(KrylovMethods.blockCG,out=-1,sym=1,PC=:jac);
		Ainv.maxIter = 1000
		Ainv.tol     = 1e-8
		Ainv.nthreads=nthreads
		Ainv.out     =-1
		Ainvt         = getBlockIterativeSolver(KrylovMethods.blockCG,out=-1,sym=1,PC=:jac);
		Ainvt.maxIter = 1000
		Ainvt.tol     = 1e-8
		Ainvt.nthreads=nthreads
		Ainvt.out     =-1

	elseif solver==3
		Ainv = getMUMPSsolver()
		Ainvt = getMUMPSsolver()
	elseif solver 		== 4
		levels      	= 4;
		numCores 	= nthreads;
		maxIter     	= 20;
		relativeTol 	= 1e-8;
		relaxType   	= "SPAI";
		relaxParam  	= 1.1;
		relaxPre 	= 2;
		relaxPost   	= 2;
		cycleType   	='V';
		coarseSolveType = "MUMPS";
		strongConnParam = 0.25;
		MG 				= getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,strongConnParam);
		Ainvt   		= getSA_AMGsolver(MG, "PCG",sym=1,out=1);
		Ainv   			= getSA_AMGsolver(MG, "PCG",sym=1,out=1);
	end

	pForp = Array{RemoteRef{Channel{Any}}}(nworkers())
	pFors = Array{RemoteRef{Channel{Any}}}(nworkers())
	Mesh2Mesh  = getInterpolationMatrix(MInv,M)'
	m          = Mesh2Mesh'*vec(m)
	glocp      = Array{RemoteRef{Channel{Any}}}(nworkers())

	for j=1:nworkers()
	println("Worker number:",j)
		# inds = unique(min(size(Q,2),(j-1)*ppw+(1:ppw)))
		pForp[j] = @spawnat workers()[j] DivSigGradParam(M,full(Q[:,1:nworkers():end]),P,[1.0],copySolver(Ainv));
		pFors[j] = @spawnat workers()[j] DivSigGradParam(M,full(Q[:,1:2]),P,[1.0],copySolver(Ainvt));
		glocp[j] = @spawnat workers()[j] sparse(I,M.nc,M.nc)
	end

	# warm up
 	wTime = @elapsed begin
	res = getData(vec(m),pFors,glocp)
	end
	@printf "%3.2f\t" wTime
	clear!(res[1])
	clear!(res[2])
	clear!(pFors)

	tTime = @elapsed begin
		getData(vec(m),pForp,glocp)
	end
	@printf "%3.2f\t" tTime

	if solver==3
		fTime = 0.0; sTime =0.0
		for p=1:length(pForp)
			pf = take!(pForp[p])
			fTime += pf.Ainv.facTime
			sTime += pf.Ainv.solveTime
		end
		fTime /= nworkers()
		sTime /= nworkers()
		@printf "%3.4f\t\t%3.4f\n" fTime sTime

		f = open(out,"a")
		str = @sprintf "%d,%d,%d,%d,%d,%d,,%3.5f,%3.5f,%3.5f,%3.5f\n" nworkers() solver n1 n2 n3 nsrc wTime tTime fTime sTime
		write(f,str )
		close(f)
	elseif solver==1 || solver==2
		nIter = 0; pcTime = 0.0; cgTime=0.0; mvTime=0.0
		for p=1:length(pForp)
			pf = take!(pForp[p])
			nIter += pf.Ainv.nIter
			pcTime += pf.Ainv.timePC
			cgTime += pf.Ainv.timeSolve
			mvTime += pf.Ainv.timeMV
			clear!(pf)
		end
		nIter /= nworkers()
		pcTime /=nworkers()
		cgTime /=nworkers()
		mvTime /=nworkers()
		@printf "%5d\t%3.4f\t%3.4f\t%3.4f\n" nIter  pcTime cgTime mvTime

		f = open(out,"a")
		str = @sprintf "%d,%d,%d,%d,%d,%d,%3.5f,%3.5f,%3.5f,%3.5f,%3.5f,%d\n" nworkers() solver n1 n2 n3 nsrc wTime tTime pcTime cgTime mvTime nIter
		write(f,str )
		close(f)
	elseif solver == 4
		nIter = 0; tsetup = 0.0; tSol = 0.0
		for p=1:length(pForp)
			pf = take!(pForp[p])
			nIter  += pf.Ainv.nIter
			tsetup += pf.Ainv.timeSetup
			tSol   += pf.Ainv.timeSolve
			clear!(pf)
		end
		#nIter /= nworkers()
		#tsetup /=nworkers()
		#tSol   /=nworkers()
		@printf "%5d\t%3.4f\t%3.4f\t\n" nIter tsetup tSol

	end

end
main()


text:
module CommunalHelperCassetteZipMover

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/CassetteZipMover" CassetteZipMover(x::Integer,
                                                                  y::Integer,
                                                                  width::Integer=Maple.defaultBlockWidth,
                                                                  height::Integer=Maple.defaultBlockHeight,
                                                                  index::Integer=0,
                                                                  tempo::Number=1.0,
                                                                  noReturn::Bool=false)

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255
)

const ropeColors = Dict{Int, Ahorn.colorTupleType}(
    1 => (194, 116, 171, 255) ./ 255,
	2 => (227, 214, 148, 255) ./ 255,
	3 => (128, 224, 141, 255) ./ 255
)

const defaultRopeColor = (110, 189, 245, 255) ./ 255
const defaultColor = (73, 170, 240, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Zip Mover ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteZipMover,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
        ),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::CassetteZipMover) = Dict{String, Any}(
    "index" => colorNames
)

Ahorn.nodeLimits(entity::CassetteZipMover) = 1, 1

Ahorn.minimumSize(entity::CassetteZipMover) = 16, 16
Ahorn.resizable(entity::CassetteZipMover) = true, true

function Ahorn.selection(entity::CassetteZipMover)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10)]
end

function getTextures(entity::CassetteZipMover)
    return "objects/cassetteblock/solid", "objects/CommunalHelper/cassetteZipMover/cog"
end
const crossSprite = "objects/CommunalHelper/cassetteMoveBlock/x"

function renderCassetteZipMover(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    block, cog = getTextures(entity)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)
    ropeColor = get(ropeColors, index, defaultRopeColor)

    # Node Rendering
    cx, cy = x + width / 2, y + height / 2
    cnx, cny = nx + width / 2, ny + height / 2
    length = sqrt((x - nx)^2 + (y - ny)^2)
    theta = atan(cny - cy, cnx - cx)

    Ahorn.Cairo.save(ctx)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)
    Ahorn.translate(ctx, cx, cy)
    Ahorn.rotate(ctx, theta)
    Ahorn.setSourceColor(ctx, ropeColor)
    # Offset for rounding errors
    Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
    Ahorn.line_to(ctx, length, 4 + (theta <= 0))
    Ahorn.move_to(ctx, 0, -4 - (theta > 0))
    Ahorn.line_to(ctx, length, -4 - (theta > 0))
    Ahorn.stroke(ctx)
    Ahorn.Cairo.restore(ctx)
    Ahorn.drawSprite(ctx, cog, cnx, cny, tint=color)

    for i in 1:tilesWidth, j in 1:tilesHeight
        tx = (i == 1) ? 0 : ((i == tilesWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tilesHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y + (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end

    if Bool(get(entity.data, "noReturn", false))
        noReturnSprite = Ahorn.getSprite(crossSprite, "Gameplay")
        Ahorn.drawImage(ctx, noReturnSprite, x + div(width - noReturnSprite.width, 2), y + div(height - noReturnSprite.height, 2), tint=color)
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover, room::Maple.Room)
    renderCassetteZipMover(ctx, entity)
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
function convolve{T, n}(A::SubArray{T, n},B::Array{T, n},dims)
    common_size=tuple(map(max,size(A),size(B))...)
    pA=zeros(Complex{T},common_size)
    pB=zeros(Complex{T},common_size)
    rangesA=[1:x for x in size(A)]
    rangesB=[1:x for x in size(B)]
    pA[rangesA...]=A
    pB[rangesB...]=B

    real(ifft!(fft!(pA,dims).*fft!(pB,dims),dims))
end

function convolve_ComplexFloat64(A,B,dims)
    common_size=tuple(map(max,size(A),size(B))...)
    pA=zeros(Complex{Float64},common_size)
    pB=zeros(Complex{Float64},common_size)
    rangesA=[1:x for x in size(A)]
    rangesB=[1:x for x in size(B)]
    pA[rangesA...]=A
    pB[rangesB...]=B

    real(ifft!(fft!(pA,dims).*fft!(pB,dims),dims))
end

# got rid of dims argument
# using real to complex fft and ifft
function convolve_Float64(A,B)
    common_size=tuple(map(max,size(A),size(B))...)
    pA=zeros(Float64,common_size)
    pB=zeros(Float64,common_size)
    rangesA=[1:x for x in size(A)]
    rangesB=[1:x for x in size(B)]
    pA[rangesA...]=A
    pB[rangesB...]=B

    irfft(rfft(pA).*rfft(pB),common_size[1])
end

function valid_convolve(A,B)
    ranges=[min(a,b):max(a,b) for (a,b) in zip(size(A),size(B))]
    convolve_Float64(A,B)[ranges...]
end

function cumsum2{T,ndim}(A::Array{T,ndim})
    # cumulative sum in two dimensions
    if ndim!=2
        throw(ArgumentError("input must be two-dimensional array"))
    end
    B = similar(A);
    # first row and column are 1D cumulative sums
    B[:,1]=cumsum(A[:,1],1);
    B[1,:]=cumsum(A[1,:],2);      # B[1,1] is redundantly computed twice
    # compute rest of matrix from recursion
    (m,n)=size(A);
    if m>1 && n>1
        for j in 2:n, i in 2:m
            B[i,j]=A[i,j]+B[i-1,j]+B[i,j-1]-B[i-1,j-1]
        end
    end
    B
end

function optimize_normxcorr2(img)
    p=plan_rfft(img,flags=FFTW.MEASURE)
    q=plan_irfft(rfft(img),flags=FFTW.MEASURE,size(img,1))
    return Void;
end

# extension:
# Params_session.jl: optimize_all_cores(params::Params)
function optimize_all_cores(img_d::Int64)
    img = rand(img_d, img_d)
	@sync begin
	  for p in 1:num_procs
	    @async begin
		remotecall_fetch(p, optimize_normxcorr2, img);
		end
	    end
	  end
    return Void;
end

"""
`NORMXCORR2` - normalized cross correlation
slide template across img, compute Pearson correlation coefficient for each
template location result has same size as MATLAB-style 'valid' convolution
efficient algorithm of J. P. Lewis
http://scribblethink.org/Work/nvisionInterface/nip.html

` correlogram = normxcorr2(template,img)`

return NaN for image patch or template with zero variance
Pearson correlation coefficient is undefined in this case

need some argument checking
e.g. sizes of template should be less than those of img
this works for arrays.  extend to Image defined in Holy's package?
"""
function normxcorr2(template,img)
    # sufficient to subtract mean from just one variable
    dt=template-mean(template)
    templatevariance=sum(dt.^2)
    numerator=valid_convolve(img,dt[end:-1:1,end:-1:1])

    ##### local statistics of img
    # zero pad image in first row and column
    # so that cumulative sums will have zeros in the same place
    (m1,m2)=size(img);
    imgpad=zeros(m1+1,m2+1); imgpad[2:end,2:end]=img;
    # define four combinations of Small and Large ranges
    (n1,n2)=size(template);
    if templatevariance==0
        return zeros(m1-n1+1,m2-n2+1)*NaN
    end
    LL=Any[1+(n1:m1),1+(n2:m2)]
    SL=LL-[n1;0]; LS=LL-[0;n2]
    SS=LL-[n1;n2]
    # sum of img and its square in template-sized neighborhoods
    s=cumsum2(imgpad)
    localsum=s[LL...]-s[SL...]-s[LS...]+s[SS...]
    s2=cumsum2(imgpad.^2)
    localsum2=s2[LL...]-s2[SL...]-s2[LS...]+s2[SS...]
    localvariance=localsum2-localsum.^2/prod(size(template))
    # localvariance is zero for image patches that are constant
    # leading to undefined Pearson correlation coefficient
    # should only be negative due to roundoff error
    localvariance[localvariance.<=0] *= NaN
    denominator=sqrt(localvariance*templatevariance)

    numerator./denominator
end


text:
@enum Suit ♣ ♦ ♥ ♠
@enum Palo Bastos Espadas Oros Copas

toString( n::Int, s::Suit ) = string( n, " ", s )
toString( n::Int, p::Palo ) = string( n, " de ", p )

println( toString( 3, ♦) )
println( toString( 7, Bastos ) )


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using DifferentialEquations
using LinearAlgebra

# Set the required model parameters for the SIRS model with two levels of transmission - Within and between households
N = 10; # Household size - Change to 10 for final analysis
betaHH = 6; # Within household transmission parameter
betaG = 1; # Population wide transmission
gamma = 1; # Rate of recovery from infection
tau = 1; # Rate of loss of protection
params = [betaHH,gamma,tau,betaG,N]; # Put all the parameters together
time = (0.0, 30.0) # Simulation time - note it defined as a float
dim = dim = convert(Int,0.5*(N+1)*(N+2)); # Number of possible configurations - works for three epidemiological classes
y0 = vec(zeros(1, dim)); # Initial condition vector
y0[end-1] = 0.00000001;
y0[end] = 0.99999999;

function hhTransitions(N,dim)
    # Function to generate transition matrices for household model
    # Input: N is the household size

    # Initialize things
    Qinf = zeros(dim,dim);
    Qrec = zeros(dim,dim);
    Qext = zeros(dim,dim);
    Qwane = zeros(dim,dim);
    dataI = Array{Int64}(zeros(dim,3))
    m = 0;
    I = Array{Int64}(zeros(N+1,N+1))

    # To help remember where to store the variables
    for ss = 0:N
        for ii = 0:(N-ss)
            m = m + 1;
            I[ss+1,ii+1] = m
        end
    end

    # Describe the epidemiological transitions

    # Counter for susceptibles
    for ss = 0:N
        # Counter for infecteds
        for ii = 0:(N-ss)
            # If susceptibles and infecteds are more than 1, then infection within the household can occur
            if (ss > 0 && ii > 0)
                Qinf[I[ss+1,ii+1],I[ss,ii+2]] = ii*ss/(N-1);
            end

            # If infecteds are more than 1, recovery can occur
            if ii > 0
                # Rate of recovery
                Qrec[I[ss+1,ii+1],I[ss+1,ii]] = ii;
            end

            # For external infection - just keep track of susceptibles
            if ss > 0
                # Rate of within household infection
                Qext[I[ss+1,ii+1],I[ss,ii+2]] = ss;
            end

            # Loss of protection hence becoming susceptible again. Possible if N-ss-ii = rr > 0
            if (N-ss-ii) > 0
                # Rate of loss of protection
                Qwane[I[ss+1,ii+1],I[ss+2,ii+1]] = N-ss-ii;
            end

            # Store the relevant indices to help identify the household configurations
            dataI[I[ss+1,ii+1],:] = [ss, ii, N-ss-ii];
        end
    end

    Qinf = Qinf - diagm(0=>vec(sum(Qinf,dims=2)));
    Qrec = Qrec - diagm(0=>vec(sum(Qrec,dims=2)));
    Qext = Qext - diagm(0=>vec(sum(Qext,dims=2)));
    Qwane = Qwane - diagm(0=>vec(sum(Qwane,dims=2)));

    # Return
    return Qinf, Qrec, Qext, Qwane, dataI
end
Qinf, Qrec, Qext, Qwane, dataI = hhTransitions(N,dim);

# +
function rateSIRS(dy_dt,y0,params,time)

    # Extract the parameters
    betaHH = params[1];
    gamma = params[2];
    tau = params[3];
    betaG = params[4];
    N = params[5];

    # Generate the transition matrices
    Qinf, Qrec, Qext, Qwane, HHconfig = hhTransitions(N,dim);

    # Combine within and external transitions
    Q = betaHH*Qinf + gamma*Qrec + tau*Qwane + betaG*Qext*((HHconfig[:,2]'*y0)/N);

    # Generate the rates of change for ODE solver
    dy_dt .= (y0'*Q)';

end

# Define the ODE problem
prob = ODEProblem(rateSIRS,y0,time,params);

# Solve
sol = solve(prob, alg=Tsit5());
# -

using Plots
#using RecursiveArrayTools

# +
Iconfig = dataI[:,2];
infProp = zeros(length(sol.t),1);

# Prepare the plots
for i = 1:length(sol.t)
    infProp[i,1] = sol[:,i]'*Iconfig/N;
end

# -

# Total infectious in the population
plot(sol.t, infProp,
    color="blue",xlabel="Time",
    ylabel="Proportion infectious",
    label=["Mean infection"],
    ylims=[0, 1]);

# Household profile at endemic prevalence
# Prepare the plots
using Plots
yprop = zeros(N+1,length(sol.t))
for j = 1:length(sol.t)
    for i = 1:N+1
        index = findall(Iconfig.==i-1);
        yprop[i,j] = sum(sol[index,j])
    end
end
plot(-0.5:1:10.5,[yprop[:,length(sol.t)];0]);

# Household profile at peak prevalence
# Prepare the plots
x = findall(infProp.==maximum(infProp))
plot(-0.5:1:10.5,[yprop[:,x];0]);


text:
# Natural Round
round2(x; digits = 0) = sign(x) * floor( abs(x) * 10.0^digits + 0.5 ) / (10.0^digits)

# Cumulative distribution function
function CDF(x, F::EG)
    if x < minimum(F.table.scale) return 0 end
    if x > maximum(F.table.scale) return 1 end
    F.table.cumprob[F.table.scale .== x][1]
end

# Percentile rank function
function PRF(x, F::EG)
    if x < (minimum(F.table.scale) - F.interval/2.0) return 0.0 end
    if x ≥ (maximum(F.table.scale) + F.interval/2.0) return 100.0 end
    x′ = round2(x)
    Fx1 = CDF(x′-F.interval, F)
    Fx = CDF(x′, F)
    P = 100*(Fx1+(x-x′+F.interval/2.0)*(Fx-Fx1))[1]
    return P
end

# Percentile Function
function p_search_descend(P, F::EG, offset)
    x = nothing;iter = length(F.table.scale)
    while x === nothing
        iter -= 1
        x =  100CDF(F.table.scale[iter], F) > P ? nothing : F.table.scale[iter+offset]
    end
    return x
end
function p_search_ascend(P, F::EG, offset)
    x = nothing;iter = 0
    while x === nothing
        iter += 1
        x = 100CDF(F.table.scale[iter], F) < P ? nothing : iter == 1 ? 0.0 : F.table.scale[iter+offset]
    end
    return x
end

# The inverse of the percentile rank function (Percentile function)
function PFu(P, F::EG)
    if P ≥ 100.0 return (maximum(F.table.scale) + .5) end
    xu = P > 50.0 ? p_search_descend(P, F, 1) : p_search_ascend(P, F, 0)
    x = (P/100 - CDF(xu-F.interval, F)) / (CDF(xu, F) - CDF(xu-F.interval, F))
    return isinf(x) || isnan(x) ? xu -F.interval/2.0 : x + xu -F.interval/2.0
end
function PFl(P, F::EG)
    if P ≤ 0.0 return -.5 end
    xl = P > 50.0 ? p_search_descend(P, F, 0) : p_search_ascend(P, F, -1)
    x = (P/100 - CDF(xl, F)) / (CDF(xl+F.interval, F) - CDF(xl, F))
    return isinf(x) || isnan(x) ? xl + F.interval/2.0 : x + xl + F.interval/2.0
end

# equipercentile equating
mutable struct SGEquateResult <: SGEquateMethod
    method
    table
    estimates
    data
end
"""
    Equipercentile(X::EG, Y::EG; case = :middle)

Equipercentile Equating under equivalent (single) group design. The smoothed frequency can be used.
A table returned by the function contains 2 columns. The first one stands for the base scale score in the test X.
The second one is equated score from test Y, which is equipercentile score on the base scale.

`case` represents which equating case use. Pass to the symbols below.

- `:upper` Calculating scores correspond to arbitrary percentile P, use the smallest integer score with a cumulative percent that is greater than P.
- `:lower` (default) Contrary to above case, use the largest integer score with a cumulative percent that is less than P
- `:both` (Not for practice) Show both case of upper and lower.
- `:middle` Use midpoint case upper between lower.

In my example cases, it seems that equate package in R uses `lower` case to equate.
"""
function Equipercentile(X::EG, Y::EG; case = :lower)
    scaleX = X.table.scale
    eYxu = zeros(Float64, length(scaleX))
    eYxl = zeros(Float64, length(scaleX))
    for (i,v) in enumerate(scaleX)
        P = PRF(v, X)
        eYxu[i] = PFu(P, Y)
        eYxl[i] = PFl(P, Y)
    end
    if case == :upper
        eYx = eYxu
    elseif case == :lower
        eYx = eYxl
    elseif case == :both
        eYx = string.(eYxu, "_", eYxl)
    elseif case == :middle
        eYx = (eYxu .+ eYxl) ./ 2.0
    end
    tbl = DataFrame(scaleX = scaleX, eYx = eYx)
    return SGEquateResult(
        Symbol("Equipercentile($(case))"),
        tbl,
        nothing,
        (X = X, Y = Y)
    )
end
# linear equating
"""
    Linear(X::EG, Y::EG)

Linear equating under the equivalent group desing.
This method, equate to match first 2 moments, is so simple to comprehend the equating result.
"""
function Linear(X::EG, Y::EG)
    μX = mean(X.raw); σX = std(X.raw)
    μY = mean(Y.raw); σY = std(Y.raw)
    slope = σY/σX; intercept = μY - slope*μX
    lYx = @. X.table.scale * slope + intercept
    tbl = DataFrame(scaleX = X.table.scale, lYx = lYx)
    return SGEquateResult(
        :Linear,
        tbl,
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end

# mean equating
"""
    Mean(X::EG, Y::EG)

Mean equating under the equivalent group desing.
This method, equate to match only first moments = mean, is so simple to comprehend the equating result.
"""
function Mean(X::EG, Y::EG)
    μX = mean(X.raw)
    μY = mean(Y.raw)
    slope = 1.0; intercept = μY - μX
    lYx = @. X.table.scale * slope + intercept
    tbl = DataFrame(scaleX = X.table.scale, lYx = lYx)
    return SGEquateResult(
        :Mean,
        tbl,
        (slope = slope, intercept = intercept),
        (X = X, Y = Y)
    )
end


text:
# testEmadoka.jl
using Test

include("Facto.jl")
include("playEmadoka.jl")

@testset "Madoka test" begin
 @testset "basic test" begin
  @test solveQ([who],[[who is mahosyojo]],MagicalWorld) == [[:Madoka],[:Homura],[:Kyoko],[:Mami],[:Sayaka]]
  @test solveQ([who],[[who have soulgem]],MagicalWorld) == [[:Homura],[:Sayaka],[:Kyoko],[:Mami]]
  @test solveQ([who],[[who have soulgem],[who love Kamijo]],MagicalWorld) == [[:Sayaka]]
  @test solveQ([who,what],[[what become witch]],MagicalWorld) == [[:who,:mahosyojo]]
 end

 ##
 @testset "reflection rule" begin
  @test solveQ([:who1,:who2],[[:who1 love :who2],[:who2 love :who1]], MagicalWorld)==[[:Hitomi,:Kamijo],[:Kamijo,:Hitomi]]

 # [] means no model valid
  @test solveQ([:who1,:who2],[[:who1 love :who2],[:who2 love :who1]], MagicalWorld2)==[]

 end
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function discretizeliver!(system::CVSystem)
    system.liver.P = zeros(system.solverparams.numsteps+1,4);
    system.liver.V = zeros(system.solverparams.numsteps+1,4);
    system.liver.Q = zeros(system.solverparams.numsteps+1,4);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
@testset "arrays" begin
    @test arraytype(ones(3, 5)) === Array
    if CUDA.has_cuda_gpu()
        @test arraytype(CuArray(ones(3, 5))) === CuArray
    end

    @test [1, 3, 1, 2, 1, 2] == Bennu.numbercontiguous([10, 29, 10, 23, 10, 23])
    TAs = ((Float64,  Array), (BigFloat, Array))
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end
    for (T, A) in TAs
        a = fieldarray(undef, T, A, ())
        @test size(a) == ()

        a = fieldarray(undef, T, A, (3,))
        @test size(a) == (3,)
        @test eltype(a) == T
        @test Tullio.storage_type(a) <: A

        a = fieldarray(undef, T, A, (3,4))
        @test size(a) == (3,4)
        @test eltype(a) == T
        @test Tullio.storage_type(a) <: A

        a = fieldarray(undef, (a=SVector{2, T}, b=T), A, (3,4))
        @test size(a) == (3,4)
        @test eltype(a) == NamedTuple{(:a, :b), Tuple{SVector{2, T}, T}}
        @test Tullio.storage_type(a) <: A
    end
end


text:
@testset "DeprecatedTest" begin
    include("alldifferent.jl")
    include("antidomain.jl")
    include("binpacking.jl")
    include("count.jl")
    include("countdistinct.jl")
    include("differentfrom.jl")
    include("domain.jl")
    include("element.jl")
    include("equivalence.jl")
    include("equivalencenot.jl")
    include("ifthenelse.jl")
    include("implication.jl")
    include("indicator.jl")
    include("inverse.jl")
    include("lexicographically.jl")
    include("minimumdistance.jl")
    include("reification.jl")
    include("strictly.jl")
    include("strictly_lexicographically.jl")
    include("truefalse.jl")
end


text:
using NeuralQuantum, QuantumOpticsBase, Statistics
using Test

N = 4
g = 0.7
V = 2.0

hilb = HomogeneousSpin(N)

H = LocalOperator(hilb)
for i=1:N
    global H += g/2.0 * sigmax(hilb, i)
    global H  += V/4.0 * sigmaz(hilb, i) * sigmaz(hilb, mod(i, N)+1)
end

Sx = LocalOperator(hilb)
Sy = LocalOperator(hilb)
Sz = LocalOperator(hilb)
for i=1:N
    global Sx += sigmax(hilb, i)/N
    global Sy += sigmay(hilb, i)/N
    global Sz += sigmaz(hilb, i)/N
end

prob  = NeuralQuantum.HamiltonianGSEnergyProblem(H)

net  = RBM( N, 2)

sampl = MetropolisSampler(LocalRule(), 300, N, burn=30)
#sampl = ExactSampler(1000)
algo  = SR(ϵ=(0.001), use_iterative=true)
algo  = Gradient()

is = BatchedSampler(net, sampl, prob, algo; batch_sz=32)
res, grad = sample!(is)

# exact value
Hmat = Matrix(H)
psi  = ket(net, hilb)

E_ex = psi'*Hmat*psi / (psi'psi)

# precise sampling
sampl_ex = ExactSampler(1000)
algo     = SR()
is = BatchedSampler(net, sampl_ex, prob, algo; batch_sz=32)

res, grad = sample!(is)

@test abs(res.mean - E_ex) < res.error

#


text:
"""
    read_dot_file(file)

Import graph from DOT-Format and store it in file `SimpleWeightedGraph` or `SimpleWeightedDiGraph``.
Return AttributeDict for graph layouts.
ToDo: ERROR-Handling if not a suitable DOT-File is not implemented

#### Arguments
- `file::AbstractString`: the filename of dot-file (i.e. "graph.dot")
"""
function read_dot_file(filename::AbstractString)

    # read_file
    #my_graph = read_graph(filename)

    # open, read & preprocessing
    my_graph = preprocessing(filename)

    # use ParserCombinator.Parsers.DOT to convert to attributes
    graphs_parser_combinator = ParserCombinator.Parsers.DOT.parse_dot(my_graph)

    # use first graph (toDo for more graphs)  !!! PROBLEM
    if length(graphs_parser_combinator) == 1
        graph_parser_combinator = graphs_parser_combinator[1]
    else
        graph_parser_combinator = graphs_parser_combinator[1]
    end

    attrs = init_Attributes!(graph_parser_combinator)

    # set! attributes
    set_attributes!(attrs, graph_parser_combinator)

    # get AbstractSimpleWeightedGraph
    g = get_AbstractSimpleWeightedGraph(attrs)

    n = nv(g)
    large_graph = 200
    if n > large_graph
        attrs = mod_attr_large_network!(attrs)
    end

    return g, attrs
end

function openfile(filename)
    f = open(filename, "r")
    lines = readlines(f)
    close(f)
    lines
end

function preprocessing(filename)
    lines = openfile(filename)
    new_lines = []
    temp = ""
    multi_line_option = false
    subgraph = false
    for line in lines

        # some small corrections:
        line = replace(line, "\t" => "")
        line = replace(line, " ]" => "]")
        line = replace(line, "[ " => "[")

        # delete comments and empty lines:
        if !isnothing(findfirst("//", lstrip(line))) || !isnothing(findfirst("/*", lstrip(line))) || !isnothing(findfirst("*", lstrip(line))) || isempty(line)
            continue
        end

        if !(multi_line_option)

            # multi subgraph lines:
            if !isnothing(findfirst("{", lstrip(line)))
                idx = findfirst("{", lstrip(line))
                if idx.stop < 2
                    subgraph = true
                    multi_line_option = true
                end
            end

            # in case of edge line:
            # if !(isnothing(findfirst("--", line))) || !(isnothing(findfirst("->", line)))
            #     push!(new_lines, line)
            #     continue
            # end

            optionend = findfirst("]", line)

            # if multi-line syntax:
            if !isnothing(findfirst("[", line)) && isnothing(optionend)
                temp = temp * line
                multi_line_option = true
                continue
            end

            # one-line sub-graph "node"
            if !isnothing(optionend)

                if isnothing(findfirst("{", line[1:optionend.start]))

                    if (isnothing(findfirst("--", line))) && (isnothing(findfirst("->", line)))

                        if !(isempty(lstrip(line[optionend.stop+1:end])))

                            if isnothing(findfirst("{", line[optionend.stop:end]))

                                if !(lstrip(line[(optionend.stop+1):end]) == ";") || (isempty(lstrip(line[(optionend.stop+1):end])))
                                    str = "subgraph {\n " * line * " \n}"
                                    push!(new_lines, str)
                                    continue
                                end
                            else
                                idx = findfirst("{", line[optionend.stop:end])
                                str = insert_at!(line, "\n", (optionend.stop + idx.start - 1))
                                push!(new_lines, str)
                                continue
                            end
                        end
                    end
                end
            end

            # add line to new_lines
            push!(new_lines, line)
        else

            if (subgraph == false)
                temp = temp * " " * line
                if !isnothing(findfirst("]", line))
                    temp = replace(temp, "\t" => "")
                    temp = replace(temp, " ]" => "]")
                    temp = replace(temp, "[ " => "[")
                    push!(new_lines, temp)
                    multi_line_option = false
                    temp = ""
                end
            else
                push!(new_lines, line)
                if !isnothing(findfirst("}", lstrip(line)))
                    subgraph = false
                    multi_line_option = false
                end
            end
        end
    end

    # #DEBUG:
    # for line in new_lines
    #     @show line
    # end

    return Base.join(new_lines, "\n")
end

function read_graph(filename)
    f = open(filename, "r")
    s = read(f, String)
    close(f)
    return s
end

function init_Attributes!(g)
    if isdefined(g, :directed)
        node_options = [PlotGraphviz.Property("directed", g.directed)]
    end

    # nodes:
    nd = ParserCombinator.Parsers.DOT.nodes(g)
    nodes = []
    count = 1
    for key in nd.dict
        push!(nodes, PlotGraphviz.gvNode(count, String(key[1]),
            [Property("label", check_value(String(key[1])))]))
        count += 1
    end

    # edges:
    eg = ParserCombinator.Parsers.DOT.edges(g)
    edges = []
    for tup in eg.dict
        to = get_id(nodes, String(tup[1][1]))
        from = get_id(nodes, String(tup[1][2]))
        push!(edges, PlotGraphviz.gvEdge(to, from, Properties()))
        if g.directed == false
            push!(edges, PlotGraphviz.gvEdge(from, to, Properties()))
        end
    end

    # subgraphs
    subgraphs = (gvSubGraph)[]

    return GraphvizAttributes(node_options, Properties(), Properties(), Properties(), subgraphs, nodes, edges)

end

function set_attributes!(attrs, g)

    for stm in g.stmts

        if stm isa ParserCombinator.Parsers.DOT.Attribute
            set!(attrs.graph_options, String(stm.name.id), check_value(String(stm.value.id)))
        elseif stm isa ParserCombinator.Parsers.DOT.GraphAttributes
            for attr in stm.attrs
                set!(attrs.graph_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.NodeAttributes
            for attr in stm.attrs
                set!(attrs.node_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.EdgeAttributes
            for attr in stm.attrs
                set!(attrs.edge_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.Edge

            if !(isempty(stm.attrs))

                for attr in stm.attrs
                    set!(attrs.edges, get_id(attrs.nodes, String(stm.nodes[1].id.id)), get_id(attrs.nodes, String(stm.nodes[2].id.id)),
                        Property(String(attr.name.id), check_value(String(attr.value.id))); override=true)

                    if (g.directed == false)
                        set!(attrs.edges, get_id(attrs.nodes, String(stm.nodes[2].id.id)), get_id(attrs.nodes, String(stm.nodes[1].id.id)),
                            Property(String(attr.name.id), check_value(String(attr.value.id))); override=true)
                    end

                end
            end
        elseif stm isa ParserCombinator.Parsers.DOT.Node
            if !(isempty(stm.attrs))
                for attr in stm.attrs
                    set!(attrs.nodes, get_id(attrs.nodes, String(stm.id.id.id)),
                        Property(String(attr.name.id), check_value(String(attr.value.id))))
                end
            end
        elseif stm isa ParserCombinator.Parsers.DOT.SubGraph
            if !isnothing(stm.id) # use only cluster as subgraphs.
                push!(attrs.subgraphs, gvSubGraph(String(stm.id.id)))
                set_subgraph!(attrs.subgraphs[end], stm, attrs.nodes, g.directed)
            else
                # all other subgraphs are lazy syntax
                attributes = Properties()
                for attr_sub in stm.stmts
                    if attr_sub isa ParserCombinator.Parsers.DOT.NodeAttributes
                        for attr in attr_sub.attrs
                            set!(attributes, String(attr.name.id), check_value(String(attr.value.id)))
                        end
                    end
                end

                # for all node in subgraph
                for attr_sub in stm.stmts
                    if attr_sub isa ParserCombinator.Parsers.DOT.Node
                        for single_attr in attributes
                            set!(attrs.nodes, get_id(attrs.nodes, String(attr_sub.id.id.id)), single_attr)
                        end
                    end

                end
            end

        end
    end

end

function set_subgraph!(subs::gvSubGraph, g::ParserCombinator.Parsers.DOT.SubGraph, nodes::gvNodes, directed::Bool)
    for stm in g.stmts
        if stm isa ParserCombinator.Parsers.DOT.SubGraph
            attrs3 = []
            for attr_sub in stm.stmts

                if attr_sub isa ParserCombinator.Parsers.DOT.NodeAttributes
                    for attr3 in attr_sub.attrs
                        set!(subs.node_options, String(attr3.name.id), check_value(String(attr3.value.id)))
                        push!(attrs3, attr3)
                    end

                elseif attr_sub isa ParserCombinator.Parsers.DOT.Node
                    push!(subs.nodes, gvNode(get_id(nodes, String(attr_sub.id.id.id)), String(attr_sub.id.id.id), Properties()))
                    for attr in attrs3
                        set!(subs.nodes, get_id(nodes, String(attr_sub.id.id.id)),
                            Property(String(attr.name.id), check_value(String(attr.value.id))))
                    end
                end
            end
        elseif stm isa ParserCombinator.Parsers.DOT.Node
            push!(subs.nodes, gvNode(get_id(nodes, String(stm.id.id.id)), String(stm.id.id.id), Properties()))
            for attr in stm.attrs
                set!(subs.nodes, get_id(nodes, String(stm.id.id.id)),
                    Property(String(attr.name.id), check_value(String(attr.value.id))))
                set!(subs.node_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.Attribute
            set!(subs.graph_options, String(stm.name.id), check_value(String(stm.value.id)))
        elseif stm isa ParserCombinator.Parsers.DOT.NodeAttributes
            for attr in stm.attrs
                set!(subs.node_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.EdgeAttributes
            for attr in stm.attrs
                set!(subs.edge_options, String(attr.name.id), check_value(String(attr.value.id)))
            end
        elseif stm isa ParserCombinator.Parsers.DOT.Edge
            for i = 1:(length(stm.nodes)-1)
                push!(subs.edges, gvEdge(get_id(nodes, String(stm.nodes[i].id.id)), get_id(nodes, String(stm.nodes[i+1].id.id)), Properties()))

                if (directed == false)
                    push!(subs.edges, gvEdge(get_id(nodes, String(stm.nodes[i+1].id.id)), get_id(nodes, String(stm.nodes[i].id.id)), Properties()))
                end

                if !(isempty(stm.attrs))
                    for attr in stm.attrs
                        set!(subs.edges, get_id(nodes, String(stm.nodes[i].id.id)), get_id(nodes, String(stm.nodes[i+1].id.id)),
                            Property(String(attr.name.id), check_value(String(attr.value.id))))
                        if (directed == false)
                            set!(subs.edges, get_id(nodes, String(stm.nodes[i+1].id.id)), get_id(nodes, String(stm.nodes[i].id.id)),
                                Property(String(attr.name.id), check_value(String(attr.value.id))))
                        end
                    end
                end
            end
        end

    end

end

function get_AbstractSimpleWeightedGraph(attrs::GraphvizAttributes)
    nodes = attrs.nodes
    ndim = length(nodes)
    adj = zeros(ndim, ndim)

    directed = val(attrs.plot_options, "directed")

    for e in attrs.edges

        weight = val(attrs.edges, e.from, e.to, "xlabel")
        if !(isempty(weight))
            adj[e.from, e.to] = parse(Float64, weight)
            (directed == false) ? adj[e.to, e.from] = parse(Float64, weight) : nothing
        else
            adj[e.from, e.to] = 1
            (directed == false) ? adj[e.to, e.from] = 1 : nothing
        end
    end

    if directed
        return SimpleWeightedDiGraph(adj)
    else
        return SimpleWeightedGraph(adj)
    end
end

const special_characters = ["+", "-", "&", " "]
function check_value(value::T) where {T}
    if value isa String
        if (contains(value, "\"") == false)
            value = "\"" * value * "\""
        end
        # if isempty(value)
        #     value = "\" " * "\""
        # else
        #     for e in special_characters
        #         if (contains(value, e) == true)
        #             value = "\"" * value * "\""
        #             break
        #         end
        #     end
        # end
    end
    return value
end

function insert_at!(str::String, insert::String, at::Int)
    len = length(str)
    str_f = ""
    str_e = ""
    new_string = ""

    if (at > 0) && (at <= len)
        str_f = str[1:at-1]
        str_e = str[at:end]
        new_string = str_f * insert * str_e
    end

    return new_string
end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz