text:
export howell_form2

# we assume that B has at least as many rows as columns
function howell_form2(A::nmod_mat)
  #A = deepcopy(B)
  n = A.r
  m = A.c
  if n < m
    A = vcat(A, zero_matrix(base_ring(A), m-n, m))
  end
  for j in 1:m
    for i in j+1:n
      #print("Computing xgcd of $(_raw_getindex(A,j,j)), $(_raw_getindex(A,i,j))\n")
      g,s,t,u,v = _xxgcd(_raw_getindex(A,j,j),_raw_getindex(A,i,j),A._n)
      #println("$g $s $t $u $v ")
      for k in 1:m
        t1 = s*_raw_getindex(A,j,k) + t*_raw_getindex(A, i, k)
        t2 = u*_raw_getindex(A,j,k) + v*_raw_getindex(A, i, k)
        t1 = ccall((:n_mod2_preinv, :libflint), Int, (Int, Int, UInt), t1, A._n, A._ninv)
        t2 = ccall((:n_mod2_preinv, :libflint), Int, (Int, Int, UInt), t2, A._n, A._ninv)
        ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &A, j - 1, k - 1, t1)
        ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &A, i - 1, k - 1, t2)
      end
    end
    #println(A)
  end
  #println("I have a triangular matrix:")
  #println(A)
  T = zero_matrix(base_ring(A), 1, A.c)
  for j in 1:m
    if _raw_getindex(A,j,j) != 0
      #println("nonzero case")
      u = _unit(_raw_getindex(A,j,j), A._n)
      for k in 1:m
        t1 = mod(u*_raw_getindex(A,j,k), A._n)
        ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &A, j - 1, k - 1, t1)
      end
      #println("Multiplied row $j with $u and obtained \n$A\n")
      for i in 1:j-1
        #println("ich nehme zeile $j und reduziere damit zeile $i")
        #println("zeile $i: $(window(A, i:i, 1:m))")
        #println("zeile $j: $(window(A, j:j, 1:m))")
        #println("A[$i,$j]: $(_raw_getindex(A,i,j)) A[$j,$j]: $(_raw_getindex(A,j,j))")
        #println("A[$i,$j]: $(A[i,j]) A[$j,$j]: $(A[j,j])")
        q = UInt(prem(- Int(div(_raw_getindex(A, i, j), _raw_getindex(A, j, j))), Int(A._n)))
        #println("quotient ist $q")
        for k in j:m
          #println("adjusting row $i, divinding by $(_raw_getindex(A, j, j))")
          #println("I am replacing A[$i,$k] with")
          #println("A[$i,$j]: $(_raw_getindex(A,i,j)) A[$j,$j]: $(_raw_getindex(A,j,j))")
          #println("quotient is $(UInt(prem(- Int(div(_raw_getindex(A, i, j), _raw_getindex(A, j, j))), Int(A._n))))")
          #print(_raw_getindex(A, i, k) + UInt(prem(- Int(div(_raw_getindex(A, i, j), _raw_getindex(A, j, j))), Int(A._n))) *_raw_getindex(A, j, k))
          #println()
          _raw_setindex(A, i, k, mod(_raw_getindex(A, i, k) + q *_raw_getindex(A, j, k), A._n))
        end
      end
      #print("computing the annihilator for $j ... ")
      l = mod(div(A._n,gcd(_raw_getindex(A,j,j),A._n)),A._n)
      #println("annihilator is $l")
      if l == 0
        continue
      end
      for k in 1:m
        _raw_setindex(T, 1, k, mod(l*_raw_getindex(A,j,k), A._n))
      end
    else
      for k in 1:m
        _raw_setindex(T, 1, k, _raw_getindex(A,j,k))
      end
    end
    #println("T looks like: $T")
    #println("Augmented matrix looks like \n$A\n$T\n")
    for i in j+1:m
      g,s,t,u,v = _xxgcd(_raw_getindex(A,i,i),_raw_getindex(T,1,i),A._n)
      for k in 1:m
        t1 = s*_raw_getindex(A,i,k) + t*_raw_getindex(T, 1, k)
        t2 = u*_raw_getindex(A,i,k) + v*_raw_getindex(T, 1, k)
        t1 = ccall((:n_mod2_preinv, :libflint), Int, (Int, Int, UInt), t1, A._n, A._ninv)
        t2 = ccall((:n_mod2_preinv, :libflint), Int, (Int, Int, UInt), t2, A._n, A._ninv)
        ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &A, i - 1, k - 1, t1)
        ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &T, 0, k - 1, t2)
      end
    end
    #println("After reducing with T the augmented matrix looks like \n$A\n$T\n")

  end
  #return A
end

function _raw_setindex(A::nmod_mat, i::Int, j::Int, x::UInt)
  ccall((:nmod_mat_set_entry, :libflint), Void, (Ptr{nmod_mat}, Int, Int, UInt), &A, i - 1, j - 1, x)
end

function _stab(a::UInt, b::UInt, N::UInt)
  g = gcd(a,b,N)
  _, c = ppio(div(N,g),div(a,g))
  #c = prem(c,N)
  return c
end

function _unit(a::UInt, N::UInt)
  g,s,_ = gcdx(a,N)
  if g == 1
    return rem(s,N)
  else
    l = div(N,g)
    d = _stab(s,l,N)
    return rem(s + d*l,N)
  end
end

function _xxgcd(a::UInt, b::UInt, N::UInt)
  g,s,t = gcdx(Int(a),Int(b))
  if g == 0
    return (UInt(0), UInt(1), UInt(0), UInt(0), UInt(1))
  end
  u = UInt(prem(-div(Int(b),g),Int(N)))
  g = UInt(g)
  v = div(a,g)
  return (g,UInt(prem(s,Int(N))),UInt(prem(t,Int(N))),u,v)
end

################################################################################
#
#  positive remainder
#
################################################################################

function prem(a::Integer, m::T) where T<:Integer
  b = a % m
  if b < 0
    return m+b
  else
    return b
  end
end

function prem(a::fmpz, m::T) where T
  return prem(BigInt(a), m)
end

###############################################################################
#
#  Howell form for Generic.Mat{fmpz}
#
###############################################################################

function _xxgcd(a::fmpz, b::fmpz, N::fmpz)
  g,s,t = gcdx(a,b)
  if g == 0
    return (fmpz(0), fmpz(1), fmpz(0), fmpz(0), fmpz(1))
  end
  u = mod(-div(b,g),N)
  v = div(a,g)
  return (g,mod(s,N),mod(t,N),u,v)
end

function _unit(a::fmpz, N::fmpz)
  g,s,_ = gcdx(a,N)
  if g == 1
    return mod(s,N)
  else
    l = div(N,g)
    d = _stab(s,l,N)
    return mod(s + d*l,N)
  end
end

function _stab(a::fmpz, b::fmpz, N::fmpz)
  g = gcd(gcd(a,b),N)
  _, c = ppio(div(N,g),div(a,g))
  return c
end

function howell_form(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
  B=deepcopy(A)
  if rows(B)<cols(B)
    B=vcat(B, zero_matrix(base_ring(B), cols(B)-rows(B), cols(B)))
  end
  howell_form!(B)
  return B
end

#
#  for the in-place function, the number of rows must be at least equal to the number of columns
#

function howell_form!(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})

  R=base_ring(A)
  A1=lift(A)
  ccall((:fmpz_mat_howell_form_mod, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz}), A1, modulus(R))
  A=MatrixSpace(R, rows(A), cols(A))(A1)
  return A

  R=base_ring(A)
  n=R.modulus

  #
  #  Get an upper triangular matrix
  #

  for j=1:cols(A)
    for i=j+1:cols(A)
      g,s,t,u,v = _xxgcd(A[j,j].data,A[i,j].data,n)
      for k in 1:cols(A)
        t1 = s* A[j,k] + t* A[i,k]
        t2 = u* A[j,k] + v* A[i,k]
        A[j,k] = t1
        A[i,k] = t2
      end
    end
  end

  #
  #  Multiply the rows by annihlator of the pivot element and reduce
  #

  T = zero_matrix(R, 1, cols(A))
  for j in 1:cols(A)
    if A[j,j] != 0
       u = _unit(A[j,j].data, n)
      for k in 1:cols(A)
        A[j,k]=u*A[j,k]
      end
      for i in 1:j-1
        q = div(A[i,j].data, A[j,j].data)
        for k in j:cols(A)
          A[i,k]-= q*A[j,k]
        end
      end
      l = div(n,gcd(A[j,j].data,n))
      if l == 0
        continue
      end
      for k in j:cols(A)
       T[1,k]=l*A[j,k]
      end
    else
      for k in j:cols(A)
        T[1, k]= A[j,k]
      end
    end
    for i in j:cols(A)
      g,s,t,u,v = _xxgcd(A[i,i].data, T[1,i].data,n)
      for k in 1:cols(A)
        t1 = s*A[i,k] + t*T[1,k]
        t2 = u*A[i,k] + v*T[1,k]
        A[i,k]= t1
        T[1,k]=t2
      end
    end
  end
  if iszero_row(A,rows(A))
    for i=1:cols(A)
      A[rows(A),i]=T[1,i]
    end
  end
end

function triangularize!(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
  R=base_ring(A)
  n=R.modulus

  #
  #  Get an upper triangular matrix
  #

  for j=1:cols(A)
    for i=j+1:cols(A)
      g,s,t,u,v = _xxgcd(A[j,j].data,A[i,j].data,n)
      for k in 1:cols(A)
        t1 = s* A[j,k] + t* A[i,k]
        t2 = u* A[j,k] + v* A[i,k]
        A[j,k] = t1
        A[i,k] = t2
      end
    end
  end
end

function triangularize(A::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
  B= triangularize!(deepcopy(A))
  return B
end

function Base.nullspace(M::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
  R = base_ring(M)
  #
  #  If the modulus is prime, the second part of the Howell form computation is useless and I can directly call the rref
  #  but I have to test if the modulus is prime. What is better?
  #

  N = zero_matrix(R, rows(M)+cols(M), rows(M)+cols(M))
  for i=1:rows(M)
    for j=1:cols(M)
      N[j,i]=M[i,j]
    end
  end
  for i=1:cols(M)
    N[i, i+rows(M)]=1
  end
  N=howell_form!(N)

  if gcd(prod([N[i,i] for i=1:rows(N)]).data,R.modulus)==1
    return MatrixSpace(R,cols(M),0, false)(), 0
  end
  nr = 1
  while nr <= rows(N) && !iszero_row(N, nr)
    nr += 1
  end
  nr -= 1
  h = sub(N, 1:nr, 1:rows(M))
  for i=1:rows(h)
    if iszero_row(h, i)
      k = sub(N, i:rows(h), rows(M)+1:cols(N))
      return k', rows(k)
    end
  end
  return MatrixSpace(R,cols(M),0, false)(0), 0
end


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

# Krawczyk method, following Tucker

# const D = differentiate

doc"""Returns two intervals, the first being a point within the
interval x such that the interval corresponding to the derivative of f there
does not contain zero, and the second is the inverse of its derivative"""
function guarded_derivative_midpoint{T}(f::Function, f_prime::Function, x::Interval{T})

    α = convert(T, 0.46875)   # close to 0.5, but exactly representable as a floating point
    m = Interval( guarded_mid(x) )

    C = inv(f_prime(m))

    # Check that 0 is not in C; if so, consider another point rather than m
    i = 0
    while zero(T) ∈ C
        m = Interval( α*x.lo + (one(T)-α)*x.hi )
        C = inv(f_prime(m))

        i += 1
        α /= 2
        i > 20 && error("""Error in guarded_deriv_midpoint:
            the derivative of the function seems too flat""")
    end

    return m, C
end

function K{T}(f::Function, f_prime::Function, x::Interval{T})
    m, C = guarded_derivative_midpoint(f, f_prime, x)
    deriv = f_prime(x)
    Kx = m - C*f(m) + (one(T) - C*deriv) * (x - m)
    Kx
end

function krawczyk_refine{T}(f::Function, f_prime::Function, x::Interval{T};
    tolerance=eps(one(T)), debug=false)

    debug && (print("Entering krawczyk_refine:"); @show x)

    while diam(x) > tolerance  # avoid problem with tiny floating-point numbers if 0 is a root
        Kx = K(f, f_prime, x)
        debug && @show(x, Kx)
        Kx = Kx ∩ x
        Kx == x && break
        x = Kx
    end

    return [Root(x, :unique)]
end

function krawczyk{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
    tolerance=eps(one(T)), debug=false, maxlevel=30)

    debug && (print("Entering krawczyk:"); @show(level); @show(x))

    # Maximum level of bisection
    level >= maxlevel && return [Root(x, :unknown)]

    isempty(x) && return Root{T}[]  # [(x, :none)]
    Kx = K(f, f_prime, x) ∩ x

    isempty(Kx ∩ x) && return Root{T}[]  # [(x, :none)]

    if Kx ⪽ x
        debug && (print("Refining "); @show(x))
        return krawczyk_refine(f, f_prime, Kx, tolerance=tolerance, debug=debug)
    end

    m = guarded_mid(x)

    debug && @show(x,m)

    isthin(x) && return [Root(x, :unknown)]

    # bisecting
    roots = vcat(
        krawczyk(f, f_prime, Interval(x.lo, m), level+1,
                 tolerance=tolerance, debug=debug, maxlevel=maxlevel),
        krawczyk(f, f_prime, Interval(m, x.hi), level+1,
                 tolerance=tolerance, debug=debug, maxlevel=maxlevel)
        )

    roots = clean_roots(f, roots)

    return roots
end

# use automatic differentiation if no derivative function given
krawczyk{T}(f::Function,x::Interval{T}; args...) =
    krawczyk(f, x->D(f,x), x; args...)

# krawczyk for vector of intervals:
krawczyk{T}(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) =
    vcat([krawczyk(f, f_prime, @interval(x); args...) for x in xx]...)

krawczyk{T}(f::Function,  xx::Vector{Interval{T}}, level; args...) =
    krawczyk(f, x->D(f,x), xx; args...)

krawczyk{T}(f::Function,  xx::Vector{Root{T}}; args...) =
    krawczyk(f, x->D(f,x), [x.interval for x in xx]; args...)

ext: jl


text:
using PackageCompiler

create_app(
    ".",
    "build",
    #precompile_execution_file="precompile.jl",
    precompile_statements_file="traced.jl",
    force=true
)


text:
#
export expKernel, rotexpKernel
export computeK # should probably not be exported
export uncoupledMKernel, GIMKernel, rotKernel, KeplerKernel

export keplerKernel

abstract type kerneltype end
# abstract type abstractkronkernel <: kerneltype  end
# abstract type abstracmatrixkernel <: kerneltype end

# temporary
abstract type scalarkernel <: kerneltype  end
abstract type matrixkernel <: kerneltype end

abstract type MatrixKernel <: kerneltype end
abstract type KronMatrixKernel <: kerneltype  end

struct uncoupledMKernel{N, K<:NTuple{N, Kernel}, Qt <: NTuple{N, AbstractArray{<:Real, 2}} } <: KronMatrixKernel
    kernels::K
    Q::Qt
end

function uncoupledMKernel(K::KT, Q::QT) where {KT<:Kernel, QT <: AbstractArray{<:Real, 2} }
    return uncoupledMKernel((K,), (Q,))
end

function uncoupledMKernel(w::Array{T, 1}) where T<:Real
    K = (pskernel(w),)
    Q = (Diagonal(ones(length(w)-1)),)
    return uncoupledMKernel(K,Q)
end

# evaluate
function (mker::uncoupledMKernel)(z1,z2)
    K = mker.kernels
    Q = mker.Q
    return mapreduce(x -> x[1](z1,z2) * x[2], +, zip(K,Q))
end

# reparametrize
function (mker::uncoupledMKernel)(w::Array{<:Real, 1})
    K = mker.kernels

    ls = length.(GaussianProcessODEs.getparam.(K))
    endIdx = cumsum(ls)
    startIdx = endIdx .- (ls .- 1)
    ra = range.(startIdx, endIdx, step = 1)
    indw = getindex.(Ref(w), ra)

    K = map( (k,w) -> k(w), K, indw)
    typeof(mker)(K, mker.Q)
end

###
# Core structure that is part of each component of the structured matrix kernel
###
struct GIMKernel{GKP, K <: KernelFunctions.Kernel } <: MatrixKernel
    kernel::K
    groupaction::Function
    # may need to be changed for more general group action?
    parameterinterval::Tuple{Float64,Float64}
    gkparams::GKP

    GIMKernel{T,K}(ker::K, grpa, psp, gkparams::T) where {T,K} = new(ker, grpa, psp, gkparams)

end

# constructor
function GIMKernel(ker::K, grpa, parameterinterval, N::Int) where K
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, parameterinterval...))
    GIMKernel{typeof(gkparams), K}(ker, grpa, parameterinterval, gkparams)
end

function GIMKernel(ker::K, grpa, parameterinterval, N::Nothing) where K
    GIMKernel{typeof(N), K}(ker, grpa, parameterinterval, N)
end

# evaluation
function (ik::GIMKernel)(z1, z2)
        grpa = ik.groupaction
        gkp = ik.gkparams
        ker = ik.kernel

        base = gkp.weights .* map(x->ker(z1, grpa(x) * z2), gkp.x)
        return sum(base .* grpa.(gkp.x))
        # return sum(ker.(Ref(z1), grpa.(gkp.x) .* [z2]) .* f.(gkp.x) .* gkp.weights)
end

###
# Specialized implementation for rotational equivariant kernel
struct rotKernel{K <: KernelFunctions.Kernel, GKP }  <: MatrixKernel
    kernel::K
    gkparams::GKP
end

# constructor
function rotKernel(ker::K; N::Int) where K
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, 0, 2pi))
    rotKernel{K, typeof(gkparams)}(ker, gkparams)
end

# group action
function rot(phi)
    c = cos(phi)
    s = sin(phi)
    [c -s; s c]
end

# evaluation
function (rk::rotKernel)(z1, z2)
    gkp = rk.gkparams
    ker = rk.kernel

    base = map(x -> ker(z1, rot(x) * z2), gkp.x) .* gkp.weights
    costerm = sum(base .* cos.(gkp.x) )
    if z1 == z2
        sinterm = 0.
    else
        sinterm = sum(base .* sin.(gkp.x) )
    end
    [costerm -sinterm; sinterm costerm]
end

# reparametize
# there is numerical issues
function (mker::rotKernel)(w::Array{T, 1}) where {T <: Real}
    l = vcat(w[1:2], w[2])
    typeof(mker)(mker.kernel(l), mker.gkparams)
end

function getparam(mker::rotKernel) #where K <: MatrixKernel
    w = getparam(mker.kernel)
    return w[[1,2]]
end

###
# Specialized implementation for kepler system kernel
struct KeplerKernel{K <: KernelFunctions.Kernel, GKP }  <: MatrixKernel
    kernel::K
    gkparams::GKP
end

# constructor
function KeplerKernel(ker::K; N::Int) where K
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, 0, 2pi))
    KeplerKernel{K, typeof(gkparams)}(ker, gkparams)
end

function KeplerKernel(ker::K, auto) where K
    gkparams = nothing
    KeplerKernel{K, typeof(gkparams)}(ker, gkparams)
end

# group action
keplerrot(phi, x) = vcat(rot(phi)*x[1:2], rot(phi)*x[3:4])

# evaluation
function _computescossin(z1, z2, ker, gkp)
    base = map(x -> ker(z1, keplerrot(x, z2) ), gkp.x) .* gkp.weights
    costerm = sum(base .* cos.(gkp.x) )
    if z1 == z2
        sinterm = 0.
    else
        sinterm = sum(base .* sin.(gkp.x) )
    end
    return costerm, sinterm
end

function _computescossin(z1, z2, ker, ::Nothing)
    costerm, _ = quadgk(phi -> ker(z1, keplerrot(phi, z2) )* cos(phi), 0, 2*pi, rtol = 1e-8)
    if z1 == z2
        sinterm = 0.
    else
        sinterm, _ = quadgk(phi -> ker(z1, keplerrot(phi, z2) )* sin(phi), 0, 2*pi, rtol = 1e-8)
    end
    return costerm, sinterm
end

function (rk::KeplerKernel)(z1, z2, )
    ker = rk.kernel
    gkp = rk.gkparams

    costerm, sinterm = _computescossin(z1,z2, ker,gkp)

    return [costerm -sinterm 0 0;
        sinterm costerm 0 0;
        0 0 costerm -sinterm;
        0 0 sinterm costerm]
end

# reparametize
# there is numerical issues
function (mker::KeplerKernel)(w::Array{T, 1}) where {T <: Real}
    # l = vcat(w[1:2], w[2], w[2], w[2])
    l = vcat(w[1:2], w[2], w[3], w[3])
    typeof(mker)(mker.kernel(l), mker.gkparams)
end

function getparam(mker::KeplerKernel) #where K <: MatrixKernel
    w = getparam(mker.kernel)
    return w[[1,2,4]]
    # return w[[1,2]]
end

### Old stuff

struct expKernel <: scalarkernel
    param::Array{Float64, 1}
end
function kernelfunctionf(z1, z2, ker::expKernel)
    d = length(z1)
    w = ker.param[2:d+1] .^2  #better without square , very odd, likely issue with forward differences
    return ker.param[1]^2 * exp(-1/2 * wsqeuclidean(z1, z2, 1.0 ./w))
end

struct rotexpKernel <: matrixkernel
    param::Array{Float64, 1}
    phis::Array{Float64, 1}
    weights::Array{Float64, 1}

    function rotexpKernel(param; N = 30)
        phis, weights = gauss(N, 0, 2*pi)
        new(param, phis, weights)
    end
end

function rotexpIntegrand(phi, x1, x2, w)
    GaussianProcessODEs.kernelfunctionf(x1, rot(phi)* x2, expKernel(w))
end

function kernelfunctionf(z1, z2, ker::rotexpKernel)
    d = length(z1)
    w = ker.param
    # w = vcat(ker.param[1], ker.param[2:d+1] .^2)
    # costerm, _ = quadgk(phi -> rotexpIntegrand(phi, z1, z2, w)* cos(phi), 0, 2*pi, rtol = 1e-6)
    # sinterm, _ = quadgk(phi -> rotexpIntegrand(phi, z1, z2, w)* sin(phi), 0, 2*pi, rtol = 1e-6)

    phi = ker.phis
    weights = ker.weights

    base = rotexpIntegrand.(phi, Ref(z1), Ref(z2), Ref(w)) .*weights

    costerm = sum(base .* cos.(phi) )
    sinterm = sum(base .* sin.(phi) )

    K = [costerm -sinterm; sinterm costerm]
end

# very ad-hoc, needs much better implementation
struct keplerKernel <: matrixkernel
    param::Array{Float64, 1}
    phis::Array{Float64, 1}
    weights::Array{Float64, 1}

    function keplerKernel(param; N = 30)
        phis, weights = gauss(N, 0, 2*pi)
        new(param, phis, weights)
    end
end

function keplrotexpIntegrand(phi, x1, x2, w)
    x2r = vcat(rot(phi)*x2[1:2], rot(phi)*x2[3:4])
    GaussianProcessODEs.kernelfunctionf(x1, x2r, expKernel(w))
end
function kernelfunctionf(z1, z2, ker::keplerKernel)
    d = length(z1)
    w = ker.param
    # w = vcat(ker.param[1], ker.param[2:d+1] .^2)

    phi = ker.phis
    weights = ker.weights

    base1 = keplrotexpIntegrand.(phi, Ref(z1), Ref(z2), Ref(w))
    # base2 = rotexpIntegrand.(phi, Ref(z1[3:4]), Ref(z2[3:4]), Ref(w))

    costerm1 = sum(base1 .* cos.(phi) .*weights)
    sinterm1 = sum(base1 .* sin.(phi) .*weights)

    # costerm2 = sum(base2 .* cos.(phi) .*weights)
    # sinterm2 = sum(base2 .* sin.(phi) .*weights)

    K = [costerm1 -sinterm1 0 0;
         sinterm1 costerm1 0 0;
         0 0 costerm1 -sinterm1;
         0 0 sinterm1 costerm1]
    # K = [0 0 costerm2 -sinterm2;
    #      0 0 sinterm2 costerm2;
    #      costerm1 -sinterm1 0 0;
    #      sinterm1 costerm1 0 0]
end

function dist(z1, z2, ker)
    d = length(z1)
    w = ker.param[2:d+1] .^2
    L = diagm( ( 1.0 ./ w))
    x = z1 - z2
    ker.param[1]^2 * exp(-1/2 * x'*L*x )
end

function derivativekernelfunctionf(z1, z2, ker::expKernel)
    d = length(z1)
    return ker.param[1]^2 * exp(-1/2 * weuclidean(z1, z2, 1.0 ./ ker.param[2:d+1])^2) * ( (z1.-z2)./ ker.param[2:d+1])
end

# function matrixkernelfunction(z1::Array{Float64, 1}, z2::Array{Float64, 1}, ker::T) where T <: scalarkernel
#     Id = Matrix{Float64}(LinearAlgebra.I, length(z1), length(z2))
#     return kernelfunctionf(z1, z2, ker) * Id
# end


text:
"""
An abstract type for implementing coating correction algorithms.

    NeXLCore.transmission(zaf::CoatingCorrection, xray::CharXRay)
"""
abstract type CoatingCorrection end

NeXLCore.minproperties(::Type{<:CoatingCorrection}) = ( )  # None are necessary but :Coating is optional

"""
    Coating

Implements a simple multi-layer coating correction.
"""
struct Coating <: CoatingCorrection
    layers::Vector{Film}
    Coating(coatings::AbstractVector{Film}) = new(coatings)
end

coatingcorrection(::Type{Coating}, film::Film) = Coating([film])
coatingcorrection(::Type{Coating}, films::AbstractVector{Film}) = Coating(films)

NeXLCore.minproperties(::Type{Coating}) = ( :Coating )  # :Coating = Film() | Film[]

"""
    carboncoating(nm)

Constructs a carbon coating of the specified thickness (in nanometers).
"""
carboncoating(nm) = Film(pure(n"C"), nm * 1.0e-7)

"""
    NeXLCore.transmission(zaf::Coating, xray::CharXRay, toa)

Calculate the transmission fraction for the specified X-ray through the coating
in the direction of the detector.
"""
NeXLCore.transmission(cc::Coating, xray::CharXRay, θtoa::AbstractFloat) =
    mapreduce(lyr -> transmission(lyr, xray, θtoa), *, cc.layers, init = 1.0)

Base.show(io::IO, coating::Coating) = length(coating.layers) == 1 ? #
    (thickness(coating.layers[1]) <= 0.0 ? Base.show(io, "No coating") : Base.show(io, coating.layers[1])) : #
    Base.show(io, coating.layers)

"""
    NullCoating

No coating (100%) transmission.
"""
struct NullCoating <: CoatingCorrection
    NullCoating(films::AbstractVector{Film}) = new()
    NullCoating(film::Film) = new()
    NullCoating() = new()
end

coatingcorrection(::Type{<:CoatingCorrection}, mss::Missing) = NullCoating()
coatingcorrection(::Type{NullCoating}, film::Film) = NullCoating()
coatingcorrection(::Type{NullCoating}, film::AbstractVector{Film}) = NullCoating()

"""
    NeXLCore.transmission(zaf::NullCoating, xray::CharXRay, toa)

Calculate the transmission fraction for the specified X-ray through no coating (1.0).
"""
NeXLCore.transmission(nc::NullCoating, xray::CharXRay, θtoa::AbstractFloat) = 1.0

Base.show(io::IO, nc::NullCoating) = print(io, "No coating.")

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
abstract type TurbineGov <: GeneratorComponent end

@def TG_ports begin
      state_input = [:ω]
      inner_input = Vector{Int64}()
end

"""
Parameters of a fixed Turbine Governor that returns a fixed mechanical torque
given by the product of P_ref*efficiency

# Conmutable structor
```julia
TGFixed(efficiency)
```

# Arguments
* `efficiency`::Float64 : Efficiency factor that multiplies P_ref

"""
mutable struct TGFixed <: TurbineGov
      efficiency::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function TGFixed(efficiency::Float64)
            n_states = 0
            states = Vector{Symbol}()
            @TG_ports
            new(efficiency, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of a Turbine Governor Type I.

# Conmutable structor
```julia
TGTypeI(R, Ts, Tc, T3, T4, T5, P_min, P_max)
```

# Arguments
* `R`::Float64 : Droop parameter
* `Ts`::Float64 : Governor time constant
* `Tc`::Float64 : Servo time constant
* `T3`::Float64 : Transient gain time constant
* `T4`::Float64 : Power fraction time constant
* `T5`::Float64 : Reheat time constant
* `P_min`::Float64 : Min Power into the Governor
* `P_max`::Float64 : Max Power into the Governor

"""
mutable struct TGTypeI <: TurbineGov
      R::Float64
      Ts::Float64
      Tc::Float64
      T3::Float64
      T4::Float64
      T5::Float64
      P_min::Float64
      P_max::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function TGTypeI(R::Float64,
                       Ts::Float64,
                       Tc::Float64,
                       T3::Float64,
                       T4::Float64,
                       T5::Float64,
                       P_min::Float64,
                       P_max::Float64)

            n_states = 3
            states = [:x_g1, :x_g2, :x_g3]
            @TG_ports
            new(R, Ts, Tc, T3, T4, T5, P_min,
                P_max, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of a Turbine Governor Type II.

# Conmutable structor
```julia
TGTypeI(R, T1, T2, τ_min, τ_max)
```

# Arguments
* `R`::Float64 : Droop parameter
* `T1`::Float64 : Transient gain time constant
* `T2`::Float64 : Governor time constant
* `τ_min`::Float64 : Min turbine torque output
* `τ_max`::Float64 : Max turbine torque output

"""
mutable struct TGTypeII <: TurbineGov
      R::Float64
      T1::Float64
      T2::Float64
      τ_min::Float64
      τ_max::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function TGTypeII(R::Float64,
                       T1::Float64,
                       T2::Float64,
                       τ_min::Float64,
                       τ_max::Float64)

            n_states = 1
            states = [:xg]
            @TG_ports
            new(R, T1, T2, τ_min, τ_max, n_states, states, Ports(state_input, inner_input))
      end
end


text:
module SymmetricSparseMatricesBenchmarks

using BenchmarkTools

begin include("SymmetricSparseMatrix/benchmarks.jl") end

end # module


text:
# N Queens Problem.

# (c) E. Adrian Garro S. Costa Rica Institute of Technology.

type NQueens
    queens_quant::Int
    solutions::Array{Array{Int}}
    solutions_num::Int
end

"""
`permutations(xs)`
Permutations without repeated elements.
"""
function permutations(xs)
    if xs == []
        produce(Int[])
    else
        for x in xs
            temp = copy(xs)
            filter!(e -> e != x, temp)
            for ys in @task permutations(temp)
                produce(union(Int[x], ys))
            end
        end
    end
end

function set_queens_quant!(problem::NQueens, queens_quant::Int)
    if queens_quant > 3
        problem.queens_quant = queens_quant
    else
        error("Invalid number for this problem.")
    end
end

function attack(queens_cols::Array{Int})
    queens_quant = length(queens_cols)
    queens_rows = 1:queens_quant
    queens_pos = zip(queens_rows, queens_cols)
    sum_of_pos = [sum(pos) for pos in queens_pos]
    negative_queens_cols = [-1 * queen_col for queen_col in queens_cols]
    queens_pos = zip(queens_rows, negative_queens_cols)
    diff_of_pos = [sum(pos) for pos in queens_pos]
    # all results are different?
    in_diag45 = length(sum_of_pos) > length(Set(sum_of_pos))
    in_diag135 = length(diff_of_pos) > length(Set(diff_of_pos))
    there_attack = in_diag45 || in_diag135
end

function not_attack(queens_cols::Array{Int})
    !attack(queens_cols)
end

function solve!(problem::NQueens)
    if problem.queens_quant != 0
        all_possible_queens_cols = @task permutations(
            collect(1:problem.queens_quant)
        )
        problem.solutions = collect(all_possible_queens_cols)
        filter!(not_attack, problem.solutions)
        problem.solutions_num = length(problem.solutions)
    else
       error("Number of queens not available.")
    end
end

function print_solution(problem::NQueens, solution_num::Int)
    if solution_num in collect(1:problem.solutions_num)
        solution = problem.solutions[solution_num]
        function create_box(queen_pos::Int)
            queen_pos -= 1
            string(
                "\n|",
                "   |" ^ queen_pos,
                " X |",
                "   |" ^ (problem.queens_quant - (queen_pos + 1))
            )
        end
        boxes = [create_box(queens_cols) for queens_cols in solution]
        chain = string("\n+", "---+"  ^ problem.queens_quant)
        chessboard = string(string(chain, join(boxes, chain)), chain)
        println(string("\nSolution #", solution_num, "\n"))
        println(chessboard)
    else
        error("Solution not exist.")
    end
end

max_stars_repo_name: adriangarro/nqueens

max_issues_repo_name: adriangarro/nqueens

max_forks_repo_name: adriangarro/nqueens


text:
export OnePlusEvo, evaluate!, populate!

"""

OnePlusEvo implements a 1+λ evolutionary algorithm. To do this, a new subtype of
AbstractEvolution is created, OnePlusEvo, and the following functions are
defined for this type:

evaluate
populate

the function generation can also be defined. Without specification, the default
generation function (which is void) is applied.
"""
mutable struct OnePlusEvo{T} <: AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{T}
    fitness::Function
    gen::Int
end

function OnePlusEvo{T}(cfg::NamedTuple, fitness::Function;
                      logfile=string("logs/", cfg.id, ".csv")) where T
    logger = CambrianLogger(logfile)
    population = initialize(T, cfg)
    OnePlusEvo(cfg, logger, population, fitness, 0)
end

function oneplus_populate!(e::AbstractEvolution)
    p1 = max_selection(e.population)
    e.population[1] = p1
    for i in 2:prod(e.config.population.size)
        e.population[i] = mutate(p1)
    end
end

evaluate!(e::OnePlusEvo) = fitness_evaluate!(e, e.fitness)
populate!(e::OnePlusEvo) = oneplus_populate!(e)


text:
"""
    hadi1994(multivariateData)

Perform Hadi (1994) algorithm for a given multivariate data.

# Arguments
- `multivariateData::Array{Float64, 2}`: Multivariate data.

# Description
Algorithm starts with an initial subset and enlarges the subset to
obtain robust covariance matrix and location estimates. This algorithm
is an extension of `hadi1992`.

# Output
- `["outliers"]`: Array of indices of outliers
- `["critical.chi.squared"]`: Threshold value for determining being an outlier
- `["rth.robust.distance"]`: rth robust distance, where (r+1)th robust distance is the first one that exceeds the threshold.

# Examples
```julia-repl
julia> multidata = hcat(hbk.x1, hbk.x2, hbk.x3);

julia> hadi1994(multidata)
Dict{Any,Any} with 3 entries:
  "outliers"              => [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
  "critical.chi.squared" => 7.81473
  "rth.robust.distance"   => 5.04541
```

# Reference
Hadi, Ali S. "A modification of a method for the dedection of outliers in multivariate samples"
Journal of the Royal Statistical Society: Series B (Methodological) 56.2 (1994): 393-396.
"""
function hadi1994(multivariateData::Array{Float64,2}; alpha=0.05)
    n, p = size(multivariateData)
    h = Int(round((n + p + 1.0) / 2.0))
    cnp = (1 + ((2) / (n - 1 - 3p)) + ((p + 1) / (n - p)))^2
    chi_50_quantile = quantile(Chisq(p), 0.50)
    critical_quantile = quantile(Chisq(p), 1 - alpha)
    allindices = collect(1:n)

    meds = coordinatwisemedians(multivariateData)
    Sm = (1.0 / (n - 1.0)) * (multivariateData .- meds')' * (multivariateData .- meds')
    mah0 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=meds, covmatrix=Sm))
    ordering_indices_mah0 = sortperm(mah0)
    best_indices_mah0 = ordering_indices_mah0[1:h]
    starting_data = multivariateData[best_indices_mah0, :]

    Cv = coordinatwisemedians(starting_data)
    Sv = (1.0 / (h - 1.0)) * (starting_data .- Cv')' * (starting_data .- Cv')
    mah1 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=Cv, covmatrix=Sv))
    ordering_indices_mah1 = sortperm(mah1)

    r = p + 1
    basic_subset_indices = []
    basic_subset = []
    sorted_mah1 = []
    Cb = []
    Sb = []
    cfactor = 0
    isFullRank = false

    while r < n

        isFullRank = false

        while !isFullRank
            basic_subset_indices = ordering_indices_mah1[1:r]
            basic_subset = multivariateData[basic_subset_indices, :]
            Cb = applyColumns(mean, basic_subset)
            Sb = cov(basic_subset)
            cfactor = cnp * sqrt(sort(mah1)[h]) / chi_50_quantile
            r += 1

            if (det(cfactor * Sb) == 0 && r < n)
                isFullRank = false
            else
                isFullRank = true
            end
        end

        mah1 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=Cb, covmatrix=(cfactor * Sb)))

        ordering_indices_mah1 = sortperm(mah1)
        basic_subset_indices = ordering_indices_mah1[1:r]

        sorted_mah1 = sort(mah1)
        if sorted_mah1[r] >= critical_quantile
            break
        end
    end

    outlierset = setdiff(allindices, basic_subset_indices)

    result = Dict()
    result["outliers"] = sort(outlierset)
    result["critical.chi.squared"] = critical_quantile
    result["rth.robust.distance"] = sorted_mah1[r - 1]
    return result
end

max_stars_repo_name: jbytecode/LinRegOutliers

max_issues_repo_name: jbytecode/LinRegOutliers

max_forks_repo_name: jbytecode/LinRegOutliers


text:
module Foo

using Requires
using Base.Test

beforeflag = false
afterflag = false

@require JSON global beforeflag = true

@test !beforeflag

using JSON

@test beforeflag

@require JSON global afterflag = true

@test afterflag

end


text:
# This file is a part of PropDicts.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package PropDicts" begin

include("test_varsubst.jl")
include("test_propdict.jl")

end # testset


text:
using ComoniconOptions
using ComoniconBuilder
using Test
using ComoniconBuilder: install, install_script, install_completion,
    cmd_script, contains_path, contains_fpath, write_path, write_fpath, detect_rcfile,
    install_env_path

home_dir = mktempdir()
test_dir = pkgdir(ComoniconBuilder, "test")
usr_dir = joinpath(test_dir, "usr")

options = ComoniconOptions.Comonicon(
    name="test",
    install=ComoniconOptions.Install(
        path=usr_dir,
        completion=false,
    ),
    sysimg=nothing,
)

module Foo
using ComoniconTestUtils
const CASTED_COMMANDS = Dict{String, Any}(
    "main" => rand_command()
)
end

@testset "cmd script" begin
    sysimg = joinpath(options.install.path, "lib", ComoniconBuilder.sysimg("test"))
    entry_script = cmd_script(
        Foo,
        joinpath(options.install.path, "bin", "test.jl");
        project=joinpath(options.install.path, "test"),
        compile=:min,
        sysimg
    )

    @test occursin("#!/bin/sh", entry_script)
    @test occursin(ComoniconBuilder.default_exename(), entry_script)
    @test occursin("-O2", entry_script)
    @test occursin("-startup-file=no", entry_script)
    @test occursin("--compile=min", entry_script)
    @test occursin("-J$sysimg", entry_script)
    @test occursin(joinpath(options.install.path, "bin", "test.jl"), entry_script)
end

@testset "detect_shell" begin
    withenv("SHELL"=>"zsh") do
        @test ComoniconBuilder.detect_shell() == ENV["SHELL"]
    end
end

@testset "contains path/write path" begin
    withenv("PATH"=>nothing, "FPATH"=>nothing) do
        dir = mktempdir()
        rcfile = joinpath(dir, ".bashrc")
        touch(rcfile)
        @test contains_path(rcfile, usr_dir, Base.EnvDict()) == false

        write_path(rcfile, usr_dir)
        @test contains_path(rcfile, usr_dir) == true

        @test contains_fpath(rcfile, usr_dir, Base.EnvDict()) == false

        write_fpath(rcfile, usr_dir)
        @test contains_fpath(rcfile, usr_dir) == true
    end

    withenv("PATH"=>joinpath(usr_dir, "bin"), "FPATH"=>joinpath(usr_dir, "completions")) do
        @test contains_path("rcfile", usr_dir) == true
        @test contains_fpath("rcfile", usr_dir) == true
    end
end

@testset "detect rcfile" begin
    @test detect_rcfile("zsh", home_dir) == joinpath(home_dir, ".zshrc")
    @test withenv("ZDOTDIR"=>"test") do
        detect_rcfile("zsh", home_dir)
    end == joinpath("test", ".zshrc")

    @test detect_rcfile("bash", home_dir) == joinpath(home_dir, ".bash_profile")
    touch(joinpath(home_dir, ".bashrc"))
    @test detect_rcfile("bash", home_dir) == joinpath(home_dir, ".bashrc")
end

@testset "install env path" begin
    install_env_path(usr_dir, "zsh", home_dir, ENV, true)

    withenv("PATH"=>nothing, "FPATH"=>nothing) do
        @test contains_path(joinpath(home_dir, ".zshrc"), usr_dir)
        @test contains_fpath(joinpath(home_dir, ".zshrc"), usr_dir)
    end
end

@testset "install completion" begin
    install_completion(Foo, options, "zsh")
    @test isfile(joinpath(usr_dir, "completions", "_test"))
    rm(usr_dir; recursive=true)
end

@testset "install script" begin
    install_script(Foo, options; project=joinpath(options.install.path, "test"))
    @test isfile(joinpath(usr_dir, "bin", "test"))
    @test isfile(joinpath(usr_dir, "bin", "test.jl"))
    rm(usr_dir; recursive=true)
end

using Pkg
Pkg.activate(joinpath(test_dir, "project"))
Pkg.develop(path=pkgdir(ComoniconBuilder))
Pkg.test()


text:
abstract type SVAE end

normalizecolumns(m::AbstractArray{T, 2}) where {T<:Number} = m ./ sqrt.(sum(m .^ 2, dims = 1) .+ eps(T))
function normalizecolumns!(m::AbstractArray{T, 2}) where {T<:Number}
	m .= normalizecolumns(m)
end

"""
	vmfentropy(m, κ)

	Entropy of Von Mises-Fisher distribution
"""
vmfentropy(m, κ) = .-κ .* besselix(m / 2, κ) ./ besselix(m / 2 - 1, κ) .- ((m ./ 2 .- 1) .* log.(κ) .- (m ./ 2) .* log(2π) .- (κ .+ log.(besselix(m / 2 - 1, κ))))

"""
	huentropy(m)

	Entropy of Hyperspherical Uniform distribution
"""
huentropy(m) = m / 2 * log(π) + log(2) - lgamma(m / 2)

"""
	kldiv(model::SVAE, κ)

	KL divergence between Von Mises-Fisher and Hyperspherical Uniform distributions
"""
kldiv(model::SVAE, κ) = .- vmfentropy(model.zdim, κ) .+ model.hue

log_normal(x) = - sum((x .^ 2), dims = 1) ./ 2 .- size(x, 1) .* log(2π) ./ 2
log_normal(x, μ) = log_normal(x - μ)
log_normal(x,μ, σ2::AbstractArray{T}) where {T<:Number} = - sum((x - μ) .^ 2 ./ σ2 .+ log.(σ2 .* 2π), dims = 1) / 2

# Likelihood estimation of a sample x under VMF with given parameters taken from https://pdfs.semanticscholar.org/2b5b/724fb175f592c1ff919cc61499adb26996b1.pdf
# normalizing constant for density function of VMF
c(p, κ) = κ ^ (p / 2 - 1) / ((2π) ^ (p / 2) * besseli(p / 2 - 1, κ))

# log likelihood of one sample under the VMF dist with given parameters
log_vmf_c(x, μ, κ) = κ * μ' * x .+ log(c(length(μ), κ))
log_vmf_wo_c(x, μ, κ) = κ * μ' * x

pairwisecos(x, y) = max.(1 .- (x' * y), 0) # This is a bit of a hack to avoid the distance being negative due to float
pairwisecos(x) = pairwisecos(x, x)

k_imq(x, y, c) = sum( c./ (c .+ pairwisecos(x, y))) / (size(x, 2) * size(y, 2))
k_imq(x::T, c) where {T <: AbstractMatrix} = sum(c ./ (c .+ pairwisecos(x))) / (size(x, 2) * (size(x, 2) - 1))
k_imq(x::T, c) where {T <: AbstractVector} = zero(eltype(x))

mmd_imq(x,y,c) = k_imq(x,c) + k_imq(y,c) - 2 * k_imq(x,y,c)

function samplehsuniform(size...)
	v = randn(size...)
	v = normalizecolumns(v)
end

function log_pxexpectedz(m::SVAE, x)
	xgivenz = m.g(zparams(m, x)[1])
	Flux.Tracker.data(log_normal(x, xgivenz))
end

cart3dtoangular(x) = [acos(1.0 * x[3] / norm(x)), 1.0 * atan(x[2] / x[1])]

function cart2spherical(x)
	scoords = vcat([acot(x[i] / norm(x[(i + 1):end])) for i in 1:(length(x) - 2)]...)
	scoords = vcat(scoords, 2acot((x[end - 1] + norm(x[(end-1):end])) / x[end]))
end

jacobian_encoder(m::SVAE, x) = Flux.Tracker.jacobian(a -> Chain(m.q, m.μzfromhidden)(a), x)

function log_pz_jacobian_encoder(m::SVAE, x)
	if size(x, 2) > 1
		xs = [x[:, i] for i in 1:size(x, 2)]
		return map(x -> log_pz_jacobian_encoder_singleinstance(m, x), xs)
	else
		return log_pz_jacobian_encoder_singleinstance(m, x)
	end
end

function log_pz_jacobian_encoder_singleinstance(m::SVAE, x)
	@assert size(x, 2) == 1
	s = svd(jacobian_encoder(m, x).data)
	d = reduce(+, log.(abs.(s.S))) * 2
	d + log_pz(m, x)
end

jacobian_decoder(m::SVAE, z) = Flux.Tracker.jacobian(a -> m.g(a), z)

function log_pz_jacobian_decoder(m::SVAE, z)
	if size(z, 2) > 1
		zs = [z[:, i] for i in 1:size(z, 2)]
		return map(z -> log_pz_jacobian_decoder_singleinstance(m, z), zs)
	else
		return log_pz_jacobian_decoder_singleinstance(m, z)
	end
end

function log_pz_jacobian_decoder_singleinstance(m::SVAE, z)
	@assert size(z, 2) == 1
	s = svd(jacobian_decoder(m, z).data)
	d = reduce(+, log.(abs.(s.S))) * 2
	-d + log_pz_from_z(m, z)
end

# function log_pz_jacobian(m::SVAE, x, z)
# 	if size(x, 2) > 1
# 		xs = [x[:, i] for i in 1:size(x, 2)]
# 		zs = [z[:, i] for i in 1:size(z, 2)]
# 		return map((x, z) -> log_pz_jacobian_singleinstance(m, x, z), xs, zs)
# 	else
# 		return log_pz_jacobian_singleinstance(m, x, z)
# 	end
# end

# function log_pz_jacobian_singleinstance(m::SVAE, x, z)
# 	@assert size(x, 2) == 1
# 	s = svd(jacobian(m, x).data)
# 	d = reduce(+, log.(abs.(s.S)))
# 	d + log(pz_from_z(m, z))
# end

"""
	infer(m::SVAE, x)

	infer latent variables and sample output x
"""
function infer(m::SVAE, x)
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	xgivenz = m.g(z)
	return xgivenz
end

"""
	zparams(model::SVAE, x)

	Computes μ and κ from the hidden layer
"""

function zparams(model::SVAE, x)
	hidden = model.q(x)
	# return model.μzfromhidden(hidden), max.(model.κzfromhidden(hidden), 100)
	return model.μzfromhidden(hidden), model.κzfromhidden(hidden)
end

"""
	zfromx(m::SVAE, x)

	convenience function that returns latent layer based on the input `x`
"""
zfromx(m::SVAE, x) = samplez(m, zparams(m, x)...)

"""
	samplez(m::SVAE, μz, κz)

	samples z layer based on its parameters
"""
function samplez(m::SVAE, μz, κz)
	ω = sampleω(m, κz)
	normal = Normal()
	v = Adapt.adapt(eltype(Flux.Tracker.data(κz)), rand(normal, size(μz, 1) - 1, size(μz, 2)))
	v = normalizecolumns(v)
	z = householderrotation(vcat(ω, sqrt.(1 .- ω .^ 2 .+ eps(Float32)) .* v), μz)
	return z
end

"""
	householderrotation(zprime, μ)

	uses Householder reflection to rotate the `zprime` vector according to mapping of e1 to μ
"""
function householderrotation(zprime, μ)
	e1 = similar(μ) .= 0
	e1[1, :] .= 1
	u = e1 .- μ
	normalizedu = normalizecolumns(u)
	return zprime .- 2 .* sum(zprime .* normalizedu, dims = 1) .* normalizedu
end

function sampleω(model::SVAE, κ)
	m = model.zdim
	c = @. sqrt(4κ ^ 2 + (m - 1) ^ 2)
	b = @. (-2κ + c) / (m - 1)
	a = @. (m - 1 + 2κ + c) / 4
	d = @. (4 * a * b) / (1 + b) - (m - 1) * log(m - 1)
	ω = rejectionsampling(m, a, b, d, κ)
	return ω
end

function rejectionsampling(m, a, b, d, κ)
	beta = Beta((m - 1) / 2, (m - 1) / 2)
	T = eltype(Flux.Tracker.data(a))
	ϵ, u = Adapt.adapt(T, rand(beta, size(a)...)), rand(T, size(a)...)

	accepted = isaccepted(ϵ, u, m, Flux.data(a), Flux.Tracker.data(b), Flux.data(d))
	it = 0
	while (!all(accepted)) & (it < 10000)
		mask = .! accepted
		ϵ[mask] = Adapt.adapt(T, rand(beta, sum(mask)))
		u[mask] = rand(T, sum(mask))
		accepted[mask] = isaccepted(mask, ϵ, u, m, Flux.data(a), Flux.data(b), Flux.data(d))
		it += 1
	end
	if it >= 10000
		println("Warning - sampler was stopped by 10000 iterations - it did not accept the sample!")
		println(m)
		println(a)
		println(b)
		println(d)
		println(κ)
		println(accepted)
	end
	return @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
end

# function rejectionsampling(m, a, b, d)
# 	beta = Beta((m - 1) / 2, (m - 1) / 2)
# 	T = eltype(Flux.Tracker.data(a))
# 	ϵ, u = Adapt.adapt(T, rand(beta, size(a)...)), rand(T, size(a)...)
#
# 	accepted = isaccepted(ϵ, u, m, Flux.data(a), Flux.Tracker.data(b), Flux.data(d))
# 	while !all(accepted)
# 		mask = .! accepted
# 		ϵ[mask] = Adapt.adapt(T, rand(beta, sum(mask)))
# 		u[mask] = rand(T, sum(mask))
# 		ia = isaccepted(mask, ϵ, u, m, Flux.data(a), Flux.data(b), Flux.data(d))
# 		accepted[mask] = ia
# 	end
# 	return @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
# end

isaccepted(mask, ϵ, u, m:: Int, a, b, d) = isaccepted(ϵ[mask], u[mask], m, a[mask], b[mask], d[mask]);
function isaccepted(ϵ, u, m:: Int, a, b, d)
	ω = @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
	t = @. 2 * a * b / (1 - (1 - b) * ϵ)
	@. (m - 1) * log(t) - t + d >= log(u)
end


text:
using Test, AlgoritmiStruttureDati.Laboratorio3
a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

@testset "sumBetweenElements.jl" begin
    @test sumBetweenElements(a, 0, 9) == 45
    @test sumBetweenElements(a, 4, 5) == 9
end


text:
module Jacobian

using PDESolver
using ODLCommonTools
using SummationByParts
using PETSc2
using Utils

# residual_evaluation.jl
export physicsRhs, assembleResidual

# sparse.jl
export getBlockSparsityCounts, INVISCID, VISCOUS, COLORING, VISCOUSTIGHT

# jacobian.jl - main functions
export physicsJac

# jacobian.jl - explicit Jacobian calculation
export _AssembleElementData, assembleElement, assembleInterface,
       assembleInterfaceVisc, assembleSharedFace, assembleBoundary,
       assembleBoundaryFull, assembleInterfaceFull, assembleSharedFaceFull,
       NullAssembleElementData, calcJacCol

# jacobian_diag.jl
export DiagJac, diagMatVec, AssembleDiagJacData, NullDiagJac,
       NullAssembleDiagJacData

# experimental stuff
export filterDiagJac, removeUnstableModes!, findStablePerturbation!

include("residual_evaluation.jl")
include("sparse.jl")
include("jacobian.jl")
include("jacobian_diag.jl")

end  # module


text:
module MyTools

using InteractiveUtils

export where

function all_loaded_modules()
	modules = []
	for mod in Base.loaded_modules_array()
		for name in names(mod, all=true)
			#println("name = ", name)
			if ! isdefined(mod, name)
				continue
			end
			obj = getfield(mod, name)
			if isa(obj, Module)
				push!(modules, obj)
			end
		end
	end
	return modules
end

"""
Returns an array of all loaded modules in which the symbol x is defined.
"""
function where(x::Symbol)
	modules = []
 	for mod in all_loaded_modules()
		if isdefined(mod, x)
			push!(modules, mod)
		end
	end
	return modules
end

function subtypes(t::DataType; recurse::Bool = false)
    #println("MyTools.subtypes: recurse = ", recurse)
    stypes = InteractiveUtils.subtypes(t)
    isempty(stypes) && return stypes
    if recurse
        for stype in stypes
            #println("stype = ", stype)
            if isa(stype, DataType)
                append!(stypes, subtypes(stype, recurse = true))
            end
        end
    end
    return unique(stypes)
end

function discover_interface(t::DataType)
    stypes = MyTools.subtypes(t, recurse=true)
    countmap = Dict{String,Int}()
    for stype in stypes
        stypemethods = methodswith(stype)
        for m in stypemethods
            mname = String(m.name)
            countmap[mname] = get(countmap, mname, 0) + 1
        end
    end

    return sort(collect(countmap), by = x->last(x), rev=true)
end

end # module

max_stars_repo_name: tgflynn/MyTools

max_issues_repo_name: tgflynn/MyTools

max_forks_repo_name: tgflynn/MyTools


text:
function rand_discr(p; m = 1)
    """
        rand_discr - discrete random generator

          y = rand_discr(p, n);

          y is a random vector of length n drawn from
          a variable X such that
              p(i) = Prob( X=i )

          Copyright (c) 2004 Gabriel Peyré
    """

    # makes sure it sums to 1
    p = p/sum(p)

    n = length(p)
    coin = rand(m)
    cumprob = [0.]
    append!(cumprob, cumsum(p))
    sample = zeros(m)

    for j in 1:n
        ind = findall((coin .> cumprob[j]) & (coin .<= cumprob[j+1]))
        sample[ind] = j
    end
    return sample
end


text:
# Draw a series of graphs as svg files

using LightGraphs
using GraphPlot
using PETScBinaryIO
using Compose
using Colors
using HDF5
using MatrixMarket
using ArgParse
using PerceptualColourMaps

s = ArgParseSettings()
@add_arg_table s begin
    "--eig"
        action = :store_true
        help = "color nodes by eigenvalue"
    "INPUT.petsc"
        help = "input graph"
        required = true
    "OUTPUT.pdf"
        help = "output filename"
        required = true
end
parsed_args = parse_args(ARGS, s)

function write_chaco(filename, g)
    open(filename, "w") do f
        write(f, "$(nv(g)) $(ne(g))\n")
        for i = 1:nv(g)
            for n in neighbors(g, i)
                write(f, "$n ")
            end
            write(f, "\n")
        end
    end
end

function map_colormap(vals, clm)
    min = minimum(vals)
    max = maximum(vals)
    map(x -> clm[Int(floor(1+(x - min)/(max-min)*(length(clm)-1)))], vals)
end

function compute_eig(A)
    mmwrite("tmp.mtx", spdiagm(vec(sum(A,1))) - A)
    run(`build/bin/ligmg_eigensolve tmp.mtx tmp.h5 -eps_smallest_magnitude -st_pc_type ligmg -st_ksp_type preonly -eps_monitor -st_type sinvert -eps_target 0 -eps_nev 1`)
    d = h5read("tmp.h5", "/")
    rm("tmp.h5")
    rm("tmp.mtx")
    d["eigen_vector_0"]
end

function layout_graph(A)
    write_chaco("tmp.graph", A)
    x = readdlm(read(`../ogdf_layout/build/ogdf_layout tmp.graph`))
    rm("tmp.graph")
    x[:,1], x[:,2]
end

println("Reading")
infile = parsed_args["INPUT.petsc"]
A = readPETSc(infile)
A = (A + A.')/2
if all(x->x==0, sum(A,1))
    A = spdiagm(diag(A)) - A
end
g = Graph(A)

plot = if parsed_args["eig"]
    println("Computing smallest eigenvalue")
    eig = compute_eig(A)
    nodefill = map_colormap(eig, cmap("D1"))
    gplot(g; layout=layout_graph, nodefillc=nodefill, NODESIZE=1/sqrt(nv(g)), edgestrokec=RGBA(0,0,0,0.01))
else
    gplot(g; layout=layout_graph, nodefillc=colorant"purple", EDGELINEWIDTH=15/sqrt(nv(g)), NODESIZE=0.5/sqrt(nv(g)), edgestrokec=RGBA(0,0,0,0.3))
end

outfile = parsed_args["OUTPUT.pdf"]
println("Saving to $outfile")
draw(PNG(outfile, 2400px, 2400px, ), plot)

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
# Aggregate expression.

mutable struct AggregateNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    name::Symbol
    distinct::Bool
    args::Vector{SQLNode}
    filter::Union{SQLNode, Nothing}

    AggregateNode(;
                  over = nothing,
                  name::Union{Symbol, AbstractString},
                  distinct = false,
                  args = SQLNode[],
                  filter = nothing) =
        new(over, Symbol(name), distinct, args, filter)
end

AggregateNode(name; over = nothing, distinct = false, args = SQLNode[], filter = nothing) =
    AggregateNode(over = over, name = name, distinct = distinct, args = args, filter = filter)

AggregateNode(name, args...; over = nothing, distinct = false, filter = nothing) =
    AggregateNode(over = over, name = name, distinct = distinct, args = SQLNode[args...], filter = filter)

"""
    Agg(; over = nothing, name, distinct = false, args = [], filter = nothing)
    Agg(name; over = nothing, distinct = false, args = [], filter = nothing)
    Agg(name, args...; over = nothing, distinct = false, filter = nothing)

An application of an aggregate function.

# Example

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Group(Get.year_of_birth) |>
           Select(Get.year_of_birth, Agg.count());

julia> print(render(q))
SELECT
  "person_1"."year_of_birth",
  COUNT(*) AS "count"
FROM "person" AS "person_1"
GROUP BY "person_1"."year_of_birth"
```

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Group() |>
           Select(Agg.count(distinct = true, Get.year_of_birth));

julia> print(render(q))
SELECT COUNT(DISTINCT "person_1"."year_of_birth") AS "count"
FROM "person" AS "person_1"
```

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id]);

julia> visit_occurrence =
           SQLTable(:visit_occurrence, columns = [:visit_occurrence_id, :person_id, :visit_start_date]);

julia> q = From(person) |>
           LeftJoin(:visit_group => From(visit_occurrence) |>
                                    Group(Get.person_id),
                    on = (Get.person_id .== Get.visit_group.person_id)) |>
           Select(Get.person_id,
                  :max_visit_start_date =>
                      Get.visit_group |> Agg.max(Get.visit_start_date));

julia> print(render(q))
SELECT
  "person_1"."person_id",
  "visit_group_1"."max" AS "max_visit_start_date"
FROM "person" AS "person_1"
LEFT JOIN (
  SELECT
    "visit_occurrence_1"."person_id",
    MAX("visit_occurrence_1"."visit_start_date") AS "max"
  FROM "visit_occurrence" AS "visit_occurrence_1"
  GROUP BY "visit_occurrence_1"."person_id"
) AS "visit_group_1" ON ("person_1"."person_id" = "visit_group_1"."person_id")
```
"""
Agg(args...; kws...) =
    AggregateNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Agg), pats::Vector{Any}) =
    dissect(scr, AggregateNode, pats)

function PrettyPrinting.quoteof(n::AggregateNode, ctx::QuoteContext)
    ex = Expr(:call,
              Expr(:., nameof(Agg),
                   QuoteNode(Base.isidentifier(n.name) ? n.name : string(n.name))))
    if n.distinct
        push!(ex.args, Expr(:kw, :distinct, n.distinct))
    end
    append!(ex.args, quoteof(n.args, ctx))
    if n.filter !== nothing
        push!(ex.args, Expr(:kw, :filter, quoteof(n.filter, ctx)))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::AggregateNode) =
    n.name

rebase(n::AggregateNode, n′) =
    AggregateNode(over = rebase(n.over, n′),
                  name = n.name,
                  distinct = n.distinct,
                  args = n.args,
                  filter = n.filter)

# Notation for making aggregate nodes.

struct AggClosure
    name::Symbol
end

Base.show(io::IO, f::AggClosure) =
    print(io, Expr(:., nameof(Agg),
                       QuoteNode(Base.isidentifier(f.name) ? f.name : string(f.name))))

Base.getproperty(::typeof(Agg), name::Symbol) =
    AggClosure(name)

Base.getproperty(::typeof(Agg), name::AbstractString) =
    AggClosure(Symbol(name))

(f::AggClosure)(args...; over = nothing, distinct = false, filter = nothing) =
    Agg(over = over, name = f.name, distinct = distinct, args = SQLNode[args...], filter = filter)

(f::AggClosure)(; over = nothing, distinct = false, args = SQLNode[], filter = nothing) =
    Agg(over = over, name = f.name, distinct = distinct, args = args, filter = filter)


text:
using HDF5

@defstruct HDF5DataLayer Layer (
  name :: AbstractString = "hdf5-data",
  (source :: AbstractString = "", source != ""),
  (batch_size :: Int = 0, batch_size > 0),
  (tops :: Vector{Symbol} = Symbol[:data,:label], length(tops) > 0),
  shuffle :: Bool = false,
  transformers :: Vector = [],
)
@characterize_layer(HDF5DataLayer,
  is_source => true
)

type HDF5DataLayerState <: LayerState
  layer :: HDF5DataLayer
  blobs :: Vector{Blob}
  epoch :: Int
  trans :: Vector{Vector{DataTransformerState}}

  sources        :: Vector{AbstractString}
  dsets          :: Vector{Any}
  curr_source    :: Int
  curr_hdf5_file :: HDF5File
  curr_index     :: Int
  shuffle_idx    :: Vector{Int}
  shuffle_src    :: Vector{Int}

  HDF5DataLayerState(backend::Backend, layer::HDF5DataLayer) = begin
    state = new(layer)

    sources = open(layer.source, "r") do s
      map(strip, filter(l -> !isspace(l), readlines(s)))
    end
    @assert(length(sources) > 0)
    state.sources = sources

    if layer.shuffle
      state.shuffle_src = randperm(length(state.sources))
    else
      state.shuffle_src = collect(1:length(state.sources))
    end

    state.epoch = 0
    state.curr_source = 1
    state.curr_hdf5_file = h5open(sources[state.shuffle_src[1]], "r")

    state.dsets = [state.curr_hdf5_file[string(x)] for x in layer.tops]
    if layer.shuffle
      # we only use mmap when random shuffle is required, as mmap is slower
      # than sequential batch read. see hdf5-read benchmark.
      state.dsets = map(readmmap, state.dsets)
      state.shuffle_idx = randperm(size(state.dsets[1])[end])
    else
      state.shuffle_idx = Int[]
    end

    state.blobs = Array(Blob, length(layer.tops))
    state.trans = Array(Vector{DataTransformerState}, length(layer.tops))
    transformers = convert(Vector{@compat(Tuple{Symbol, DataTransformerType})}, layer.transformers)
    for i = 1:length(state.blobs)
      dims = size(state.dsets[i])

      dims = tuple(dims[1:end-1]..., layer.batch_size)

      dset = state.dsets[i]
      state.blobs[i] = make_blob(backend, eltype(dset), dims)

      state.trans[i] = [setup(backend, convert(DataTransformerType, t), state.blobs[i])
          for (k,t) in filter(kt -> kt[1] == layer.tops[i], transformers)]
    end
    state.curr_index = 1

    return state
  end
end

function setup(backend::Backend, layer::HDF5DataLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  @assert length(inputs) == 0
  state = HDF5DataLayerState(backend, layer)
  return state
end
function shutdown(backend::Backend, state::HDF5DataLayerState)
  map(destroy, state.blobs)
  map(ts -> map(t -> shutdown(backend, t), ts), state.trans)
  close(state.curr_hdf5_file)
end

function forward(backend::Backend, state::HDF5DataLayerState, inputs::Vector{Blob})
  n_done = 0
  while n_done < state.layer.batch_size
    n_remain = size(state.dsets[1])[end] - state.curr_index + 1
    if n_remain == 0
      close(state.curr_hdf5_file)
      state.curr_source = state.curr_source % length(state.sources) + 1
      state.curr_hdf5_file = h5open(state.sources[state.shuffle_src[state.curr_source]], "r")
      state.dsets = [state.curr_hdf5_file[string(x)] for x in state.layer.tops]

      if state.layer.shuffle
        state.dsets = map(readmmap, state.dsets)
        state.shuffle_idx = randperm(size(state.dsets[1])[end])
      end

      state.curr_index = 1
      n_remain = size(state.dsets[1])[end]
    end

    n1 = min(state.layer.batch_size-n_done, n_remain)
    if n1 > 0
      for i = 1:length(state.blobs)
        idx = map(x -> 1:x, size(state.blobs[i])[1:end-1])
        dset = state.dsets[i]
        if state.layer.shuffle
          the_data = dset[idx..., state.shuffle_idx[state.curr_index:state.curr_index+n1-1]]
        else
          the_data = dset[idx..., state.curr_index:state.curr_index+n1-1]
        end
        set_blob_data(the_data, state.blobs[i], n_done+1)
      end
    end
    state.curr_index += n1
    n_done += n1

    # update epoch
    if state.curr_index > size(state.dsets[1])[end] &&
        state.curr_source == length(state.sources)
      state.epoch += 1
    end
  end

  for i = 1:length(state.blobs)
    for j = 1:length(state.trans[i])
      forward(backend, state.trans[i][j], state.blobs[i])
    end
  end
end

function set_blob_data(data::Array, blob::CPUBlob, blob_idx::Int)
  n_fea = get_fea_size(blob)
  idx_start = (blob_idx-1)*n_fea
  blob.data[idx_start+1:idx_start+length(data)] = data
end

function backward(backend::Backend, state::HDF5DataLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end


text:
# https://adventofcode.com/2021/day/3
using AdventOfCode
using Statistics
input = (x->parse.(Int,split(x,""))).(readlines("2021/data/day_3.txt"))
function parsebin(binvec::BitArray)
    sum(binvec.*2 .^((length(binvec)-1:-1:0)))
end
parsebin(vec::Vector{<:Int})=parsebin(Bool.(vec))
function part_1(input)
    gammavec= >(0.5).(mean(input))
    parsebin(gammavec)*
        parsebin(.!gammavec)
end

function sift(vectors, meancond, sumvec, digit=1)
    len=length(vectors)
    if len==1 return(parsebin(vectors[1])) end
    filter!(x->x[digit]==meancond(sumvec[digit]/len),vectors)
    sumvec.=sum(vectors)
    sift(vectors,meancond,sumvec,digit+1)
end
function part_2(input)
    sift(copy(input),≥(0.5),sum(input))*
        sift(copy(input),<(0.5),sum(input))
end


text:
#####
##### generate boundary conditions for each nutrient field
#####

# struct to combine the left and right boundary conditions in one direction
mutable struct BoundaryConditions
    west::Union{Nothing, Number, AbstractArray}
    east::Union{Nothing, Number, AbstractArray}
    north::Union{Nothing, Number, AbstractArray}
    south::Union{Nothing, Number, AbstractArray}
    top::Union{Nothing, Number, AbstractArray}
    bottom::Union{Nothing, Number, AbstractArray}
end

# generate the default boundary conditions for one field
function default_bcs()
    bcs = BoundaryConditions(nothing, nothing, nothing, nothing, nothing, nothing)
    return bcs
end

"""
    set_bc!(model, tracer::Symbol, pos::Symbol, bc_value::Union{Number, AbstractArray})
Set the boundary condition of `tracer` on `pos` with `bc_value`.
"""
function set_bc!(model, tracer::Symbol, pos::Symbol, bc_value::Union{Number, AbstractArray})
    @assert tracer in nut_names

    bc_value_d = bc_value
    if isa(bc_value, AbstractArray)
        bc_value_d = vc_value |> array_type(model.arch)
    end
    setproperty!(model.nutrients[tracer].bc, pos, bc_value_d)
end

# get boundary condition at each grid point
@inline getbc(bc::Number, i, j, t) = bc
@inline getbc(bc::AbstractArray{Float64,2}, i, j, t) = bc[i,j]
@inline getbc(bc::AbstractArray{Float64,3}, i, j, t) = bc[i,j,t]

# validate boundary conditions, check if the grid information is compatible with nutrient field
function validate_bc(bc, bc_size, nΔT)
    if typeof(bc) <: AbstractArray{Float64,2}
        if size(bc) == bc_size
            return nothing
        else
            throw(ArgumentError("BC west: grid mismatch, size(bc) must equal to $(bc_size) for a constant flux boundary condition."))
        end
    end
    if typeof(bc) <: AbstractArray{Float64,3}
        if size(bc) == (bc_size..., nΔT)
            return nothing
        else
            throw(ArgumentError("BC west: grid mismatch, size(bc) must equal to $((bc_size..., nΔT)) for a time-dependent flux boundary condition."))
        end
    end
end

function validate_bcs(nut, grid, nΔT)
    for name in nut_names
        bc_west    = nut[name].bc.west
        bc_east    = nut[name].bc.east
        bc_south   = nut[name].bc.south
        bc_north   = nut[name].bc.north
        bc_bottom  = nut[name].bc.bottom
        bc_top     = nut[name].bc.top

        bc_size_x = (grid.Ny+2*grid.Hy, grid.Nz+2*grid.Hz)
        bc_size_y = (grid.Nx+2*grid.Hx, grid.Nz+2*grid.Hz)
        bc_size_z = (grid.Nx+2*grid.Hx, grid.Ny+2*grid.Hy)

        validate_bc(bc_west, bc_size_x, nΔT)
        validate_bc(bc_east, bc_size_x, nΔT)
        validate_bc(bc_south, bc_size_y, nΔT)
        validate_bc(bc_north, bc_size_y, nΔT)
        validate_bc(bc_bottom, bc_size_z, nΔT)
        validate_bc(bc_top, bc_size_z, nΔT)
    end
end


text:
include("required.jl")
include("functions.jl")

###  Load data, Clean and convert data to right types ###
df = CSV.read("WoodEtal_Append1_v2.csv")
edgelist_T = load_data(df, 1)
edgelist_Q = load_data(df, 2)
edgelist_A = load_data(df, 3)
edgelist_L = load_data(df, 4)
edgelist_S = load_data(df, 5)

# Rename IDS to start at 1
keys_T = renameIDs(edgelist_T)
keys_Q = renameIDs(edgelist_Q)
keys_A = renameIDs(edgelist_A)
keys_L = renameIDs(edgelist_L)
keys_S = renameIDs(edgelist_S)

# Compute adjacency matrix
M_T = adjacencyMatrix(keys_T)
M_Q = adjacencyMatrix(keys_Q)
M_A = adjacencyMatrix(keys_A)
M_L = adjacencyMatrix(keys_L)
M_S = adjacencyMatrix(keys_S)

# Compute LP/BRIM
Q_T, Qr_T, B_T = computeModularity(M_T)
Q_Q, Qr_Q, B_Q = computeModularity(M_Q)
Q_A, Qr_A, B_A = computeModularity(M_A)
Q_L, Qr_L, B_L = computeModularity(M_L)
Q_S, Qr_S, B_S = computeModularity(M_S)

# Compute spectral radius
ρ_T = computeNestedness(B_T)
ρ_Q = computeNestedness(B_Q)
ρ_A = computeNestedness(B_A)
ρ_L = computeNestedness(B_L)
ρ_S = computeNestedness(B_S)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
import Base.show

mutable struct Neo4jInstance
  connection::Neo4j.Connection
  graph::Neo4j.Graph
end

mutable struct CloudGraphsDFG{T <: AbstractParams} <: AbstractDFG
    neo4jInstance::Neo4jInstance
    description::String
    userId::String
    robotId::String
    sessionId::String
    encodePackedTypeFunc
    getPackedTypeFunc
    decodePackedTypeFunc
    rebuildFactorMetadata!
    labelDict::Dict{Symbol, Int64}
    variableCache::Dict{Symbol, DFGVariable}
    factorCache::Dict{Symbol, DFGFactor}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
    useCache::Bool
end

function show(io::IO, c::CloudGraphsDFG)
    println(io, "CloudGraphsDFG:")
    println(io, " - Neo4J instance: $(c.neo4jInstance.connection.host)")
    println(io, " - Session: $(c.userId):$(c.robotId):$(c.sessionId)")
    println(io, " - Caching: $(c.useCache)")
end


text:
##  a bare bones Test-Suite for some critical functions
using Test, Distributions, Random, LinearAlgebra

include("../src/JS_SAA_main.jl")
include("consts.jl")

@testset "All Tests" begin

#Tests shrink/crossval other simpler helpers
@testset "Vector Functions Baby Newsvendor" begin
	Random.seed!(8675309)
	K = 100
	ps = rand(K)
	alpha = .1
	s = .9
	p0 = .5
	Nhat = 20

	@test isapprox(JS.exp_nv_obj(ps, p0, alpha, Nhat, s),
					0.4670605847641166)
	@test isapprox(JS.exp_nv_fullInfo(ps, s),
					0.4472649361852286)
end

@testset "bestAlpha Functions" begin
	Random.seed!(8675309)
	K = 100
	ps = rand(K)
	Nhats = rand(Poisson(20), K)
	alpha_grid = range(0, stop=80, length=80)
	ss = rand(K)/4 .+ .75
	p0 = .5
	mhats = JS.nv_sim_path(ps, Nhats)
	alphaOR, jstar, out = JS.nv_oracle_alpha(mhats, ps, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 0.4323864805448839)
	@test isapprox(alphaOR, 2.0253164556962027)

	alphaAP, jstar, out = JS.nv_apriori_alpha(ps, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 0.4338668678652777)
	@test isapprox(alphaAP, 2.025316455696202)

	alphaLOO, jstar, out = JS.nv_loo_alpha(mhats, p0, alpha_grid, Nhats, ss)
	@test isapprox(out[jstar], 7.926847617453273)
	@test isapprox(alphaLOO, 1.0126582278481013)

	p0_grid = range(0, stop=1, length=20)
	p0OR, alphaOR, jstar = JS.nv_oracle_both(mhats, ps, p0_grid, alpha_grid, Nhats, ss)
	@test isapprox(p0OR, 0.21052631578947367)
	@test isapprox(alphaOR, 14.177215189873417)
	@test jstar == 15

	p0LOO, alphaLOO, jstar = JS.nv_loo_both(mhats, p0_grid, alpha_grid, Nhats, ss)
	@test isapprox(p0LOO, 0.15789473684210525)
	@test isapprox(alphaLOO, 11.139240506329115)
	@test jstar == 12

end

#Needs a function to test gen Newsvendors

end #all tests

max_stars_repo_name: CausalML/DataPooling

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
module TimingDiagrams

using ..Plots: plot, plot!, annotate!, Shape, text
using LinearAlgebra

export timingdiagram, node!, tx!

struct Node
  name::String
  pos::Vector{Float64}
  ypos::Float64
end

const nodes = Node[]
const speed = Ref(1.0)
const duration = Ref(0.0)
const timings = Float64[]

const theme = Dict{Symbol,Any}(
  :txheight => 10.0,
  :fontsize => 8,
  :ypad => 5.0,
  :tpad => 0.01,
  :txcolor => :blue,
  :rxcolor => :green,
  :linecolor => :lightgray,
  :nodecolor => :black
)

speed!(c) = (speed[] = c)
duration!(d) = (duration[] = d)

function timingdiagram(f; kwargs...)
  plt = plot(; border=:none, legend=false, yflip=true, kwargs...)
  empty!(nodes)
  empty!(timings)
  f()
  for n ∈ nodes
    plot!(plt, [minimum(timings) - theme[:tpad], maximum(timings) + theme[:tpad]], [n.ypos, n.ypos]; color=theme[:nodecolor])
    annotate!(plt, -theme[:tpad], n.ypos, text(n.name * " ", theme[:fontsize], :right))
  end
  empty!(nodes)
  empty!(timings)
  plt
end

function node!(name, pos; ypos=norm(pos))
  push!(nodes, Node(name, pos isa Real ? [Float64(pos)] : pos, ypos))
  length(nodes)
end

function tx!(t, from, to=setdiff(eachindex(nodes),from); duration=duration[], txlabel=nothing, txtiminglabel=nothing, rxlabels=nothing, rxtiminglabels=nothing)
  rxtimes = [t + norm(nodes[i].pos - nodes[from].pos) / speed[] for i ∈ to]
  duration > 0.0 && rect!(t, nodes[from].ypos, duration, theme[:txheight]; color=theme[:txcolor])
  txlabel === nothing || annotate!(t + duration/2, nodes[from].ypos - theme[:txheight] - theme[:ypad], text(txlabel, theme[:fontsize], :bottom))
  txtiminglabel === nothing || annotate!(t, nodes[from].ypos + theme[:ypad], text(txtiminglabel, theme[:fontsize], :top))
  for (i, rx) ∈ enumerate(to)
    plot!([t, rxtimes[i]], [nodes[from].ypos, nodes[rx].ypos]; color=theme[:linecolor])
    duration > 0.0 && rect!(rxtimes[i], nodes[rx].ypos, duration, theme[:txheight]; color=theme[:rxcolor])
    rxlabels === nothing || annotate!(rxtimes[i] + duration/2, nodes[rx].ypos - theme[:txheight] - theme[:ypad], text(rxlabels[i], theme[:fontsize], :bottom))
    rxtiminglabels === nothing || annotate!(rxtimes[i], nodes[rx].ypos + theme[:ypad], text(rxtiminglabels[i], theme[:fontsize], :top))
    end
  push!(timings, t)
  push!(timings, minimum(rxtimes))
  push!(timings, maximum(rxtimes) + duration)
  rxtimes
end

function rect!(t, y, dt, dy; color)
  plot!(Shape([t, t+dt, t+dt, t], [y, y, y-dy, y-dy]); color)
end

end # module


text:
struct BoxDomain{N} <: Domain{N}
    size::NTuple{N,Float64} # in m
    γ::Float64 # global admittance (loss at boundary)
end

function BoxDomain(size::NTuple{N,Number}; gamma=0) where N
    BoxDomain{N}(Float64.(size), Float64(gamma))
end

function BoxDomain(size...; gamma=0)
    BoxDomain(Float64.(size), Float64(gamma))
end

Base.size(box::BoxDomain) = box.size
gridsize(domain::Domain, dx) = ceil.(Int, size(domain) ./ dx)

mutable struct BoxState{N,T<:AbstractFloat,D,A<:AbstractArray{T,N},Q<:AbstractArray} <: State{N,T,D}
    box::D
    previous::A
    current::A
    q::Q  # domain indicator grid
    t::T  # accumulated time
    iter::Int # accumulated simulation steps
end

function state_init(f0, backend::Backend, box::BoxDomain{N}, sim::Simulator{N,T}) where {N,T}
    dims = gridsize(box, sim.wave.dx)
    A = zeros(T, dims)
    for I in CartesianIndices(ntuple(i->2:(dims[i]-1), N))
        A[I] = T(f0(((I.I.-1).*sim.wave.dx)...))
    end
    Ψ₀ = move_backend(A, backend)
    Ψ₁ = move_backend(copy(A), backend)
    qhost = zeros(Int8, dims)
    for I in CartesianIndices(ntuple(i->2:(size(qhost)[i]-1), N))
        qhost[I] = Int8(1)
    end
    q = move_backend(qhost, backend)
    BoxState{N,T,typeof(box),typeof(Ψ₀),typeof(q)}(box,Ψ₀,Ψ₁,q,zero(T),0)
end

function state_init(backend::Backend, box::BoxDomain{N}, sim::Simulator{N,T}) where {N,T}
    dims = gridsize(box, sim.wave.dx)
    Ψ₀ = move_backend(zeros(T, dims), backend)
    Ψ₁ = move_backend(zeros(T, dims), backend)
    qhost = zeros(Int8, dims)
    for I in CartesianIndices(ntuple(i->3:(size(qhost)[i]-1), N))
        qhost[I] = Int8(1)
    end
    q = move_backend(qhost, backend)
    BoxState{N,T,typeof(box),typeof(Ψ₀),typeof(q)}(box,Ψ₀,Ψ₁,q,zero(T),0)
end


text:
__precompile__()

module Agtor

using Parameters
using DataFrames, Dates, YAML, CSV
using Distributed

import OrderedCollections: LittleDict
import DataStructures: OrderedDict, DefaultDict

include("AgBase/Units.jl")
include("AgBase/Component.jl")
include("AgBase/Parameter.jl")
include("AgBase/io.jl")
include("AgBase/properties.jl")
include("AgBase/ResultsHandler.jl")
include("Infrastructure/Infrastructure.jl")
include("Infrastructure/Irrigation.jl")
include("Infrastructure/Pump.jl")
include("Climate.jl")
include("WaterSource.jl")
include("Crop.jl")
include("Field.jl")
include("Zone.jl")
include("Manager.jl")
include("RigidManager.jl")
include("Basin.jl")
include("Model.jl")

AgUnion = Union{Int64, Float64, Agtor.AgParameter}

const component_sep = "└──"
const component_name_sep = "-"
const component_field_sep = "~"

export @def, AgUnion

export load_yaml, load_spec, load_climate, create
export generate_agparams, create, reset!
export AgComponent, AgParameter, Infrastructure, Irrigation, Pump
export WaterSource, Crop, FarmField, CropField, FarmZone
export Manager, BaseManager, RigidManager
export Climate, Basin, run_timestep!, subtotal_costs, total_costs, update_available_water!
export collect_results, min_max, extract_values, param_info, extract_spec, add_prefix!
export set_params!, extract_agparams, collect_agparams!, collect_agparams, update_model!
export water_used_by_source
export get_data_for_timestep, aggregate_field_logs, run_model
export run_scenarios!, save_results!, save_state!
export collate_results, collate_results!, component_relationship
export select_results, scenario_stats, results_df

end # module

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
export pi_axis_labels, set_pi_axis_labels

function pi_axis_labels(min, max, step = 1//2)
    first, last = round.([min, max] / (pi * step))
    locs = Float64[]
    labels = String[]
    for k in first:last
        factor = Integer(k) * step
        push!(locs, factor * π)

        sign = factor < 0 ? "-" : ""

        abs_num = abs(numerator(factor))
        num_string =
            if abs_num == 0
                "0"
            elseif abs_num == 1
                "\\pi"
            else
                "$abs_num \\pi"
            end

        if denominator(factor) == 1
            push!(labels, "\$ $sign $num_string \$")
        else
            push!(labels, "\$ $sign \\frac{$num_string}{$(denominator(factor))} \$")
        end
    end

    locs, labels
end

"""
Places ticks and labels on an axis in given multiples of π.
"""
function set_pi_axis_labels(axis, min, max, step = 1//2)
    locs, labels = pi_axis_labels(min, max, step)
    axis.set_ticks(locs)
    axis.set_ticklabels(labels)
end


text:
using ModelConstructors, FileIO, Random, SMC

### Estimate a single factor CAPM model
# R_{it} = α_i + β_i R_{Mt} + ϵ_{it}, i = 1,...,N; t = 1,...,T
# where R_{Mt} is the excess return on a market index in time period t,
# and ϵ_{it} is an i.i.d. normally distributed mean zero shock with variance σ_i^2

### Construct a generic model and populate it with parameters
capm = GenericModel()
fn = dirname(@__FILE__)
capm <= Setting(:dataroot, "$(fn)/../save/input_data/")
capm <= Setting(:saveroot, "$(fn)/../save/")

capm <= parameter(:α1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:β1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:σ1, 1., (1e-5, 1e5), (1e-5, 1e5), SquareRoot(), Uniform(0, 1e3),
                  fixed = false)
capm <= parameter(:α2, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:β2, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:σ2, 1., (1e-5, 1e5), (1e-5, 1e5), SquareRoot(), Uniform(0, 1e3),
                  fixed = false)
capm <= parameter(:α3, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:β3, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                  fixed = false)
capm <= parameter(:σ3, 1., (1e-5, 1e5), (1e-5, 1e5), SquareRoot(), Uniform(0, 1e3),
                  fixed = false)

### Estimate with SMC

## Get data
N = 3 # number of asset returns
lik_data = load("../../../save/input_data/capm.jld2", "lik_data")
market_data = load("../../../save/input_data/capm.jld2", "market_data")

## Construct likelihood function:
# likelihood function is just R_{it} ∼ N(α_i + β_i R_{Mt}, σ_i)
# parameters to estimate are α_i, β_i, σ_i
# data is a time series of individual factor returns
# use S&P 500 data and a returns on a couple stocks. To get into returns,
# just take the log of the prices and regress those on each other.
function likelihood_fnct(p, d)
    # we assume the ordering of (α_i, β_i, σ_i)
    Σ = zeros(N,N)
    α = Vector{Float64}(undef,N)
    β = Vector{Float64}(undef,N)
    for i in 1:N
        α[i]   = p[i * 3 - 2]
        β[i]   = p[i * 3 - 2]
        Σ[i,i] = p[i * 3]^2
    end
    det_Σ = det(Σ)
    inv_Σ = inv(Σ)
    term1 = -N / 2 * log(2 * π) - 1 /2 * log(det_Σ)
    logprob = 0.
    errors = d .- α .- β .* market_data
    for t in 1:size(d,2)
        logprob += term1 - 1/2 * dot(errors, inv_Σ * errors)
    end
    return exp(logprob)
end

Random.seed!(1793)
println("Starting to estimate CAPM with SMC . . .")
@everywhere using SMC, OrderedCollections
smc(likelihood_fnct, capm.parameters, lik_data)


text:
function paformat(β;scalefactor, decimals, prefix, suffix)
  valstring = num2str(β, decimals; scalefactor)
  return "$(prefix)$(valstring)$(suffix)"
end

pabetaformat =(; scalefactor=1.0, decimals=2, prefix="", suffix="\\:\\:")
pabetax2format =(; scalefactor=2.0, decimals=2, prefix="", suffix="\\:\\:")
papercentageformat =(; scalefactor=100.0, decimals=2, prefix="", suffix="\\:\\:")
paseformat = (;scalefactor=1.0, decimals=2, prefix="\\mathit{(", suffix=")}")
pasex2format = (;scalefactor=1.0, decimals=2, prefix="\\mathit{(", suffix=")}")
paNformat = (;scalefactor=1.0, decimals=0, prefix="", suffix="")
par2format = (;scalefactor=100.0, decimals=2, prefix="", suffix="")

#put all of the field information here
pafields = Dict{Symbol, Any}(
  :rhs=>Dict(
    :A_WXmomentum121s=>(;
      id=:A_WXmomentum121s,
      label="\$A_0\$ (estimated) L+S",
      format=pabetax2format,
      se=(;
        id=:A_WXmomentum121s_se,
        label="",
        format=pasex2format)),
    :A_WXmomentum121s0=>(;
      id=:A_WXmomentum121s0,
      label="\$A_0\$ (\$G=1 \\: \\forall t\$) L+S",
      format=pabetax2format,
      se=(;
        id=:A_WXmomentum121s0_se,
        label="",
        format=pasex2format)),
    :intercept=>(;
      id=:intercept,
      label="intercept",
      format=pabetaformat,
      se=(;
        id=:intercept_se,
        label="",
        format=paseformat)),
    ),
    :description => Dict(
      :permno=>(;
        id=:permno,
        label="firm FE"),
      :date=>(;
        id=:date,
        label="time FE"),
      :N=>(;
        id=:N,
        label="N",
        format=paNformat),
      :R2 => (
        ;
        id=:R2,
        label = "\\midrule \$R^2\$",
        format = par2format),
      :R2delta => (
        ;
        id=:R2delta,
        label = "\$R^2 - R^2(mean)\$",
        format = par2format),
    )
  )

  #start of a table spec
function patablefull()
  Fmeasures=[
    [:intercept],
    #[],
    [:A_WXmomentum121s, :intercept],
    #[:A_WXmomentum121s],
    #[:A_WXmomentum121s],
    #[:A_WXmomentum121s],
    [:A_WXmomentum121s0, :intercept],
    #[:A_WXmomentum121s0],
    #[:A_WXmomentum121s0],
    #[:A_WXmomentum121s0],
    ]
  tablename="$(PARAM[:crspdatalabel])_pafull"
  Ffes=[
    Symbol[],
    #[:date],
    Symbol[],
    #[:date],
    #[:permno],
    #[:permno, :date],
    Symbol[],
    #[:date],
    #[:permno],
    #[:permno, :date]
    ]

  ses = [()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    ()->Vcov.cluster(:date,:permno),
    ()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    #()->Vcov.cluster(:date,:permno),
    ]
  Fvols = [:dvolscaled for i ∈ 1:length(Fmeasures)]
  Fse=:se

  Fdescriptions = [:R2, :R2delta, :N, #=:date, :permno=#]
  Frowheaders = [:A_WXmomentum121s; :A_WXmomentum121s0; :intercept; Fdescriptions;]
  alignmentstring="l | rrr" #"l | r | rrr | rrr"

  #verify integrity here for the sake of debugging
  allmeasures = unique([Fmeasures...;])
  @assert length(Fmeasures) == length(Ffes)
  @assert length(Fvols) == length(Ffes)
  @assert setdiff(allmeasures, Frowheaders) |> isempty

  #build the row list
  Frows = Symbol[]
  for Frowheader ∈ Frowheaders
    push!(Frows, Frowheader)
    @assert haskey(pafields[:rhs], Frowheader) || haskey(pafields[:description], Frowheader
      ) "field not found in pafield: Frowheader=$Frowheader"

    if haskey(pafields[:rhs], Frowheader)
      push!(Frows, Symbol(Frowheader, :_, Fse))
      @assert pafields[:rhs][Frowheader].se.id ≡ Symbol(Frowheader, :_, Fse) "
        could not find $(Symbol(Frowheader, :_, Fse)) ∈ Frowheader $(pafields[:rhs][Frowheader])"
    end
  end

  return (; Fmeasures, Frows, tablename, Ffes, ses, Fvols, Fdescriptions, Fse, alignmentstring)
end

#run it with the Gs and without, cross-sectionally as well
#include standard errors- maybe fixed effects?
function ivoltests(panel, results, ms::MeasureSpec, returns)

  #get the field names
  Fξs = ms.Fξs
  Fw::Vector{Symbol} = Fweights(ms, :Fw)
  FRtw::Vector{Symbol} = Fweights(ms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(ms, :FRLw)
  FW::Vector{Symbol} = ms.FW
  FWgroup = Fgroupedcontrols(ms)

  Fcoefs::Vector{Symbol} = findcols(results, :asset)
  Fgcoefs::Vector{Symbol} = findcols(results, :growth)
  Fcoefs0 = Fcoefs .|> (Fcoef->Symbol(Fcoef, 0))

  #construct the model parts using cpu only data types and no controls
  p = modelparts(panel, Float64, Vector{Float64}, Matrix{Float64},
    weightdata = (;Fw, FRtw, FRLw), Fvol=ms.Fvol,
    ;FW, FWgroup, ms)

  dims = p.dims
  ws = p.dat[:Fw]
  Rtws = p.dat[:FRtw]
  RLws = p.dat[:FRLw]
  W = p.dat[:FW]
  Wgroup = p.dat[:FWgroup]
  ts::Vector{Int} = p.ts
  v = p.dat[:Fvol]
  permno = p.dat[:permno]
  dates = p.dat[:date]
  cleanpanel = p.cleanpanel

  @assert (dates .== (t->results.date[t]).(ts)) |> all

  #construct simplest possible structures for computing predicted volumes
  @info "Computing raw controls for output file"
  Xvraw = XYIterControl(ws, Rtws, RLws, v, ts) #raw implies no controls or transformations
  Θ = VolumePartsIter(dims.T, dims.K, ts, Float64, Matrix{Float64}, Vector{Float64})

  A = Matrix{Float64}(results[!, Fcoefs])
  G = Matrix{Float64}(results[2:end, Fgcoefs])
  #now follow the sequence for estimating the volume, minus the actual estimation
  someones = ones(Float64, dims.K, 1)
  prodmatG = cumprodbyrow(G' |> Matrix{Float64})
  Ã = hcat(someones, prodmatG) |> Matrix{Float64}

  RHSraw::Matrix{Float64} = projectbyt(Ã, Xvraw.xsection) #no controls or transofmrations
  RHSraw0::Matrix{Float64} = projectbyt(ones(Float64, size(Ã)), Xvraw.xsection)

  #form the sample for regressions
  expanded::Matrix{Float64}, controlsused = expandcontrolsbyt(Wgroup, ts) |> (nt
    )->(nt.expanded, nt.Wgroupcolsused)
  focal = hcat(
    DataFrame(date=dates, permno=permno),
    DataFrame(RHSraw, Fcoefs),
    DataFrame(RHSraw0, Fcoefs0))
  focal[!, PARAM[:Fvol]] = v

  #construct intercepts by hand
  focal.intercept = ones(length(v))

  ############################
  patable(focal; patablefull()...)
  pax = paxsectional(focal; Fcoefs,Fcoefs0)
  results = leftjoin(results, pax, on=:date, validate=(true,true))
  results = sort!(results, :date)

  return results
end

function patable(focal;
    tablename,
    Fmeasures,
    Fvols::TFvols,
    Fse,
    ses,
    Ffes,
    Frows,
    Fdescriptions,
    verbose = PARAM[:tabpaverbose],
    tablepath = PARAM[:tablepath],
    colnames=[1:length(Fmeasures) .|> i->"($i)",],
    widthcolnames= [ones(Int,length(Fmeasures)),],
    alignmentcolnames = [
      (rc-> "c").(colnamerow) for (r,colnamerow) ∈ enumerate(colnames)],
    alignmentstring = string("l",["r" for i ∈ 1:length(Fmeasures)]...),

    #construct the info for each row
    #first extract lightly constructed field objects IN ORDER
    datafields = reduce(vcat, (reduce(vcat, Fmeasures) |> unique) .|> (
      (s)->[pafields[:rhs][s], pafields[:rhs][s][Fse]])),
    descriptionfields = Fdescriptions .|> ((s)->pafields[:description][s])

    ) where {TFvols}

  #create an index of row information. order is not preserved.
  Frowindex = Dict(f.id=>f for f ∈ [datafields; descriptionfields])

  #double check alignment
  @assert (id->Frowindex[id].id≡id).(Frows) |> all

  #extract the row labels, using the order preserved in Frowid
  descrownames = (id->Frowindex[id].label).(Frows)

  #form the regression formulai
  formulai = FormulaTerm[]
  fmlmrhs = Any[]
  for (Fmeasure, Ffe, Fvol) ∈ zip(Fmeasures, Ffes, Fvols)

    termrhs = term.([Fmeasure; fe.(Ffe)]) |> sum
    push!(fmlmrhs, Meta.parse(join((string).([Fmeasure; Ffe]),"+") * "+ 0"))

    if :intercept ∈ Fmeasure
      @assert Ffe |> isempty "Intercept detected in fixed effects spec!!"
      termrhs += term(0) # do not want an itnercept if we are doing the regressions by hand
    else
      @assert Ffe |> !isempty "No intercept detected without FE!"
    end

    push!(formulai, term(Fvol) ~ termrhs)
  end

  reg(f,err) = FixedEffectModels.reg(focal, f, err())

  #run the full regressions
  femlms = reg.(formulai, ses)

  verbose && println.(femlms)

  allFfes = reduce(vcat, Ffes) |> unique
  regressioncols = DataFrame([id=>String[] for id ∈ Frows])

  #make sure that the first column si the baseline
  @assert (length(Fmeasures[1]) == 1) && (Fmeasures[1][1] ≡ :intercept)
  r2base = r2(femlms[1])
  Ndates = length(unique(focal.date))
  for (c, femlm) ∈ enumerate(femlms)
    #get the main coefficients and SEs
    coefindex = Dict(Symbol(n) => A for (n,A) ∈ zip(coefnames(femlm), femlm.coef))
    tindex = Dict(Symbol(n, :_, Fse)=>A/σ for (n, A, σ) ∈ zip(
      coefnames(femlm), femlm.coef, vcov(femlm) |> (Σ)->diag(Σ).^0.5))
    regressioncolvals = merge(coefindex, tindex)

    #check the coefficients
    Ncoef = Fmeasures[c] |> length

    #= TOO SLOWWWWWWW
    fmlm = FMLM(focal, fmlmrhs[c], Fvols[c],
      #=withinsym=length(Ffes[c]) > 0 ? Ffes[c][1] : nothing,=# clustersyms=[:permno, :date])
    fmlmdf = DataFrame(coef=Fmeasures[c])
    fmlmdf.beta = fmlm.β[1:Ncoef]
    fmlmdf.femlmbeta = (F -> coefindex[F]).(Fmeasures[c])
    fmlmdf.femlmse = (F -> coefindex[F]/tindex[Symbol(F, :_, Fse)]).(Fmeasures[c])
    fmlmdf.robust = diag(modifiedwhiteΣ!(fmlm))[1:Ncoef].^0.5
    fmlmdf.clustered = diag(clusteredΣ!(fmlm))[1:Ncoef].^0.5
    fmlmdf.hac2way = diag(neweywestΣ!(fmlm, nrow(focal)^0.33333 |> ceil |> Int))[1:Ncoef].^0.5

    fmlm1way = FMLM(focal, fmlmrhs[c], Fvols[c],
      #=withinsym=length(Ffes[c]) > 0 ? Ffes[c][1] : nothing,=# clustersyms=[:permno])
    fmlmdf.panelhac1way = diag(neweywestpanelΣ!(fmlm1way, Ndates^0.33333 |> ceil |> Int))[1:Ncoef].^0.5

    @assert fmlmdf.beta ≈ fmlmdf.femlmbeta
    @info "$(Fvols[c]) ~ $(fmlmrhs[c])"
    println(fmlmdf)=#

    #format the values
    regressioncol = Dict(k =>paformat(regressioncolvals[k]; Frowindex[k].format...)
      for k ∈ keys(regressioncolvals))

    Fempties = setdiff(Frows,
      [Fdescriptions; Fmeasures[c]; (n->Symbol(n, :_, Fse)).(Fmeasures[c])])
    regressioncol = merge(regressioncol, Dict(F=>"" for F ∈ Fempties))

    #now do fixed effects
    regressioncol = merge(regressioncol,Dict(F=>F ∈ Ffes[c] ? "X" : "" for F ∈ allFfes))

    regressioncol[:R2] = paformat(r2(femlm); Frowindex[:R2].format...)
    regressioncol[:R2delta] = paformat(r2(femlm) - r2base; Frowindex[:R2delta].format...)
    regressioncol[:N] = paformat(
      sum(completecases(focal, [Fmeasures[c]; Ffes[c];])); Frowindex[:N].format...)

    #@info "propertynames: $(propertynames(regressioncols))"
    #@info "row keys: $(keys(regressioncol))"
    push!(regressioncols, regressioncol)
  end

  @assert (propertynames(regressioncols) .≡ Frows) |> all
  desccontent = eachcol(regressioncols) |> collect .|> Vector

  t = textable(;colnames, descrownames, desccontent,
    widthcolnames, alignmentstring, alignmentcolnames)

  writetextable(t, path=tablepath, outname="$tablename.tex")
end

function paxsectional(focal;
    Fcoefs,Fcoefs0, Fvol=PARAM[:Fvol],
    )

  Fcoefsx = Fcoefs .|> (Fcoef->Symbol(:x, Fcoef))
  Fcoefsxse = Fcoefs .|> (Fcoef->Symbol(:xse, Fcoef))
  Fcoefsx0 = Fcoefs0 .|> (Fcoef->Symbol(:x, Fcoef))
  Fcoefsx0se = Fcoefs0 .|> (Fcoef->Symbol(:xse, Fcoef))
  Fcoefsxd = Fcoefs .|> (Fcoef->Symbol(:xt, Fcoef))

  #this will hold the results
  xs = DataFrame(date=Date[])
  for F ∈ [Fcoefsx; Fcoefsxse; Fcoefsx0; Fcoefsxse]
    xs[!, F] = MFloat64[]
  end

  #some housekeeping
  intercepts = [:intercept for i ∈ 1:length(Fcoefs)]
  focal.intercept0 = deepcopy(focal.intercept)
  intercepts0 = [:intercept0 for i ∈ 1:length(Fcoefs)]
  Fvols = [Fvol for i ∈ 1:length(Fcoefs)]

  #run xsectional regressions
  function runxsectional(t)
    f = term(Fvol) ~ sum(term.(Fcoefs)) + term(:intercept) + term(0)
    femlm = FixedEffectModels.reg(t, f, Vcov.robust())

    #run the regressions
    coefindex= Dict(Symbol(:x, n, ) => 2*A for (n,A) ∈ zip(coefnames(femlm), femlm.coef))
    seindex = Dict(Symbol(:xse, n, )=> 2*σ for (n, σ) ∈ zip(
      coefnames(femlm), vcov(femlm) |> (Σ)->diag(Σ).^0.5))

    #run the 0 regressions (note the 0 is already added)
    f0 = term(Fvol) ~ sum(term.(Fcoefs0)) + term(:intercept0) + term(0)
    femlm0 = FixedEffectModels.reg(t, f0, Vcov.robust())
    coefindex0= Dict(Symbol(:x, n, ) => 2*A for (n,A) ∈ zip(coefnames(femlm0), femlm0.coef))
    seindex0 = Dict(Symbol(:xse, n, )=> 2*σ for (n, σ) ∈ zip(
      coefnames(femlm0), vcov(femlm0) |> (Σ)->diag(Σ).^0.5))

    return Ref(merge(coefindex, coefindex0, seindex, seindex0))
  end

  xfocal = groupby(focal, :date)
  xresults = combine(xfocal, zip(Fcoefs, Fcoefs0, Fvols, intercepts, intercepts0) .|> collect .|> AsTable
    .=> runxsectional .=> Fcoefsxd)

  #splat all the gathered regression results
  for i ∈ 1:length(Fcoefs)
    transform!(xresults,
      Fcoefsxd => ((d)-> (di->di[Fcoefsx[i]]).(d)) => Fcoefsx[i],
      Fcoefsxd => ((d)-> ((di->di[Fcoefsxse[i]]).(d))) => Fcoefsxse[i],
      Fcoefsxd => ((d)-> ((di->di[Fcoefsx0[i]]).(d))) => Fcoefsx0[i],
      Fcoefsxd => ((d)-> ((di->di[Fcoefsx0se[i]]).(d))) => Fcoefsx0se[i])
  end
  transform!(xresults,
     Fcoefsxd => ((d)->(di->di[:xintercept]).(d))=>:xintercept,
     Fcoefsxd => ((d)-> (di->di[:xseintercept]).(d))=>:xseintercept,
     Fcoefsxd => ((d)-> (di->di[:xintercept0]).(d))=>:xintercept0,
     Fcoefsxd => ((d)-> (di->di[:xseintercept0]).(d))=>:xseintercept0)

  select!(xresults, Not(Fcoefsxd))
  return xresults
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using OutlierDetectionTest
using Test

@testset "OutlierDetectionTest.jl" begin
    # Make sure the import worked
end


text:
function get_column_dimension(obj::RealLinearOperator)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_row_dimension(obj::RealLinearOperator)
    return jcall(obj, "getRowDimension", jint, ())
end

function operate(obj::RealLinearOperator, arg0::RealVector)
    return jcall(obj, "operate", RealVector, (RealVector,), arg0)
end


text:
const gdt_lookup = Dict{DataType, GDAL.GDALDataType}(
    UInt8 => GDAL.GDT_Byte,
    UInt16 => GDAL.GDT_UInt16,
    Int16 => GDAL.GDT_Int16,
    UInt32 => GDAL.GDT_UInt32,
    Int32 => GDAL.GDT_Int32,
    Float32 => GDAL.GDT_Float32,
    Float64 => GDAL.GDT_Float64
)

function read(fn::AbstractString)
    # GDAL specific init
    GDAL.allregister()

    isfile(fn) || error("File not found.")
    dataset = ArchGDAL.unsafe_read(fn)
    A = ArchGDAL.read(dataset)
    am = get_affine_map(dataset)

    # nodata masking
    # A = Array{Union{Missing, eltype(A)}}(A)
    mask = falses(size(A))
    for i = 1:size(A)[end]
        band = ArchGDAL.getband(dataset, i)
        maskflags = mask_flags(band)

        # All values are valid, skip masking
        if :GMF_ALL_VALID in maskflags
            @debug "No masking"
            continue
        # Mask is valid for all bands
        elseif :GMF_PER_DATASET in maskflags
            @debug "Mask for each band"
            maskband = ArchGDAL.getmaskband(band)
            m = ArchGDAL.read(maskband) .== 0
            mask[:,:,i] = m
        # Alpha layer
        elseif :GMF_ALPHA in maskflags
            @warn "Dataset has band $i with an Alpha band, which is unsupported for now."
            continue
        # Nodata values
        elseif :GMF_NODATA in maskflags
            @debug "Flag NODATA"
            nodata = get_nodata(band)
            mask[:, :, i] = A[:,:,i] .== nodata
        else
            @warn "Unknown/unsupported mask."
        end
    end

    if any(mask)
        A = Array{Union{Missing, eltype(A)}}(A)
        A[mask] .= missing
    end

    # crs
    wkt = ArchGDAL.getproj(dataset)

    # GDAL specific cleanup
    ArchGDAL.destroy(dataset)
    GDAL.destroydrivermanager()
    GeoArray(A, am, wkt)
end

function write!(fn::AbstractString, ga::GeoArray, nodata=nothing)
    GDAL.allregister()
    shortname = find_shortname(fn)
    options = String[]
    w, h, b = size(ga)
    dtype = eltype(ga)
    data = copy(ga.A)
    use_nodata = false

    # Set compression options for GeoTIFFs
    shortname == "GTiff" && (options = ["COMPRESS=DEFLATE","TILED=YES"])

    # Slice data and replace missing by nodata
    if isa(dtype, Union) && dtype.a == Missing
        dtype = dtype.b
        nodata == nothing && (nodata = typemax(dtype))
        m = ismissing.(data)
        data[m] .= nodata
        data = Array{dtype}(data)
        use_nodata = true
    end

    ArchGDAL.create(fn, shortname, width=w, height=h, nbands=b, dtype=dtype, options=options) do dataset
        for i=1:b
            band = ArchGDAL.getband(dataset, i)
            ArchGDAL.write!(band, data[:,:,i])
            use_nodata && GDAL.setrasternodatavalue(band.ptr, nodata)
        end

        # Set geotransform and crs
        gt = affine_to_geotransform(ga.f)
        GDAL.setgeotransform(dataset.ptr, gt)
        GDAL.setprojection(dataset.ptr, ga.crs)

    end
    # GDAL specific cleanup
    GDAL.destroydrivermanager()
end


text:
using Lora

function ess(x::Matrix{Float64})
    (p, T) = size(x)

    ESS = zeros(p)

    for i = 1:p
        ESS[i] = T*var(x[i,:])/Lora.mcvar_imse(squeeze(x[i,:],1))
    end

    ESS
end

function ess(chain, varname)

    theta_chain = zeros(length(chain[1][1][varname]), length(chain))
    for i = 1:length(chain)
        theta = chain[i][1][varname]

        theta_chain[:,i] = theta[:]

    end

    ess(theta_chain)

end

function ref_autocor(chain, ref_mean, ref_var)
    (p, T) = size(chain)

    centered_chain = broadcast(-,chain, ref_mean)

    max_lag = min(T-1, ceil(Integer, 100log10(T)))
    acor = zeros(p,max_lag)
    for i = 1:p
        for lag = 1:max_lag
            A = centered_chain[i, 1:T-lag]*centered_chain[i, lag+1:T]'
            @assert length(A) == 1
            acor[i,lag] = A[1]/((T-lag)*ref_var[i])
        end
    end

    acor
end

function ref_autocov(chain, ref_mean, ref_var)
    (p, T) = size(chain)

    centered_chain = broadcast(-,chain, ref_mean)

    max_lag = min(T-1, ceil(Integer, 100log10(T)))
    acov = zeros(p,max_lag)
    for i = 1:p
        for lag = 1:max_lag
            A = centered_chain[i, 1:T-lag]*centered_chain[i, lag+1:T]'
            @assert length(A) == 1
            acov[i,lag] = A[1]/(T)
        end
    end

    ref_var + 2*sum(acov,2)

end

function ref_ess_f(chain, reference_chain, f)
    c = [ f(chain[:,i][:])::Float64 for j in 1:1, i in 1:size(chain,2)]
    r = [ f(reference_chain[:,i][:])::Float64 for j in 1:1, i in 1:size(reference_chain,2)]

    ref_ess(c,r)
end

function ref_ess(chain, reference_chain)
    (p,T) = size(chain)
    ref_mean = mean(reference_chain,2)
    ref_var = var(reference_chain,2)

    rho = ref_autocor(chain, ref_mean, ref_var)
    (p, M) = size(rho)
    threshhold_acor = [ minimum([M; find(rho[i,:] .< 0.05).-1]) for i = 1:p]

    ESS = zeros(p)

    one_sT = 1 - collect(1:M)/T

    for i = 1:p
        Mt = threshhold_acor[i]
        ESS[i] = T/(1+2*dot(one_sT[1:Mt],rho[i,1:Mt][:]))
    end

    centered_chain = broadcast(-, chain, ref_mean)
    moment_2 = centered_chain.^2

    ref_centered = broadcast(-, reference_chain, ref_mean)
    ref_moment_2 = ref_centered.^2

    ref_m2_mean = mean(ref_moment_2, 2)
    ref_m2_var = var(ref_moment_2, 2)

    rho = ref_autocor(moment_2, ref_m2_mean, ref_m2_var)
    (p, M) = size(rho)
    threshhold_acor = [ minimum([M; find(rho[i,:] .< 0.05).-1]) for i = 1:p]

    ESS2 = zeros(p)

    one_sT = 1 .- collect(1:M)/T
    for i = 1:p
        Mt = threshhold_acor[i]
        ESS2[i] = T/(1+2*dot(one_sT[1:Mt],rho[i,1:Mt][:]))
    end

    ESS, ESS2
end

function compute_bias2_var(chain, true_mean, true_var, g)
    (p,T) = size(chain)
    gs_arr = [g(chain[:,i]) for i in 1:T]
    gs = [gs_arr[j][i]::Float64 for i in 1:length(gs_arr[1]), j in 1:T]

    est_mean = sum(gs,2)/T
    bias2 = (est_mean .- true_mean).^2
    bias = sqrt(bias2)

#    println("bias: $bias")
#    println("gs[1]: $(gs[1])")

    #debiased_gs = gs .- bias

    #chain_var = [Lora.mcvar_imse(squeeze(debiased_gs[i,:],1))::Float64 for i = 1:p]
    chain_var = ref_autocov(gs, true_mean + bias, true_var)

    println("chain_var: $chain_var")

    return sum(bias2), sum(chain_var)
end


text:
include(joinpath("Occurrence", "logreg_chosenVariables.jl"))
include(joinpath("Occurrence", "logreg_research.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO