text:
# ---------------------------------------------------------------------------- #
#
#   hdgSolveElas.jl
#
#   Solve linear elasticity equations (n-dimensional)
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    hdgSolve( master::Master, mesh::Mesh, material::Material,
              problem::Problem)

Solves the linear elasticity equations for n-dimensional problems.
"""
function hdgSolve( master::Master, mesh::Mesh, problem::Elas; method=:LU)

dim     = mesh.dim
nelem   = size( mesh.nodes, 3 ) # Mumber of elements in mesh
nnodes  = size( mesh.nodes, 1 ) # Number of nodes in one element
nodfac  = master.nodfac # Number of nodes on a face

unkUhat = nodfac * dim          # Total of uhat unknowns per face
nfaces  = dim + 1               # Number of faces per element

# Compute the stiffness tensor
Cstiff = problem.mat.Cstiff

# Stability parameter
τ = Cstiff

### Initialize quantities

# Equations of motion (Newton's second law)
A = fill(0.0, nnodes*dim    , nnodes*dim^2      )    # (v_{i,j}, σ^h_{ij})_{T^h} - <v_i, σ^h_{ij} n_j>_{∂T^h}
B = fill(0.0, nnodes*dim    , nnodes*dim        )    # <v_i, func(u^h_k) >_{∂T^h}
C = fill(0.0, nnodes*dim    , nodfac*nfaces*dim )    # <v_i, func(\hat{u}^h_k) >_{∂T^h}

F = fill(0.0, nnodes*dim    , 1                 )    # (v_{i}, b_{i})_{T^h}

# strain-displacement
D = fill(0.0, nnodes*dim^2  , nnodes*dim^2      )    # (w_{ij}  , ϵ^h_{ij})_{T^h}
H = fill(0.0, nnodes*dim^2  , nnodes*dim        )    # <w_{ij}  , (\hat{u}^h_i*n_j+ \hat{u}^h_j*n_i)>_{∂T^h}
J = fill(0.0, nnodes*dim^2  , nodfac*nfaces*dim )    # (w_{ij,j}, u^h_i)_{T^h} + (w_{ij,i}, u^h_j*n_i)_{T^h}

# Hooke's law
K = fill(0.0, nnodes*dim^2  , nnodes*dim^2      )    # (z_{ij}  , σ^h_{ij})_{T^h}
M = fill(0.0, nnodes*dim^2  , nnodes*dim^2      )    # (z_{ij}  , C_{ijkl}ϵ^h_{kl})_{T^h}

# Flux jumps
N = fill(0.0, nfaces*nodfac*dim, nnodes*dim^2      ) # <μ_{i} , σ^h_{ij}*n_j>_{∂T^h\∂Ω_D}
P = fill(0.0, nfaces*nodfac*dim, nnodes*dim        ) # <μ_{i} , func(u^h_k)>_{∂T^h\∂Ω_D}
Q = fill(0.0, nfaces*nodfac*dim, nodfac*nfaces*dim ) # <μ_{i} , func(\hat{u}^h_k)>_{∂T^h\∂Ω_D} or <μ_{i} , \hat{u}^h_{i})_{∂Ω_D}

# Dirichlet BC
#   BC RHS
G = fill(0.0, nfaces*nodfac*dim, 1                ) # (μ_{i} , t_i)_{∂T^h\∂Ω_D}) or (μ_{i} , \bar{u}_i)_{∂Ω_D}

# Matrix with unknowns for uhath
A_UHATH = fill(0.0, nfaces*nodfac*dim, nodfac*nfaces*dim)
B_UHATH = fill(0.0, nfaces*nodfac*dim, 1                )

# Matrices that need to be saved in order to recover u, σ and ϵ at the end
sysUUH_11_12 = fill(0.0, nnodes*dim    , nodfac*nfaces*dim, nelem )
sysUUH_11_F  = fill(0.0, nnodes*dim    , 1,                 nelem )
KMDH         = fill(0.0, nnodes*dim^2  , nnodes*dim,        nelem )
KMDJ         = fill(0.0, nnodes*dim^2  , nodfac*nfaces*dim, nelem )

# Set up column and row indices for sparse matrix
indRow = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )
indCol = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )
indUnk = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )

Rfull = fill( 0.0, size(mesh.f,1)*unkUhat ) # RHS vector

rr = 1 # iterator for unknowns in sparsity matrix

# preallocate
jcwd = fill( 0.0, size(master.∇ϕ,2), size(master.∇ϕ,2) )
∂ξ∂x = fill( 0.0, size(master.∇ϕ,2), dim^2 )
∂x∂ξ = fill( 0.0, size(master.∇ϕ,2), dim^2 )

tempf  = fill( 0.0, nodfac, nodfac )
ϕjϕn   = [Array{Float64}(nodfac,nodfac) for idx in 1:dim]
ϕjϕnl  = [Array{Float64}(nodfac,nodfac) for idx in 1:dim, idy in 1:dim]
τϕjϕnl = [Array{Float64}(nodfac,nodfac) for idx in 1:dim, idy in 1:dim]

@time for pp in 1:nelem # Loop over all elements

  # zero out all matrices
  A .*= 0.0
  B .*= 0.0
  C .*= 0.0
  F .*= 0.0
  D .*= 0.0
  H .*= 0.0
  J .*= 0.0
  K .*= 0.0
  M .*= 0.0
  N .*= 0.0
  P .*= 0.0
  Q .*= 0.0
  G .*= 0.0
  A_UHATH .*= 0.0
  B_UHATH .*= 0.0

  # Compute Jacobians
  compJacob!( master, mesh.nodes[:,:,pp], ∂ξ∂x, jcwd, ∂x∂ξ )

  pLoc = master.ϕ' * mesh.nodes[:,:,pp]

  ∇ϕc = getderbfel( master, ∂ξ∂x )

  # mass matrices
  ϕj  = master.ϕ * jcwd
  ϕjϕ = master.ϕ * jcwd * master.ϕ'

  # ------------------------- Volume integrals ------------------------------- #
  ## A (first part)
  # (v_{i,j}, σ^h_{ij})_{T^h}
  for ii in 1:dim, jj in 1:dim
    ij = (ii-1)*dim + jj
    A[1+(ii-1)*nnodes:ii*nnodes,1+(ij-1)*nnodes:ij*nnodes] = ∇ϕc[:,:,jj] * ϕj'
  end

  ## D
  # (w_{ij}  , ϵ^h_{ij})_{T^h}
  for ii in 1:dim, jj in 1:dim
    ij = (ii-1)*dim + jj
    D[1+(ij-1)*nnodes:ij*nnodes,1+(ij-1)*nnodes:ij*nnodes] = ϕjϕ
  end

  ## H
  # (w_{ij,j}, u^h_i)_{T^h} + (w_{ij,i}, u^h_j)_{T^h}
  for ii in 1:dim, jj in 1:dim
    ij = (ii-1)*dim + jj
    H[1+(ij-1)*nnodes:ij*nnodes,1+(ii-1)*nnodes:ii*nnodes] += 0.5 * ∇ϕc[:,:,jj] * ϕj'
    H[1+(ij-1)*nnodes:ij*nnodes,1+(jj-1)*nnodes:jj*nnodes] += 0.5 * ∇ϕc[:,:,ii] * ϕj'
  end

  ## K
  # (z_{ij}  , σ^h_{ij})_{T^h}
  for ii in 1:dim, jj in 1:dim
    ij = (ii-1)*dim + jj
    K[1+(ij-1)*nnodes:ij*nnodes,1+(ij-1)*nnodes:ij*nnodes] = ϕjϕ
  end

  ## M
  # (z_{ij}  , C_{ijkl}ϵ^h_{kl})_{T^h}
  for ii in 1:dim, jj in 1:dim
    ij = (ii-1)*dim + jj
    for kk in 1:dim, ll in 1:dim
      kl = (kk-1)*dim + ll
      M[1+(ij-1)*nnodes:ij*nnodes,1+(kl-1)*nnodes:kl*nnodes] -= Cstiff[ii,jj,kk,ll] * ϕjϕ
    end
  end

  ## F
  # (v_{i}, b_{i})_{T^h}
  src = problem.source(pLoc)
  for ii in 1:dim
      F[1+(ii-1)*nnodes:ii*nnodes,1] = ϕj * src[:,ii]
  end

  # -------------------------------------------------------------------------- #
  #
  # ------------------------- Boundary integrals ----------------------------- #
  for qq in 1:nfaces # loop over number of faces

    (ϕdm, pdm, nod, normal, jcwdm) = compJacobFace( mesh, master, pp, qq )

    jcwddm = diagm( jcwdm )

    faceInd = dim * (qq-1) * nodfac + ( 1:nodfac )
    # We need to premultiply with "dim" because those are the number of unknowns per face

    # mass matrices
    ϕjf  = ϕdm * jcwddm
    ϕjϕf = ϕjf * ϕdm'
    for ii in 1:dim
        ϕjϕn[ii] = ϕjf * ( ϕdm * diagm(normal[:,ii]) )'
    end
    τϕjϕnl .*= 0.0
    for ll in 1:dim, jj in 1:dim
        ϕjϕnl[ll,jj] = ϕjf * ( ϕdm * diagm(normal[:,ll].*normal[:,jj]) )'
        for kk in 1:dim, ii in 1:dim
            τϕjϕnl[ii,kk] += τ[ii,jj,kk,ll] * ϕjϕnl[ll,jj]
        end
    end

    ## A (second part)
    # -<v_i, σ^h_{ij} n_j>_{∂T^h}
    for ii in 1:dim, jj in 1:dim
      ij = (ii-1)*dim + jj
      A[(ii-1)*nnodes + nod,(ij-1)*nnodes + nod] -= ϕjϕn[jj]
    end

    ## B
    # <v_i, τ_{ijkl} u^h_k * n_l * n_j >_{∂T^h}
    for kk in 1:dim, ii in 1:dim # column-major ordering
      B[(ii-1)*nnodes + nod,(kk-1)*nnodes + nod] +=
        τϕjϕnl[ii,kk]
    end

    ## C
    # do not have to broadcast these variables,
    #   because \hat{u}^h is different for each face
    # <v_i, τ_{ijkl} \hat{u}^h_k * n_l * n_j >_{∂T^h}
    for kk in 1:dim, ii in 1:dim # column-major ordering
      C[(ii-1)*nnodes + nod, (kk-1)*nodfac + faceInd] =
          -τϕjϕnl[ii,kk]
    end

    ## J
    # <w_{ij}  , (\hat{u}^h_i*n_j + \hat{u}^h_j*n_i)>_{∂T^h}
    for ii in 1:dim, jj in 1:dim
      ij = (ii-1)*dim + jj
      J[(ij-1)*nnodes + nod, (ii-1)*nodfac + faceInd] -= 0.5 * ϕjϕn[jj]
      J[(ij-1)*nnodes + nod, (jj-1)*nodfac + faceInd] -= 0.5 * ϕjϕn[ii]
    end

    indF = abs( mesh.t2f[pp,qq] )

    if ( mesh.f[ indF, end ] < 0 )
      # At boundary

      bNo = - mesh.f[ indF, end ]

      if problem.bctype[bNo] == 1
        # Dirichlet boundary condition

        ## Q
        # (μ_{i} , \hat{u}^h_{i})_{∂Ω_D}
        for ii in 1:dim     #i
          Q[(ii-1)*nodfac + faceInd,(ii-1)*nodfac + faceInd] = ϕjϕf
        end

        # BC RHS
        ## G
        # (μ_{i} , \bar{u}_i)_{∂Ω_D}
        bcout = problem.bcfunc[bNo](pdm)
        for ii in 1:dim     #i
          G[(ii-1)*nodfac + faceInd, 1] = ϕjf * bcout[:,ii]
        end

      elseif problem.bctype[bNo] == 2
        # Neumann boundary condition

        # do not have to broadcast these variables,
        #   because μ is different for each face

        ## N
        # <μ_{i} , σ^h_{ij}*n_j>_{∂T^h\∂Ω_D}
        for ii in 1:dim, jj in 1:dim   #j
          ij = (ii-1)*dim + jj
          N[(ii-1)*nodfac + faceInd, (ij-1)*nnodes + nod] =
              ϕjϕn[jj]
        end

        ## P
        # <μ_{i} , -τ_{ijkl} u^h_k n_l n_j ) >_{∂T^h\∂Ω_D}
        for kk in 1:dim, ii in 1:dim # column-major ordering
          P[(ii-1)*nodfac + faceInd,(kk-1)*nnodes + nod] =
            -τϕjϕnl[ii,kk]
        end

        ## Q
        # <μ_{i} , τ_{ijkl} \hat{u}^h_k n_l n_j ) >_{∂T^h\∂Ω_D}
        for kk in 1:dim, ii in 1:dim # column-major ordering
          Q[(ii-1)*nodfac + faceInd,(kk-1)*nodfac + faceInd] =
            τϕjϕnl[ii,kk]
        end

        # BC RHS
        ## G
        # (μ_{i} , t_i)_{∂T^h\∂Ω_D})
        if problem.bcnorm
          # Boundary conditions are defined normal to the surface
          bcout = problem.bcfunc[bNo](pdm)
          for ii in 1:dim     #i
            temp = normal[:,ii] .* bcout[:,1]# + tL[:,ii] .* bcout[:,2] # NOTE: No tangential component, because hard to define in 3D
            G[(ii-1)*nodfac + faceInd, 1] = ϕjf * temp
          end
        else
          # Boundary conditions are defined in the general coordinate system
          bcout = problem.bcfunc[bNo](pdm)
          for ii in 1:dim     #i
            G[(ii-1)*nodfac + faceInd, 1] = ϕjf * bcout[:,ii]
          end
        end
      else
          error("hdgSolveElas:: BC type not recognized. 1 = Dirichlet, 2 = Neumann")
      end

    else
      # In interior

      # do not have to broadcast these variables,
      #  because μ is different for each face

      ## N
      # <μ_{i} , σ^h_{ij}*n_j>_{∂T^h\∂Ω_D}
      for ii in 1:dim, jj in 1:dim
        ij = (ii-1)*dim + jj
        N[(ii-1)*nodfac + faceInd, (ij-1)*nnodes + nod] =
          ϕjϕn[jj]
      end

      ## P
      # <μ_{i} , -τ_{ijkl} u^h_k n_l n_j ) >_{∂T^h\∂Ω_D}
      for kk in 1:dim, ii in 1:dim # column-major ordering
        P[(ii-1)*nodfac + faceInd,(kk-1)*nnodes + nod] =
          -τϕjϕnl[ii,kk]
      end

      ## Q
      # <μ_{i} , τ_{ijkl} \hat{u}^h_k n_l n_j ) >_{∂T^h\∂Ω_D}
      for kk in 1:dim, ii in 1:dim # column-major ordering
        Q[(ii-1)*nodfac + faceInd,(kk-1)*nodfac + faceInd] =
          τϕjϕnl[ii,kk]
      end

    end # boundary if-statement

  end # end loop over element faces
  # -------------------------------------------------------------------------- #

  # ------------------------- Form elemental quantities ---------------------- #
  invKM = K \ M
  invDH = D \ H
  invDJ = D \ J

  SYS_U_UHATH_11 = B + A * invKM * invDH
  SYS_U_UHATH_12 = C + A * invKM * invDJ
  SYS_U_UHATH_21 = P + N * invKM * invDH
  SYS_U_UHATH_22 = Q + N * invKM * invDJ

  tempMat =  SYS_U_UHATH_11 \ ( [ -SYS_U_UHATH_12 F ] )
  A_UHATH =  SYS_U_UHATH_21 * tempMat[:,1:end-1] + SYS_U_UHATH_22
  B_UHATH = -SYS_U_UHATH_21 * tempMat[:,end]     + G

  sysUUH_11_12[:,:,pp] = tempMat[:,1:end-1]
  sysUUH_11_F[:,:,pp]  = tempMat[:,end]
  KMDH[:,:,pp]         = invKM * invDH
  KMDJ[:,:,pp]         = invKM * invDJ
  # -------------------------------------------------------------------------- #

  # ------------------------- Fill up complete H and R matrices -------------- #
  for qq = 1:nfaces
    indJ = [i for i=1:nfaces]
    deleteat!(indJ,qq)

    indFall = abs.(mesh.t2f[pp,1:nfaces])
    deleteat!(indFall,qq)

    indF = abs(mesh.t2f[pp,qq])
    Rfull[1+(indF-1)*unkUhat:indF*unkUhat] += B_UHATH[1+(qq-1)*unkUhat:qq*unkUhat,1]

    indRow[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, unkUhat )
    indCol[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, 1, unkUhat )'[:]
    indUnk[rr:rr-1+unkUhat^2] = A_UHATH[ 1+(qq-1)*unkUhat:qq*unkUhat, 1+(qq-1)*unkUhat:qq*unkUhat ][:]

    rr = rr + unkUhat^2

    for tt = 1:(nfaces-1)
      indRow[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, unkUhat )
      indCol[rr:rr-1+unkUhat^2] = repmat(  1+(indFall[tt]-1)*unkUhat:indFall[tt]*unkUhat, 1, unkUhat )'[:]
      indUnk[rr:rr-1+unkUhat^2] = A_UHATH[ 1+(qq-1)*unkUhat:qq*unkUhat, 1+(indJ[tt]-1)*unkUhat:indJ[tt]*unkUhat ][:]

      rr = rr + unkUhat^2
    end
  end
  # -------------------------------------------------------------------------- #
end # end element loop

### Compute approximate trace
Hfull = sparse( indRow, indCol, indUnk, size(mesh.f,1) * unkUhat, size(mesh.f,1) * unkUhat )

uhath = [0.]
if method == :LU
    @time uhath = Hfull \ Rfull
elseif method == :GMRES
    @time uhath = IterativeSolvers.gmres( Hfull, Rfull )
elseif method == :GMRESILU
    @time fact = ILU.crout_ilu( Hfull, τ = 0.1 )
    @time uhath = IterativeSolvers.gmres( Hfull, Rfull, Pl=fact )
end
# ---------------------------------------------------------------------------- #

## Compute approximate scalar value and flux
uhathTri = fill( 0.0, unkUhat*nfaces, 1 )
uh       = fill( 0.0, nnodes,         dim,   nelem )
σh       = fill( 0.0, nnodes,         dim^2, nelem )

@time for pp in 1:nelem
    # ----------- Find uhath corresponding to this element ------------------- #
    indF = abs.( mesh.t2f[pp,1:nfaces] )
    for qq in 1:nfaces
      uhathTri[ 1+(qq-1)*unkUhat:qq*unkUhat, 1 ] = uhath[ 1+(indF[qq]-1)*unkUhat:indF[qq]*unkUhat ]
    end
    # ------------------------------------------------------------------------ #
    # ----------- Compute approximate displacement value --------------------- #
    uTemp = sysUUH_11_F[:,:,pp] + sysUUH_11_12[:,:,pp] * uhathTri
    for ii in 1:dim
        uh[:,ii,pp] = uTemp[ 1 + (ii-1)*nnodes:nnodes*ii ]
    end
    # ------------------------------------------------------------------------ #
    # ----------- Compute approximate stress field --------------------------- #
    sTemp = KMDH[:,:,pp] * uTemp + KMDJ[:,:,pp] * uhathTri
    for ii in 1:dim^2
        σh[ :,ii,pp ] = sTemp[ 1 + (ii-1)*nnodes:nnodes*ii ]
    end
    # ------------------------------------------------------------------------ #
end

return (uhath, uh, σh)

end # end hdgSolveElas


text:
module State

import ..Config

mutable struct SharedState
    config::Config.AppConfig
    healthy::Bool
end

function SharedState(config::Config.AppConfig)
    SharedState(config, true)
end

end

ext: jl

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
struct UnsignedIterator{U}
    u::U
end

Base.IteratorSize(::Type{<:UnsignedIterator}) = Base.HasShape{1}()
Base.IteratorEltype(::Type{<:UnsignedIterator}) = Base.HasEltype()

Base.eltype(::UnsignedIterator) = UInt32
Base.length(u::UnsignedIterator) = count_ones(u.u)
Base.size(u::UnsignedIterator) = (count_ones(u.u),)

# @inline function Base.iterate(u::UnsignedIterator, uu = u.u)
#     tz = trailing_zeros(uu) % UInt32
#     # tz ≥ 0x00000020 && return nothing
#     tz > 0x0000001f && return nothing
#     uu ⊻= (0x00000001 << tz)
#     tz, uu
# end
@inline function Base.iterate(u::UnsignedIterator, (i,uu) = (0x00000000,u.u))
    tz = trailing_zeros(uu) % UInt32
    tz == 0x00000020 && return nothing
    i += tz
    tz += 0x00000001
    uu >>>= tz
    (i, (i+0x00000001,uu))
end

"""
    UnsignedIteratorEarlyStop(thread_mask[, num_threads = count_ones(thread_mask)])

Iterator, returning `(i,t) = Tuple{UInt32,UInt32}`, where `i` iterates from `1,2,...,num_threads`, and `t` gives the threadids to call `ThreadingUtilities.taskpointer` with.

Unfortunately, codegen is suboptimal when used in the ergonomic `for (i,tid) ∈ thread_iterator` fashion. If you want to microoptimize,
You'd get better performance from a pattern like:
```julia
function sumk(u,l = count_ones(u) % UInt32)
    uu = ServiceSolicitation.UnsignedIteratorEarlyStop(u,l)
    s = zero(UInt32); state = ServiceSolicitation.initial_state(uu)
    while true
        iter = iterate(uu, state)
        iter === nothing && break
        (i,t),state = iter
        s += t
    end
    s
end
```

This iterator will iterate at least once; it's important to check and exit early with a single threaded version.
"""
struct UnsignedIteratorEarlyStop{U}
    u::U
    i::UInt32
end
UnsignedIteratorEarlyStop(u) = UnsignedIteratorEarlyStop(u, count_ones(u) % UInt32)
UnsignedIteratorEarlyStop(u, i) = UnsignedIteratorEarlyStop(u, i % UInt32)

mask(u::UnsignedIteratorEarlyStop) = getfield(u, :u)
Base.IteratorSize(::Type{<:UnsignedIteratorEarlyStop}) = Base.HasShape{1}()
Base.IteratorEltype(::Type{<:UnsignedIteratorEarlyStop}) = Base.HasEltype()

Base.eltype(::UnsignedIteratorEarlyStop) = Tuple{UInt32,UInt32}
Base.length(u::UnsignedIteratorEarlyStop) = getfield(u, :i)
Base.size(u::UnsignedIteratorEarlyStop) = (getfield(u, :i),)

@inline function initial_state(u::UnsignedIteratorEarlyStop)
  # LLVM should figure this out if you check?
  assume(0x00000000 ≠ u.i)
  (0x00000000,u.u)
end
@inline function iter(i, uu)
  assume(uu ≠ zero(uu))
  tz = trailing_zeros(uu) % UInt32
  tz += 0x00000001
  i += tz
  uu >>>= tz
  i, uu
end
@inline function Base.iterate(u::UnsignedIteratorEarlyStop, ((i,uu),j) = (initial_state(u),0x00000000))
  # assume(u.i ≤ 0x00000020)
  # assume(j ≤ count_ones(uu))
  # iszero(j) && return nothing
  j == u.i && return nothing
  j += 0x00000001
  i, uu = iter(i, uu)
  ((j, i), ((i, uu), j))
end
function Base.show(io::IO, u::UnsignedIteratorEarlyStop)
    l = length(u)
    s = Vector{Int}(undef, l)
    if l > 0
        s .= last.(u)
    end
    print("Thread ($l) Iterator: U", s)
end

# @inline function Base.iterate(u::UnsignedIteratorEarlyStop, (i,uu) = (0xffffffff,u.u))
#     tz = trailing_zeros(uu) % UInt32
#     tz == 0x00000020 && return nothing
#     tz += 0x00000001
#     i += tz
#     uu >>>= tz
#     (i, (i,uu))
# end


text:
function sqlite3_errmsg()
    return ccall( (:sqlite3_errmsg, libsqlite),
        Ptr{UInt8}, ()
        )
end
function sqlite3_errmsg(db::Ptr{Cvoid})
    @NULLCHECK db
    return ccall( (:sqlite3_errmsg, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid},),
        db)
end
function sqlite3_open(file::AbstractString, handle)
    return ccall( (:sqlite3_open, libsqlite),
        Cint, (Ptr{UInt8}, Ptr{Cvoid}),
        file, handle)
end
function sqlite3_close(handle::Ptr{Cvoid})
    @NULLCHECK handle
    return ccall( (:sqlite3_close, libsqlite),
        Cint, (Ptr{Cvoid},),
        handle)
end
function sqlite3_next_stmt(db::Ptr{Cvoid}, stmt::Ptr{Cvoid})
    @NULLCHECK db
    return ccall( (:sqlite3_next_stmt, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}),
        db, stmt)
end
function sqlite3_prepare_v2(handle::Ptr{Cvoid}, query::AbstractString, stmt, unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare_v2, libsqlite),
        Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
            handle, query, sizeof(query), stmt, unused)
end
function sqlite3_prepare16_v2(handle::Ptr{Cvoid}, query::AbstractString, stmt, unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare16_v2, libsqlite),
        Cint, (Ptr{Cvoid}, Ptr{UInt16}, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
        handle, query, sizeof(query), stmt, unused)
end
function sqlite3_finalize(stmt::Ptr{Cvoid})
    @NULLCHECK stmt
    return ccall( (:sqlite3_finalize, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end

# SQLITE_API int sqlite3_bind_paramter_count(sqlite3_stmt*)
function sqlite3_bind_parameter_count(stmt::Ptr{Cvoid})
    return ccall( (:sqlite3_bind_parameter_count, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end
#SQLITE_API const char* sqlite3_bind_parameter_name(sqlite3_stmt*, int)
function sqlite3_bind_parameter_name(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_bind_parameter_name, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid}, Cint),
        stmt, col)
end
# SQLITE_API int sqlite3_bind_parameter_index(sqlite3_stmt*, const char *zName);
function sqlite3_bind_parameter_index(stmt::Ptr{Cvoid}, value::AbstractString)
    return ccall( (:sqlite3_bind_parameter_index, libsqlite),
        Cint, (Ptr{Cvoid}, Ptr{UInt8}),
        stmt, value)
end
# SQLITE_API int sqlite3_bind_double(sqlite3_stmt*, int, double);
function sqlite3_bind_double(stmt::Ptr{Cvoid}, col::Int, value::Float64)
    return ccall( (:sqlite3_bind_double, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Float64),
        stmt, col, value)
end
# SQLITE_API int sqlite3_bind_int(sqlite3_stmt*, int, int);
function sqlite3_bind_int(stmt::Ptr{Cvoid}, col::Int, value::Int32)
    return ccall( (:sqlite3_bind_int, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Int32),
        stmt, col, value)
end
# SQLITE_API int sqlite3_bind_int64(sqlite3_stmt*, int, sqlite3_int64);
function sqlite3_bind_int64(stmt::Ptr{Cvoid}, col::Int, value::Int64)
    return ccall( (:sqlite3_bind_int64, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Int64),
        stmt, col, value)
end
# SQLITE_API int sqlite3_bind_null(sqlite3_stmt*, int);
function sqlite3_bind_null(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_bind_null, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        stmt, col)
end
# SQLITE_API int sqlite3_bind_text(sqlite3_stmt*, int, const char*, int n, void(*)(void*));
function sqlite3_bind_text(stmt::Ptr{Cvoid}, col::Int, value::AbstractString)
    return ccall( (:sqlite3_bind_text, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Ptr{UInt8}, Cint, Ptr{Cvoid}),
        stmt, col, value, sizeof(value), C_NULL)
end
function sqlite3_bind_text(stmt::Ptr{Cvoid}, col::Int, ptr::Ptr{UInt8}, len::Int)
    return ccall( (:sqlite3_bind_text, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Ptr{UInt8}, Cint, Ptr{Cvoid}),
        stmt, col, ptr, len, C_NULL)
end
# SQLITE_API int sqlite3_bind_text16(sqlite3_stmt*, int, const void*, int, void(*)(void*));
function sqlite3_bind_text16(stmt::Ptr{Cvoid}, col::Int, ptr::Ptr{UInt16}, len::Int)
    return ccall( (:sqlite3_bind_text16, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Ptr{UInt16}, Cint, Ptr{Cvoid}),
        stmt, col, ptr, len, C_NULL)
end

# SQLITE_API int sqlite3_bind_blob(sqlite3_stmt*, int, const void*, int n, void(*)(void*));
function sqlite3_bind_blob(stmt::Ptr{Cvoid}, col::Int, value)
    return ccall( (:sqlite3_bind_blob, libsqlite),
        Cint, (Ptr{Cvoid}, Cint, Ptr{UInt8}, Cint, Ptr{Cvoid}),
        stmt, col, value, sizeof(value), SQLITE_STATIC)
end
# SQLITE_API int sqlite3_bind_zeroblob(sqlite3_stmt*, int, int n);
# SQLITE_API int sqlite3_bind_value(sqlite3_stmt*, int, const sqlite3_value*);

# SQLITE_API int sqlite3_busy_timeout(sqlite3*, int ms);
function sqlite3_busy_timeout(db::Ptr{Cvoid}, ms)
    @NULLCHECK db
    return ccall( (:sqlite3_busy_timeout, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        db, ms)
end

function sqlite3_clear_bindings(stmt::Ptr{Cvoid})
    return ccall( (:sqlite3_clear_bindings, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end

function sqlite3_step(stmt::Ptr{Cvoid})
    return ccall( (:sqlite3_step, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end
function sqlite3_column_count(stmt::Ptr{Cvoid})
    return ccall( (:sqlite3_column_count, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end
function sqlite3_column_type(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_type, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_blob(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_blob, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_bytes(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_bytes, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_bytes16(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_bytes16, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_double(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_double, libsqlite),
        Cdouble, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_int(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_int, libsqlite),
        Cint, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_int64(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_int64, libsqlite),
        Clonglong, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_text(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_text, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_text16(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_text16, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
# function sqlite3_column_value(stmt::Ptr{Cvoid}, col::Cint)
#     return ccall( (:sqlite3_column_value, libsqlite),
#             Ptr{Cvoid}, (Ptr{Cvoid}, Cint),
#             stmt, col-1)
# end
# SQLITE_API sqlite3_value *sqlite3_column_value(sqlite3_stmt*, int iCol);
function sqlite3_reset(stmt::Ptr{Cvoid})
    return ccall( (:sqlite3_reset, libsqlite),
        Cint, (Ptr{Cvoid},),
        stmt)
end

# SQLITE_API const char *sqlite3_column_name(sqlite3_stmt*, int N);
function sqlite3_column_name(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_name, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
function sqlite3_column_name16(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_name16, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end

function sqlite3_changes(db::Ptr{Cvoid})
    @NULLCHECK db
   return ccall( (:sqlite3_changes, libsqlite),
        Cint, (Ptr{Cvoid},),
        db)
end
function sqlite3_total_changes(db::Ptr{Cvoid})
    @NULLCHECK db
   return ccall( (:sqlite3_changes, libsqlite),
        Cint, (Ptr{Cvoid},),
        db)
end
# SQLITE_API const void *sqlite3_column_name16(sqlite3_stmt*, int N);

# SQLITE_API const char *sqlite3_column_database_name(sqlite3_stmt*, int);
# SQLITE_API const void *sqlite3_column_database_name16(sqlite3_stmt*, int);
# SQLITE_API const char *sqlite3_column_table_name(sqlite3_stmt*, int);
# SQLITE_API const void *sqlite3_column_table_name16(sqlite3_stmt*, int);
# SQLITE_API const char *sqlite3_column_origin_name(sqlite3_stmt*, int);
# SQLITE_API const void *sqlite3_column_origin_name16(sqlite3_stmt*, int);

function sqlite3_column_decltype(stmt::Ptr{Cvoid}, col::Int)
    return ccall( (:sqlite3_column_decltype, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid}, Cint),
        stmt, col-1)
end
# SQLITE_API const char *sqlite3_column_decltype(sqlite3_stmt*, int);
# SQLITE_API const void *sqlite3_column_decltype16(sqlite3_stmt*, int);

# SQLITE_API int sqlite3_data_count(sqlite3_stmt *pStmt);

# SQLITE_API void sqlite3_result_double(sqlite3_context*, double);
function sqlite3_result_double(context::Ptr{Cvoid}, value::Float64)
    return ccall( (:sqlite3_result_double, libsqlite),
        Cvoid, (Ptr{Cvoid}, Float64),
        context, value)
end
# SQLITE_API void sqlite3_result_error(sqlite3_context*, const char*, int)
function sqlite3_result_error(context::Ptr{Cvoid}, msg::AbstractString)
    return ccall( (:sqlite3_result_error, libsqlite),
        Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, Cint),
        context, value, sizeof(msg)+1)
end
# SQLITE_API void sqlite3_result_error16(sqlite3_context*, const void*, int)
# SQLITE_API void sqlite3_result_int(sqlite3_context*, int);
function sqlite3_result_int(context::Ptr{Cvoid}, value::Int32)
    return ccall( (:sqlite3_result_int, libsqlite),
        Cvoid, (Ptr{Cvoid}, Int32),
        context, value)
end
# SQLITE_API void sqlite3_result_int64(sqlite3_context*, sqlite3_int64);
function sqlite3_result_int64(context::Ptr{Cvoid}, value::Int64)
    return ccall( (:sqlite3_result_int64, libsqlite),
        Cvoid, (Ptr{Cvoid}, Int64),
        context, value)
end
# SQLITE_API void sqlite3_result_null(sqlite3_context*);
function sqlite3_result_null(context::Ptr{Cvoid})
    return ccall( (:sqlite3_result_null, libsqlite),
        Cvoid, (Ptr{Cvoid},),
        context)
end
# SQLITE_API void sqlite3_result_text(sqlite3_context*, const char*, int n, void(*)(void*));
function sqlite3_result_text(context::Ptr{Cvoid}, value::AbstractString)
    return ccall( (:sqlite3_result_text, libsqlite),
        Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cvoid}),
        context, value, sizeof(value)+1, SQLITE_TRANSIENT)
end
# SQLITE_API void sqlite3_result_text16(sqlite3_context*, const void*, int, void(*)(void*));
# SQLITE_API void sqlite3_result_blob(sqlite3_context*, const void*, int n, void(*)(void*));
function sqlite3_result_blob(context::Ptr{Cvoid}, value)
    return ccall( (:sqlite3_result_blob, libsqlite),
        Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cvoid}),
        context, value, sizeof(value), SQLITE_TRANSIENT)
end
# SQLITE_API void sqlite3_result_zeroblob(sqlite3_context*, int n);
# SQLITE_API void sqlite3_result_value(sqlite3_context*, const sqlite3_value*);
# SQLITE_API void sqlite3_result_error_toobig(sqlite3_context*)
# SQLITE_API void sqlite3_result_error_nomem(sqlite3_context*)
# SQLITE_API void sqlite3_result_error_code(sqlite3_context*, int)

function sqlite3_create_function_v2(db::Ptr{Cvoid}, name::AbstractString, nargs::Integer,
                                    enc::Integer, data::Ptr{Cvoid}, func,
                                    step, final,
                                    destructor::Ptr{Cvoid})
    @NULLCHECK db
    return ccall(
        (:sqlite3_create_function_v2, libsqlite),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Cvoid},
         Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        db, name, nargs, enc, data, func, step, final, destructor)
end

# SQLITE_API void* sqlite3_aggregate_context(sqlite3_context*, int nBytes)
function sqlite3_aggregate_context(context::Ptr{Cvoid}, nbytes::Integer)
    return ccall( (:sqlite3_aggregate_context, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid}, Cint),
        context, nbytes)
end

# SQLITE_API int sqlite3_value_type(sqlite3_value*)
function sqlite3_value_type(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_type, libsqlite),
        Cint, (Ptr{Cvoid},),
        value)
end

# SQLITE_API const void* sqlite3_value_blob(sqlite3_value*)
function sqlite3_value_blob(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_blob, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid},),
        value)
end
# SQLITE_API int sqlite3_value_bytes(sqlite3_value*)
function sqlite3_value_bytes(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_bytes, libsqlite),
        Cint, (Ptr{Cvoid},),
        value)
end
# SQLITE_API int sqlite3_value_bytes16(sqlite3_value*)
function sqlite3_value_bytes16(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_bytes16, libsqlite),
        Cint, (Ptr{Cvoid},),
        value)
end
# SQLITE_API double sqlite3_value_double(sqlite3_value*)
function sqlite3_value_double(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_double, libsqlite),
        Cdouble, (Ptr{Cvoid},),
        value)
end
# SQLITE_API int sqlite3_value_int(sqlite3_value*)
function sqlite3_value_int(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_int, libsqlite),
        Cint, (Ptr{Cvoid},),
        value)
end
# SQLITE_API sqlite_int64 sqlite3_value_int64(sqlite3_value*)
function sqlite3_value_int64(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_int64, libsqlite),
        Clonglong, (Ptr{Cvoid},),
        value)
end
# SQLITE_API const unsigned char* sqlite3_value_text(sqlite3_value*)
function sqlite3_value_text(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_text, libsqlite),
        Ptr{UInt8}, (Ptr{Cvoid},),
        value)
end
# SQLITE_API const void* sqlite3_value_text16(sqlite3_value*)
function sqlite3_value_text16(value::Ptr{Cvoid})
    return ccall( (:sqlite3_value_text16, libsqlite),
        Ptr{Cvoid}, (Ptr{Cvoid},),
        value)
end
# SQLITE_API int sqlite3_value_numeric_type(sqlite3_value*)

function sqlite3_initialize()
    return ccall( (:sqlite3_initialize, libsqlite),
        Cint, (),
        )
end
function sqlite3_shutdown()
    return ccall( (:sqlite3_shutdown, libsqlite),
        Cint, (),
        )
end
function sqlite3_os_init()
    return ccall( (:sqlite3_os_init, libsqlite),
        Cint, (),
        )
end
function sqlite3_os_end()
    return ccall( (:sqlite3_os_end, libsqlite),
        Cint, (),
        )
end
function sqlite3_free_table(result::Array{AbstractString, 1})
    return ccall( (:sqlite3_free_table, libsqlite),
        Cvoid, (Ptr{Ptr{Cvoid}},),
        result)
end

# SQLITE_API const char *sqlite3_uri_parameter(const char *zFilename, const char *zParam);
# SQLITE_API int sqlite3_uri_boolean(const char *zFile, const char *zParam, int bDefault);
# SQLITE_API sqlite3_int64 sqlite3_uri_int64(const char*, const char*, sqlite3_int64);
function sqlite3_errcode(db::Ptr{Cvoid})
    @NULLCHECK db
    return ccall( (:sqlite3_errcode, libsqlite),
        Cint, (Ptr{Cvoid},),
        db)
end
function sqlite3_extended_errcode(db::Ptr{Cvoid})
    @NULLCHECK db
    return ccall( (:sqlite3_extended_errcode, libsqlite),
        Cint, (Ptr{Cvoid},),
        db)
end
# SQLITE_API int sqlite3_errcode(sqlite3 *db);
# SQLITE_API int sqlite3_extended_errcode(sqlite3 *db);
function sqlite3_errstr(ret::Cint)
    return ccall( (:sqlite3_errstr, libsqlite),
        Ptr{UInt8}, (Cint,),
        ret)
end
# SQLITE_API const char *sqlite3_errstr(int);

# SQLITE_API int sqlite3_limit(sqlite3*, int id, int newVal);

# SQLITE_API int sqlite3_stmt_readonly(sqlite3_stmt *pStmt);

# SQLITE_API int sqlite3_stmt_busy(sqlite3_stmt*);

# SQLITE_API int sqlite3_table_column_metadata(
#   sqlite3 *db,                /* Connection handle */
#   const char *zDbName,        /* Database name or NULL */
#   const char *zTableName,     /* Table name */
#   const char *zColumnName,    /* Column name */
#   char const **pzDataType,    /* OUTPUT: Declared data type */
#   char const **pzCollSeq,     /* OUTPUT: Collation sequence name */
#   int *pNotNull,              /* OUTPUT: True if NOT NULL constraint exists */
#   int *pPrimaryKey,           /* OUTPUT: True if column part of PK */
#   int *pAutoinc               /* OUTPUT: True if column is auto-increment */
# );

# SQLITE_API int sqlite3_db_status(sqlite3*, int op, int *pCur, int *pHiwtr, int resetFlg);

# Not directly used
function sqlite3_open_v2(file::AbstractString, handle, flags::Cint, vfs::AbstractString)
    return ccall( (:sqlite3_open_v2, libsqlite),
            Cint, (Ptr{UInt8}, Ptr{Cvoid}, Cint, Ptr{UInt8}),
            file, handle, flags, vfs)
end
function sqlite3_prepare(handle::Ptr{Cvoid}, query::AbstractString, stmt, unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare, libsqlite),
        Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
            handle, query, sizeof(query), stmt, unused)
end
function sqlite3_prepare16(handle::Ptr{Cvoid}, query::AbstractString, stmt, unused)
    @NULLCHECK handle
    return ccall( (:sqlite3_prepare16, libsqlite),
        Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
            handle, query, sizeof(query), stmt, unused)
end
function sqlite3_close_v2(handle::Ptr{Cvoid})
    @NULLCHECK handle
    try
        return ccall( (:sqlite3_close_v2, libsqlite),
            Cint, (Ptr{Cvoid},),
            handle)
    catch
        # Older versions of the library don't have this, abort to other close
        warn("sqlite3_close_v2 not available.")
        sqlite3_close(handle)
    end
end
function sqlite3_last_insert_rowid(handle::Ptr{Cvoid})
    @NULLCHECK handle
    return ccall( (:sqlite3_last_insert_rowid, libsqlite),
        Clong, (Ptr{Cvoid},),
        handle)
end


text:
export sample, to_sampled_dataframes

import DataFrames: DataFrame, mapcols!
import Random: default_rng

#####################
# Circuit sampling
#####################

"""
    sample(pc::ProbCircuit, num_samples)
    sample(pc::ProbCircuit, num_samples, evidences)

Sample states from the probabilistic circuit distribution. Also can do conditional sampling if evidence is given (any subset of features).
"""
function sample(pc::ProbCircuit; rng = default_rng())
    states, prs = sample(pc, 1, [missing for i=1:num_variables(pc)]...; rng)
    return states[1,:], prs[1]
end

sample(pc::ProbCircuit, num_samples; rng = default_rng(), gpu=false) =
    sample(pc, num_samples, [missing for i=1:num_variables(pc)]...; rng, gpu)

sample(pc::ProbCircuit, num_samples, inputs::Union{Bool,Missing}...;
        rng = default_rng(), gpu=false) =
    sample(pc, num_samples, collect(Union{Bool,Missing}, inputs); rng, gpu)

function sample(pc::ProbCircuit, num_samples, inputs::AbstractVector{Union{Bool,Missing}};
                rng = default_rng(), gpu=false)
    data = DataFrame(reshape(inputs, 1, :), :auto)
    data = gpu ? to_gpu(data) : data
    states, prs = sample(pc, num_samples, data; rng)
    return states[:,1,:], prs[:,1]
end

sample(circuit::ProbCircuit, num_samples, data::DataFrame; rng = default_rng()) =
    sample(same_device(ParamBitCircuit(circuit, data), data), num_samples, data; rng)

function sample(pbc::ParamBitCircuit, num_samples, data; Float = Float32, rng = default_rng())
    @assert isgpu(data) == isgpu(pbc) "ParamBitCircuit and data need to be on the same device"
    values = marginal_all(pbc, data)
    return sample_down(pbc, num_samples, data, values, rng, Float)
end

"Convert an array of samples into a vector of dataframes"
function to_sampled_dataframes(states)
    dfs = mapslices(DataFrame, states, dims = [2,3])
    map(dfs) do df
        mapcols!(c -> BitVector(c), df)
    end
    return dfs
end

# CPU code

function sample_down(pbc, num_samples, data, values::Array, rng, ::Type{Float}) where Float
    state = zeros(Bool, num_samples, num_examples(data), num_features(data))
    logprob = zeros(Float, num_samples, num_examples(data))
    Threads.@threads for (s_id, ex_id) = collect(Iterators.product(1:size(state,1), 1:size(state,2)))
        sample_rec(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), ex_id, s_id, num_nodes(pbc), values, state, logprob, rng)
    end
    return state, logprob
end

function sample_rec(nl, params, nodes, elements, ex_id, s_id, dec_id, values, state, logprob, rng)
    if isleafgate(nl, dec_id)
        if isliteralgate(nl, dec_id)
            l = literal(nl, dec_id)
            @inbounds state[s_id, ex_id, lit2var(l)] = (l > 0)
        end
    else
        edge_log_pr, prime, sub = sample_child_cpu(params, nodes, elements, ex_id, dec_id, values, rng)
        @inbounds logprob[s_id, ex_id] += edge_log_pr
        sample_rec(nl, params, nodes, elements, ex_id, s_id, prime, values, state, logprob, rng)
        sample_rec(nl, params, nodes, elements, ex_id, s_id, sub,   values, state, logprob, rng)
    end
end

function sample_child_cpu(params, nodes, elements, ex_id, dec_id, values, rng)
    @inbounds els_start = nodes[1,dec_id]
    @inbounds els_end = nodes[2,dec_id]
    threshold = log(rand(rng)) + values[ex_id, dec_id]
    cumul_prob = -Inf
    j_sampled = els_end - els_start + 1 # give all numerical error probability to the last node
    for j = els_start:els_end
        @inbounds prime = elements[2,j]
        @inbounds sub = elements[3,j]
        @inbounds pr = values[ex_id, prime] + values[ex_id, sub] + params[j]
        Δ = ifelse(cumul_prob == pr, zero(cumul_prob), abs(cumul_prob - pr))
        cumul_prob = max(cumul_prob, pr) + log1p(exp(-Δ))
        if cumul_prob > threshold
            j_sampled = j
            break
        end
    end
    @inbounds return params[j_sampled], elements[2,j_sampled], elements[3,j_sampled]
end

# GPU code

function sample_down(pbc, num_samples, data, values::CuArray, rng, ::Type{Float}) where Float
    CUDA.seed!(rand(rng, UInt))
    state = CUDA.zeros(Bool, num_samples, num_examples(data), num_features(data))
    logprob = CUDA.zeros(Float, num_samples, num_examples(data))
    stack = CUDA.zeros(Int32, num_samples, num_examples(data), num_features(data)+3)
    @inbounds stack[:,:,1] .= 1 # start with 1 dec_id in the stack
    @inbounds stack[:,:,2] .= num_nodes(pbc) # start with the pc in the stack
    num_threads =  balance_threads(num_samples, num_examples(data), 8)
    num_blocks = (ceil(Int, num_samples/num_threads[1]),
                  ceil(Int, num_examples(data)/num_threads[2]))
    CUDA.@sync  while true
        r = CUDA.rand(num_samples, num_examples(data))
        kernel = @cuda name="sample_cuda_kernel" launch=false sample_cuda_kernel(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), values, state, logprob, stack, r, Float)
        config = launch_configuration(kernel.fun)
        threads, blocks = balance_threads_2d(num_samples, num_examples(data), config.threads)
        kernel(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), values, state, logprob, stack, r, Float
            ; threads, blocks)

        all_empty(stack) && break
    end
    CUDA.unsafe_free!(values) # save the GC some effort
    return state, logprob
end

function sample_cuda_kernel(nl, params, nodes, elements, values, state, logprob, stack, r, ::Type{Float}) where Float
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for s_id = index_x:stride_x:size(state,1)
        for ex_id = index_y:stride_y:size(state,2)
            dec_id = pop_cuda!(stack, s_id, ex_id)
            if dec_id > zero(eltype(stack))
                if isleafgate(nl, dec_id)
                    if isliteralgate(nl, dec_id)
                        l = literal(nl, dec_id)
                        var = lit2var(l)
                        @inbounds state[s_id, ex_id, var] = (l > 0)
                    end
                else
                    edge_log_pr, prime, sub = sample_child_cuda(params, nodes, elements, s_id, ex_id, dec_id, values, r, Float)
                    @inbounds logprob[s_id, ex_id] += edge_log_pr
                    push_cuda!(stack, prime, s_id, ex_id)
                    push_cuda!(stack, sub, s_id, ex_id)
                end
            end
        end
    end
    return nothing
end

function sample_child_cuda(params, nodes, elements, s_id, ex_id, dec_id, values, r, ::Type{Float}) where Float
    @inbounds els_start = nodes[1,dec_id]
    @inbounds els_end = nodes[2,dec_id]
    @inbounds threshold = CUDA.log(r[s_id, ex_id]) + values[ex_id, dec_id]
    cumul_prob::Float = -Inf
    j_sampled = els_end - els_start + 1 # give all numerical error probability to the last node
    for j = els_start:els_end
        @inbounds prime = elements[2,j]
        @inbounds sub = elements[3,j]
        @inbounds pr::Float = values[ex_id, prime] + values[ex_id, sub] + params[j]
        cumul_prob = logsumexp_cuda(cumul_prob, pr)
        if cumul_prob > threshold
            j_sampled = j
            break
        end
    end
    @inbounds return params[j_sampled], elements[2,j_sampled], elements[3,j_sampled]
end


text:
const package_regex = r"package\s(\S*)[\s]*;.*"
const service_regex = r"service\s(\S*)[\s]*.*"

function write_header(io, generated_module, package, client_module_name)
    print(io, """module $(client_module_name)
    using gRPCClient

    include("$(generated_module).jl")
    using .$(package)

    import Base: show
    """)
end

function write_trailer(io, client_module_name)
    print(io, """

    end # module $(client_module_name)
    """)
end

function write_service(io, package, service, methods)
    print(io, """

    # begin service: $(package).$(service)

    export $(service)BlockingClient, $(service)Client

    struct $(service)BlockingClient
        controller::gRPCController
        channel::gRPCChannel
        stub::$(service)BlockingStub

        function $(service)BlockingClient(api_base_url::String; kwargs...)
            controller = gRPCController(; kwargs...)
            channel = gRPCChannel(api_base_url)
            stub = $(service)BlockingStub(channel)
            new(controller, channel, stub)
        end
    end

    struct $(service)Client
        controller::gRPCController
        channel::gRPCChannel
        stub::$(service)Stub

        function $(service)Client(api_base_url::String; kwargs...)
            controller = gRPCController(; kwargs...)
            channel = gRPCChannel(api_base_url)
            stub = $(service)Stub(channel)
            new(controller, channel, stub)
        end
    end

    show(io::IO, client::$(service)BlockingClient) = print(io, "$(service)BlockingClient(", client.channel.baseurl, ")")
    show(io::IO, client::$(service)Client) = print(io, "$(service)Client(", client.channel.baseurl, ")")
    """)

    for method in methods
        write_service_method(io, package, service, method)
    end

    print(io, """

    # end service: $(package).$(service)
    """)
end

typename(ch::Type{T}) where {T <: Channel} = string("Channel{", typename(eltype(ch)), "}")
typename(T) = last(split(string(T), '.'; limit=2))

function write_service_method(io, package, service, method)
    method_name = method.name
    input_type = typename(method.input_type)
    output_type = typename(method.output_type)

    print(io, """

    import .$(package): $(method_name)
    \"\"\"
        $(method_name)

    - input: $input_type
    - output: $output_type
    \"\"\"
    $(method_name)(client::$(service)BlockingClient, inp::$(input_type)) = $(method_name)(client.stub, client.controller, inp)
    $(method_name)(client::$(service)Client, inp::$(input_type), done::Function) = $(method_name)(client.stub, client.controller, inp, done)
    """)
end

function detect_services(proto::String)
    package = ""
    services = String[]

    for line in readlines(proto)
        line = strip(line)
        if startswith(line, "package")
            regexmatches = match(package_regex, line)
            if (regexmatches !== nothing) && (length(regexmatches.captures) == 1)
                package = string(first(regexmatches.captures))
            end
        elseif startswith(line, "service")
            regexmatches = match(service_regex, line)
            if (regexmatches !== nothing) && (length(regexmatches.captures) == 1)
                service = string(first(regexmatches.captures))
                push!(services, service)
            end
        end
    end
    package, services
end

function get_generated_method_table(s::String)
    T = Main
    for t in split(s, ".")
        T = Base.eval(T, Symbol(t))
    end
    T
end

"""
    generate(proto::String; outdir::String=pwd())

Generate a gRPC client from protobuf specification file.

- `proto`: Path to the protobuf specification to used.
- `outdir`: Directory to write generated code into, created if not present
    already. Existing files if any will be overwtitten.
"""
function generate(proto::String; outdir::String=pwd(), includes::Vector{String}=String[])
    if !isfile(proto)
        throw(ArgumentError("No such file - $proto"))
    end
    proto = abspath(proto)

    @info("Generating gRPC client", proto, outdir)

    # determine the package name and service name
    package, services = detect_services(proto)
    protodir = dirname(proto)
    includeflag = `-I=$protodir`
    for inc in includes
        includeflag = `$includeflag -I=$inc`
    end
    @info("Detected", package, services, includes)

    # generate protobuf services
    mkpath(outdir)
    bindir = Sys.BINDIR
    pathenv = string(ENV["PATH"], Sys.iswindows() ? ";" : ":", bindir)
    withenv("PATH"=>pathenv) do
        ProtoBuf.protoc(`$includeflag --julia_out=$outdir $proto`)
    end

    # include the generated code and detect service method names
    generated_module = first(split(package, '.'; limit=2))
    generated_module_file = joinpath(outdir, string(generated_module, ".jl"))
    Main.eval(:(include($generated_module_file)))

    # generate the gRPC client code
    client_module_name = string(titlecase(generated_module; strict=false), "Clients")
    open(joinpath(outdir, "$(client_module_name).jl"), "w") do grpcservice
        write_header(grpcservice, generated_module, package, client_module_name)
        for service in services
            methods = get_generated_method_table(string(package, "._", service, "_methods"))
            write_service(grpcservice, package, service, methods)
        end
        write_trailer(grpcservice, client_module_name)
    end

    @info("Generated", outdir)
end


text:
using NamedTupleTools
using Test

namesofvalues  = (:instrument, :madeby)
matchingvalues = ("violin", "Stradivarius")
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = [:instrument, :madeby]
matchingvalues = ["violin", "Stradivarius"]
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = (:instrument, :madeby)
matchingvalues = ["violin", "Stradivarius"]
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = [:instrument, :madeby]
matchingvalues = ("violin", "Stradivarius")
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

nt = NamedTuple{(:a, :b)}(1.0, "two")

@test fieldnames(nt) == (:a, :b)
@test fieldtypes(nt) == (Float64, String)
@test valtype(nt) == Tuple{Float64, String}
@test fieldtypes(nt) == (Float64, String)

namedtuple(:a, :b) == NamedTuple{(:a, :b),T} where T<:Tuple

ntproto1 = namedtuple(:a, :b, :c, :d)
ntproto2 = namedtuple(:a, :b)

@test ntproto1 == NamedTuple{(:a, :b, :c, :d),T} where T<:Tuple

nt1 = ntproto1(1, 2, 3, 4)
nt2 = ntproto2("one", "two")

@test isprototype(ntproto1) === true
@test isprototype(nt1) === false
@test isprototype(UnionAll) === false

@test nt2 === (a = "one", b = "two")
@test ntproto2(nt1) === (a = 1, b = 2)

@test delete(ntproto1, :a) === NamedTuple{(:b, :c, :d),T} where T<:Tuple
@test delete(ntproto1, :a, :c) === NamedTuple{(:b, :d),T} where T<:Tuple
@test delete(ntproto1, (:b, :c)) === NamedTuple{(:a, :d),T} where T<:Tuple

@test delete(nt1, :a) == (b = 2, c = 3, d = 4)
@test delete(nt1, :a, :c) == (b = 2, d = 4)
@test delete(nt1, (:a, :b, :c)) === (d = 4,)

@test merge(nt1, nt2) === (a = "one", b  = "two", c = 3, d = 4)

@test select(nt1, :a) == nt1[:a]
@test select(nt1, nt2) == (a=1,b=2)
@test select(nt1, nt2) == select(nt1, keys(nt2))

struct MyStruct
    tally::Int
    team::String
end

mystruct = MyStruct(5, "hometeam")
mynamedtuple = ntfromstruct(mystruct)

@test mynamedtuple == (tally = 5, team = "hometeam")

nt = (tally=5, team="hometeam")
ntstruct = structfromnt(MyStruct, nt)

@test ntstruct == mystruct

v = [:a => 1, :b => 2]
@test namedtuple(v) == NamedTuple{(:a, :b)}(1, 2)
v = ["a" => 1, "b" => 2]
@test namedtuple(v) == NamedTuple{(:a, :b)}(1, 2)
v = [:b => 2, :a => 1]
@test namedtuple(v) == NamedTuple{(:b, :a)}(2, 1)
v = ["b" => 2, "a" => 1]
@test namedtuple(v) == NamedTuple{(:b, :a)}(2, 1)

dict1 = Dict(:a=>1)
nt1 = (a = 1,)
dict3 = Dict(:a=>1, :b=>2//11, :c=>"three")
nt3 = (a = 1, b = 2//11, c = "three")

@test convert(Dict, nt1) == dict1
@test namedtuple(dict1) == nt1
@test convert(Dict, nt3) == dict3
@test namedtuple(dict3) == nt3

nt = (a = 1, b = 2)
a = 1; b = 2;
nt_ab = @namedtuple(a, b)
@test nt_ab == nt


text:
""" Vertical circulation associated with horizontal velocities u, v. """
@inline Γᶠᶠᶜ(i, j, k, grid, u, v) = δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, v) - δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, u)

""" Vertical vorticity associated with horizontal velocities u, v. """
@inline ζ₃ᶠᶠᶜ(i, j, k, grid, u, v) = Γᶠᶠᶜ(i, j, k, grid, u, v) / Azᶠᶠᶜ(i, j, k, grid)

#####
##### Vertical circulation at the corners of the cubed sphere needs to treated in a special manner.
##### See: https://github.com/CliMA/Oceananigans.jl/issues/1584
#####

@inline function Γᶠᶠᶜ(i, j, k, grid::ConformalCubedSphereFaceGrid, u, v)
    # South-west corner
    if i == 1 && j == 1
        return Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u)

    # South-east corner
    elseif i == grid.Nx+1 && j == 1
        return - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u)

    # North-west corner
    elseif i == 1 && j == grid.Ny+1
        return Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u)

    # North-east corner
    elseif i == grid.Nx+1 && j == grid.Ny+1
        return - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u)

    # Not a corner
    else
        return Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u)
    end
end


text:
using APCA
using Test
using Random

@testset "APCA.jl" begin
    v= [7.,5.,5.,3.,3.,3.,4.,6.]
    M= 3
    @test runAPCA(v,M) == [6.,6.,3.5,3.5,3.5,3.5,5.,5.]

    Random.seed!(1234)
    v = vcat(rand(10).+3,rand(10).+2,rand(10).-1)
    M = 4
    @test runAPCA(v,M) == [3.1124992225766124, 3.1124992225766124, 3.1124992225766124, 3.1124992225766124,
    3.1124992225766124, 3.1124992225766124, 3.1124992225766124, 3.1124992225766124,
    3.1124992225766124, 3.1124992225766124, 3.1124992225766124, 3.1124992225766124, 3.1124992225766124,
    3.1124992225766124, 3.1124992225766124, 3.1124992225766124, 2.456760310559139, 2.456760310559139, 2.456760310559139,
    2.456760310559139, -0.08728818214843481, -0.08728818214843481, -0.08728818214843481, -0.08728818214843481, -0.7055353573422534,
    -0.7055353573422534, -0.7055353573422534, -0.7055353573422534, -0.7055353573422534, -0.7055353573422534]

end


text:
## This file provides a wrapper to Convex.jl for solving nonnegative
## least-squares problems. In general, the active-set methods seem
## to outperform Convex.jl, but this function can be used to directly
## compare the two approaches.

using Convex
using SCS
using ECOS

"""
x = convex_nnls(A, b, solver; kwargs...)

Solves the nonnegative least-squares problem using Convex.jl and the solver
specified by "solver".
"""
function convex_nnls(
	A::Matrix{Float64},
	B::Matrix{Float64};
	variant::Symbol = :SCS,
	kwargs...)

	# Primal optimization variables
	X = Variable(size(A,2),size(B,2))

	# declare problem
	problem  = minimize(sumsquares(B - (A*X)), [X >= 0])

	# solve problem
	if variant == :SCS
		solve!(problem,SCSSolver(kwargs))
	elseif variant == :ECOS
		solve!(problem,ECOSSolver(kwargs))
	else
		error("Specified solver :",variant," not recognized")
	end

	# check solution
	if problem.status != :Optimal
		warn("Solving for H problem status is :",problem.status)
	end

	# return solution
	return X.value
end


text:
#DECORRELATION SEARCH PATTERNS
#random permutation of single sources until model with log likelihood>contour found or iterates limit reached. will always produce at least one weight shift per iteration for weight_shift_freq>0, more than this or length changes depend on the supplied probabilities. one length change per iterate
function permute_source(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer=length(m.sources), weight_shift_freq::AbstractFloat=.5, length_change_freq::AbstractFloat=1.0, weight_shift_dist::Distributions.ContinuousUnivariateDistribution=Weibull(1.5,.1))
#weight_shift_dist is given in decimal probability values- converted to log space in permute_source_lengths!
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources);
    flags=deepcopy(m.flags); flags[1]="PS from $(m.name)"

    a, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources);
        s = rand(1:S)
        clean=Vector{Bool}(trues(O))
        clean[m.mix_matrix[:,s]].=false #all obs with source are dirty

        weight_shift_freq > 0 && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits))
        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, m.mix_matrix, new_log_Li, flags), time()-instrct_start) : (return consolidate_srcs(cons_idxs, new_sources, m.mix_matrix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits), time()-instrct_start)
end

function permute_mix(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer=10, mix_move_range::UnitRange=Int(ceil(.001*length(m.mix_matrix))):length(m.mix_matrix))
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_mix=falses(size(m.mix_matrix))
    flags=deepcopy(m.flags); flags[1]="PM from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"),flags))

    a, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)
    dirty=false

    while (new_log_Li <= contour || !dirty) && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        mix_moves=rand(mix_move_range)
        mix_moves > length(m.mix_matrix) && (mix_moves = length(m.mix_matrix))

        new_mix, clean = mix_matrix_decorrelate(m.mix_matrix, mix_moves) #generate a decorrelated candidate mix
        c_log_li, c_cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #calculate the model with the candidate mix
        positive_indices=c_cache.>(cache) #obtain any obs indices that have greater probability than we started with
        clean=Vector{Bool}(trues(O)-positive_indices)
        if any(positive_indices) #if there are any such indices
            new_log_Li, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #calculate the new model
            dirty=true
        end
        iterate += 1
    end

    return ICA_PWM_model("candidate",m.sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li, flags)
end

function perm_src_fit_mix(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer=length(m.sources), weight_shift_freq::AbstractFloat=.25, length_change_freq::AbstractFloat=.5, weight_shift_dist::Distributions.ContinuousUnivariateDistribution=Weibull(1.0,.1))
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    flags=deepcopy(m.flags); flags[1]="PSFM from $(m.name)"

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix); tm_one=deepcopy(m.mix_matrix);
        clean=Vector{Bool}(trues(O))
        s = rand(1:S);
        clean[new_mix[:,s]].=false #all obs starting with source are dirty

        new_mix[:,s].=false;tm_one[:,s].=true

        weight_shift_freq > 0 && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits))

        l,zero_cache=IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true)
        l,one_cache=IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, tm_one, true, true)

        fit_mix=one_cache.>=zero_cache

        new_mix[:,s]=fit_mix

        clean[fit_mix].=false #all obs ending with source are dirty

        new_log_Li = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, false, zero_cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li, flags)) : (return consolidate_srcs(cons_idxs, new_sources, new_mix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function fit_mix(m::ICA_PWM_model, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, exclude_src::Integer=0)
    instrct_start=time(); T,O = size(observations); T=T-1; S = length(m.sources)
    new_mix=deepcopy(m.mix_matrix); test_mix=falses(size(m.mix_matrix))
    flags=deepcopy(m.flags); flags[1]="FM from $(m.name)"; push!(flags,"nofit")

    l,zero_cache=IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, test_mix, true, true)

    for (s,source) in enumerate(m.sources)
        if s!=exclude_src
            test_mix=falses(size(m.mix_matrix))
            test_mix[:,s].=true

            l,src_cache=IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, test_mix, true, true)

            fit_mix=src_cache.>=zero_cache
            new_mix[:,s]=fit_mix
        end
    end

    new_mix==m.mix_matrix ? new_log_Li=-Inf : new_log_Li = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, new_mix, true)
    return ICA_PWM_model("candidate",m.sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li, flags)
end

function random_decorrelate(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer=length(m.sources), weight_shift_freq::AbstractFloat=.1, length_change_freq::AbstractFloat=.3, weight_shift_dist::Distributions.ContinuousUnivariateDistribution=Weibull(1.5,.1), mix_move_range::UnitRange=1:Int(ceil(size(m.mix_matrix,1)*.1)))
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    flags=deepcopy(m.flags); flags[1]="RD from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"),flags))

    a, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))
        s = rand(1:length(m.sources))
        clean[new_mix[:,s]].=false #all obs starting with source are dirty
        weight_shift_freq > 0 && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits))
        new_mix[:,s]=mixvec_decorrelate(new_mix[:,s],rand(mix_move_range))
        clean[new_mix[:,s]].=false #all obs ending with source are dirty
        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true, cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li, ["RD from $(m.name)"])) : (return consolidate_srcs(cons_idxs, new_sources, new_mix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function distance_merge(models::Vector{nnlearn.Model_Record}, m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    flags=deepcopy(m.flags); flags[1]="DM from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"), flags))

    a, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        merger_m = deserialize(rand(models).path) #randomly select a model to merge
        s = rand(1:S) #randomly select a source to merge
        s > m.informed_sources ? #if the source is on an uninformative mix prior, the merger model source will be selected by mixvector similarity
            merge_s=most_dissimilar(new_mix,merger_m.mix_matrix) : merge_s=s

        clean[new_mix[:,s]].=false #mark dirty any obs that start with the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source
        new_mix[:,s] = merger_m.mix_matrix[:,merge_s]
        clean[new_mix[:,s]].=false #mark dirty any obs that end with the source

        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #assess likelihood
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, new_mix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function distance_merge(remote_models::Vector{nnlearn.Model_Record}, m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    flags=deepcopy(m.flags); flags[1]="DM from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"),flags))

    a, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))
        merger_m = remotecall_fetch(deserialize, 1, rand(remote_models).path) #randomly select a model to merge
        s = rand(1:S) #randomly select a source to merge
        s > m.informed_sources ? #if the source is on an uninformative prior, the merger model source will be selected by mixvector dissimilarity
            merge_s=most_dissimilar(new_mix,merger_m.mix_matrix) : merge_s=s

        clean[new_mix[:,s]].=false #mark dirty any obs that start with the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source
        new_mix[:,s] = merger_m.mix_matrix[:,merge_s]
        clean[new_mix[:,s]].=false #mark dirty any obs that end with the source

        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #assess likelihood
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, new_mix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function similarity_merge(models::Vector{nnlearn.Model_Record}, m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources)
    flags=deepcopy(m.flags); flags[1]="SM from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"),flags))

    a, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        merger_m = deserialize(rand(models).path) #randomly select a model to merge
        s = rand(1:S) #randomly select a source in the model to merge
        s > m.informed_sources ? #if the source is on an uninformative prior, the merger model source will be selected by mixvector similarity
        merge_s=most_similar(m.mix_matrix[:,s],merger_m.mix_matrix) : merge_s=s

        clean[m.mix_matrix[:,s]].=false #mark dirty any obs that have the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source

        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, m.mix_matrix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, m.mix_matrix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function similarity_merge(remote_models::Vector{nnlearn.Model_Record}, m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources)
    flags=deepcopy(m.flags); flags[1]="SM from $(m.name)"
    "nofit" in flags && deleteat!(flags, findfirst(isequal("nofit"),flags))

    a, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)
    clean=Vector{Bool}(trues(O))

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        merger_m = remotecall_fetch(deserialize, 1, rand(remote_models).path) #randomly select a model to merge
        s = rand(1:S) #randomly select a source in the model to merge
        s > m.informed_sources ? #if the source is on an uninformative prior, the merger model source will be selected by mixvector similarity
        merge_s=most_similar(m.mix_matrix[:,s],merger_m.mix_matrix) : merge_s=s

        clean[m.mix_matrix[:,s]].=false #mark dirty any obs that have the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source

        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood

        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, m.mix_matrix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, m.mix_matrix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

function reinit_src(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, mix_prior::Tuple{BitMatrix,AbstractFloat}, iterates::Integer, inform_mix_prior::Bool=true)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    flags=deepcopy(m.flags); flags[1]="RS from $(m.name)"

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix); tm_one=deepcopy(m.mix_matrix);
        clean=Vector{Bool}(trues(O))
        s = rand(1:S);
        clean[new_mix[:,s]].=false #all obs starting with source are dirty

        new_mix[:,s].=false;tm_one[:,s].=true

        new_sources[s] = init_logPWM_sources([source_priors[s]], m.source_length_limits)[1] #reinitialize the source

        l,zero_cache=IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, true)
        l,one_cache=IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, tm_one, true, true)

        fit_mix=one_cache.>=zero_cache

        new_mix[:,s]=fit_mix

        clean[fit_mix].=false #all obs ending with source are dirty

        new_log_Li = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, new_mix, true, false, zero_cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, new_mix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, new_mix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

Set()

function erode_model(m::ICA_PWM_model, contour::AbstractFloat, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, iterates::Integer, info_thresh::AbstractFloat=1.)
    instrct_start=time(); new_log_Li=-Inf;  iterate = 1
    T,O = size(observations); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources)
    flags=deepcopy(m.flags); flags[1]="EM from $(m.name)"

    erosion_sources=Set{Integer}()
    for (s,src) in enumerate(m.sources)
        pwm,pi=src
        if size(pwm,1)>m.source_length_limits[1]
            infovec=get_pwm_info(pwm)
            any(info-><(info, info_thresh),infovec) && push!(erosion_sources,s)
        end
    end

    length(erosion_sources)==0 && return perm_src_fit_mix(m, contour, observations,obs_lengths,bg_scores,source_priors, iterates)#if we got a model we cant erode bail out to PSFM

    a, cache = IPM_likelihood(m.sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && length(erosion_sources) > 0 #until we produce a model more likely than the lh contour or there are no more sources to erode
        clean=Vector{Bool}(trues(O))
        s=pop!(erosion_sources)

        new_sources[s]=erode_source(new_sources[s], m.source_length_limits, info_thresh)
        clean[m.mix_matrix[:,s]].=false

        new_log_Li, cache = IPM_likelihood(new_sources, observations, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_model("candidate",new_sources, m.informed_sources, m.source_length_limits, m.mix_matrix, new_log_Li,flags)) : (return consolidate_srcs(cons_idxs, new_sources, m.mix_matrix, observations, obs_lengths, bg_scores, source_priors, m.informed_sources, m.source_length_limits))
end

##ORTHOGONALITY HELPER
function consolidate_srcs(con_idxs::Vector{Vector{Integer}}, sources::Vector{Tuple{Matrix{AbstractFloat},Integer}}, mix::BitMatrix, observations::Matrix{Integer}, obs_lengths::Vector{Integer}, bg_scores::Matrix{AbstractFloat}, source_priors::Vector{Vector{Dirichlet{AbstractFloat}}}, informed_sources::Integer, source_length_limits::UnitRange)
    new_sources=deepcopy(sources);new_mix=deepcopy(mix)
    consrc=0
    for (s,convec) in enumerate(con_idxs)
        if length(convec)>0
            consrc=s
            for cons in convec
                new_sources[cons]=init_logPWM_sources([source_priors[cons]], source_length_limits)[1]
                new_mix[:,s]=[any([new_mix[i,s],mix[i,cons]]) for i in 1:size(mix,1)] #consolidate the new_mix source with the consolidated mix
            end
            break #dont do this for more than one consolidation at a time
        end
    end

    return fit_mix(ICA_PWM_model("consolidate", new_sources, informed_sources, source_length_limits, new_mix, -Inf, [""]), observations, obs_lengths, bg_scores, consrc)
end

function consolidate_check(sources::Vector{Tuple{Matrix{AbstractFloat},Integer}}; thresh=.035)
    pass=true
    con_idxs=Vector{Vector{Integer}}()
    for (s1,src1) in enumerate(sources)
        s1_idxs=Vector{Integer}()
        for (s2,src2) in enumerate(sources)
            if !(s1==s2)
                pwm1=src1[1]; pwm2=src2[1]
                if -3<=(size(pwm1,1)-size(pwm2,1))<=3
                    if pwm_distance(pwm1,pwm2)<thresh
                        push!(s1_idxs,s2)
                        pass=false
                    end
                end
            end
        end
        push!(con_idxs,s1_idxs)
    end
    return pass, con_idxs
end

                function pwm_distance(pwm1,pwm2)
                    minwml=min(size(pwm1,1),size(pwm2,1))
                    return sum([euclidean(exp.(pwm1[pos,:]), exp.(pwm2[pos,:])) for pos in 1:minwml])/minwml
                end

##BASIC UTILITY FUNCTIONS
#SOURCE PERMUTATION
function permute_source_weights(source::Tuple{Matrix{AbstractFloat},Integer}, shift_freq::AbstractFloat, PWM_shift_dist::Distributions.ContinuousUnivariateDistribution)
    dirty=false; source_length=size(source[1],1)
    new_source=deepcopy(source)

    for source_pos in 1:source_length
        if rand() <= shift_freq
            pos_WM = exp.(source[1][source_pos,:]) #leaving logspace, get the wm at that position
            new_source[1][source_pos,:] = log.(wm_shift(pos_WM, PWM_shift_dist)) #accumulate probabilty at a randomly selected base, reassign in logspace and carry on
            !dirty && (dirty=true)
        end
    end

    if !dirty #if no positions were shifted, pick one and shift
        rand_pos=rand(1:source_length)
        pos_WM = exp.(source[1][rand_pos,:])
        new_source[1][rand_pos,:]=log.(wm_shift(pos_WM, PWM_shift_dist))
    end

    return new_source
end

                function wm_shift(pos_WM::Vector{AbstractFloat}, PWM_shift_dist::Distributions.ContinuousUnivariateDistribution)
                    base_to_shift = rand(1:4) #pick a base to accumulate probability
                    permute_sign = rand(-1:2:1)
                    shift_size = rand(PWM_shift_dist)
                    new_wm=zeros(4)

                    for base in 1:4 #ACGT
                        if base == base_to_shift
                            new_wm[base] =
                            clamp(0, #no lower than 0 prob
                            (pos_WM[base]          #selected PWM posn
                            + permute_sign * shift_size), #randomly permuted by size param
                            1) #no higher than prob 1
                        else
                            size_frac = shift_size / 3 #other bases shifted in the opposite direction by 1/3 the shift accumulated at the base to permute
                            new_wm[base] =
                            clamp(0,
                            (pos_WM[base]
                            - permute_sign * size_frac),
                            1)
                        end
                    end
                    new_wm = new_wm ./ sum(new_wm) #renormalise to sum 1 - necessary in case of clamping at 0 or 1
                    !isprobvec(new_wm) && throw(DomainError(new_wm, "Bad weight vector generated in wm_shift!")) #throw assertion exception if the position WM is invalid
                    return new_wm
                end

function permute_source_length(source::Tuple{Matrix{AbstractFloat},Integer}, prior::Vector{Dirichlet{AbstractFloat}}, length_limits::UnitRange{Integer}, permute_range::UnitRange{Integer}=1:3, uninformative::Dirichlet=Dirichlet([.25,.25,.25,.25]))
    source_PWM, prior_idx = source
    source_length = size(source_PWM,1)

    permute_sign, permute_length = get_length_params(source_length, length_limits, permute_range)

    permute_sign==1 ? permute_pos = rand(1:source_length+1) :
        permute_pos=rand(1:source_length-permute_length)

    if permute_sign == 1 #if we're to add positions to the PWM
        ins_WM=zeros(permute_length,4)
        for pos in 1:permute_length
            prior_position=permute_pos+prior_idx
            prior_position<1 || prior_position>length(prior) ?
                ins_WM[pos,:] = log.(transpose(rand(uninformative))) :
                ins_WM[pos,:] = log.(transpose(rand(prior[prior_position])))
                !isprobvec(exp.(ins_WM[pos,:])) && throw(DomainError(ins_WM, "Bad weight vector generated in permute_source_length!"))
        end
        upstream_source=source_PWM[1:permute_pos-1,:]
        downstream_source=source_PWM[permute_pos:end,:]
        source_PWM=vcat(upstream_source,ins_WM,downstream_source)
        permute_pos==1 && (prior_idx-=permute_length)
    else #if we're to remove positions
        upstream_source=source_PWM[1:permute_pos-1,:]
        downstream_source=source_PWM[permute_pos+permute_length:end,:]
        source_PWM=vcat(upstream_source,downstream_source)
        permute_pos==1 && (prior_idx+=permute_length)
    end

    return (source_PWM, prior_idx) #return a new source
end

                function get_length_params(source_length::Integer, length_limits::UnitRange{Integer}, permute_range::UnitRange{Integer})
                    extendable = length_limits[end]-source_length
                    contractable =  source_length-length_limits[1]

                    if extendable == 0 && contractable > 0
                        permute_sign=-1
                    elseif contractable == 0 && extendable > 0
                        permute_sign=1
                    else
                        permute_sign = rand(-1:2:1)
                    end

                    permute_sign==1 && extendable<permute_range[end] && (permute_range=permute_range[1]:extendable)
                    permute_sign==-1 && contractable<permute_range[end] && (permute_range=permute_range[1]:contractable)
                    permute_length = rand(permute_range)

                    return permute_sign, permute_length
                end

function erode_source(source::Tuple{Matrix{AbstractFloat},Integer},length_limits::UnitRange{Integer},info_thresh)
    pwm,prior_idx=source
    infovec=get_pwm_info(pwm)
    start_idx,end_idx=get_erosion_idxs(infovec, info_thresh, length_limits)

    return new_source=(pwm[start_idx:end_idx,:], prior_idx+start_idx-1)
end

    function get_pwm_info(pwm::Matrix{AbstractFloat}; logsw::Bool=true)
        wml=size(pwm,1)
        infovec=zeros(wml)
        for pos in 1:wml
            logsw ? wvec=deepcopy(exp.(pwm[pos,:])) : wvec=deepcopy(pwm[pos,:])
            !isprobvec(wvec) && throw(DomainError(wvec, "Bad wvec in get_pwm_info -Original sources must be in logspace!!"))
            wvec.+=10^-99
            infscore = (2.0 + sum([x*log(2,x) for x in wvec]))
            infovec[pos]=infscore
        end
        return infovec
    end

    function get_erosion_idxs(infovec::Vector{AbstractFloat}, info_thresh::AbstractFloat, length_limits::UnitRange{Integer})
        srcl=length(infovec)
        contractable =  srcl-length_limits[1]
        contractable <=0 && throw(DomainError(contractable, "erode_source passed a source at its lower length limit!"))
        centeridx=findmax(infovec)[2]

        start_idx=findprev(info-><(info,info_thresh),infovec,centeridx)
        start_idx===nothing ? (start_idx=1) : (start_idx+=1)
        end_idx=findnext(info-><(info, info_thresh),infovec,centeridx)
        end_idx===nothing ? (end_idx=srcl) : (end_idx-=1)

        pos_to_erode=srcl-(end_idx-start_idx)
        if pos_to_erode > contractable
            pos_to_restore = pos_to_erode-contractable
            while pos_to_restore>0
                end_die=rand()
                if end_die <= .5
                    start_idx>1 && (pos_to_restore-=1; start_idx-=1)
                else
                    end_idx<srcl && (pos_to_restore-=1; end_idx+=1)
                end
            end
        end

        return start_idx, end_idx
    end

#MIX MATRIX FUNCTIONS
function mixvec_decorrelate(mix::BitVector, moves::Integer)
    new_mix=deepcopy(mix)
    idxs_to_flip=rand(1:length(mix), moves)
    new_mix[idxs_to_flip] .= .!mix[idxs_to_flip]
    return new_mix
end

function mix_matrix_decorrelate(mix::BitMatrix, moves::Integer)
    clean=Vector{Bool}(trues(size(mix,1)))
    new_mix=deepcopy(mix)
    indices_to_flip = rand(CartesianIndices(mix), moves)
    new_mix[indices_to_flip] .= .!mix[indices_to_flip]
    clean[unique([idx[1] for idx in indices_to_flip])] .= false #mark all obs that had flipped indices dirty
    return new_mix, clean
end

function most_dissimilar(mix1, mix2)
    S1=size(mix1,2);S2=size(mix2,2)
    dist_mat=zeros(S1,S2)
    for s1 in 1:S1, s2 in 1:S2
        dist_mat[s1,s2]=sum(mix1[:,s1].==mix2[:,s2])
    end
    scores=vec(sum(dist_mat,dims=1))
    return findmin(scores)[2]
end

function most_similar(src_mixvec, target_mixmat)
    src_sim = [sum(src_mixvec.==target_mixmat[:,s]) for s in 1:size(target_mixmat,2)] #compose array of elementwise equality comparisons between mixvectors and sum to score
    merge_s=findmax(src_sim)[2] #source from merger model will be the one with the highest equality comparison score
end


text:
module TroposphereWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Class = @jimport java.lang.Class
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const DateTimeComponents = @jimport org.orekit.time.DateTimeComponents
const DiscreteTroposphericModel = @jimport org.orekit.models.earth.troposphere.DiscreteTroposphericModel
const Enum = @jimport java.lang.Enum
const EstimatedTroposphericModel = @jimport org.orekit.models.earth.troposphere.EstimatedTroposphericModel
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FixedTroposphericDelay = @jimport org.orekit.models.earth.troposphere.FixedTroposphericDelay
const GlobalMappingFunctionModel = @jimport org.orekit.models.earth.troposphere.GlobalMappingFunctionModel
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const MappingFunction = @jimport org.orekit.models.earth.troposphere.MappingFunction
const MariniMurrayModel = @jimport org.orekit.models.earth.troposphere.MariniMurrayModel
const MendesPavlisModel = @jimport org.orekit.models.earth.troposphere.MendesPavlisModel
const NavigableSet = @jimport java.util.NavigableSet
const NiellMappingFunctionModel = @jimport org.orekit.models.earth.troposphere.NiellMappingFunctionModel
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SaastamoinenModel = @jimport org.orekit.models.earth.troposphere.SaastamoinenModel
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeSpanEstimatedTroposphericModel = @jimport org.orekit.models.earth.troposphere.TimeSpanEstimatedTroposphericModel
const ViennaModelCoefficientsLoader = @jimport org.orekit.models.earth.troposphere.ViennaModelCoefficientsLoader
const ViennaModelType = @jimport org.orekit.models.earth.troposphere.ViennaModelType
const ViennaOneModel = @jimport org.orekit.models.earth.troposphere.ViennaOneModel
const ViennaThreeModel = @jimport org.orekit.models.earth.troposphere.ViennaThreeModel

export DiscreteTroposphericModel
export EstimatedTroposphericModel
export FixedTroposphericDelay
export GlobalMappingFunctionModel
export MappingFunction
export MariniMurrayModel
export MendesPavlisModel
export NiellMappingFunctionModel
export SaastamoinenModel
export TimeSpanEstimatedTroposphericModel
export ViennaModelCoefficientsLoader
export ViennaModelType
export ViennaOneModel
export ViennaThreeModel
export add_tropospheric_model_valid_after
export add_tropospheric_model_valid_before
export compare_to
export compute_zenith_delay
export describe_constable
export equals
export extract_parameters
export get_a
export get_class
export get_declaring_class
export get_default_model
export get_low_elevation_threshold
export get_parameters
export get_parameters_drivers
export get_standard_model
export get_supported_names
export get_transitions
export get_tropospheric_model
export get_zenith_delay
export hash_code
export load_data
export load_vienna_coefficients
export mapping_factors
export name
export notify
export notify_all
export ordinal
export path_delay
export set_low_elevation_threshold
export still_accepts_data
export to_string
export value_of
export values
export wait

include("discrete_tropospheric_model.jl")
include("estimated_tropospheric_model.jl")
include("fixed_tropospheric_delay.jl")
include("global_mapping_function_model.jl")
include("mapping_function.jl")
include("marini_murray_model.jl")
include("mendes_pavlis_model.jl")
include("niell_mapping_function_model.jl")
include("saastamoinen_model.jl")
include("time_span_estimated_tropospheric_model.jl")
include("vienna_model_coefficients_loader.jl")
include("vienna_model_type.jl")
include("vienna_one_model.jl")
include("vienna_three_model.jl")

end


text:
module fileLog_wind

include( "LogLstLib.jl")

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
#####
##### `norm`
#####

function frule((_, Δx), ::typeof(norm), x)
    y = norm(x)
    return y, _norm2_forward(x, Δx, norm(x))
end

function frule((_, ẋ), ::typeof(norm), x::Number, p::Real)
    Δx = unthunk(ẋ)
    y = norm(x, p)
    ∂y = if iszero(Δx) || iszero(p)
        zero(real(x)) * zero(real(Δx))
    else
        signx = x isa Real ? sign(x) : x * pinv(y)
        _realconjtimes(signx, Δx)
    end
    return y, ∂y
end

function rrule(::typeof(norm), x::AbstractArray{<:Number}, p::Real)
    y = LinearAlgebra.norm(x, p)
    function norm_pullback_p(ȳ)
        Δy = unthunk(ȳ)
        ∂x = InplaceableThunk(
            # out-of-place versions
            @thunk(if isempty(x) || p == 0
                zero.(x) .* (zero(y) * zero(real(Δy)))
            elseif p == 2
                _norm2_back(x, y, Δy)
            elseif p == 1
                _norm1_back(x, y, Δy)
            elseif p == Inf
                _normInf_back(x, y, Δy)
            elseif p == -Inf
                _normInf_back(x, y, Δy)
            else
                _normp_back_x(x, p, y, Δy)
            end)
            , # in-place versions -- can be fixed when actually useful?
            dx -> if isempty(x) || p == 0
                dx
            elseif p == 2
                _norm2_back!(dx, x, y, Δy)
            elseif p == 1
                _norm1_back!(dx, x, y, Δy)
            elseif p == Inf
                dx .+= _normInf_back(x, y, Δy)  # not really in-place! could perhaps be improved
            elseif p == -Inf
                dx .+= _normInf_back(x, y, Δy)
            else
                dx .+= _normp_back_x(x, p, y, Δy)
            end
            )
        ∂p = @thunk _normp_back_p(x, p, y, Δy)
        return (NoTangent(), ∂x, ∂p)
    end
    norm_pullback_p(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent())
    return y, norm_pullback_p
end

function rrule(::typeof(norm), x::AbstractArray{<:Number})
    y = LinearAlgebra.norm(x)
    function norm_pullback_2(ȳ)
        Δy = unthunk(ȳ)
        ∂x = InplaceableThunk(
            @thunk(if isempty(x)
                zero.(x) .* (zero(y) * zero(real(Δy)))
            else
                _norm2_back(x, y, Δy)
            end)
            ,
            dx -> if isempty(x)
                dx
            else
                _norm2_back!(dx, x, y, Δy)
            end
            )
        return (NoTangent(), ∂x)
    end
    norm_pullback_2(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, norm_pullback_2
end

function rrule(::typeof(norm), x::LinearAlgebra.AdjOrTransAbsVec{<:Number}, p::Real)
    y, inner_pullback = rrule(norm, parent(x), p)
    function norm_pullback(Δy)
        (∂self, ∂x′, ∂p) = inner_pullback(Δy)
        fdual = x isa Transpose ? transpose : adjoint
        ∂x = @thunk fdual(unthunk(∂x′))
        return (∂self, ∂x, ∂p)
    end
    return y, norm_pullback
end

function rrule(::typeof(norm), x::Number, p::Real)
    y = norm(x, p)
    function norm_pullback(ȳ)
        Δy = unthunk(ȳ)
        ∂x = if iszero(Δy) || iszero(p)
            zero(x) * zero(real(Δy))
        else
            signx = x isa Real ? sign(x) : x * pinv(y)
            signx * real(Δy)
        end
        return (NoTangent(), ∂x, ZeroTangent())
    end
    norm_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent())
    return y, norm_pullback
end

#####
##### `normp`
#####

function rrule(::typeof(LinearAlgebra.normp), x::AbstractArray{<:Number}, p)
    y = LinearAlgebra.normp(x, p)
    function normp_pullback(ȳ)
        Δy = unthunk(ȳ)
        ∂x = @thunk _normp_back_x(x, p, y, Δy)
        ∂p = @thunk _normp_back_p(x, p, y, Δy)
        return (NoTangent(), ∂x, ∂p)
    end
    normp_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent())
    return y, normp_pullback
end

function _normp_back_x(x, p, y, Δy)
    c = real(Δy) / y
    ∂x = map(x) do xi
        a = norm(xi)
        ∂xi = xi * ((a / y)^(p - 2) * c)
        return ifelse(isfinite(∂xi), ∂xi, zero(∂xi))
    end
    return ∂x
end

function _normp_back_x(x::WithSomeZeros, p, y, Δy) # Diagonal, UpperTriangular, etc.
    c = real(Δy) / y
    ∂x_data = map(parent(x)) do xi
        a = norm(xi)
        ∂xi = xi * ((a / y)^(p - 2) * c)
        return ifelse(isfinite(∂xi), ∂xi, zero(∂xi))
    end
    return withsomezeros_rewrap(x, ∂x_data)
end

function _normp_back_p(x, p, y, Δy)
    y > 0 && isfinite(y) && !iszero(p) || return zero(real(Δy)) * zero(y) / one(p)
    s = sum(x) do xi
        a = norm(xi)
        c = (a / y)^(p - 1) * a * log(a)
        return ifelse(isfinite(c), c, zero(c))
    end
    ∂p = real(Δy) * (s - y * log(y)) / p
    return ∂p
end

#####
##### `normMinusInf`/`normInf`
#####

function rrule(::typeof(LinearAlgebra.normMinusInf), x::AbstractArray{<:Number})
    y = LinearAlgebra.normMinusInf(x)
    normMinusInf_pullback(Δy) = (NoTangent(), _normInf_back(x, y, Δy))
    normMinusInf_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, normMinusInf_pullback
end

function rrule(::typeof(LinearAlgebra.normInf), x::AbstractArray{<:Number})
    y = LinearAlgebra.normInf(x)
    normInf_pullback(Δy) = (NoTangent(), _normInf_back(x, y, Δy))
    normInf_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, normInf_pullback
end

function _normInf_back(x, y, Δy)
    Δu = real(Δy)
    T = typeof(zero(float(eltype(x))) * zero(Δu))
    ∂x = fill!(similar(x, T), 0)
    # if multiple `xi`s have the exact same norm, then they must have been identically
    # produced, e.g. with `fill`. So we set only one to be non-zero.
    # we choose last index to match the `frule`.
    yind = findlast(xi -> norm(xi) == y, x)
    yind === nothing && throw(ArgumentError("y is not the correct norm of x"))
    @inbounds ∂x[yind] = sign(x[yind]) * Δu
    return ∂x
end

#####
##### `norm1`
#####

function rrule(::typeof(LinearAlgebra.norm1), x::AbstractArray{<:Number})
    y = LinearAlgebra.norm1(x)
    norm1_pullback(Δy) = (NoTangent(), InplaceableThunk(
        @thunk(_norm1_back(x, y, Δy)),
        dx -> _norm1_back!(dx, x, y, Δy),
    ))
    norm1_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, norm1_pullback
end

function _norm1_back(x, y, Δy)
    ∂x = sign.(x) .* real(Δy)
    return ∂x
end
function _norm1_back(x::WithSomeZeros, y, Δy)
    ∂x_data = sign.(parent(x)) .* real(Δy)
    return withsomezeros_rewrap(x, ∂x_data)
end
function _norm1_back!(∂x, x, y, Δy)
    ∂x .+= sign.(x) .* real(Δy)
    return ∂x
end

#####
##### `norm2`
#####

function frule((_, Δx), ::typeof(LinearAlgebra.norm2), x)
    y = LinearAlgebra.norm2(x)
    return y, _norm2_forward(x, Δx, y)
end

function rrule(::typeof(LinearAlgebra.norm2), x::AbstractArray{<:Number})
    y = LinearAlgebra.norm2(x)
    norm2_pullback(Δy) = (NoTangent(), InplaceableThunk(
        @thunk(_norm2_back(x, y, Δy)),
        dx -> _norm2_back!(dx, x, y, Δy),
    ))
    norm2_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, norm2_pullback
end

function _norm2_forward(x, Δx, y)
    ∂y = real(dot(x, Δx)) * pinv(y)
    return ∂y
end
function _norm2_back(x, y, Δy)
    ∂x = x .* (real(Δy) * pinv(y))
    return ∂x
end
function _norm2_back(x::WithSomeZeros, y, ȳ)
    Δy = unthunk(ȳ)
    T = typeof(one(eltype(x)) / one(real(eltype(Δy))))
    ∂x_data = parent(x) .* (real(Δy) * pinv(y))
    return withsomezeros_rewrap(x, ∂x_data)
end
function _norm2_back!(∂x, x, y, Δy)
    ∂x .+= x .* (real(Δy) * pinv(y))
    return ∂x  # must return after mutating
end

#####
##### `normalize`
#####

function rrule(::typeof(normalize), x::AbstractVector{<:Number}, p::Real)
    nrm, inner_pullback = rrule(norm, x, p)
    Ty = typeof(first(x) / nrm)
    y = copyto!(similar(x, Ty), x)
    LinearAlgebra.__normalize!(y, nrm)
    function normalize_pullback(Δy)
        invnrm = pinv(nrm)
        ∂nrm = -dot(y, Δy) * invnrm
        (_, ∂xnorm, ∂p) = inner_pullback(∂nrm)
        ∂x = @thunk unthunk(∂xnorm) .+ Δy .* invnrm
        return (NoTangent(), ∂x, ∂p)
    end
    normalize_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent())
    return y, normalize_pullback
end

function rrule(::typeof(normalize), x::AbstractVector{<:Number})
    nrm = LinearAlgebra.norm2(x)
    Ty = typeof(first(x) / nrm)
    y = copyto!(similar(x, Ty), x)
    LinearAlgebra.__normalize!(y, nrm)
    function normalize_pullback(ȳ)
        Δy = unthunk(ȳ)
        ∂x = (Δy .- real(dot(y, Δy)) .* y) .* pinv(nrm)
        return (NoTangent(), ∂x)
    end
    normalize_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent())
    return y, normalize_pullback
end