text:
###############################################################################
#
#   fmpz_abs_series.jl : Power series over flint fmpz integers
#
###############################################################################

export fmpz_abs_series, FmpzAbsSeriesRing, PowerSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fmpz_abs_series)
   if iszero(a)
      return deepcopy(a)    # 0 + O(x^n)
   end
   prec = length(a) - 1
   prec < 0 && throw(DomainError())
   z = fmpz_abs_series(Array{fmpz}(0), 0, prec)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FmpzAbsSeriesRing}) = fmpz_abs_series

parent_type(::Type{fmpz_abs_series}) = FmpzAbsSeriesRing

base_ring(R::FmpzAbsSeriesRing) = R.base_ring

var(a::FmpzAbsSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FmpzAbsSeriesRing) = R.prec_max

function normalise(a::fmpz_abs_series, len::Int)
   if len > 0
      c = fmpz()
      ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_abs_series}, Int), c, a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_abs_series}, Int), c, a, len - 1)
      end
   end

   return len
end

function length(x::fmpz_abs_series)
   return ccall((:fmpz_poly_length, :libflint), Int, (Ref{fmpz_abs_series},), x)
end

precision(x::fmpz_abs_series) = x.prec

function coeff(x::fmpz_abs_series, n::Int)
   if n < 0
      return fmpz(0)
   end
   z = fmpz()
   ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_abs_series}, Int), z, x, n)
   return z
end

zero(R::FmpzAbsSeriesRing) = R(0)

one(R::FmpzAbsSeriesRing) = R(1)

function gen(R::FmpzAbsSeriesRing)
   z = fmpz_abs_series([fmpz(0), fmpz(1)], 2, max_precision(R))
   z.parent = R
   return z
end

function deepcopy_internal(a::fmpz_abs_series, dict::ObjectIdDict)
   z = fmpz_abs_series(a)
   z.prec = a.prec
   z.parent = parent(a)
   return z
end

function isgen(a::fmpz_abs_series)
   return precision(a) == 0 || ccall((:fmpz_poly_is_x, :libflint), Bool,
                            (Ref{fmpz_abs_series},), a)
end

iszero(a::fmpz_abs_series) = length(a) == 0

isunit(a::fmpz_abs_series) = valuation(a) == 0 && isunit(coeff(a, 0))

function isone(a::fmpz_abs_series)
   return precision(a) == 0 || ccall((:fmpz_poly_is_one, :libflint), Bool,
                                (Ref{fmpz_abs_series},), a)
end

# todo: write an fmpz_poly_valuation
function valuation(a::fmpz_abs_series)
   for i = 1:length(a)
      if !iszero(coeff(a, i - 1))
         return i - 1
      end
   end
   return precision(a)
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FmpzAbsSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fmpz_abs_series}) = show_minus_one(GenRes{fmpz})

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fmpz_abs_series)
   z = parent(x)()
   ccall((:fmpz_poly_neg, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}),
               z, x)
   z.prec = x.prec
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpz_abs_series, b::fmpz_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function -(a::fmpz_abs_series, b::fmpz_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fmpz_poly_sub_series, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function *(a::fmpz_abs_series, b::fmpz_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(a)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   z = parent(a)()
   z.prec = prec

   if lena == 0 || lenb == 0
      return z
   end

   lenz = min(lena + lenb - 1, prec)

   ccall((:fmpz_poly_mullow, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, a, b, lenz)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpz_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fmpz_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, y, x)
   return z
end

*(x::fmpz_abs_series, y::Int) = y * x

function *(x::fmpz, y::fmpz_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fmpz_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz}),
               z, y, x)
   return z
end

*(x::fmpz_abs_series, y::fmpz) = y * x

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpz_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   z.prec = x.prec + len
   ccall((:fmpz_poly_shift_left, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, x, len)
   return z
end

function shift_right(x::fmpz_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   if len >= xlen
      z.prec = max(0, x.prec - len)
   else
      z.prec = x.prec - len
      ccall((:fmpz_poly_shift_right, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, x, len)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fmpz_abs_series, prec::Int)
   prec < 0 && throw(DomainError())
   if x.prec <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   ccall((:fmpz_poly_set_trunc, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, x, prec)
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpz_abs_series, b::Int)
   b < 0 && throw(DomainError())
   if precision(a) > 0 && isgen(a) && b > 0
      return shift_left(a, b - 1)
   elseif length(a) == 1
      return parent(a)([coeff(a, 0)^b], 1, a.prec)
   elseif b == 0
      z = one(parent(a))
      set_prec!(z, precision(a))
      return z
   else
      z = parent(a)()
      z.prec = a.prec + (b - 1)*valuation(a)
      z.prec = min(z.prec, max_precision(parent(a)))
      ccall((:fmpz_poly_pow_trunc, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int, Int),
               z, a, b, z.prec)
   end
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fmpz_abs_series, y::fmpz_abs_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)

   n = max(length(x), length(y))
   n = min(n, prec)

   return Bool(ccall((:fmpz_poly_equal_trunc, :libflint), Cint,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               x, y, n))
end

function isequal(x::fmpz_abs_series, y::fmpz_abs_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || length(x) != length(y)
      return false
   end
   return Bool(ccall((:fmpz_poly_equal, :libflint), Cint,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               x, y, length(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpz_abs_series, y::fmpz)
   if length(x) > 1
      return false
   elseif length(x) == 1
      z = fmpz()
      ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
                       (Ref{fmpz}, Ref{fmpz_abs_series}, Int), z, x, 0)
      return ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}, Int), z, y, 0)
   else
      return precision(x) == 0 || iszero(y)
   end
end

==(x::fmpz, y::fmpz_abs_series) = y == x

==(x::fmpz_abs_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::fmpz_abs_series) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_abs_series, y::fmpz_abs_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   v2 = valuation(y)
   v1 = valuation(x)
   if v2 != 0
      if v1 >= v2
         x = shift_right(x, v2)
         y = shift_right(y, v2)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec, y.prec - v2 + v1)
   z = parent(x)()
   z.prec = prec
   ccall((:fmpz_poly_div_series, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, x, y, prec)
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz_abs_series, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fmpz_poly_scalar_divexact_si, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, x, y)
   return z
end

function divexact(x::fmpz_abs_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fmpz_poly_scalar_divexact_fmpz, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz}),
               z, x, y)
   return z
end

divexact(x::fmpz_abs_series, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpz_abs_series)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ccall((:fmpz_poly_inv_series, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function setcoeff!(z::fmpz_abs_series, n::Int, x::fmpz)
   ccall((:fmpz_poly_set_coeff_fmpz, :libflint), Void,
                (Ref{fmpz_abs_series}, Int, Ref{fmpz}),
               z, n, x)
   return z
end

function mul!(z::fmpz_abs_series, a::fmpz_abs_series, b::fmpz_abs_series)
   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(z)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = min(lena + lenb - 1, prec)
   if lenz < 0
      lenz = 0
   end

   z.prec = prec
   ccall((:fmpz_poly_mullow, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function addeq!(a::fmpz_abs_series, b::fmpz_abs_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   a.prec = prec
   ccall((:fmpz_poly_add_series, :libflint), Void,
                (Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Ref{fmpz_abs_series}, Int),
               a, a, b, lenz)
   return a
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fmpz_abs_series}, ::Type{T}) where {T <: Integer} = fmpz_abs_series

promote_rule(::Type{fmpz_abs_series}, ::Type{fmpz}) = fmpz_abs_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FmpzAbsSeriesRing)()
   z = fmpz_abs_series()
   z.prec = a.prec_max
   z.parent = a
   return z
end

function (a::FmpzAbsSeriesRing)(b::Integer)
   if b == 0
      z = fmpz_abs_series()
      z.prec = a.prec_max
   else
      z = fmpz_abs_series([fmpz(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FmpzAbsSeriesRing)(b::fmpz)
   if iszero(b)
      z = fmpz_abs_series()
      z.prec = a.prec_max
   else
      z = fmpz_abs_series([b], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FmpzAbsSeriesRing)(b::fmpz_abs_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FmpzAbsSeriesRing)(b::Array{fmpz, 1}, len::Int, prec::Int)
   z = fmpz_abs_series(b, len, prec)
   z.parent = a
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

function PowerSeriesRing(R::FlintIntegerRing, prec::Int, s::AbstractString;  model=:capped_relative, cached = true)
   S = Symbol(s)

   if model == :capped_relative
      parent_obj = FmpzRelSeriesRing(prec, S, cached)
   elseif model == :capped_absolute
      parent_obj = FmpzAbsSeriesRing(prec, S, cached)
   end

   return parent_obj, gen(parent_obj)
end


text:
function tt_normal_model_test(t::Trace, data::Vector{T}, params...) where T
    loc = sample(t, "loc", Normal())
    scale = sample(t, "scale", LogNormal())
    obs = Vector{Float64}(undef, length(data))
    for (i, d) in enumerate(data)
        obs[i] = observe(t, "obs $i", Normal(loc, scale), d)
    end
    obs
end

@testset "typed trace 1" begin
    data = rand(5)
    # to compile
    _ = likelihood_weighting(tt_normal_model_test, data)
    _ = likelihood_weighting(tt_normal_model_test, (String, Float64), data)
    # to time
    @time ut_res = likelihood_weighting(tt_normal_model_test, data; nsamples=100000)
    @time t_res = likelihood_weighting(tt_normal_model_test, (String, Float64), data; nsamples=100000)
end

function is_proposal(t, data, loc_guess, std_guess)
    propose(t, "loc", Normal(loc_guess, std_guess * 2.0))
    propose(t, "scale", truncated(Normal(std_guess, std_guess * 2.0), 0.1, Inf))
end

@testset "typed trace 2" begin
    data = rand(5)
    loc_guess = mean(data)
    std_guess = std(data)
    types = (String, Float64)
    _ = importance_sampling(tt_normal_model_test, is_proposal; params = (data, loc_guess, std_guess))
    _ = importance_sampling(tt_normal_model_test, is_proposal, types; params = (data, loc_guess, std_guess))
    @time importance_sampling(tt_normal_model_test, is_proposal; params = (data, loc_guess, std_guess), nsamples=100000)
    @time importance_sampling(tt_normal_model_test, is_proposal, types; params = (data, loc_guess, std_guess), nsamples=100000)
end

tt_loc_proposal(t0, t1, params...) = propose(t1, "loc", Normal(t0["loc"].value, 0.5))
tt_scale_proposal(t0, t1, params...) = propose(t1, "scale", LogNormal(log(t0["scale"].value), 0.5))

@testset "typed trace 3" begin
data = rand(5)
# to compile
proposals = [tt_loc_proposal, tt_scale_proposal]
types = (String, Float64)
_ = mh(tt_normal_model_test, proposals; params = (data,), inverse_verbosity=Inf)
_ = mh(tt_normal_model_test, proposals, types; params = (data,), inverse_verbosity=Inf)
# to time
@time ut_res = mh(tt_normal_model_test, proposals; params = (data,), burn=1000, thin=100, num_iterations=50000, inverse_verbosity=Inf)
@time t_res = mh(tt_normal_model_test, proposals, types; params = (data,), burn=1000, thin=100, num_iterations=50000, inverse_verbosity=Inf)
end

ext: jl


text:
# EXCLUDE FROM TESTING
using KernelAbstractions
using  CUDAapi
if CUDAapi.has_cuda_gpu()
    using CUDAdrv
    using CuArrays
    CuArrays.allowscalar(false)
else
    exit()
end
using MPI
using Test

device(A) = typeof(A) <: Array ? CPU() : CUDA()

@kernel function kernel!(c, a, b)
  I = @index(Global)
  @inbounds a_val = a[I]
  @inbounds b_val = b[I]
  @inbounds c_val = c[I]

  for j in 1:33
    a_val = a_val * b_val + c_val
    b_val = a_val * b_val + c_val
    c_val = a_val * b_val + c_val
  end

  @inbounds c[I] = a_val * b_val + c_val
end

function mpiyield()
    MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, MPI.COMM_WORLD)
    yield()
end

function __Irecv!(request, buf, src, tag, comm)
    request[1] = MPI.Irecv!(buf, src, tag, comm)
end

function __Isend!(request, buf, dst, tag, comm)
    request[1] = MPI.Isend(buf, dst, tag, comm)
end

function __testall!(requests; dependencies=nothing)
    done = false
    while !done
        done, _ = MPI.Testall!(requests)
        yield()
    end
end

function exchange!(d_send_buf, h_send_buf, d_recv_buf, h_recv_buf, src_rank,
                   dst_rank, send_request, recv_request, comm;
                   dependencies=nothing)

    event = Event(__Irecv!, recv_request, h_recv_buf, src_rank, 666, comm;
                  dependencies = dependencies, progress=mpiyield)
    event = Event(__testall!, recv_request; dependencies = event, progress=mpiyield)

    recv_event = async_copy!(device(d_recv_buf), d_recv_buf, h_recv_buf;
                             dependencies = event, progress=mpiyield)

    event = async_copy!(device(d_send_buf), h_send_buf, d_send_buf;
                        dependencies = dependencies, progress=mpiyield)

    event = Event(__Isend!, send_request, h_send_buf, dst_rank, 666, comm;
                  dependencies = dependencies, progress=mpiyield)
    send_event = Event(__testall!, send_request; dependencies = event, progress=mpiyield)

    return recv_event, send_event
end

function init(T, M, comm)
    d_send_buf = CuArrays.zeros(T, M)
    d_recv_buf = CuArrays.zeros(T, M)
    h_send_buf = zeros(T, M)
    h_recv_buf = zeros(T, M)

    fill!(d_send_buf, MPI.Comm_rank(comm))
    fill!(d_recv_buf, -1)
    fill!(h_send_buf, -1)
    fill!(h_recv_buf, -1)

    send_request = fill(MPI.REQUEST_NULL, 1)
    recv_request = fill(MPI.REQUEST_NULL, 1)

    return (d_send_buf = d_send_buf,
            d_recv_buf = d_recv_buf,
            h_send_buf = h_send_buf,
            h_recv_buf = h_recv_buf,
            send_request = send_request,
            recv_request = recv_request,
            comm = MPI.Comm_dup(comm))
end

function main()
    if !MPI.Initialized()
        MPI.Init()
    end
    comm = MPI.COMM_WORLD
    MPI.Barrier(comm)

    dst_rank = mod(MPI.Comm_rank(comm)+1, MPI.Comm_size(comm))
    src_rank = mod(MPI.Comm_rank(comm)-1, MPI.Comm_size(comm))

    T = Float64
    M = 10
    N = 1024^2

    a = CuArrays.rand(T, N)
    b = CuArrays.rand(T, N)
    c = CuArrays.rand(T, N)

    e1 = init(T, M, comm)
    e2 = init(T, M, comm)

    pc_kernel! = kernel!(CUDA(), 256)
    KernelAbstractions.precompile(pc_kernel!, c, a, b; ndrange=(length(c),))

    synchronize()

    event = Event(CUDA())

    @info "Here 1"
    recv_e1, send_e1 = exchange!(e1.d_send_buf, e1.h_send_buf,
                                 e1.d_recv_buf, e1.h_recv_buf,
                                 src_rank, dst_rank,
                                 e1.send_request, e1.recv_request,
                                 e1.comm; dependencies=event)

    @info "Here 3"
    event = kernel!(CUDA(), 256)(c, a, b, ndrange=length(c), dependencies = event)

    @info "Here 4"
    wait(recv_e1)

    @info "Here 5"
    @test all(e1.d_recv_buf .== src_rank)

    @info "Here 6"
    wait(send_e1)

    @info "Here 7"
    wait(event)
end

main()


text:
import Test: @test, @testset, @test_broken
import SpecialFunctions: besselj
import Cubature: hcubature
import LinearAlgebra: dot, norm
import Luna: Modes, Capillary, Grid
import Luna.PhysData: c, roomtemp, ref_index_fun, ε_0, μ_0
import Luna.PhysData: wlfreq

@testset "Capillary" begin
@testset "loss" begin
    λ = 800e-9
    ω = wlfreq(λ)
    a = 125e-6
    m = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced)
    @test isapprox(Capillary.losslength(m, ω), 7.0593180702769, rtol=1e-5)
    @test isapprox(Capillary.dB_per_m(m, ω), 0.6152074146252722, rtol=1e-5)
    @test Capillary.dB_per_m(m, ω) ≈ 8*Capillary.dB_per_m(Capillary.MarcatiliMode(2a, :He, 1.0, model=:reduced), ω)
end

@testset "normalisation" begin
    # Copied definitions from Modes.jl to force manual calculation of Aeff and N
    function Aeff(m; z=0)
        em = Modes.absE(m, z=z)
        dl = Modes.dimlimits(m, z=z)
        # Numerator
        function Aeff_num(xs)
            e = em(xs)
            dl[1] == :polar ? xs[1]*e^2 : e^2
        end
        val, err = hcubature(Aeff_num, dl[2], dl[3])
        num = val^2
        # Denominator
        function Aeff_den(xs)
            e = em(xs)
            dl[1] == :polar ? xs[1]*e^4 : e^4
        end
        den, err = hcubature(Aeff_den, dl[2], dl[3])
        return num / den
    end
    function N(m; z=0)
        f = Modes.field(m, z=z)
        dl = Modes.dimlimits(m, z=z)
        function Nfunc(xs)
            E = f(xs)
            ret = sqrt(ε_0/μ_0)*dot(E, E)
            dl[1] == :polar ? xs[1]*ret : ret
        end
        val, err = hcubature(Nfunc, dl[2], dl[3])
        0.5*abs(val)
    end

    a = 125e-6
    for n=1:10
        m = Capillary.MarcatiliMode(a, :He, 1.0, n=n)
        @test Modes.N(m) ≈ N(m)
        @test Modes.Aeff(m) ≈ Aeff(m)
    end
    m = Capillary.MarcatiliMode(a, :He, 1.0, n=0, kind=:TE)
    @test Modes.N(m) ≈ N(m)
    @test Modes.Aeff(m) ≈ Aeff(m)
    m = Capillary.MarcatiliMode(a, :He, 1.0, n=0, kind=:TM)
    @test Modes.N(m) ≈ N(m)
    @test Modes.Aeff(m) ≈ Aeff(m)

    a0 = a
    aL = a/2
    L = 1
    afun = let a0=a0, aL=aL, L=L
        afun(z) = a0 + (aL-a0)*z/L
    end
    m = Capillary.MarcatiliMode(afun, :He, 1, loss=false, model=:full)
    @test Modes.N(m) ≈ N(m)
    @test Modes.N(m, z=L/2) ≈ N(m, z=L/2)
    @test Modes.N(m, z=L) ≈ N(m, z=L)
    @test Modes.Aeff(m) ≈ Aeff(m)
    @test Modes.Aeff(m, z=L/2) ≈ Aeff(m, z=L/2)
    @test Modes.Aeff(m, z=L) ≈ Aeff(m, z=L)
end

@testset "β, α" begin
    a = 125e-6
    m = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced)
    λ = 1e-9 .* collect(range(70, stop=7300, length=128))
    ω = wlfreq.(λ)
    @test all(isfinite.(Capillary.β.(m, ω)))
    @test all(isreal.(Capillary.β.(m, ω)))
    @test all(isreal.(Capillary.β.(Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced), ω)))
    @test all(isreal.(Capillary.β.(Capillary.MarcatiliMode(a, :He, 10.0, model=:reduced), ω)))
    @test all(isreal.(Capillary.β.(Capillary.MarcatiliMode(a, :He, 50.0, model=:reduced), ω)))
    @test all(isfinite.(Capillary.α.(m, ω)))
    @test all(isreal.(Capillary.α.(m, ω)))
end

@testset "ZDW/Aeff" begin
    @test abs(1e9*Capillary.zdw(Capillary.MarcatiliMode(125e-6, :He, 0.4, model=:reduced)) - 379) < 1
    @test abs(1e9*Capillary.zdw(Capillary.MarcatiliMode(75e-6, :He, 5.9, model=:reduced)) - 562) < 1
    @test Capillary.Aeff(Capillary.MarcatiliMode(75e-6, :He, 1.0, model=:reduced)) ≈ 8.42157534886545e-09
end

@testset "dispersion" begin
    # tests based on symbolic results in symbolic_marcatilli.py
    m = Capillary.MarcatiliMode(50e-6, :Ar, 2.0, model=:reduced, T=294.0)
    ω = wlfreq(800e-9)
    @test isapprox(Capillary.β(m, ω), 7857866.63899973, rtol=1e-15)
    @test isapprox(Capillary.dispersion(m, 1, ω), 3.33744405855838e-9, rtol=1e-13)
    @test isapprox(Capillary.dispersion(m, 2, ω), -1.68187063100098e-29, rtol=3e-7)
    @test isapprox(Capillary.dispersion(m, 3, ω), 8.44025458282146e-44, rtol=2e-6)
    @test isapprox(Capillary.dispersion(m, 4, ω), -1.13289638668388e-58, rtol=8e-6)
    @test isapprox(Capillary.dispersion(m, 5, ω), 2.45046358889987e-73, rtol=2e-4)
    @test isapprox(Capillary.zdw(m), 7.289431065526978e-07, rtol=1e-7)

    rfg = ref_index_fun(:Ar, 2.0, 294.0)
    coren = (ω; z) -> rfg(wlfreq(ω))
    cladn = (ω; z) -> 1.45
    m = Capillary.MarcatiliMode(50e-6, 1, 1, :HE, 0.0, coren, cladn)
    @test isapprox(Capillary.β(m, ω), 7857865.68069111, rtol=1e-15)
    @test isapprox(Capillary.dispersion(m, 1, ω), 3.33744357233436e-9, rtol=1e-13)
    @test isapprox(Capillary.dispersion(m, 2, ω), -1.68088637707608e-29, rtol=7e-7)
    @test isapprox(Capillary.dispersion(m, 3, ω), 8.43560292665704e-44, rtol=1e-6)
    @test isapprox(Capillary.dispersion(m, 4, ω), -1.13224034494362e-58, rtol=1e-5)
    @test isapprox(Capillary.dispersion(m, 5, ω), 2.44896931597213e-73, rtol=1e-4)
    @test isapprox(Capillary.zdw(m), 7.289459086128214e-07, rtol=1e-8)
    @test isapprox(Capillary.α(m, ω), 2.21505826015050, rtol=1e-14)

    rfg = ref_index_fun(:Ar, 2.0, 294.0)
    coren = (ω; z) -> rfg(wlfreq(ω))
    cladn = (ω; z) -> 0.036759+im*5.5698
    m = Capillary.MarcatiliMode(50e-6, 1, 1, :HE, 0.0, coren, cladn)
    @test isapprox(Capillary.β(m, ω), 7857863.68326066, rtol=1e-15)
    @test isapprox(Capillary.dispersion(m, 1, ω), 3.33744527275634e-9, rtol=1e-13)
    @test isapprox(Capillary.dispersion(m, 2, ω), -1.89802237417366e-29, rtol=2e-7)
    @test isapprox(Capillary.dispersion(m, 3, ω), 8.80530032930868e-44, rtol=7e-7)
    @test isapprox(Capillary.dispersion(m, 4, ω), -1.21092167154788e-58, rtol=1e-5)
    @test isapprox(Capillary.dispersion(m, 5, ω), 2.64991510536236e-73, rtol=5e-5)
    @test isapprox(Capillary.zdw(m), 7.225347947615157e-07, rtol=2e-8)
    @test isapprox(Capillary.α(m, ω), 0.0290115706883820, rtol=1e-14)
    @test Capillary.Aeff(Capillary.MarcatiliMode(75e-6, :He, 1.0)) ≈ 8.42157534886545e-09
end

@testset "to_space" begin
    ms = (Capillary.MarcatiliMode(125e-6, :He, 1.0),
          Capillary.MarcatiliMode(125e-6, :He, 1.0, m=2, ϕ=π/2))
    components = :xy
    xs = (10e-6, π/7)
    ts = Modes.ToSpace(ms, components=components)
    Emω = Array{ComplexF64,2}(undef, 8192, 2)
    fill!(Emω, 0.3+0.5im)
    Erω1 = Modes.to_space(Emω, xs, ts)
    Erω! = copy(Erω1)
    Modes.to_space!(Erω!, Emω, xs, ts)
    Erω2 = Modes.to_space(Emω, xs, ms, components=components)
    @test all(Erω2 .== Erω!)
    @test all(Erω1 .== Erω!)
end

end


text:
function SeisPatchProcess(in,out,param=Dict())
    # read data from disk, split into multidimensional overlapping patches,
    # apply processes, merge patches with tapers, and write to disk.
    # Processing of patches can be done in parallel by running the code using
    # (for example): julia -p 2 script_name.jl
    #
    # Important Notice: You must make declare global variables f and f_param on
    # every processor. You can do this in your main function by typing (for
    # example):
    # @everywhere global f_param = ["style"=>"mxmyhxhy", "Niter"=> 100,
    #                               "alpha"=>1,"fmax"=>80]
    # @everywhere global f = [SeisPOCS]
    #
    # f is an array of functions that have the following syntax:
    # d2, h2 = f(d1,h1,f_param), where
    # param is a dictionary (Dict) of parameters for the function.
    #
    # note that param should contain parameters for the patching and unpatching
    # operations.
    #
    # to execute the code type:  julia -p 4 my_code.jl where 4 can be replaced
    # with the number of processors you wish to use.

    patch_names = SeisPatch(in,out,param)
    a = pmap(MyProcess,patch_names)
    SeisUnPatch(patch_names,out,param)
    #for ipatch = 1 : length(patch_names)
    #	SeisRemove(patch_names[ipatch])
    #end

end

function MyProcess(filename)
    d1,h1 = SeisRead(filename)
    for ifunc = 1 : length(f)
	func = f[ifunc]
	d2,h2 = func(d1,h1,f_param)
	d1 = copy(d2)
	h1 = copy(h2)
    end
    SeisWrite(filename,d1,h1)
    return(1)
end


text:
using Gurobi

function LP_init(nlp::Class_CUTEst, pars::Class_parameters, timer::class_advanced_timer)
    start_advanced_timer(timer, "INIT/x")
    x0 = suggested_starting_point(nlp)

    #@show norm(x0)

    if pars.init.start_satisfying_bounds
      x = projection_onto_bounds_ipopt_style( nlp, pars, x0 )
    end

    pause_advanced_timer(timer, "INIT/x")

    a, J, g = eval_init(nlp, pars, timer, x)

    s = a - 2 * minimum(a)
    @assert(all(s .>= 0.0))

    # make bounds exact.
    bis = bound_indicies(nlp)
    s[bis] = a[bis]
    @assert(all(s .>= 0.0))

    ncon = length(a)
    nvar = length(g)

    m = Model(solver=GurobiSolver())
    @variable(m, yVar[1:ncon] >= 0.0)
    @variable(m, muVar >= 0.0)

    y0 = 10.0 * ones(length(a)) # TEMP

    theta = 1e-4 # parameter
    penalty = theta
    scale = mean(abs.(J))
    sc_J = J / scale
    sc_g = g / scale
    #@show sc_J
    @show mean(abs.(J))

    if false
        @objective(m, Min, (sc_J' * yVar - sc_g)' * (sc_J' * yVar - sc_g) + penalty * (yVar - y0)' * (yVar - y0))
        status = solve(m)
    else
        @variable(m, z1[1:nvar] >= 0.0)
        @variable(m, z2[1:nvar] >= 0.0)
        @variable(m, t >= 0.0)
        @objective(m, Min, sum(z1[i]^2 for i = 1:nvar) + penalty * sum((yVar - y0).^2)) # + sum(log(yVar[j]) for j=1:ncon))

        if false
            for j = 1:ncon
                @constraint(m, yVar[j] <= muVar * 10.0 / s[j])
                @constraint(m, yVar[j] >= muVar * 0.1 / s[j])
            end
        end

        for i = 1:nvar
            @constraint(m, sum(sc_J[i,j] * yVar[j] for j = 1:nvar) - sc_g[i] <= z1[i])
            @constraint(m, sc_g[i] - sum(sc_J[i,j] * yVar[j] for j = 1:nvar) <= z1[i])
            #@constraint(m, sum(sc_J[i,j] * yVar[j] for j = 1:nvar) - sc_g[i] <= z2[i])
            #@constraint(m, sc_g[i] - sum(sc_J[i,j] * yVar[j] for j = 1:nvar) <= z2[i])
        end
        status = solve(m)
    end

    #mu = getvalue(muVar)
    y = getvalue(yVar)
    #@show getvalue(muVar)
    mu = - minimum(a) * mean(y) #/ length(y)
    #mu = mean(s .* y)
    obj = dot(J' * y - g, J' * y - g) + penalty * dot(y - y0, y - y0)
    @show maximum(y), getobjectivevalue(m), penalty, obj
    buffer = 2.0
    y_c = mu ./ s;
    @show norm(y_c,Inf), mu, minimum(s)
    y = min.( y_c / (pars.ls.comp_feas * buffer), max.(y, pars.ls.comp_feas * y_c * buffer)) # project onto complementarity constraints

    #@show mu, theta
    #error("")

    start_advanced_timer(timer, "INIT/construct_class")
    init_point = Class_point(x, y, s, mu)
    check_for_nan(init_point)
    init_it = Class_iterate(init_point, nlp, Class_local_info(), timer, pars);
    pause_advanced_timer(timer, "INIT/construct_class")

    return init_it
end


text:
function readIDs( xs::XLSX.Worksheet, idCol::Int, nRecords::Int, mpSim::MPsim )

    if idCol == 0
        return string.( "Init", 1:nRecords ), falses( nRecords )
    end  # if idCol == 0

    ids = xs[XLSX.CellRange( 2, idCol, nRecords + 1, idCol )][:]

    let
        isMissing = isa.( ids, Missing )
        nMissing = count( isMissing )

        if nMissing > 0
            ids[ isMissing ] = string.( "Init", 1:nMissing )

            if mpSim.showInfo
                @info string( "Generated ids for ", nMissing,
                    " record(s)." )
            end  # if mpSim.showInfo
        end  # if nMissing > 0
    end  # let

    ids = Vector{String}( strip.( string.( ids ) ) )
    isDuplicate = flagDuplicates( ids )
    nDuplicate = count( isDuplicate )

    if nDuplicate > 0
        ids = ids[.!isDuplicate]

        if mpSim.showInfo
            @info string( "Removed ", nDuplicate, " duplicate record(s)." )
        end  # if mpSim.showInfo
    end  # if nDuplicate > 0

    return ids, isDuplicate

end  # readIDs( xs, idCol, nRecords, mpSim )

function flagDuplicates( ids::Vector{String} )

    uniqueSet = Set{String}()
    duplicateSet = Set{String}()

    for id in ids
        if id ∉ uniqueSet
            push!( uniqueSet, id )
        elseif id ∉ duplicateSet
            push!( duplicateSet, id )
        end  # if id ∉ uniqueSet
    end  # for id in ids

    return map( id -> id ∈ duplicateSet, ids )

end  # flagDuplicates( ids )

function readTimes( xs::XLSX.Worksheet, timeCol::Int, refDate::Date,
    dateType::Symbol, nRecords::Int, isBad::Vector{Bool}, mpSim::MPsim,
    timeFactor::Float64=1.0 )

    if timeCol == 0
        return zeros( Float64, nRecords - count( isBad ) )
    end  # if tenureCol == 0

    times = xs[XLSX.CellRange( 2, timeCol, nRecords + 1, timeCol )][:]
    times = times[.!isBad]

    let
        isMissing = map( times ) do timeVal
            return ( timeVal isa Missing ) || (
                ( timeVal isa AbstractString ) &&
                ( tryparse( Float64, timeVal ) isa Nothing ) )
        end  # map( times ) do ltt

        nMissing = count( isMissing )

        if nMissing > 0
            times[ isMissing ] .= 0.0

            if mpSim.showInfo
                @info string( "Missing time values for ", nMissing,
                    " record(s), set them to 0.0." )
            end  # if mpSim.showInfo
        end  # if nMissing > 0

        computeTimes!( times, refDate, dateType )
        times = max.( times, 0.0 )
    end  # let

    return times * timeFactor

end  # readTimes( xs, timeCol, refDate, dateType, nRecords, isBad, mpSim,
     #   timeFactor )

function readNodes( xs::XLSX.Worksheet, nodeCol::Int, nRecords::Int,
    isBad::Vector{Bool}, mpSim::MPsim )

    if nodeCol == 0
        return fill( "NULL", nRecords - count( isBad ) )
    end  # if tenureCol == 0

    nodes = xs[XLSX.CellRange( 2, nodeCol, nRecords + 1, nodeCol )][:]
    nodes = nodes[.!isBad]

    let
        isMissing = isa.( nodes, Missing )
        nMissing = count( isMissing )

        if nMissing > 0
            nodes[ isMissing ] .= "NULL"

            if mpSim.showInfo
                @info string( "Missing initial nodes for ", nMissing,
                    " record(s), set them to NULL." )
            end  # if mpSim.showInfo
        end  # if nMissing > 0
    end  # let

    return Vector{String}( strip.( string.( nodes ) ) )

end # readNodes( xs, nodeCol, nRecords, isBad, mpSim )

function readData( xs::XLSX.Worksheet, ids::Vector{String},
    lastTransTimes::Vector{Float64}, tenures::Vector{Float64},
    ages::Vector{Float64}, nodes::Vector{String}, colsToImport::Vector{Int},
    nRecords::Int, isBad::Vector{Bool}, mpSim::MPsim )

    attributeNames = Vector{String}( strip.( string.( getindex.(
        Ref(xs), 1, colsToImport ) ) ) )
    personnelData = map( colNr -> xs[XLSX.CellRange( 2, colNr,
        nRecords + 1, colNr )][:], colsToImport )
    personnelData = hcat( personnelData... )[.!isBad, :]
    nRecords -= count( isBad )

    isBad = falses( nRecords )

    for ii in eachindex( colsToImport )
        isBad .|= isa.( personnelData[:, ii], Missing )
    end  # for ii in eachindex( personnelData )

    personnelData = Array{String}( strip.( string.(
        personnelData[.!isBad, :] ) ) )
    nBad = count( isBad )

    if nBad > 0
        ids = ids[.!isBad]
        lastTransTimes = lastTransTimes[.!isBad]
        tenures = tenures[.!isBad]
        ages = ages[.!isBad]
        nodes = nodes[.!isBad]

        if mpSim.showInfo
            @info string( "Removed ", nBad,
                " record(s) with missing data." )
        end  # # if mpSim.showInfo
    end  # if nBad > 0

    tenures = max.( tenures, lastTransTimes )
    ages = max.( ages, tenures )
    personnelData = hcat( personnelData, ids, lastTransTimes, tenures, ages,
        nodes )

    return personnelData, attributeNames

end  # readData( xs, ids, lastTransTimes, tenures, ages, colsToImport,
     #   nRecords, isBad, mpSim )

function validateAttributeData( personnelData::Array{Any,2},
    attributeNames::Vector{String}, mpSim::MPsim )

    isBad = falses( size( personnelData, 1 ) )

    for ii in eachindex( attributeNames )
        # No need to validate anything if the attribute isn't defined.
        if !haskey( mpSim.attributeList, attributeNames[ii] )
            continue
        end  # if !haskey( mpSim.attributeList, attributeNames[ii] )

        attribute = mpSim.attributeList[attributeNames[ii]]
        isBad .|= map( val -> val ∉ attribute.possibleValues,
            personnelData[:, ii] )
    end  # for ii in eachindex( attributeNames )

    nBad = count( isBad )

    if nBad > 0
        personnelData = personnelData[.!isBad, :]

        if mpSim.showInfo
            @info string( "Removed ", nBad,
                " record(s) with inconsistent data." )
        end  # if mpSim.showInfo
    end  # if nBad > 0

    return personnelData

end  # validateAttributeData( personnelData, attributeNames, mpSim )

function generateSimAttributes!( mpSim::MPsim, personnelData::Array{Any,2}, attributeNames::Vector{String}, generateAttribute::BitArray{1} )

    attributes = Vector{Attribute}()

    for ii in filter( ii -> generateAttribute[ii],
        eachindex( generateAttribute ) )
        name = attributeNames[ii]
        attribute = Attribute( name )
        setPossibleAttributeValues!( attribute,
            Vector{String}( unique( personnelData[:, ii] ) ) )
        push!( attributes, attribute )
    end  # for ii in filter( ... )

    if isempty( attributes )
        return
    end  # if isempty( attributes )

    addSimulationAttribute!( mpSim, attributes... )

end  # generateSimAttributes!( mpSim, personnelData, attributeNames,
     #   generateAttribute )

function generateSimNodes!( mpSim::MPsim, personnelData::Array{Any,2},
    attributeNames::Vector{String} )

    nodes = personnelData[:, end]
    baseNodes = Vector{BaseNode}()

    for name in filter( name -> ( name != "NULL" ) &&
        ( !haskey( mpSim.baseNodeList, name ) ), unique( nodes ) )
        recordsOfNode = findall( nodes .== name )
        node = BaseNode( name )

        for ii in eachindex( attributeNames )
            attrName = attributeNames[ii]
            possibleVals =  unique( personnelData[recordsOfNode, ii] )

            if length( possibleVals ) == 1
                addNodeRequirement!( node, attrName, possibleVals[1] )
            end  # if length( possibleVals ) == 1
        end  # for ii in eachindex( attributeNames )

        push!( baseNodes, node )
    end  # for name in filter( ... )

    if isempty( baseNodes )
        return
    end  # if isempty( baseNodes )

    addSimulationBaseNode!( mpSim, baseNodes... )

end  # generateSimNodes!( mpSim, personnelData, attributeNames )

function inferBaseNodes( mpSim::MPsim, personnelData::Array{Any,2},
    attributeNames::Vector{String} )

    isNull = personnelData[:, end] .== "NULL"

    if !any( isNull )
        return
    end  # if !any( isNull )

    nodes = collect( keys( mpSim.baseNodeList ) )
    attributeVals = map( nodes ) do name
        node = mpSim.baseNodeList[name]
        attrVals = get.( Ref(node.requirements), attributeNames, nothing )
    end  # map( nodes ) do name

    for ii in findall( isNull )
        possibleNodes = map( eachindex( nodes ) ) do jj
            return all( isa.( attributeVals[jj], Nothing ) .|
                ( attributeVals[jj] .== personnelData[ii, 1:(end-5)] ) )
        end  # map( eachindex( nodes ) ) do jj

        possibleNodes = nodes[possibleNodes]

        if length( possibleNodes ) == 1
            personnelData[ii, end] = possibleNodes[1]
        end  # if length( possibleNodes ) == 1
    end  # for ii in findall( isNull )

end  # inferBaseNodes( mpSim, personnelData, attributeNames )

function addSnapshotAttrition( mpSim::MPsim, personnelData::Array{Any,2},
    nRecords::Int )

    attritionTimes = fill( +Inf, nRecords )

    for node in keys( mpSim.baseNodeList )
        nodeInds = personnelData[:, end] .== node
        attrition = mpSim.baseNodeList[node].attrition
        attritionTimes[nodeInds] = generateTimeToAttrition(
            mpSim.attritionSchemes[attrition],
            Vector{Float64}( personnelData[nodeInds, end-3] ) )
    end  # for node in keys( mpSim.baseNodeList )

    attritionTimes -= personnelData[:, end-3]

    return hcat( personnelData, attritionTimes )

end  # addSnapshotAttrition( mpSim, personnelData )

function saveSnapshotToDatabase( mpSim::MPsim, personnelData::Array{Any,2},
    attributeNames::Vector{String} )

    resetSimulation( mpSim, false )
    nRecords = size( personnelData, 1 )
    nAttributes = length( attributeNames )
    ageAtRecruitment = personnelData[:, end-2] - personnelData[:, end-3]
    nodes = personnelData[:, end-1]
    isValidNode = nodes .!= "NULL"
    DBInterface.execute( mpSim.simDB, "BEGIN TRANSACTION" )

    # The SQLite command that fills the personnel database.
    personnelCmd = map( 1:nRecords ) do ii
        node = nodes[ii]
        return string( "\n    ('",
            join( personnelData[ii, 1:(nAttributes + 1)], "', '" ), "', ",
            -personnelData[ii, end-3], ", ", ageAtRecruitment[ii], ", ",
            isValidNode[ii] ? string( "'", node, "'" ) : node, ", ",
            -personnelData[ii, end-4], ", ",
            personnelData[ii, end] == +Inf ? "NULL" : personnelData[ii, end],
            ", 'active')" )
    end  # map( 1:nRecords ) do ii

    personnelCmd = string( "INSERT INTO `", mpSim.persDBname, "` (`",
        join( attributeNames, "`, `" ), "`, `", mpSim.idKey,
        "`, timeEntered, ageAtRecruitment, currentNode, inNodeSince, expectedAttritionTime, status) VALUES",
        join( personnelCmd, "," ) )
    DBInterface.execute( mpSim.simDB, personnelCmd )

    # The SQLite command that fills the transition database.
    if any( isValidNode )
        transitionCmd = map( findall( isValidNode ) ) do ii
            return string( "\n    ('", personnelData[ii, end-5], "', ",
                - personnelData[ii, end-4], ", 'Init', '", nodes[ii], "')" )
        end  # map( findall( isValidNode ) ) do ii

        transitionCmd = string( "INSERT INTO `", mpSim.transDBname, "` (`",
            mpSim.idKey, "`, timeIndex, transition, targetNode) VALUES",
            join( transitionCmd, "," ) )
        DBInterface.execute( mpSim.simDB, transitionCmd )
    end  # if any( isValidNode )

    # The SQLite command that fills the attribute history database.
    historyCmd = map( eachindex( attributeNames ) ) do ii
        return join( string.( "\n    ('", personnelData[:, end-5], "', ",
            -personnelData[:, end-4], ", '", attributeNames[ii], "', '",
            personnelData[:, ii], "')" ), "," )
    end  # map( eachindex( attributeNames ) ) do ii

    historyCmd = string( "INSERT INTO `", mpSim.histDBname, "` (`", mpSim.idKey,
        "`, timeIndex, attribute, value) VALUES", join( historyCmd, "," ) )
    DBInterface.execute( mpSim.simDB, historyCmd )
    DBInterface.execute( mpSim.simDB, "COMMIT" )

end  # saveSnapshotToDatabase( mpSim, personnelData, attributeNames )

function saveNodeTime!( mpSim::MPsim, personnelData::Array{Any,2} )

    for node in keys( mpSim.baseNodeList )
        isIDinNode = personnelData[:, end-1] .== node
        setindex!.( Ref(mpSim.baseNodeList[node].inNodeSince),
            -personnelData[isIDinNode, end-4],
            personnelData[isIDinNode, end-5] )
    end  # for node in keys( mpSim.baseNodeList )

end  # saveNodeTime!( mpSim, personnelData )

function dateDiff( date1::Date, date2::Date, dateType::Symbol )
    dateType === :days && return (date1 - date2).value
    dateType === :years && return year(date1) - year(date2) +
        ( dayofyear(date1) - dayofyear(date2) ) / daysinyear(date1)

    ddiff = ( year(date1) - year(date2) ) * 12 + ( month(date1) - month(date2) )
    ddiff + ( day(date1) - day(date2) ) / daysinmonth(date1)
end  # dateDiff( date1, date2, dateType )

function computeTimes!( timeData::Vector, refDate::Date, dateType::Symbol )
    dateInds = findall(isa.( timeData, Date ))
    timeData[dateInds] = dateDiff.( refDate, timeData[dateInds],
        dateType )
end  # computeTimes!( timeData, refDate, dateType )


text:
import Arpack,LinearAlgebra#, NLsolve
#TODO: Degeneracy
# function householder_update(f,a=1)
#     order=length(f)-1
#     a=1.0/a
#
#     if order==1   #aka Newton's method
#         dz=-f[1]/(a*f[2])
#     elseif order==2 #aka Halley's method
#         dz=-2*f[1]*f[2]/((a+1)*f[2]^2-f[1]*f[3])
#     elseif order==3
#         dz=-3*f[1]*((a+1)*f[2]^2-f[1]*f[3]) / ((a^2 + 3*a + 2) *f[2]^3 - 3* (a + 1)* f[1]* f[2] *f[3] + f[1]^2* f[4])
#     elseif order==4
#         dz=-4*f[1]*((a^2 + 3*a + 2) *f[2]^3 - 3* (a + 1)* f[1]* f[2] *f[3] + f[1]^2* f[4]) / (-6* (a^2 + 3 *a + 2) *f[1]* f[2]^2* f[3] + (a^3 + 6 *a^2 + 11* a + 6)* f[2]^4 + f[1]^2 *(3* (a + 1)*f[3]^2 - f[1]* f[5]) + 4* (a + 1)*f[1]^2* f[4] *f[2])
#     else
#         dz=-5*f[1]*(-6* (a^2 + 3 *a + 2) *f[1]* f[2]^2* f[3] + (a^3 + 6 *a^2 + 11* a + 6)* f[2]^4 + f[1]^2 *(3* (a + 1)*f[3]^2 - f[1]* f[5]) + 4* (a + 1)*f[1]^2* f[4] *f[2]) / ((a + 1)*(a + 2)*(a + 3)*(a + 4)* f[2]^5 + 10* (a + 1)*(a + 2)*f[1]^2*f[4]*f[2]^2-10*(a + 1)*(a + 2)*(a + 3)*f[1]*f[2]^3*f[3]+f[1]^3*(f[1]*f[6] - 10 *(a + 1) *f[4]* f[3]) + 5 *(a + 1) *f[1]^2*f[2]* (3* (a + 2)* f[3]^2 - f[1]*f[5]))
#     end
#     return dz
# end

function householder_update(f)
    order=length(f)-1
    if order==1
        dz=-f[1]/f[2]
    elseif order==2
        dz=-f[1]*f[2]/ (f[2]^2-0.5*f[1]*f[3])
    elseif order==3
        dz=-(6*f[1]*f[2]^2-3*f[1]^2*f[3]) / (6*f[2]^3-6*f[1]*f[2]*f[3]+f[1]^2*f[4])
    elseif order == 4
        dz=-(4 *f[1] *(6* f[2]^3 - 6 *f[1] *f[2]* f[3] + f[1]^2* f[4]))/(24 *f[2]^4 - 36* f[1] *f[2]^2 *f[3] + 6*f[1]^2* f[3]^2 + 8*f[1]^2* f[2]* f[4] - f[1]^3* f[5])
    else
        dz=(5 *f[1] *(24 *f[2]^4 - 36* f[1] *f[2]^2* f[3] + 6 *f[1]^2* f[3]^2 + 8* f[1]^2* f[2]* f[4] - f[1]^3* f[5]))/(-120* f[2]^5 + 240* f[1]* f[2]^3* f[3] - 60* f[1]^2* f[2]^2* f[4] + 10* f[1]^2* f[2]* (-9* f[3]^2 + f[1]* f[5]) + f[1]^3* (20* f[3]* f[4] - f[1]* f[6]))
    end
    return dz
end
"""
    sol::Solution, n, flag = householder(L::LinearOperatorFamily, z; <keyword arguments>)

Use a Householder method to iteratively find an eigenpar of `L`, starting the the iteration from `z`.

# Arguments
- `L::LinearOperatorFamily`: Definition of the nonlinear eigenvalue problem.
- `z`: Initial guess for the eigenvalue.
- `maxiter::Integer=10`: Maximum number of iterations.
- `tol=0`: Absolute tolerance to trigger the stopping of the iteration. If the difference of two consecutive iterates is `abs(z0-z1)<tol` the iteration is aborted.
- `relax=1`: relaxation parameter
- `lam_tol=Inf`: tolerance for the auxiliary eigenvalue to test convergence. The default is infinity, so there is effectively no test.
- `order::Integer=1`: Order of the Householder method, Maximum is 5
- `nev::Integer=1`: Number of Eigenvalues to be searched for in intermediate ARPACK calls.
- `v0::Vector`: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with ones.
- `v0_adj::Vector`: Initial vector for Krylov subspace generation in ARPACK calls. If not provided the vector is initialized with `v0`.
- `output::Bool`: Toggle printing online information.

# Returns
- `sol::Solution`
- `n::Integer`: Number of perforemed iterations
- `flag::Integer`: flag reporting the success of the method. Most important values are `1`: method converged, `0`: convergence might be slow, `-1`:maximum number of iteration has been reached. For other error codes see the source code.

# Notes
Housholder methods are a generalization of Newton's method. If order=1 the Housholder method is identical to Newton's method. The solver then reduces to the "generalized Rayleigh Quotient iteration" presented in [1]. If no relaxation is used (`relax == 1`), the convergence rate is of `order+1`. With relaxation (`relax != 1`) the convergence rate is 1.
Thus, with a higher order less iterations will be necessary. However, the computational time must not necessarily improve nor do the convergence properties. Anyway, if the method converges, the error in the eigenvalue is bounded above by `tol`. For more details on the solver, see the thesis [2].

# References
[1] P. Lancaster, A Generalised Rayleigh Quotient Iteration for Lambda-Matrices,Arch. Rational Mech Anal., 1961, 8, p. 309-322, https://doi.org/10.1007/BF00277446

[2] G.A. Mensah, Efficient Computation of Thermoacoustic Modes, Ph.D. Thesis, TU Berlin, 2019

See also: [`beyn`](@ref)
"""
function householder(L,z;maxiter=10,tol=0.,relax=1.,lam_tol=Inf,order=1,nev=1,v0=[],v0_adj=[],output=true)
    if output
        println("Launching Householder...")
        println("Iter    Res:     dz:     z:")
        println("----------------------------------")
    end
    z0=complex(Inf)
    lam=float(Inf)

    n=0
    active=L.active #IDEA: better integrate activity into householder
    mode=L.mode #IDEA: same for mode
    if v0==[]
        v0=ones(ComplexF64,size(L(0))[1])
    end
    if v0_adj==[]
        v0_adj=conj.(v0)#ones(ComplexF64,size(L(0))[1])
    end

    flag=1
    M=-L.terms[end].coeff
    try
        while  abs(z-z0)>tol   && n<maxiter #&& abs(lam)>lam_tol
            if output; println(n,"\t\t",abs(lam),"\t",abs(z-z0),"\t", z );flush(stdout); end
            z0=z
            L.params[L.eigval]=z
            L.params[L.auxval]=0
            A=L(z)
            lam,v = Arpack.eigs(A,M,nev=nev, sigma = 0,v0=v0)
            lam_adj,v_adj = Arpack.eigs(A',M', nev=nev, sigma = 0,v0=v0_adj)
            #TODO: consider constdouton
            #TODO: multiple eigenvalues
            indexing=sortperm(lam, by=abs)
            lam=lam[indexing]
            v=v[:,indexing]
            indexing=sortperm(lam_adj, by=abs)
            lam_adj=lam_adj[indexing]
            v_adj=v_adj[:,indexing]
            delta_z =[]
            L.active=[L.auxval,L.eigval]
            for i in 1:nev
                L.params[L.auxval]=lam[i]
                sol=Solution(L.params,v[:,i],v_adj[:,i],L.auxval)
                perturb!(sol,L,L.eigval,order,mode=:householder)
                f=[factorial(idx-1)*coeff for (idx,coeff) in enumerate(sol.eigval_pert[Symbol("$(string(L.eigval))/Taylor")])]
                dz=householder_update(f) #TODO: implement multiplicity
                #print(z+dz)
                push!(delta_z,dz)
            end
            L.active=[L.eigval]
            indexing=sortperm(delta_z, by=abs)
            lam=lam[indexing[1]]
            L.params[L.auxval]=lam #TODO remove this from the loop body
            z=z+relax*delta_z[indexing[1]]
            v0=(1-relax)*v0+relax*v[:,indexing[1]]
            v0_adj=(1-relax)*v0_adj+relax*v_adj[:,indexing[1]]
            n+=1
        end

    catch excp
        if output
            println("Error occured:")
            println(excp)
            println(typeof(excp))
            println("...aborted Householder!")
        end

        flag=-2
        if typeof(excp) <: Arpack.ARPACKException
           flag=-4
           if excp==Arpack.ARPACKException(-9999)
               flag=-9999
           end
       elseif excp== LinearAlgebra.SingularException(0) #This means that the solution is so good that L(z) cannot be LU factorized...TODO: implement other strategy into perturb
           flag=-6
           L.params[L.eigval]=z
        end
    end
    if flag==1
        L.params[L.eigval]=z
        if output;  println(n,"\t\t",abs(lam),"\t",abs(z-z0),"\t", z ); end

        if n>=maxiter
            flag=-1
            if output; println("Warning: Maximum number of iterations has been reached!");end

        elseif abs(lam)<=lam_tol
            flag=1
            if output; println("Solution has converged!"); end
        elseif abs(z-z0)<=tol
            flag=0
            if output; println("Warning: Slow convergence!"); end
        elseif isnan(z)
            flag=-5
            if output; println("Warning: computer arithmetics problem. Eigenvalue is NaN"); end
        else
            if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
            flag=-3
            println(z)
        end

        if output
            println("...finished Householder!")
            println("#####################")
            println(" Householder results ")
            println("#####################")
            println("Number of steps: ",n)
            println("Last step parameter variation:",abs(z0-z))
            println("Auxiliary eigenvalue λ residual (rhs):", abs(lam))
            println("Eigenvalue:",z)
            println("Eigenvalue/(2*pi):",z/2/pi)
        end
    end
    L.active=active
    L.mode=mode
    #normalization
    v0/=sqrt(v0'*M*v0)
    v0_adj/=conj(v0_adj'*L(L.params[L.eigval],1)*v0)
    return Solution(L.params,v0,v0_adj,L.eigval), n, flag
end

##PAde solver
function poly_roots(p)
  N=length(p)-1
  C=zeros(ComplexF64,N,N)
  for i=2:N
    C[i,i-1]=1
  end
  C[:,N].=-p[1:N]./p[N+1]
  return LinearAlgebra.eigvals(C)
end

function padesolve(L,z;maxiter=10,tol=0.,relax=1.,lam_tol=Inf,order=1,nev=1,v0=[],v0_adj=[],output=true,num_order=1)
    if output
        println("Launching Pade solver...")
        println("Iter    Res:     dz:     z:")
        println("----------------------------------")
    end
    z0=complex(Inf)
    lam=float(Inf)
    lam0=float(Inf)

    n=0
    active=L.active #IDEA: better integrate activity into householder
    mode=L.mode #IDEA: same for mode
    if v0==[]
        v0=ones(ComplexF64,size(L(0))[1])
    end
    if v0_adj==[]
        v0_adj=conj.(v0)#ones(ComplexF64,size(L(0))[1])
    end

    flag=1
    M=-L.terms[end].coeff
    try
        while  abs(z-z0)>tol   && n<maxiter #&& abs(lam)>lam_tol
            if output; println(n,"\t\t",abs(lam),"\t",abs(z-z0),"\t", z );flush(stdout); end

            L.params[L.eigval]=z
            L.params[L.auxval]=0
            A=L(z)
            lam,v = Arpack.eigs(A,M,nev=nev, sigma = 0,v0=v0)
            lam_adj,v_adj = Arpack.eigs(A',M', nev=nev, sigma = 0,v0=v0_adj)
            #TODO: consider constdouton
            #TODO: multiple eigenvalues
            indexing=sortperm(lam, by=abs)
            lam=lam[indexing]
            v=v[:,indexing]
            indexing=sortperm(lam_adj, by=abs)
            lam_adj=lam_adj[indexing]
            v_adj=v_adj[:,indexing]
            delta_z =[]
            back_delta_z=[]
            L.active=[L.auxval,L.eigval]
            #println("#############")
            #println("lam0:$lam0 ")
            for i in 1:nev
                L.params[L.auxval]=lam[i]
                sol=Solution(L.params,v[:,i],v_adj[:,i],L.auxval)
                perturb!(sol,L,L.eigval,order,mode=:householder)
                coeffs=sol.eigval_pert[Symbol("$(L.eigval)/Taylor")]
                num,den=pade(coeffs,num_order,order-num_order)

                #forward calculation (has issues with multi-valuedness)
                roots=poly_roots(num)
                indexing=sortperm(roots, by=abs)
                dz=roots[indexing[1]]
                push!(delta_z,dz)
                #poles=sort(poly_roots(den),by=abs)
                #poles=poles[1]
                #println(">>>$i<<<")
                #println("$coeffs")
                #println("residue:$(LinearAlgebra.norm((A-lam[i]*M)*v[:,i])) and $(LinearAlgebra.norm((A'-lam_adj[i]*M')*v_adj[:,i]))")
                #println("poles: $(poles+z) r:$(abs(poles))")
                #backward check(for solving multi-valuedness problems by continuity)
                if z0!=Inf
                    back_lam=polyval(num,z0-z)/polyval(den,z0-z)
                    #println("back:$back_lam")
                    #println("root: $(z+dz)")
                    #estm_lam=polyval(num,dz)/polyval(den,dz)
                    #println("estm. lam: $estm_lam")
                    back_lam=lam0-back_lam
                    push!(back_delta_z,back_lam)
                end
            end
            L.active=[L.eigval]
            if z0!=Inf
                indexing=sortperm(back_delta_z, by=abs)
            else
                indexing=sortperm(delta_z, by=abs)
            end
            lam=lam[indexing[1]]
            L.params[L.auxval]=lam #TODO remove this from the loop body
            z0=z
            lam0=lam
            z=z+relax*delta_z[indexing[1]]
            v0=(1-relax)*v0+relax*v[:,indexing[1]]
            v0_adj=(1-relax)*v0_adj+relax*v_adj[:,indexing[1]]
            n+=1
        end

    catch excp
        if output
            println("Error occured:")
            println(excp)
            println(typeof(excp))
            println("...aborted Householder!")
        end

        flag=-2
        if typeof(excp) <: Arpack.ARPACKException
           flag=-4
           if excp==Arpack.ARPACKException(-9999)
               flag=-9999
           end
       elseif excp== LinearAlgebra.SingularException(0) #This means that the solution is so good that L(z) cannot be LU factorized...TODO: implement other strategy into perturb
           flag=-6
           L.params[L.eigval]=z
        end
    end
    if flag==1
        L.params[L.eigval]=z
        if output;  println(n,"\t\t",abs(lam),"\t",abs(z-z0),"\t", z ); end

        if n>=maxiter
            flag=-1
            if output; println("Warning: Maximum number of iterations has been reached!");end

        elseif abs(lam)<=lam_tol
            flag=1
            if output; println("Solution has converged!"); end
        elseif abs(z-z0)<=tol
            flag=0
            if output; println("Warning: Slow convergence!"); end
        elseif isnan(z)
            flag=-5
            if output; println("Warning: computer arithmetics problem. Eigenvalue is NaN"); end
        else
            if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
            flag=-3
            println(z)
        end

        if output
            println("...finished Householder!")
            println("#####################")
            println(" Householder results ")
            println("#####################")
            println("Number of steps: ",n)
            println("Last step parameter variation:",abs(z0-z))
            println("Auxiliary eigenvalue λ residual (rhs):", abs(lam))
            println("Eigenvalue:",z)
            println("Eigenvalue/(2*pi):",z/2/pi)
        end
    end
    L.active=active
    L.mode=mode
    #normalization
    v0/=sqrt(v0'*M*v0)
    v0_adj/=conj(v0_adj'*L(L.params[L.eigval],1)*v0)
    return Solution(L.params,v0,v0_adj,L.eigval), n, flag
end


text:
module BenchCollect

using BenchmarkTools
using Transducers

suite = BenchmarkGroup()

suite["identity-float"] =
    @benchmarkable(collect(Map(identity), xs), setup = (xs = randn(10_000)))

suite["identity-union"] = @benchmarkable(
    collect(Map(identity), xs),
    setup = (xs = [abs(x) < 1 ? x : missing for x in randn(10_000)]),
)
suite["filter-missing"] = @benchmarkable(
    collect(Map(!ismissing), xs),
    setup = (xs = [abs(x) < 1 ? x : missing for x in randn(10_000)]),
)

end  # module
BenchCollect.suite


text:
module TestParseUInts

import GraphIdx.Io: parse_uints
using Test

function record_ints(io::IO)
    args = Vector{Tuple{UInt,Bool}}()
    parse_uints(io) do u::UInt, b::Bool
        push!(args, (u, b))
    end
    return args
end

@testset "parse_uints                    " begin
    @test record_ints(IOBuffer("")) == []
    @test record_ints(IOBuffer("0")) == [(0, true)]
    @test record_ints(IOBuffer("0120 3")) == [(120, false), (3, true)]
    @test record_ints(IOBuffer("0120 \n 3")) == [(120, true), (3, true)]
    @test record_ints(IOBuffer(", 0120 3")) == [(120, false), (3, true)]
    @test record_ints(IOBuffer(", 0120 , 3")) == [(120, false), (3, true)]
    @test record_ints(IOBuffer(", 0120 ,\n 3")) == [(120, true), (3, true)]
    @test record_ints(IOBuffer(", 0120 ,\n 3_000")) == [(120, true), (3000, true)]
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
"""
    `streakpoke(r::DataFrames.DataFrameRow)`
    Given poke information as a DataFrame row compose a box whose y position is
    equal to the trial number, and its x extension is given by the poke duration
"""
function streakpoke(r::DataFrames.DataFrameRow)
    y_start = r.Streak - 0.25
    y_finish = r.Streak + 0.25
    x_start = r.In
    x_finish = r.Out
    Shape([(x_start,y_start),(x_start,y_finish),(x_finish,y_finish),(x_finish, y_start),(x_start,y_start)])
end

function outcome_col(r::DataFrames.DataFrameRow)
    if r.Reward
        fcol = :green
    elseif !r.Reward
        fcol = :grey
    end
    if r.Side != r.SideHigh
        fcol = :white
    end
    return fcol
end

function side_annotate!(p,r::DataFrames.DataFrameRow)
    if typeof(r.Side) <: Real
        side =  r.Side == 0 ? "R" : "L"
    else
        side = r.Side
    end
    annotate!(p,[(-0.5, r.Streak, Plots.text(side, 4, :center))])
end

"""
    `session_plot!(p,r::DataFrames.DataFrameRow)`
     Given poke information as a DataFrame row, plots a box on an existing
     figure with its positioned and color given by the trial number, poke time
     and outcome.
"""
function session_plot!(p,r::DataFrames.DataFrameRow)
    bp = streakpoke(r)
    fcol = outcome_col(r)
    plot!(p,bp, fillcolor = plot_color(fcol, 0.6),linecolor = :black,linewidth = 1)
    side_annotate!(p,r)
    p
end

"""
    `session_plot(df::DataFrames.AbstractDataFrame)`
     Given DataFrame containing all pokes in a session loops each row, to plot
     a box per each box.
"""

function session_plot(df::DataFrames.AbstractDataFrame)
    xyfont = font(18, "Bookman Light")
    legfont = font(14, "Bookman Light")
    xprop = ("Time (seconds)", (-1,16), xyfont)
    yprop = ("Trials", xyfont,false, false)
    plt = plot(legend = false, xaxis = xprop, yaxis = yprop, yticks = false)
    for r in eachrow(df)
        FLPDevelopment.session_plot!(plt, r)
    end
    plt
end


text:
# Visualization
function visualize!(vis, model::Hopper, q; Δt = 0.1, fixed_camera=true)

    r_body = model.body_radius
    r_foot = model.foot_radius
    r_leg = 0.5 * r_foot

    default_background!(vis)

    setobject!(vis["body"], GeometryBasics.Sphere(Point3f0(0),
        convert(Float32, model.body_radius)),
        MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    setobject!(vis["foot"], GeometryBasics.Sphere(Point3f0(0),
        convert(Float32, model.foot_radius)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))

    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics_foot(model, q[t])[1], 0.0, kinematics_foot(model, q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics_foot(model, q_tmp)[1], 0.0, kinematics_foot(model, q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics_foot(model, q_tmp)[1], 0.0, kinematics_foot(model, q_tmp)[2]]

        z_shift = [0.0; 0.0; 0.0]

        MeshCat.atframe(anim, t) do
            MeshCat.settransform!(vis["body"], MeshCat.compose(MeshCat.Translation(p_body + z_shift), MeshCat.LinearMap(Rotations.RotY(0.0))))
            MeshCat.settransform!(vis["foot"], MeshCat.compose(MeshCat.Translation(p_foot + z_shift), MeshCat.LinearMap(Rotations.RotY(0.0))))

            for i = 1:n_leg
                MeshCat.settransform!(vis["leg$i"], MeshCat.compose(MeshCat.Translation(p_leg[i] + z_shift), MeshCat.LinearMap(Rotations.RotY(0.0))))
            end
        end
    end

    if fixed_camera
        MeshCat.settransform!(vis["/Cameras/default"],
        MeshCat.compose(MeshCat.Translation(0.0, -50.0, -1.0), MeshCat.LinearMap(Rotations.RotZ(-pi / 2.0))))
        setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)
    end

    MeshCat.setanimation!(vis, anim)

    return true
end


text:
#!/usr/local/bin/julia
#
# メカ女子将棋システム (C) 2017 メカ女子将棋部☆
#
# メカ女子将棋部：
#
# 竹部　さゆり（女流三段）様
#   メカさゆりん
# 渡辺　弥生（女流初段）様
#   メカみおたん
# 酒井　美由紀（裁縫顧問）様
# 　メカみゅーん
# 辻　理絵子（女子大学院生）様
#   メカりえぽん
# 木村　健（メカウーサーメカ担当、実装責任者）
#   メカきむりん、プロジェクトリーダー、小六女子w
#
#
# 使い方
# cd Julia
# ./Main.jl stdio           # 標準入出力を使う場合
# ./Main.jl tcp 192.168.0.1 # tcp思考サーバを立ち上げる場合
#

const MECHA_JYOSHI_SHOGI = 1
const MECHAJYO_VERSION = "Monkey Mechajyo"

const THINK_TIME_IN_NS = UInt64(1000000000*14)
const THINK_BYOYOMI_IN_NS = UInt64(1000000000*8)
const BEGIN_BYOYOMI_IN_NS = UInt64(1000000000*60)

srand(1234)

#using IProfile

include("BoardConsts.jl")
include("Board.jl")
include("Move.jl")
include("GameStatus.jl")
include("OldGenMove.jl")
include("BitBoard.jl")
include("fvbin.jl")
include("Eval.jl")
include("EvalBonanza.jl")
include("GenMove.jl")
include("Search.jl")
include("PVS.jl")

# おまじない
function setupIO()
    gs = InitGS()
    InitTables(gs)
    #BBTest2(gs)
    #BBTestForMoveGeneration(gs)
    return gs
end

function producer()
    if length(ARGS) == 1
        if ARGS[1] == "stdio"
            sock  = STDOUT
            stdin = STDIN
            ret = main(stdin,sock)
        end
    end
    if (length(ARGS) == 2) && (ARGS[1] == "tcp")
        addr = ARGS[2]
        println("establish server (", addr, ") port: ", "4091")
        begin
            #server = listen(getaddrinfo("127.0.0.1"), parse(UInt32,"4091"))
            server = listen(getaddrinfo(addr), parse(UInt32,"4091"))
            while true
                println("waiting for connection...")
                sock  = accept(server)
                println("establish connection!")
                while true
                    ret = main(sock,sock)
                    break
                end
                close(sock)
                break
            end
        end
    end
end

function main(stdin, sock)
    sfenHirate = "lnsgkgsnl/1r5b1/ppppppppp/9/9/9/PPPPPPPPP/1B5R1/LNSGKGSNL w - 1"
    li::Array{String,1} = ["" for x = 1:100]::Array{String,1}
    li2::Array{String,1} = ["" for x = 1:100]::Array{String,1}
    lilen = 0
    lilen2 = 0
    side::Int = SENTE # default
    count::Int = 0
    count2::Int = 0
    gs = setupIO()

    out = [Move(0,0,0,0,0,0) for n = 1:300000]
    history = [0 for x = 1:NumSQ, y = 1:NumSQ]
    #GenMoveTest(gs)
    ###t = Task(producer)
    while true
        st = readline(stdin) #consume(t)
        # st = chomp(st)
        # println(sock,"consumed $st")
        if st == "quit" || st == "exit"
            break
        elseif st == "usi"
            println(sock,"id name Mecha Jyoshi Shogi ",MECHAJYO_VERSION)
            println(sock,"id author Sayuri TAKEBE, Mio WATANABE, Miyuki SAKAI, Rieko TSUJI and Takeshi KIMURA")
            println(sock,"option name BookFile type string default $(gs.bookfile)")
            println(sock,"option name UseBook type check default $(gs.usebook)")
            println(sock,"usiok");
        elseif st == "isready"
            srand(time_ns())
            println(sock,"readyok")
        elseif startswith(st,"setoption")
            if startswith(st,"setoption name USI_Ponder value true")
                gs.canponder = true
            elseif startswith(st,"setoption name USI_Ponder value false")
                gs.canponder = false
            elseif startswith(st,"setoption name USI_Hash value ")
                gs.hashsize = parse(UInt32,st[length("setoption name USI_Hash value "):end])
            elseif startswith(st,"setoption name BookFile value ")
                gs.bookfile = st[length("setoption name BookFile value "):end]
            else
            end
        elseif startswith(st,"usinewgame")
            # do nothing
        elseif startswith(st,"position startpos")
            li = split(st)
            count = 0
            count2 = 0
            if li[2] == "startpos" # 平手
                if size(li,1) == 2 # 初手
                    gs.board = InitSFEN(sfenHirate, gs.board)
                    count = 0
                    side = SENTE
                else
                    if size(li,1) == 4
                        side = GOTE
                    end
                    gs.board = InitSFEN(sfenHirate, gs.board)
                    count = 0

                    for x = 4:size(li,1)
                        count = 0
                        if iseven(x) # 先手
                            count2 = generateMoves(gs.board, out, SENTE, count, gs)
                            idx = findIndex(out,li[x],SENTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println(sock,"Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println(sock,"Warning: invalid moves:", li[x])
                            else
                                #println(sock,"index = ", idx)
                                makeMove(gs.board,
                                idx,
                                out,
                                SENTE)
                                #println(sock,"move=",move2USIString(out[idx]))
                            end
                        else # 後手
                            count2 = generateMoves(gs.board, out, GOTE, count, gs)
                            idx = findIndex(out,li[x],GOTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println(sock,"Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println(sock,"Warning: invalid moves:", li[x])
                            else
                                #println(sock,"index = ", idx)
                                makeMove(gs.board,
                                idx,
                                out,
                                GOTE)
                                #println(sock,"move=",move2USIString(out[idx]))
                            end
                        end

                    end
                end
            end
            DisplayBoard(gs.board)
        elseif startswith(st,"position sfen")
            li = split(st)
            count = 0
            count2 = 0
            sfen = join(li[3:6]," ")
            #println(sock,"sfen=",sfen)
            if li[2] == "sfen" # 平手
                if size(li,1) == 6 # 初手
                    gs.board = InitSFEN(sfen, gs.board)
                    count = 0
                    side = SENTE
                else
                    if size(li,1) == 8
                        side = GOTE
                    end
                    gs.board = InitSFEN(sfen, gs.board)
                    count = 0

                    for x = 8:size(li,1)
                        count = 0
                        if iseven(x) # 先手
                            count2 = generateMoves(gs.board, out, SENTE, count, gs)
                            idx = findIndex(out,li[x],SENTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println(sock,"Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                printlpn("Warning: invalid moves:", li[x])
                            else
                                #println(sock,"index = ", idx)
                                makeMove(gs.board,
                                idx,
                                out,
                                SENTE)
                                #println(sock,"move=",move2USIString(out[idx]))
                            end
                        else # 後手
                            count2 = generateMoves(gs.board, out, GOTE, count, gs)
                            idx = findIndex(out,li[x],GOTE,count+1,count2)
                            if idx == -1
                                idx = 1
                                println(sock,"Warning: cannot find move:", li[x])
                            elseif idx == 0 # resign
                                # do nothing
                                println(sock,"Warning: invalid moves:", li[x])
                            else
                                #println(sock,"index = ", idx)
                                makeMove(gs.board,
                                idx,
                                out,
                                GOTE)
                                #println(sock,"move=",move2USIString(out[idx]))
                            end
                        end
                    end
                end
            end
            DisplayBoard(gs.board)
        elseif startswith(st,"go")
            li2 = split(st)
            start, goal, isByoyomi, thinkTime = parseGo(li2,side)
            count2 = generateMoves(gs.board, out, side, 0, gs)
            count3 = generateBB(gs.board, out, side, 0, gs)
            if count2 == 0
                println(sock,"bestmove resign")
            else
                # chose random moves
                #randomIndex::Int = rand(UInt32) % (count2)
                Index::Int = -1
                #m::Move = think(side,gs)
                gs.maxThinkingTime = thinkTime
                m::Move = thinkASP(side,gs,sock)
                #println(sock,"move=",m)
                for q = 1:count2
                    if out[q].move == m.move
                        Index = q
                        break
                    end
                end
                if Index == -1
                    println(sock,"bestmove resign")
                else
                    makeMove(gs.board,Index,out,side)
                    if in_check( side, gs.board)
                        # println("check!")
                    end
                    println(sock,"bestmove ",move2USIString(out[Index]))
                    DisplayBoard(gs.board)
                end
            end
        elseif startswith(st,"gameover")
            # do nothing
        else
            println("debug COMMAND NOT FOUND($st)")
        end
        #println(sock,"COMMAND: $st")
    end
    return "quit"
end

function parseGo(list,side)
    i = 1
    btime = 0
    wtime = 0
    byoyomi = 0
    winc = 0
    binc = 0
    for x in list
        if btime == -1
            btime = parse(UInt64,x*"000000")
        end
        if wtime == -1
            wtime = parse(UInt64,x*"000000")
        end
        if byoyomi == -1
            byoyomi = parse(UInt64,x*"000000")
        end
        if winc == -1
            winc = parse(UInt64,x*"000000")
        end
        if binc == -1
            binc = parse(UInt64,x*"000000")
        end

        if x == "btime"
            btime = -1
        elseif x == "wtime"
            wtime = -1
        elseif x == "byoyomi"
            byoyomi = -1
        elseif x == "winc"
            winc = -1
        elseif x == "binc"
            binc = -1
        end
        i = i + 1
    end

    isByoyomi = (side == SENTE) ? (btime < BEGIN_BYOYOMI_IN_NS) : (wtime < BEGIN_BYOYOMI_IN_NS)

    remainTime = (side == SENTE) ? btime : wtime

    println("remain time = ", remainTime)

    minByoyomi = 2 # 2 secs
    byoyomiThinkingTime = 10
    thinkTime = 0
    if isByoyomi
        if (winc > 0)&&(binc > 0) # fisher rule
            thinkTime = Int64(1000000000 * byoyomiThinkingTime)
        else
            if byoyomi < Int64(1000000000 * minByoyomi)
                thinkTime = Int64(UInt64(byoyomi >>> 1))
            else
                thinkTime = byoyomi - Int64(1000000000 * minByoyomi)
            end
        end
    else
        thinkTime = THINK_TIME_IN_NS
    end
    start = time_ns() # sampling time value
    goal = start + thinkTime
    tObj = (start, goal, isByoyomi, thinkTime)
    # println(tObj)
    return tObj
end

producer()
#@iprofile report

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
import FFTW.Frequencies

"""The Mel spectrogram"""
struct MelSpectrogram{T, F <: Union{Frequencies,AbstractRange}} <: DSP.Periodograms.TFR{T}
    power::Matrix{T}
    mels::F
    time::StepRangeLen{Float64}
end
heatmap(tfr::MelSpectrogram) = log10(tfr.power)
DSP.freq(tfr::MelSpectrogram) = tfr.mels
Base.time(tfr::MelSpectrogram) = tfr.time

"""A TFR subtype where each column represents an MFCC vector"""
struct MFCC{T, F <: Union{Frequencies,AbstractRange}} <: DSP.Periodograms.TFR{T}
    mfcc::Matrix{T}
    number::F
    time::StepRangeLen{Float64}
end
heatmap(tfr::MFCC) = tfr.mfcc
DSP.freq(tfr::MFCC) = tfr.number
Base.time(tfr::MFCC) = tfr.time

""""""
function hz_to_mel(frequencies::Union{F, AbstractArray{F}}) where {F <: Real}
    f_min = 0f0
    f_sp = 200f0 / 3

    mels = collect((frequencies - f_min) / f_sp)

    min_log_hz = 1000f0
    min_log_mel = (min_log_hz - f_min) / f_sp
    logstep = log(6.4f0) / 27f0

    @inbounds for i = 1:length(mels)
        if frequencies[i] >= min_log_hz
            mels[i] = min_log_mel + log(frequencies[i] / min_log_hz) / logstep
        end
    end

    mels
end

""""""
function mel_to_hz(mels::Union{F, AbstractArray{F}}) where {F <: Real}
    f_min = 0f0
    f_sp = 200f0 / 3
    frequencies = collect(f_min + f_sp * mels)

    min_log_hz = 1000f0
    min_log_mel = (min_log_hz - f_min) / f_sp
    logstep = log(6.4f0) / 27f0

    @inbounds for i = 1:length(frequencies)
        if mels[i] >= min_log_mel
            frequencies[i] = min_log_hz * exp(logstep * (mels[i] - min_log_mel))
        end
    end

    frequencies
end

""""""
function mel_frequencies(nmels::Int = 128, fmin::Real = 0.0f0, fmax::Real = 11025f0)
    min_mel = hz_to_mel(fmin)[1]
    max_mel = hz_to_mel(fmax)[1]

    mels = linspace(min_mel, max_mel, nmels)
    mel_to_hz(mels)
end

""""""
function mel(samplerate::Real, nfft::Int, nmels::Int = 128, fmin::Real = 0f0, fmax::Real = samplerate/2f0)
    weights = zeros(Float32, nmels, (nfft >> 1) + 1)
    fftfreqs = fft_frequencies(samplerate, nfft)
    melfreqs = mel_frequencies(nmels + 2, fmin, fmax)
    enorm = 2f0 ./ (melfreqs[3:end] - melfreqs[1:nmels])

    for i in 1:nmels
        lower = (fftfreqs - melfreqs[i]) / (melfreqs[i+1] - melfreqs[i])
        upper = (melfreqs[i+2] - fftfreqs) / (melfreqs[i+2] - melfreqs[i+1])

        weights[i, :] = max(0, min(lower, upper)) * enorm[i]
    end

    weights
end

function melspectrogram(audio::SampleBuf{T, 1}, windowsize::Int = 1024, hopsize::Int = windowsize >> 2;
                        nmels::Int = 128, fmin::Real = 0f0, fmax::Real = audio.samplerate.val / 2f0) where T
    samplerate = audio.samplerate.val
    nfft = DSP.nextfastfft(windowsize)
    S = spectrogram(audio, windowsize, hopsize).power
    data = mel(samplerate, nfft, nmels, fmin, fmax) * S
    nframes = size(data, 2)
    MelSpectrogram(data, linspace(hz_to_mel(fmin)[1], hz_to_mel(fmax)[1], nmels), (0.0:nframes-1) * hopsize / samplerate)
end

function mfcc(audio::SampleBuf{T, 1}, windowsize::Int = 1024, hopsize::Int = windowsize >> 2;
              nmfcc::Int = 20, nmels::Int = 128, fmin::Real = 0f0, fmax::Real = audio.samplerate.val / 2f0) where T
    if nmfcc >= nmels
        error("number of mfcc components should be less than the number of mel frequency bins")
    end

    M = melspectrogram(audio, windowsize, hopsize; nmels = nmels, fmin = fmin, fmax = fmax)
    mfcc = dct(nmfcc, nmels) * power(M)

    for frame in 1:size(mfcc, 2)
        mfcc[:, frame] /= norm(mfcc[:, frame])
    end

    MFCC(mfcc, 1:nmfcc, time(M))
end


text:
using Base.Test
# spec2d.jl
import WaveletScattering: Spec2D, checkspec,
    default_max_aspectratio, default_n_orientations
# waveletclass.jl
import WaveletScattering: MexicanHat, Morlet

# Spec2D constructor
spec = Spec2D()
@test isa(spec, Spec2D)

# checkspec
@test_throws ErrorException checkspec(
    Spec2D(class=Morlet(), max_aspectratio = 0.5))
@test_throws ErrorException checkspec(
    Spec2D(class=Morlet(), n_orientations=1))
@test_throws ErrorException checkspec(
    Spec2D(class=MexicanHat(), max_aspectratio=2.0))
@test_throws ErrorException checkspec(
    Spec2D(class=MexicanHat(), max_aspectratio=0.5))
@test_throws ErrorException checkspec(
    Spec2D(class=MexicanHat(), n_orientations=2))
@test_throws ErrorException checkspec(
    Spec2D(class=MexicanHat(), n_orientations=0))

# default_max_aspectratio
@test default_max_aspectratio(Morlet(), nothing) ≈ 2.0
@test default_max_aspectratio(MexicanHat(), nothing) ≈ 1.0
@test default_max_aspectratio(Morlet(), 3.0) ≈ 3.0

# default_n_orientations
@test default_n_orientations(Morlet(), nothing) == 4
@test default_n_orientations(MexicanHat(), nothing) == 1
@test default_n_orientations(Morlet(), 3) == 3


text:
include("SV.jl")
posterior, chain = main()
# Extract the parameter and plot the posterior
σe_hat = [i[1][1] for i in posterior]
#σu_hat = 2.0*log.(σu_hat)
dstats(σe_hat)
p1 = npdensity(σe_hat)
# Extract the parameter and plot the posterior
ρhat = [i[2][1] for i in posterior]
dstats(ρhat)
p2 = npdensity(ρhat)
# Extract the parameter and plot the posterior
σu_hat = [i[3][1] for i in posterior]
dstats(σu_hat)
p3 = npdensity(σu_hat)
plot(p1, p2, p3)
plot!(legend=false)
gui()
#savefig("nuts.svg")

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
"""
    IndependentMOGP(fs)

A multi-output GP with independent outputs where output `i` is modelled by the
single-output GP fs[i].

# Arguments:
- fs: a vector of `p` single-output GPs where `p` is the dimension of the output.
"""
struct IndependentMOGP{Tfs<:Vector{<:AbstractGP}} <: AbstractGP
    fs::Tfs
end

"""
    independent_mogp(fs)

Returns an IndependentMOGP given a list of single output GPs `fs`.

```jldoctest
julia> ind_mogp1 = independent_mogp([GP(KernelFunctions.SEKernel())]);

julia> ind_mogp2 = IndependentMOGP([GP(KernelFunctions.SEKernel())]);

julia> typeof(ind_mogp1) == typeof(ind_mogp2)
true

julia> ind_mogp1.fs == ind_mogp2.fs
true
```
"""
independent_mogp(fs::Vector{<:AbstractGP}) = IndependentMOGP(fs)

"""
    finite_gps(fx)

Returns a list of of the finite GPs for all latent processes, given a finite
IndependentMOGP and *isotopic inputs*.
"""
function finite_gps(fx::FiniteGP{<:IndependentMOGP, <:MOInputIsotopicByOutputs}, σ²::Real)
    return [f(fx.x.x, σ²) for f in fx.f.fs]
end

const IsotropicByOutputsFiniteIndependentMOGP = FiniteGP{
    <:IndependentMOGP,<:MOInputIsotopicByOutputs,<:Diagonal{<:Real,<:Fill}
}

# Optimisations for MOInputIsotopicByOutputs.

# See AbstractGPs.jl API docs.
function AbstractGPs.mean(f::IndependentMOGP, x::MOInputIsotopicByOutputs)
    return reduce(vcat, map(f -> mean(f, x.x), f.fs))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.var(f::IndependentMOGP, x::MOInputIsotopicByOutputs)
    return reduce(vcat, map(f -> var(f, x.x), f.fs))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.cov(f::IndependentMOGP, x::MOInputIsotopicByOutputs)
    Cs = map(f -> cov(f, x.x), f.fs)
    return Matrix(BlockDiagonal(Cs))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.cov(
    f::IndependentMOGP, x::MOInputIsotopicByOutputs, y::MOInputIsotopicByOutputs
)
    Cs = map(f -> cov(f, x.x, y.x), f.fs)
    return Matrix(BlockDiagonal(Cs))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.logpdf(
    ft::IsotropicByOutputsFiniteIndependentMOGP, y::AbstractVector{<:Real}
)
    finiteGPs = finite_gps(ft, ft.Σy[1])
    ys = collect(eachcol(reshape(y, (length(ft.x.x), :))))
    return sum(map(logpdf, finiteGPs, ys))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.rand(rng::AbstractRNG, ft::IsotropicByOutputsFiniteIndependentMOGP)
    finiteGPs = finite_gps(ft, ft.Σy[1])
    return reduce(vcat, map(fx -> rand(rng, fx), finiteGPs))
end

# See AbstractGPs.jl API docs.
AbstractGPs.rand(ft::FiniteGP{<:IndependentMOGP}) = rand(Random.GLOBAL_RNG, ft)

# See AbstractGPs.jl API docs.
function AbstractGPs.rand(
    rng::AbstractRNG, ft::IsotropicByOutputsFiniteIndependentMOGP, N::Int
)
    return reduce(hcat, [rand(rng, ft) for _ in 1:N])
end

# See AbstractGPs.jl API docs.
AbstractGPs.rand(ft::FiniteGP{<:IndependentMOGP}, N::Int) = rand(Random.GLOBAL_RNG, ft, N)

# See AbstractGPs.jl API docs.
function Distributions._rand!(
    rng::AbstractRNG,
    fx::FiniteGP{<:IndependentMOGP},
    y::AbstractVecOrMat{<:Real}
)
    N = size(y, 2)
    if N == 1
        y .= AbstractGPs.rand(rng, fx)
    else
        y .= AbstractGPs.rand(rng, fx, N)
    end
end

"""
Posterior implementation for isotopic inputs, given diagonal Σy (OILMM).
See AbstractGPs.jl API docs.
"""
function AbstractGPs.posterior(
    ft::IsotropicByOutputsFiniteIndependentMOGP, y::AbstractVector{<:Real}
)
    finiteGPs = finite_gps(ft, ft.Σy[1])
    ys = collect(eachcol(reshape(y, (length(ft.x.x), :))))
    ind_posts = [AbstractGPs.posterior(fx, y_i) for (fx, y_i) in zip(finiteGPs, ys)]
    return independent_mogp(ind_posts)
end

# AbstractGPs APIs implementations for MOInputIsotopicByFeatures.

const IsotropicByFeaturesFiniteIndependentMOGP = FiniteGP{
    <:IndependentMOGP,<:MOInputIsotopicByFeatures,<:Diagonal{<:Real,<:Fill}
}

# Indices which, when applied to a vector ordered by features, will order it by outputs.
function indices_which_reorder_features_to_outputs(
    x::Union{MOInputIsotopicByOutputs,MOInputIsotopicByFeatures}
)
    return collect(vec(reshape(1:length(x), length(x.x), x.out_dim)'))
end

reorder_by_outputs(x::MOInputIsotopicByFeatures) = MOInputIsotopicByOutputs(x.x, x.out_dim)

function reorder_by_outputs(Σy::Diagonal{<:Real}, x::MOInputIsotopicByFeatures)
    return Diagonal(Σy.diag[indices_which_reorder_features_to_outputs(x)])
end

reorder_by_outputs(Σy::Diagonal{<:Real,<:Fill}, x::MOInputIsotopicByFeatures) = Σy

function reorder_by_outputs(
    fx::FiniteGP{<:IndependentMOGP,<:MOInputIsotopicByFeatures,<:Diagonal{<:Real}}
)
    return FiniteGP(fx.f, reorder_by_outputs(fx.x), reorder_by_outputs(fx.Σy, fx.x))
end

@non_differentiable indices_which_reorder_features_to_outputs(::Any)
@non_differentiable reorder_by_outputs(::Any)

function finite_gps(fx::FiniteGP{<:IndependentMOGP,<:MOInputIsotopicByFeatures}, σ²::Real)
    return [f(fx.x.x, σ²) for f in fx.f.fs]
end

function AbstractGPs.mean(f::IndependentMOGP, x::MOInputIsotopicByFeatures)
    x_by_outputs = reorder_by_outputs(x)
    mean_by_outputs = mean(f, x_by_outputs)
    return mean_by_outputs[indices_which_reorder_features_to_outputs(x_by_outputs)]
end

function AbstractGPs.var(f::IndependentMOGP, x::MOInputIsotopicByFeatures)
    x_by_outputs = reorder_by_outputs(x)
    var_by_outputs = var(f, x_by_outputs)
    return var_by_outputs[indices_which_reorder_features_to_outputs(x_by_outputs)]
end

function AbstractGPs.cov(f::IndependentMOGP, x::MOInputIsotopicByFeatures)
    x_by_outputs = reorder_by_outputs(x)
    C_by_outputs = cov(f, x_by_outputs)
    idx = indices_which_reorder_features_to_outputs(x_by_outputs)
    return C_by_outputs[idx, idx]
end

function AbstractGPs.cov(
    f::IndependentMOGP, x::MOInputIsotopicByFeatures, y::MOInputIsotopicByFeatures
)
    x_by_outputs = reorder_by_outputs(x)
    y_by_outputs = reorder_by_outputs(y)
    C_by_outputs = cov(f, x_by_outputs, y_by_outputs)
    idx_x = indices_which_reorder_features_to_outputs(x_by_outputs)
    idx_y = indices_which_reorder_features_to_outputs(y_by_outputs)
    return C_by_outputs[idx_x, idx_y]
end

function AbstractGPs.cov(
    f::IndependentMOGP, x::MOInputIsotopicByFeatures, y::MOInputIsotopicByOutputs
)
    x_by_outputs = reorder_by_outputs(x)
    C_by_outputs = cov(f, x_by_outputs, y)
    idx_x = indices_which_reorder_features_to_outputs(x_by_outputs)
    return C_by_outputs[idx_x, :]
end

function AbstractGPs.cov(
    f::IndependentMOGP, x::MOInputIsotopicByOutputs, y::MOInputIsotopicByFeatures
)
    y_by_outputs = reorder_by_outputs(y)
    C_by_outputs = cov(f, x, y_by_outputs)
    idx_y = indices_which_reorder_features_to_outputs(y_by_outputs)
    return C_by_outputs[:, idx_y]
end

function AbstractGPs.rand(rng::AbstractRNG, ft::IsotropicByFeaturesFiniteIndependentMOGP)
    finiteGPs = finite_gps(ft, ft.Σy[1])
    return vec(reduce(hcat, map(fx -> rand(rng, fx), finiteGPs))')
end

function AbstractGPs.logpdf(
    ft::FiniteGP{<:IndependentMOGP,<:MOInputIsotopicByFeatures,<:Diagonal{<:Real}},
    y::AbstractVector{<:Real},
)
    return logpdf(
        reorder_by_outputs(ft), y[indices_which_reorder_features_to_outputs(ft.x)]
    )
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries
using Test

@testset "Test inspired by Fateman (takes a few seconds)" begin
    x, y, z, w = set_variables(Int128, "x", numvars=4, order=40)

    function fateman2(degree::Int)
        T = Int128
        oneH = HomogeneousPolynomial([one(T)], 0)
        # s = 1 + x + y + z + w
        s = TaylorN(
            [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree)
        s = s^degree
        # s is converted to order 2*ndeg
        s = TaylorN(s.coeffs, 2*degree)
        return s^2 + s
    end

    function fateman3(degree::Int)
        s = x + y + z + w + 1
        s = s^degree
        s * (s+1)
    end

    f2 = fateman2(20)
    f3 = fateman3(20)
    c = getcoeff(f2,[1,6,7,20])
    @test c == 128358585324486316800
    @test getcoeff(f3,[1,6,7,20]) == c
end


text:
using SatSolver
using Test

include("utils.jl")

@testset "SatSolver.jl" begin
  include("formula.jl")
  include("nnf_parser.jl")
  include("io.jl")
  include("dpll.jl")
end


text:
##############################################################################################################
##############################################################################################################
########################################         MODEL         ###############################################
##############################################################################################################
##############################################################################################################
# abstract type PredictionModel end
abstract type Model end

abstract type PredictionModelParameters end

mutable struct InferenceConfiguration
	convergence_criterion::Float64
	min_iter::Int64
	max_iter::Int64
	target_dir::String
	fold_num::Int64
	do_gridsearch::Bool
	do_cross_validation::Bool
	compute_wpc_index::Bool
	subsume_pathways::Bool
	do_variance_filtering::Bool

	InferenceConfiguration() = new(1e-3, 5, 200, "results/", 1, false, false, false, true, true)
end

mutable struct PredictionModel{M<:Model}
	# model_type::Type{<:Model}
	data::Experiment
	precomputations::Dict{String, Any}

	PredictionModel(::Type{M}, data::Experiment) where {M <: Model} = new{M}(data, Dict{String, Any}())
end

struct ModelConfiguration
	parameters::Dict{String, Any}
	function ModelConfiguration()
		new(Dict{String, Any}())
	end
end


text:
__precompile__()

module App

using JLD
using Images
using ImageMagick

function run()
  include("run.jl")
end

end

function main()
  App.run()
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL