text:
#=
  Cookie problem.

  From Think Bayes, page 3
  """
  Suppose there are two bowls of cookies.
  Bowl 1 contains 30 vanilla cookies and 10 chocolate cookies.
  Bowl 2 contains 20 of each.

  Now suppose you choose one of the bowls at random and, without looking,
  select a cookie at random. The cookie is vanilla.
  What is the probability that it came from Bowl 1?
  """

  Distributions of variable bowl (num:0)
  1.00000 =>    5991  (0.599100)
  2.00000 =>    4009  (0.400900)

  Cf  ~/webppl/cookie_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function cookie_problem()
    vanilla = 1
    chocolate = 2

    bowl1 = 1
    bowl2 = 2

    bowl ~ DiscreteUniform(bowl1,bowl2)
    cookie ~ (bowl == bowl1) ? Categorical(simplex([30,10])) : # [vanilla,chocolate]
                               Categorical(simplex([20,20]))

    true ~ Dirac(cookie == vanilla)

end

model = cookie_problem()
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :bowl)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using DeepForest, Flux, ProgressMeter, Plots, StatsPlots
using StatsBase: mode

x = rand(Float32, 2, 1000);
x[1, :] .*= 2*π;
x[1, :] .-= π;
x[2, :] .*= 3;
x[2, :] .-= 1.5;
y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

model = Forest(25, 2, 10, 2, 1.0)

forest_train!(1000, model, x, y)

ŷ = predict(model, x)
scatter(x[1, :], x[2, :], color=ŷ, legend=false)
plot!((x) -> sin(x), xlims=(-π, π), color=:green, linewidth=5)

imp = importance(model, x)
plot_importance(imp, x)

ext: jl


text:
using ForwardDiff, NeighbourLists

# slaterkoster.jl
#
# Collect all generic stuff for Slater-Koster-type Tight-binding
# models (which is 99.99% of non-selfconsistent TB models)
#

"""
`indexblock`:
a little auxiliary function to compute indices of Slater Koster orbitals,
this is returned as an SVector, i.e. it is generated on the stack so that no
heap memory is allocated.
"""
indexblock{IO,NORB}(n::Integer, H::SKHamiltonian{IO,NORB}) =
   SVector{NORB, Int}( ((n-1)*NORB+1):(n*NORB) )

function indexblock{T <: Integer}(Iat::AbstractVector{T}, H::SKHamiltonian)
   out = T[]
   for n in Iat
      append!(out, indexblock(n, H))
   end
   return out
end

############################################################
### Hamiltonian entries

######## s-orbital model

function sk!{IO}(H::SKHamiltonian{IO, 1}, U, bonds, out)
   setindex!(out, bonds[1], 1)
   return out
end

function sk_d!{IO}(H::SKHamiltonian{IO, 1}, r, R, b, db, dH_nm)
   # dH_nm is 3 x 1 x 1 so we can just index it linearly    (NORB = 1)
   for a = 1:3
      dH_nm[a] = db[1] * R[a] / r
   end
   return dH_nm
end

######## sp-orbital model

sk!{IO}(H::SKHamiltonian{IO, 4}, U, bonds, out) = _sk4!(U, bonds, out)
sk_d!{IO}(H::SKHamiltonian{IO, 4}, r, R, b, db, dout) = _sk4_d!(R/r, r, b, db, dout)

######## spd-orbital model

sk!{IO}(H::SKHamiltonian{IO, 9}, U, bonds, out) = _sk9!(U, bonds, out)
sk_d!{IO}(H::SKHamiltonian{IO, 9}, r, R, b, db, dout) = _sk9_d!(R/r, r, b, db, dout)

# ======================================================================
# prototypes for the functions needed to assemble the hamiltonian
# and its derivatives

# TODO: documentation, what is this called?

"""
`hop(H::SKHamiltonian, r, i)`: where `r` is real, `i` integer, this
should return the
"""
@protofun hop(::SKHamiltonian, ::Any, ::Any)

# TODO: this is probably type-unstable!
hop_d(H::SKHamiltonian, r, i) = ForwardDiff.derivative(s -> hop(H,s,i), r)

function hop!(H::SKHamiltonian, r, bonds)
   for i = 1:nbonds(H)
      bonds[i] = hop(H, r, i)
   end
   return bonds
end

# TODO: hop_d should return hop and hop_d, since we need both!!!!
function hop_d!(H::SKHamiltonian, r, b, db)
   for i = 1:nbonds(H)
      b[i] = hop(H, r, i)
      db[i] = hop_d(H, r, i)
   end
   return b, db
end

@protofun overlap(::SKHamiltonian, ::Real, ::Integer)

overlap_d(H::SKHamiltonian, r::Real, i) = ForwardDiff.derivative(s->overlap(H,s,i), r)

function overlap!(H::SKHamiltonian, r, bonds)
   for i = 1:nbonds(H)
      bonds[i] = overlap(H, r, i)
   end
   return bonds
end

function overlap_d!(H::SKHamiltonian, r, b, db)
   for i = 1:nbonds(H)
      b[i] = overlap(H, r, i)
      db[i] = overlap_d(H, r, i)
   end
   return b, db
end

# rewrite onsite! as diagonals???
# but we don't provide AD for these since they depend on too many variables,
# so AD will necessarily be inefficient.

function onsite! end
function onsite_grad! end

# =========================  Intermediate Hamiltonian Type

const SKBlock{NORB} = SMatrix{NORB, NORB, Float64}

"""
a triplet sparse matrix kind of thing that stores pre-computed
hamiltonian blocks, and can efficiently generate k-dependent Hamiltonians,
either sparse of full.

### Methods:
* `full(A::SparseSKH, k)`: returns full `H`, `M` for given `k`-vector.
* `full!(out, A, k)`: same as `full`, but in-place

### Methods to be implemented:
* `sparse(A::SparseSKH, k)`: same but sparse
* `collect(A::SparseSKH, k)`: chooses full (small systems) or sparse (large
                           systems) based on some simple heuristic

### Fields:
* `H` : the hamiltonian used to construct it
* `at` : the associated atoms object
* `i, j` : row and column indices
* `first` : `first[n]` is the index in `i, j` for which `i[idx] == n`
* `vH` : Hamiltonian blocks
* `vM` : overlap blocks
* `Rcell` : each hamiltonian block is associated with an e^{i k ⋅ S} multiplier
            from the Bloch transform; this S is stored in Rcell.
"""
immutable SparseSKH{HT, TV}  # v0.6: require that TV <: SKBlock{NORB}
   H::HT
   at::AbstractAtoms
   i::Vector{Int32}
   j::Vector{Int32}
   first::Vector{Int32}
   vH::Vector{TV}
   vM::Vector{TV}
   Rcell::Vector{JVecF}
end

Base.length(skh::SparseSKH) = length(skh.i)

function SparseSKH{ISORTH, NORB}(H::SKHamiltonian{ISORTH, NORB}, at::AbstractAtoms)
   if has_data(at, :SKH)
      return get_data(at, :SKH)
   end

   # here is a little hack that will turn an abstract type into a concrete type
   SKB = typeof(zero(SKBlock{NORB}))

   # get back positions since we still can't do without to get the cell-shifts
   # TODO: >>> maybe implement an alternative in JuLIP?
   X = positions(at)

   # ALLOCATIONS: TODO: move to a different method?
   # neighbourlist
   nlist = neighbourlist(at, cutoff(H))
   #      off-diagonal + diagonal
   nnz = npairs(nlist) + length(at)
   # allocate space for ordered triplet format
   i = zeros(Int32, nnz)
   j = zeros(Int32, nnz)
   first = zeros(Int32, length(at))
   vH = zeros(SKB, nnz)
   vM = ISORTH ? Vector{SKB}() : zeros(SKB, nnz)
   Rcell = zeros(JVecF, nnz)
   # index into the triplet format
   idx = 0

   # allocate space to assemble the hamiltonian blocks, we use MMatrix
   # here, but is this really necessary? Matrix should do fine?
   H_nm = zero(MMatrix{NORB, NORB, Float64})
   M_nm = zero(MMatrix{NORB, NORB, Float64})
   bonds = zeros(nbonds(H))     # temporary array for storing the potentials

   # loop through sites
   for (n, neigs, r, R) in sites(nlist)
      first[n] = idx+1          # where do the triplet entries for atom n start?

      # add the diagonal/on-site entries
      onsite!(H, r, R, H_nm)
      idx += 1
      i[idx], j[idx], vH[idx] = n, n, SKB(H_nm)
      # on-site overlap matrix block (only if the model is NONORTHOGONAL)
      if !ISORTH
         overlap!(H, M_nm)
         vM[idx] = SKB(M_nm)
      end

      # compute the Rvector corresponding to this block
      Rcell[idx] = zero(JVecF)

      # loop through the neighbours of the current atom (i.e. the bonds)
      for m = 1:length(neigs)
         U = R[m]/r[m]
         # hamiltonian block
         sk!(H, U, hop!(H, r[m], bonds), H_nm)
         idx += 1
         i[idx], j[idx], vH[idx] = n, neigs[m], SKB(H_nm)
         # overlap block (only if the model is NONORTHOGONAL)
         if !ISORTH
            sk!(H, U, overlap!(H, r[m], bonds), M_nm)
            vM[idx] = SKB(M_nm)
         end
         # compute the Rcell vector for these blocks
         Rcell[idx] = R[m] - (X[neigs[m]] - X[n])
      end
   end

   skh = SparseSKH(H, at, i, j, first, vH, vM, Rcell)
   set_transient!(at, :SKH, skh)
   return skh
end

_alloc_full(skh::SparseSKH) = _alloc_full(skh.H, skh.at)

_alloc_full(H::SKHamiltonian{NONORTHOGONAL}, at::AbstractAtoms) =
   Matrix{Complex128}(ndofs(H, at), ndofs(H, at)), Matrix{Complex128}(ndofs(H, at), ndofs(H, at))

_alloc_full(H::SKHamiltonian{ORTHOGONAL}, at::AbstractAtoms) =
   Matrix{Complex128}(ndofs(H, at), ndofs(H, at)), Matrix{Complex128}(0, 0)

Base.full(H::SparseSKH, k::AbstractVector = zero(JVecF)) =
   full!(_alloc_full(H), H, k)

full!(out, H::SparseSKH, k::AbstractVector = zero(JVecF)) =
   _full!(out[1], out[2], H, k, H.H)

function _full!{NORB}(Hout, Mout, skh, k, H::SKHamiltonian{NONORTHOGONAL, NORB})
   fill!(Hout, 0.0)
   fill!(Mout, 0.0)
   k = JVecF(k)
   for (i, j, H_ij, M_ij, S) in zip(skh.i, skh.j, skh.vH, skh.vM, skh.Rcell)
      eikR = exp( im * dot(k, S) )
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      @inbounds for a = 1:NORB, b = 1:NORB
         Hout[Ii[a], Ij[b]] += H_ij[a, b] * eikR
         Mout[Ii[a], Ij[b]] += M_ij[a, b] * eikR
      end
   end
   Hout[diagind(Hout)] = real(Hout[diagind(Hout)])
   Mout[diagind(Mout)] = real(Mout[diagind(Mout)])
   return Hout, Mout
end

function _full!{NORB}(Hout, _Mout_, skh, k, H::SKHamiltonian{ORTHOGONAL, NORB})
   fill!(Hout, 0.0)
   k = JVecF(k)
   for (i, j, H_ij, S) in zip(skh.i, skh.j, skh.vH, skh.Rcell)
      if i == 0 || j == 0
         error("unexplained i or j = 0")
      end
      eikR = exp( im * dot(k, S) )
      Ii, Ij = indexblock(i, H), indexblock(j, H)
      @inbounds for a = 1:NORB, b = 1:NORB
         Hout[Ii[a], Ij[b]] += H_ij[a, b] * eikR
      end
   end
   Hout[diagind(Hout)] = real(Hout[diagind(Hout)])
   return Hout, I
end

function Base.nnz(skh::SparseSKH)
   norb = norbitals(skh.H)
   return length(skh) * norb^2 + length(skh.at) * norb^2
end

# append to triplet format: version 1 for H and M (non-orth TB)
function _append!{NORB}(H::SKHamiltonian{NONORTHOGONAL, NORB},
                  It, Jt, Ht, Mt, In, Im, H_nm, M_nm, exp_i_kR, idx)
   @inbounds for i = 1:NORB, j = 1:NORB
      idx += 1
      It[idx] = In[i]
      Jt[idx] = Im[j]
      Ht[idx] = H_nm[i,j] * exp_i_kR
      Mt[idx] = M_nm[i,j] * exp_i_kR
   end
   return idx
end

# append to triplet format: version 2 for H only (orthogonal TB)
function _append!{NORB}(H::SKHamiltonian{ORTHOGONAL, NORB},
                  It, Jt, Ht, _Mt_, In, Im, H_nm, _Mnm_, exp_i_kR, idx)
   @inbounds for i = 1:NORB, j = 1:NORB
      idx += 1
      It[idx] = In[i]
      Jt[idx] = Im[j]
      Ht[idx] = H_nm[i,j] * exp_i_kR
   end
   return idx
end

function alloc_sparse(skh::SparseSKH)
   nnzest = nnz(skh)
   It = zeros(Int32, nnzest)
   Jt = zeros(Int32, nnzest)
   Ht = zeros(nnzest)
   Mt = zeros(nnzest)
   return It, Jt, Ht, Mt
end

#  * exp_i_kR = complex multiplier needed for BZ integration
#  * we could use cell * S instead of R[m] - (X[neigs[m]] - X[n])
#       but this would actually be less efficient, and less clear to read
function _assemble!{ISORTH, NORB}(H::SKHamiltonian{ISORTH, NORB},
                                  k, It, Jt, Ht, Mt, skh)
   N = ndofs(H, skh.at)
   M_nm = @SMatrix zeros(NORB, NORB)  # empty M_nm for ORTHOGONAL case
   idx = 0  # initialise index into triplet format
   for t = 1:length(skh)
      n, m, H_nm, S = skh.i[t], skh.j[t], skh.vH[t], skh.Rcell[t]
      if !ISORTH; M_nm = skh.vM[t]; end
      In, Im = indexblock(n, skh.H), indexblock(m, skh.H)
      exp_i_kR = exp( im * dot(k, S) )
      idx = _append!(H, It, Jt, Ht, Mt, In, Im, H_nm, M_nm, exp_i_kR, idx)
   end
   It, Jt, Ht, Mt = It[1:idx], Jt[1:idx], Ht[1:idx], Mt[1:idx]

   # convert M, H into Sparse CCS and return
   return ISORTH ? (sparse(It, Jt, Ht, N, N), I) :
                   (sparse(It, Jt, Ht, N, N), sparse(It, Jt, Mt, N, N))
end

Base.sparse(skh::SparseSKH, k = JVecF(0.0,0.0,0.0)) =
   _assemble!(skh.H, k, alloc_sparse(skh)..., skh)

function best(skh, k)
   const FULLSPARSE_CROSSOVER = 0.4  # based on testing with sp NRLTB model for Si
   nnz_full = (length(skh.at) * norbitals(skh.H))^2
   nnz_sparse = nnz(skh)
   if nnz_sparse / nnz_full  <= FULLSPARSE_CROSSOVER
      return sparse(skh, k)
   else
      return full(skh, k)
   end
end

############################################################
### Hamiltonian Evaluation

evaluate(H::SKHamiltonian, at::AbstractAtoms, k::AbstractVector; T=best) =
      T(SparseSKH(H, at), k)

# ================== Intermediate Storage for dHij/dRa elements

"""
a sparse-matrix kind of thing that stores pre-computed hamiltonian
derivative blocks
"""
mutable struct SparseSKHgrad{HT, TV}
   H::HT
   at::AbstractAtoms
   i::Vector{Int32}
   j::Vector{Int32}
   first::Vector{Int32}
   dH::Vector{TV}    # derivatives of hopping blocks
   dOS::Vector{TV}   # derivatives of on-site blocks
   dM::Vector{TV}    # derivatives of off-diagonal overlap blocks
   Rcell::Vector{JVecF}
end

# TODO
# this is a work-around until we move to v0.6; at that point
# StaticArrays changes and we can do
#    const SKBlockGrad{NORB} = SArray{Tuple{3,NORB,NORB}, ...}
# instead. But there is still the problem with D and L, which we need to
# figure out as well!
#
SKBlockGradType{IO, NORB}(H::SKHamiltonian{IO, NORB}) = typeof(@SArray zeros(3, NORB, NORB))

function SparseSKHgrad{ISORTH, NORB}(H::SKHamiltonian{ISORTH, NORB}, at::AbstractAtoms)
   # if the array has previously been computed, just return it
   if has_data(at, :SKBg)
      return get_data(at, :SKBg)
   end
   # (if not then generate it from scratch >>> rest of this function)

   # entry-type
   SKBg = SKBlockGradType(H)

   # TODO: >>> maybe implement an alternative in JuLIP?
   X = positions(at)

   # ALLOCATIONS: TODO: move to a different method?
   # neighbourlist
   nlist = neighbourlist(at, cutoff(H))
   #      off-diagonal + diagonal
   nnz = npairs(nlist)
   # allocate space for ordered triplet format
   i = zeros(Int32, nnz)
   j = zeros(Int32, nnz)
   first = zeros(Int32, length(at))
   dH = zeros(SKBg, nnz)
   dOS = zeros(SKBg, nnz)
   dM = ISORTH ? Vector{SKBg}() : zeros(SKBg, nnz)
   Rcell = zeros(JVecF, nnz)
   # index into the triplet format
   idx = 0

   # allocate work arrays
   maxneigs = maximum(length(jj) for (_1, jj, _2, _3) in sites(nlist))
   const dH_nn = zeros(3, NORB, NORB, maxneigs)
   const dH_nm = zeros(3, NORB, NORB)
   const dM_nm = zeros(3, NORB, NORB)
   const bonds = zeros(nbonds(H))
   const dbonds = zeros(nbonds(H))

   for (n, neigs, r, R) in sites(nlist)
      first[n] = idx+1

      # onsite derivative; copied into the new format in the following loop
      onsite_grad!(H, r, R, dH_nn)

      # loop through neighbours
      for i_n = 1:length(neigs)
         m = neigs[i_n]

         # compute ∂H_nm/∂y_n (hopping terms) and ∂M_nm/∂y_n
         hop_d!(H, r[i_n], bonds, dbonds)
         sk_d!(H, r[i_n], R[i_n], bonds, dbonds, dH_nm)
         if !ISORTH
            overlap_d!(H, r[i_n], bonds, dbonds)
            sk_d!(H, r[i_n], R[i_n], bonds, dbonds, dM_nm)
         end

         # store in new format
         idx += 1
         i[idx], j[idx], dH[idx] = n, m, SKBg(dH_nm)
         if !ISORTH; dM[idx] = SKBg(dM_nm); end
         dOS[idx] = SKBg(view(dH_nn, :, :, :, i_n))
         Rcell[idx] = R[i_n] - (X[m] - X[n])
      end
   end

   skhg = SparseSKHgrad(H, at, i, j, first, dH, dOS, dM, Rcell)
   set_transient!(at, :SKBg, skhg)
   return skhg
end


text:
# A collection of models for which the mean-of-means for the posterior should
# be same.
@model function gdemo1(x=10 * ones(2), ::Type{TV}=Vector{Float64}) where {TV}
    # `dot_assume` and `observe`
    m = TV(undef, length(x))
    m .~ Normal()
    return x ~ MvNormal(m, 0.5)
end

@model function gdemo2(x=10 * ones(2), ::Type{TV}=Vector{Float64}) where {TV}
    # `assume` with indexing and `observe`
    m = TV(undef, length(x))
    for i in eachindex(m)
        m[i] ~ Normal()
    end
    return x ~ MvNormal(m, 0.5)
end

@model function gdemo3(x=10 * ones(2))
    # Multivariate `assume` and `observe`
    m ~ MvNormal(length(x), 1.0)
    return x ~ MvNormal(m, 0.5)
end

@model function gdemo4(x=10 * ones(2), ::Type{TV}=Vector{Float64}) where {TV}
    # `dot_assume` and `observe` with indexing
    m = TV(undef, length(x))
    m .~ Normal()
    for i in eachindex(x)
        x[i] ~ Normal(m[i], 0.5)
    end
end

# Using vector of `length` 1 here so the posterior of `m` is the same
# as the others.
@model function gdemo5(x=10 * ones(1))
    # `assume` and `dot_observe`
    m ~ Normal()
    return x .~ Normal(m, 0.5)
end

@model function gdemo6(::Type{TV}=Vector{Float64}) where {TV}
    # `assume` and literal `observe`
    m ~ MvNormal(2, 1.0)
    return [10.0, 10.0] ~ MvNormal(m, 0.5)
end

@model function gdemo7(::Type{TV}=Vector{Float64}) where {TV}
    # `dot_assume` and literal `observe` with indexing
    m = TV(undef, 2)
    m .~ Normal()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end
end

@model function gdemo8(::Type{TV}=Vector{Float64}) where {TV}
    # `assume` and literal `dot_observe`
    m ~ Normal()
    return [10.0] .~ Normal(m, 0.5)
end

@model function _prior_dot_assume(::Type{TV}=Vector{Float64}) where {TV}
    m = TV(undef, 2)
    m .~ Normal()

    return m
end

@model function gdemo9()
    # Submodel prior
    m = @submodel _prior_dot_assume()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end
end

@model function _likelihood_dot_observe(m, x)
    return x ~ MvNormal(m, 0.5)
end

@model function gdemo10(x=10 * ones(2), ::Type{TV}=Vector{Float64}) where {TV}
    m = TV(undef, length(x))
    m .~ Normal()

    # Submodel likelihood
    @submodel _likelihood_dot_observe(m, x)
end

@model function gdemo11(x=10 * ones(2, 1), ::Type{TV}=Vector{Float64}) where {TV}
    m = TV(undef, length(x))
    m .~ Normal()

    # Dotted observe for `Matrix`.
    return x .~ MvNormal(m, 0.5)
end

const gdemo_models = (
    gdemo1(),
    gdemo2(),
    gdemo3(),
    gdemo4(),
    gdemo5(),
    gdemo6(),
    gdemo7(),
    gdemo8(),
    gdemo9(),
    gdemo10(),
    gdemo11(),
)

@testset "loglikelihoods.jl" begin
    for m in gdemo_models
        vi = VarInfo(m)

        vns = vi.metadata.m.vns
        if length(vns) == 1 && length(vi[vns[1]]) == 1
            # Only have one latent variable.
            DynamicPPL.setval!(vi, [1.0], ["m"])
        else
            DynamicPPL.setval!(vi, [1.0, 1.0], ["m[1]", "m[2]"])
        end

        lls = pointwise_loglikelihoods(m, vi)

        if isempty(lls)
            # One of the models with literal observations, so we just skip.
            continue
        end

        loglikelihood = if length(keys(lls)) == 1 && length(m.args.x) == 1
            # Only have one observation, so we need to double it
            # for comparison with other models.
            2 * sum(lls[first(keys(lls))])
        else
            sum(sum, values(lls))
        end

        @test loglikelihood ≈ -324.45158270528947
    end
end


text:
using Documenter
using State

makedocs(
    sitename = "State.jl",
    format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [State],
    pages = [
             "Home" => "index.md",
             "API" => "api.md"
            ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/Goysa2/State.jl")#


text:
_qr(::AbstractBandedLayout, ax, A) = _banded_qr(ax, A)
_banded_qr(_, A) = qr!(BandedMatrix{float(eltype(A))}(A, (bandwidth(A,1),bandwidth(A,1)+bandwidth(A,2))))

qr(A::Tridiagonal{T}) where T = qr!(BandedMatrix{float(T)}(A, (1,2)))
_qr!(::AbstractBandedLayout, _, A) = banded_qr!(A)

function _banded_qr!(R::AbstractMatrix{T}, τ) where T
    m,n=size(R)
    _banded_qr!(R, τ, min(m - 1 + !(T<:Real), n))
end

function _banded_qr!(R::AbstractMatrix, τ, ncols)
    D = bandeddata(R)
    l,u = bandwidths(R)
    ν = l+u+1
    m,n=size(R)

    for k = 1:ncols
        x = view(D,u+1:min(ν,m-k+u+1), k)
        τk = reflector!(x)
        τ[k] = τk
        for j = 1:min(u,n-k)
            reflectorApply!(x, τk, view(D, u+1-j:min(ν-j,m-k-j+u+1), k+j:k+j))
        end
    end
    R, τ
end

function banded_qr!(R::AbstractMatrix{T}, τ) where T
   _banded_qr!(R, τ)
   QR(R, convert(Vector{T},τ)) # TODO: remove convert
end

banded_qr!(R::AbstractMatrix{T}) where T = banded_qr!(R, zeros(T, min(size(R)...)))

function banded_qr_lmul!(A, B)
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = bandeddata(Afactors)
    @inbounds begin
        for k = min(mA,nA):-1:1
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:min(k+l,mB)
                    vBj += conj(D[i-k+u+1,k])*B[i,j]
                end
                vBj = A.τ[k]*vBj
                B[k,j] -= vBj
                for i = k+1:min(k+l,mB)
                    B[i,j] -= D[i-k+u+1,k]*vBj
                end
            end
        end
    end
    B
end

function banded_qr_lmul!(adjA::Adjoint, B)
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = bandeddata(Afactors)
    @inbounds begin
        for j = 1:nB
            cs = colsupport(B,j)
            for k = max(1,minimum(cs)):min(mA,nA,maximum(cs)+l,length(A.τ))
                vBj = B[k,j]
                for i = k+1:min(k+l,mB)
                    vBj += conj(D[i-k+u+1,k])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = k+1:min(k+l,mB)
                    B[i,j] -= D[i-k+u+1,k]*vBj
                end
            end
        end
    end
    B
end

function banded_qr_rmul!(A, Q)
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    l,u = bandwidths(Qfactors)
    D = bandeddata(Qfactors)
    @inbounds begin
        for k = 1:min(mQ,nQ)
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:min(k+l,mQ)
                    vAi += A[i,j]*D[j-k+u+1,k]
                end
                vAi = vAi*Q.τ[k]
                A[i,k] -= vAi
                for j = k+1:min(k+l,nA)
                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
                end
            end
        end
    end
    A
end
function banded_qr_rmul!(A, adjQ::Adjoint)
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    l,u = bandwidths(Qfactors)
    D = bandeddata(Qfactors)
    @inbounds begin
        for k = min(mQ,nQ):-1:1
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:min(k+l,mQ)
                    vAi += A[i,j]*D[j-k+u+1,k]
                end
                vAi = vAi*conj(Q.τ[k])
                A[i,k] -= vAi
                for j = k+1:min(k+l,nA)
                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
                end
            end
        end
    end
    A
end

banded_lmul!(A::QRPackedQ, B::AbstractVecOrMat) = banded_qr_lmul!(A, B)
banded_lmul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) = banded_qr_lmul!(adjA, B)
banded_rmul!(A::AbstractMatrix, Q::QRPackedQ) = banded_qr_rmul!(A, Q)
banded_rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ}) = banded_qr_rmul!(A, adjQ)

lmul!(A::QRPackedQ{<:Any,<:AbstractBandedMatrix}, B::AbstractVecOrMat) = banded_lmul!(A,B)
lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}, B::AbstractVecOrMat) = banded_lmul!(adjA,B)
lmul!(A::QRPackedQ{<:Any,BandedSubBandedMatrix{T,C,R,I1,I2}}, B::AbstractVecOrMat) where {T,C,R,I1<:AbstractUnitRange,I2<:AbstractUnitRange} =
    banded_lmul!(A,B)
lmul!(adjA::Adjoint{T,<:QRPackedQ{T,<:BandedSubBandedMatrix{T,C,R,I1,I2,t}}}, B::AbstractVecOrMat) where {T,C,R,I1<:AbstractUnitRange,I2<:AbstractUnitRange,t} =
    banded_lmul!(adjA,B)
# rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}) = banded_rmul!(A, adjA)
# rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}) = banded_rmul!(A, adjA)
rmul!(A::StridedVecOrMat{T}, Q::QRPackedQ{T,B}) where {T<:BlasFloat,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, Q)
rmul!(A::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,QRPackedQ{T,B}}) where {T<:BlasComplex,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, adjQ)
rmul!(A::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,QRPackedQ{T,B}}) where {T<:BlasReal,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, adjQ)

function _banded_widerect_ldiv!(A::QR{T}, B) where T
    error("Not implemented")
end
function _banded_longrect_ldiv!(A::QR, B)
    m, n = size(A)
    R = A.factors
    lmul!(adjoint(A.Q), B)
    B̃ = view(B, 1:n, :)
    B̃ .= Ldiv(UpperTriangular(view(R, 1:n, 1:n)), B̃)
    B
end
function _banded_square_ldiv!(A::QR, B)
    R = A.factors
    lmul!(adjoint(A.Q), B)
    B .= Ldiv(UpperTriangular(R), B)
    B
end

const UpperLayoutMatrix{T} = UpperTriangular{T,<:LayoutMatrix{T}}

for  BTyp in (:AbstractBandedMatrix, :BandedSubBandedMatrix), Typ in (:StridedVector, :StridedMatrix, :AbstractVecOrMat, :UpperLayoutMatrix )
    @eval function ldiv!(A::QR{T,<:$BTyp}, B::$Typ{T}) where T
        m, n = size(A)
        if m == n
            _banded_square_ldiv!(A, B)
        elseif n > m
            _banded_widerect_ldiv!(A, B)
        else
            _banded_longrect_ldiv!(A, B)
        end
    end
end


text:
const COLORS = [:blue, :red, :green, :yellow, :cyan, :magenta]

function paddedprint(io, s...; color=:default, pad=[])
    @nospecialize
    for (c, p) in pad
        printstyled(io, p, color=c)
    end
    printstyled(io, s..., color=color)
end

_printkeys(ch) = ["" for _ in eachindex(ch)]
_printkeys(ch::Union{NamedTuple, Dict, OrderedDict}) = ["$k: " for k in keys(ch)]
_printkeys(ch::PairVec) = ["$k: " for (k,v) in ch]

function _printtree(io::IO, n, C, d, p, e, trav, htrunc, vtrunc, comments)
    @nospecialize
    c = isa(NodeType(n), LeafNode) ? :default : C[1+d%length(C)]
    nr = sprint(nodeshow, n, context=io)
    paddedprint(io, nr * (trav ? ' ' * "[\"$(stringify(e))\"]" : ""), color=c)
    if comments
        nc = sprint(nodecommshow, n, context=io)
        paddedprint(io, isempty(nc) ? "" : " ", nc, color=:light_black)
    end

    CH = printchildren(n)
    PK = _printkeys(CH)
    CH = _iter(CH)
    nch = length(CH)

    gap = " " ^ clamp(length(nr)-1, 0, 2)

    function _printchild(i, ch, pk, l)
        println(io)
        line = (i == nch ? '└' : '├') * '─' ^ (2 + l-length(pk))
        paddedprint(io, gap * line * ' ' * pk, color=c, pad=p)
        ns = gap * (i == nch ? ' ' : '│') * ' ' ^ (3+l)
        _printtree(io, ch, C, d+1, [p; (c, ns)], e * encode(i, nch), trav, htrunc, vtrunc, comments)
    end

    if nch > 0
        if d+1 >= htrunc
            println(io)
            paddedprint(io, gap * '⋮', color=c, pad=p)
        elseif nch > vtrunc
            l2 = vtrunc < 2 ? 0 : (vtrunc < 4 ? 1 : 2)
            l1 = vtrunc - l2
            int1, int2 = 1:l1, nch-l2+1:nch
            l = reduce(max, vcat(length.(PK)[int1], length.(PK)[int2]), init=0)
            coll = enumerate(zip(CH, PK)) |> collect

            for (i, (ch, pk)) in coll[int1]
                _printchild(i, ch, pk, l)
            end

            println(io)
            paddedprint(io, gap * '⋮' , color=c, pad=p)

            for (i, (ch, pk)) in coll[int2]
                _printchild(i, ch, pk, l)
            end
        else
            l = maximum(length.(PK))
            for (i, (ch, pk)) in enumerate(zip(CH, PK))
                _printchild(i, ch, pk, l)
            end
        end
    end
end

printtree(n; kwargs...) = printtree(stdout, n; kwargs...)
function printtree(io::IO, n; trav::Bool=false, htrunc::Real=Inf, vtrunc::Real=Inf,
                        breakline::Bool=true, comments::Bool=true)
    @nospecialize
    @assert htrunc ≥ 0 "htrunc must be ≥ 0"
    @assert vtrunc ≥ 0 "vtrunc must be ≥ 0"
    _printtree(io, n, COLORS, 0, [], "", trav, htrunc, vtrunc, comments)
    breakline && println(io)
end


text: include("chains.jl")


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

using BAT
using Test
using Statistics, StatsBase, Distributions
using HypothesisTests

@testset "Multimodal Tests" begin
    cm = MixtureModel([Distributions.TDist(1), Distributions.TDist(1)])
    prod = Distributions.Product([Distributions.TDist(1), Distributions.TDist(1)])

    #Instantiates multimodal Student T with Mixture Model and Product
    @test @inferred(BAT.MultimodalStudentT(cm, 0.1, 2, prod)) isa BAT.MultimodalStudentT
    #Instantiates multimodal Student T with μ, σ and n
    @test @inferred(BAT.MultimodalStudentT(μ=60, σ=1., ν = 3, n = 2)) isa BAT.MultimodalStudentT
    #Error if n = 1
    @test_throws ArgumentError  BAT.MultimodalStudentT(μ=60, σ=1., ν = 3, n = 1)

    #If n = 1, student t = cauchy
    mt = BAT.MultimodalStudentT(μ = 1., σ = 0.1, ν = 1, n=4)

    @test @inferred(Base.size(mt)) == (4,)
    @test @inferred(Base.length(mt)) == 4
    @test @inferred(Base.eltype(mt)) == Float64
    @test all(isnan, Statistics.var(mt))
    @test all(isnan, Statistics.mean(mt))

    params = @inferred(StatsBase.params(mt))

    #Parameters
    @test params[1] == [-1., 1., 0., 0.]
    @test params[2] == [0.1, 0.1, 0.1, 0.1]
    @test params[3] == 4

    #logpdf
    @test @inferred(Distributions._logpdf(mt, [-1., 0., 1., 2.])) == -11.283438151658

    #If μ = 0, d = 1 the Distribution should be Cauchy-like in every dimension
    mmc = BAT.MultimodalStudentT(μ = 0., σ = 0.1, ν = 1, n=2)
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[1,:], Cauchy(0., 0.1))#First dimension
    @test pvalue(ks_test) > 0.05
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[2,:], Cauchy(0., 0.1))#Second dimension
    @test pvalue(ks_test) > 0.05

    #If μ = 0, ν > 1 the Distribution should be Student t-like in every dimension
    mmc = BAT.MultimodalStudentT(μ = 0., σ = 0.1, ν = 3, n=4)
    tdist = LocationScale(0, 0.1, TDist(3))
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[1,:], tdist)
    @test pvalue(ks_test) > 0.05
    ks_test = HypothesisTests.ExactOneSampleKSTest(rand(mmc, 10^6)[2,:], tdist)
    @test pvalue(ks_test) > 0.05

    @test ~any(isnan, Statistics.var(mmc))
    @test ~any(isnan, Statistics.mean(mmc))

end


text:
# Write a program that outputs all possibilities to put + or - or nothing
# between the numbers 1,2,…,9 (in this order) such that the result is 100. For
# example 1 + 2 + 3 - 4 + 5 + 6 + 78 + 9 = 100.
using Combinatorics

num = [1,2,3,4,5,6,7,8,9]

function comb(x::Vector{Int64}, s::Int64)
    k = [s,s+1]
    t = parse(Int64, join(x[k]))
    return t, k
end

function  eq_100(x::Vector{Int64})
	co = collect(combinations(1:8))

	for j in co
		for i in 1:8
			c,n = comb(x,i)
			nu = insert!(filter(x -> x ∉ n  , x),n[1],c)
			su = sum(.-(nu[j])) + sum(.+(filter(x -> x ∉ j  , nu)))
			if su == 100
				println(  map(x -> string.("-",x),nu[j]), string.(filter(x -> x ∉ j  , nu)) )
			end
		end
	end

end

eq_100(num)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
module LMDB_Env
    using LMDB
    using Test

    const dbname = "testdb"

    # Open environemnt
    env = create()
    @test env.handle != C_NULL
    @test env[:Readers] == 126
    @test env[:KeySize] == 511
    @test env[:Flags] == 0

    # Manipulate flags
    @test !isflagset(env[:Flags], Cuint(LMDB.NOSYNC))
    set!(env, LMDB.NOSYNC)
    @test isflagset(env[:Flags], Cuint(LMDB.NOSYNC))
    unset!(env, LMDB.NOSYNC)
    @test !isflagset(env[:Flags], Cuint(LMDB.NOSYNC))

    # Parameters
    @test (env[:Readers] = 100) == 100
    @test (env[:MapSize] = 1000^2) == 1000^2
    @test (env[:DBs] = 10) == 10
    @test env[:Readers] == 100

    # open db
    mkdir(dbname)
    try
        ret = open(env, dbname)
        @test ret[1] == 0

        # Close environment
        close(env)
        @test !isopen(env)

        # do block
        create() do env
            set!(env, LMDB.NOSYNC)
            open(env, dbname)
            @test isopen(env)
        end
    finally
        rm(dbname, recursive=true)
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Cairo

const FONT_NAME = "Sans"
const FONT_SIZE = 16.0
const LINE_HEIGHT = 20.0
const IMAGE_WIDTH = 300.0
const IMAGE_HEIGHT = 600.0

# ==================================================================================================
"""
    struct Line{T}
        words::Vector{SubString{T}}
        slant::Int32
        weight::Int32
    end

Construct a "line" of text.  Handles italics an bold via the constructor.

    Line(text::T; italic::Bool = false, bold::Bool = false) where T <: AbstractString
"""
struct Line{T}
    words::Vector{SubString{T}}
    slant::Int32
    weight::Int32
    Line(text::T; italic::Bool = false, bold::Bool = false) where T <: AbstractString = new{T}(
        split(text),
        italic ? Cairo.FONT_SLANT_ITALIC : Cairo.FONT_SLANT_NORMAL,
        bold ? Cairo.FONT_WEIGHT_BOLD : Cairo.FONT_WEIGHT_NORMAL
    )
end

const DEFAULT_TEXT = [
    Line("The quick brown fox jumped over the lazy dog."),
    Line("Suzy sells sea shells by the sea shore.", italic=true),
    Line("The sixth sick sheik's sixth sheep's sick.", bold=true)
]

# ==================================================================================================
"""
    wrap_words(
        cr::CairoContext,
        margin::Float64,
        start::Float64,
        words
    )::Float64

Writes the list of words to the image wrapping if they will extend off the right side.  Returns the
next Y offset to start writing at to avoid overwriting the existing text.
"""
function wrap_words(
            cr::CairoContext,
            margin::Float64,
            start::Float64,
            words
        )::Float64
    local y = 0
    local w = 1

    # -----------------------------------------------------------------------------------------
    # Process all the words.
    while w < length(words)
        local good = "$(words[w])"
        local line = ""

        # -----------------------------------------------------------------------------------------
        # Keep increasing the line length until we overflow.
        while w < length(words)
            w       = w + 1
            line    = "$good $(words[w])"
            local e = text_extents(cr, line)

            if e[3] > IMAGE_WIDTH - (2.0 * margin)
                w -= 1
                break
            end

            good = line
        end

        # -----------------------------------------------------------------------------------------
        # Write out the line.
        move_to(cr, margin, start + Float64(y) * LINE_HEIGHT)
        show_text(cr, good)
        y += 1
    end

    if start + Float64(y) * LINE_HEIGHT > IMAGE_HEIGHT
        @error "Ran off the bottom of the image."
    end

    return start + Float64(y) * LINE_HEIGHT
end

# ==================================================================================================
"""
    draw_image(text = DEFAULT_TEXT)::Vector{UInt8}

Creates an image with the specified text as a PNG.  The text should be a list of Line objects to
write to the image.
"""
function draw_image(text = DEFAULT_TEXT)::Vector{UInt8}
    local c = CairoImageSurface(
        IMAGE_WIDTH,
        IMAGE_HEIGHT,
        Cairo.FORMAT_RGB24
        )
    local cr = CairoContext(c)

    # Create a white image.
    save(cr)
        set_source_rgb(cr, 1.0, 1.0, 1.0)
        rectangle(cr, 0.0, 0.0, IMAGE_WIDTH, IMAGE_HEIGHT)
        fill(cr)
    restore(cr)

    save(cr)
        set_source_rgb(cr, 0.0, 0.0, 0.0);
        set_font_size(cr, FONT_SIZE)

        local y = 0.0
        for line in text
            y += LINE_HEIGHT
            select_font_face(cr, "Sans", line.slant, line.weight)
            y = wrap_words(cr, 10.0, y, line.words)
        end
    restore(cr)

    local buffer = IOBuffer()
    write_to_png(c, buffer)
    return take!(buffer)
end

# =================================================================================================
"""
    pix_with(text = DEFAULT_TEXT)::Pix

Load a Pix image into memory.
"""
function pix_with(text = DEFAULT_TEXT)::Pix
    local buffer = draw_image(text)
    return pix_read_png(buffer)
end

# =================================================================================================
"""
    safe_tmp_file()::String

Create a temporary file we can write to.
"""
function safe_tmp_file()::String
    local path, io = mktemp()
    close(io)
    return path
end

# =================================================================================================
"""
    save_image(data)::String

Save an image to a temporary file.  Returns the path of the temporary file.
"""
function save_image(data)::String
    local path, io = mktemp()
    write(io, data)
    close(io)
    return path
end

# =================================================================================================
"""
    bmp_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a BMP image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function bmp_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_bmp(pix)
end

# =================================================================================================
"""
    gif_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a GIF image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function gif_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_gif(pix)
end

# =================================================================================================
"""
    jp2k_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a J2K image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function jp2k_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_jp2k(pix)
end

# =================================================================================================
"""
    jpeg_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a JPG image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function jpeg_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_jpeg(pix)
end

# =================================================================================================
"""
    png_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a PNG image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function png_with(text = DEFAULT_TEXT)::Vector{UInt8}
    return draw_image(text)
end

# =================================================================================================
"""
    pnm_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a PNM image with the specified text in memory.  Returns the image data. The text should be a
list of Line objects.
"""
function pnm_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_pnm(pix)
end

# =================================================================================================
"""
    spix_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a SPIX image with the specified text in memory.  Returns the image data. The text should be
a list of Line objects.
"""
function spix_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_spix(pix)
end

# =================================================================================================
"""
    tiff_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a TIFF image with the specified text in memory.  Returns the image data. The text should be
a list of Line objects.
"""
function tiff_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_tiff(pix)
end

# =================================================================================================
"""
    webp_with(text = DEFAULT_TEXT)::Vector{UInt8}

Create a WEBP image with the specified text in memory.  Returns the image data. The text should be
a list of Line objects.
"""
function webp_with(text = DEFAULT_TEXT)::Vector{UInt8}
    local pix = pix_with(text)
    return pix_write_webp(pix)
end

# =================================================================================================
"""
    bmp_file_with(text = DEFAULT_TEXT)::String

Create a BMP file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function bmp_file_with(text = DEFAULT_TEXT)::String
    local data = bmp_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    gif_file_with(text = DEFAULT_TEXT)::String

Create a GIF file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function gif_file_with(text = DEFAULT_TEXT)::String
    local data = gif_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    jp2k_file_with(text = DEFAULT_TEXT)::String

Create a J2K file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function jp2k_file_with(text = DEFAULT_TEXT)::String
    local data = jp2k_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    jpeg_file_with(text = DEFAULT_TEXT)::String

Create a JPG file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function jpeg_file_with(text = DEFAULT_TEXT)::String
    local data = jpeg_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    png_file_with(text = DEFAULT_TEXT)::String

Create a PNG file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function png_file_with(text = DEFAULT_TEXT)::String
    local data = png_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    pnm_file_with(text = DEFAULT_TEXT)::String

Create a PNM file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function pnm_file_with(text = DEFAULT_TEXT)::String
    local data = pnm_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    spix_file_with(text = DEFAULT_TEXT)::String

Create a SPIX file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function spix_file_with(text = DEFAULT_TEXT)::String
    local data = spix_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    tiff_file_with(text = DEFAULT_TEXT)::String

Create a TIFF file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function tiff_file_with(text = DEFAULT_TEXT)::String
    local data = tiff_with(text)
    return save_image(data)
end

# =================================================================================================
"""
    webp_file_with(text = DEFAULT_TEXT)::String

Create a WEBP file with the specified text.  Returns the path to the file that was created.  The
text should be a list of Line objects.
"""
function webp_file_with(text = DEFAULT_TEXT)::String
    local data = webp_with(text)
    return save_image(data)
end


text:
include("isValidPair.jl");
function isValidSetCard(set, newCard)
    setLength = length(set);
    if (setLength < 1)
        return true;
    end
    if (in(newCard, set))
        return false;
    end
    for j = 1:setLength
        valid = isValidPair(set[j], newCard)
        if (!valid)
            return false;
        end
    end
    return true;
end

max_stars_repo_name: driesverhees/dobble

max_issues_repo_name: driesverhees/dobble

max_forks_repo_name: driesverhees/dobble


text:
#Link to ODBC Driver Manager (system-dependent)
let
    global odbc_dm
    local lib
    succeeded=false
    @linux_only lib_choices = ["libodbc", "libodbc.so", "libodbc.so.1", "libodbc.so.2", "libodbc.so.3"]
	@windows_only lib_choices = ["odbc32"]
	@osx_only lib_choices = ["libiodbc","libiodbc.dylib","libiodbc.1.dylib","libiodbc.2.dylib","libiodbc.3.dylib"]
    try
    	if isdefined(odbc_dm) insert!(lib_choices,1,strip(odbc_dm)) end
    end
    for lib in lib_choices
        try
            dlopen(lib)
            succeeded=true
            break
        end
    end
    if !succeeded error("ODBC library not found") end
    @eval const odbc_dm = $lib
end

#Translation of sqltypes.h; C typealiases for SQL functions
#http://msdn.microsoft.com/en-us/library/windows/desktop/ms716298(v=vs.85).aspx
#http://msdn.microsoft.com/en-us/library/windows/desktop/aa383751(v=vs.85).aspx
typealias SQLCHAR       Uint8
typealias SQLSCHAR      Uint8
typealias SQLVARCHAR    Uint8
typealias SQLDECIMAL    Cdouble
typealias SQLNUMERIC    Cdouble
typealias SQLDOUBLE     Cdouble
typealias SQLFLOAT      Cdouble
typealias SQLREAL       Cfloat
typealias SQLINTEGER    Cint
typealias SQLUINTEGER   Cuint
typealias SQLSMALLINT   Cshort
typealias SQLUSMALLINT  Cushort
typealias SQLTIME       Cuchar
typealias SQLTIMESTAMP  Cuchar

if contains(odbc_dm,"iodbc")
	typealias SQLWCHAR Uint32
else
	typealias SQLWCHAR Uint16
end
if WORD_SIZE == 64
	typealias SQLLEN        Int64
	typealias SQLULEN       Uint64
	typealias SQLSETPOSIROW Uint64
else
	typealias SQLLEN        SQLINTEGER
	typealias SQLULEN       SQLUINTEGER
	typealias SQLSETPOSIROW SQLUSMALLINT
end

typealias SQLROWCOUNT   SQLULEN
typealias SQLROWSETSIZE SQLULEN
typealias SQLTRANSID    SQLULEN
typealias SQLROWOFFSET  SQLLEN
typealias SQLPOINTER    Ptr{Void}
typealias SQLRETURN     SQLSMALLINT
typealias SQLHANDLE     Ptr{Void}
typealias SQLHENV       SQLHANDLE
typealias SQLHDBC       SQLHANDLE
typealias SQLHSTMT      SQLHANDLE
typealias SQLHDESC      SQLHANDLE
typealias ULONG         Cuint
typealias PULONG        Ptr{ULONG}
typealias USHORT        Cushort
typealias PUSHORT       Ptr{USHORT}
typealias UCHAR         Cuchar
typealias PUCHAR        Ptr{Cuchar}
typealias PSZ           Ptr{Cchar}
typealias SCHAR         Cchar
typealias SDWORD        Cint
typealias SWORD         Cshort
typealias UDWORD        Cuint
typealias UWORD         Cushort
typealias SLONG         Cint
typealias SSHORT        Cshort
typealias SDOUBLE       Cdouble
typealias LDOUBLE       Cdouble
typealias SFLOAT        Cfloat
typealias PTR           Ptr{Void}
typealias HENV          Ptr{Void}
typealias HDBC          Ptr{Void}
typealias HSTMT         Ptr{Void}
typealias RETCODE       Cshort
typealias SQLHWND       Ptr{Void}

#################

# provide lowercase conversion functions for all types
# e.g., sqlchar(x) = convert(SQLCHAR, x)
for t in [:SQLCHAR, :SQLSCHAR, :SQLWCHAR, :SQLDATE, :SQLDECIMAL,
		  :SQLDOUBLE, :SQLFLOAT, :SQLINTEGER, :SQLUINTEGER,
		  :SQLSMALLINT, :SQLUSMALLINT, :SQLLEN, :SQLULEN,
		  :SQLSETPOSIROW, :SQLROWCOUNT, :SQLROWSETSIZE, :SQLTRANSID,
		  :SQLROWOFFSET, :SQLNUMERIC, :SQLPOINTER, :SQLREAL, :SQLTIME,
		  :SQLTIMESTAMP, :SQLVARCHAR, :SQLRETURN, :SQLHANDLE,
		  :SQLHENV, :SQLHDBC, :SQLHSTMT, :SQLHDESC, :ULONG, :PULONG,
		  :USHORT, :PUSHORT, :UCHAR, :PUCHAR, :PSZ, :SCHAR, :SDWORD,
		  :SWORD, :UDWORD, :UWORD, :SLONG, :SSHORT, :SDOUBLE,
		  :LDOUBLE, :SFLOAT, :PTR, :HENV, :HDBC, :HSTMT, :RETCODE,
		  :SQLHWND]
	fn = symbol(lowercase(string(t)))
	@eval $fn(x) = convert($t, x)
end

# Data Type Mappings
# SQL data types are returned in resultset metadata calls (ODBCMetadata)
# C data types are used in SQLBindCols (ODBCFetch) to allocate column memory; the driver manager converts from the SQL type to this C type in memory
# Julia types indicate how julia should read the returned C data type memory from the previous step

# Data Type Status: Pretty good, I think we're at 95% support, really only missing native date, time, interval, and GUID types (currently just read in as strings)

# SQL Data Type 					C Data Type 						Julia Type
# ---------------------------------------------------------------------------------
# SQL_CHAR 							SQL_C_CHAR 							Uint8
# SQL_VARCHAR 						SQL_C_CHAR 							Uint8
# SQL_LONGVARCHAR 					SQL_C_CHAR 							Uint8
# SQL_WCHAR 						SQL_C_WCHAR 						Uint16
# SQL_WVARCHAR 						SQL_C_WCHAR 						Uint16
# SQL_WLONGVARCHAR 					SQL_C_WCHAR 						Uint16
# SQL_DECIMAL 						SQL_C_DOUBLE 						Float64
# SQL_NUMERIC 						SQL_C_DOUBLE 						Float64
# SQL_SMALLINT 						SQL_C_SHORT 						Int16
# SQL_INTEGER 						SQL_C_LONG 							Int32
# SQL_REAL 							SQL_C_FLOAT 						Float64
# SQL_FLOAT 						SQL_C_DOUBLE 						Float64
# SQL_DOUBLE 						SQL_C_DOUBLE 						Float64
# SQL_BIT 							SQL_C_BIT 							Int8
# SQL_TINYINT 						SQL_C_TINYINT 						Int8
# SQL_BIGINT 						SQL_C_BIGINT 						Int64
# SQL_BINARY 						SQL_C_BINARY 						Uint8
# SQL_VARBINARY 					SQL_C_BINARY 						Uint8
# SQL_LONGVARBINARY 				SQL_C_BINARY 						Uint8
# SQL_TYPE_DATE 					SQL_C_TYPE_DATE 					SQLDate
# SQL_TYPE_TIME 					SQL_C_TYPE_TIME 					SQLTime
# SQL_TYPE_TIMESTAMP 				SQL_C_TYPE_TIMESTAMP 				SQLTimestamp
# SQL_INTERVAL_MONTH				SQL_C_INTERVAL_MONTH 				Uint8
# SQL_INTERVAL_YEAR					SQL_C_INTERVAL_YEAR 				Uint8
# SQL_INTERVAL_YEAR_TO_MONTH		SQL_C_INTERVAL_YEAR_TO_MONTH 		Uint8
# SQL_INTERVAL_DAY					SQL_C_INTERVAL_DAY 					Uint8
# SQL_INTERVAL_HOUR					SQL_C_INTERVAL_HOUR 				Uint8
# SQL_INTERVAL_MINUTE				SQL_C_INTERVAL_MINUTE 				Uint8
# SQL_INTERVAL_SECOND				SQL_C_INTERVAL_SECOND 				Uint8
# SQL_INTERVAL_DAY_TO_HOUR			SQL_C_INTERVAL_DAY_TO_HOUR 			Uint8
# SQL_INTERVAL_DAY_TO_MINUTE		SQL_C_INTERVAL_DAY_TO_MINUTE 		Uint8
# SQL_INTERVAL_DAY_TO_SECOND		SQL_C_INTERVAL_DAY_TO_SECOND 		Uint8
# SQL_INTERVAL_HOUR_TO_MINUTE		SQL_C_INTERVAL_HOUR_TO_MINUTE 		Uint8
# SQL_INTERVAL_HOUR_TO_SECOND		SQL_C_INTERVAL_HOUR_TO_SECOND 		Uint8
# SQL_INTERVAL_MINUTE_TO_SECOND		SQL_C_INTERVAL_MINUTE_TO_SECOND 	Uint8
# SQL_GUID							SQL_C_GUID 							Uint8

#SQL Data Type Definitions
const SQL_CHAR          = int16(  1) # Character string of fixed string length n.
const SQL_VARCHAR       = int16( 12) # Variable-length character string with a maximum string length n.
const SQL_LONGVARCHAR   = int16( -1) # Variable length character data. Maximum length is data source–dependent.
const SQL_WCHAR         = int16( -8) # Unicode character string of fixed string length n
const SQL_WVARCHAR      = int16( -9) # Unicode variable-length character string with a maximum string length n
const SQL_WLONGVARCHAR  = int16(-10) # Unicode variable-length character data. Maximum length is data source–dependent
const SQL_DECIMAL       = int16(  3) #
const SQL_NUMERIC       = int16(  2)
const SQL_SMALLINT      = int16(  5) # Exact numeric value with precision 5 and scale 0 (signed: –32,768 <= n <= 32,767, unsigned: 0 <= n <= 65,535)
const SQL_INTEGER       = int16(  4) # Exact numeric value with precision 10 and scale 0 (signed: –2[31] <= n <= 2[31] – 1, unsigned: 0 <= n <= 2[32] – 1)
const SQL_REAL          = int16(  7) # Signed, approximate, numeric value with a binary precision 24 (zero or absolute value 10[–38] to 10[38]).
const SQL_FLOAT         = int16(  6) # Signed, approximate, numeric value with a binary precision of at least p. (The maximum precision is driver-defined.)
const SQL_DOUBLE        = int16(  8) # Signed, approximate, numeric value with a binary precision 53 (zero or absolute value 10[–308] to 10[308]).
const SQL_BIT           = int16( -7) # Single bit binary data.
const SQL_TINYINT       = int16( -6) # Exact numeric value with precision 3 and scale 0 (signed: –128 <= n <= 127, unsigned: 0 <= n <= 255)
const SQL_BIGINT        = int16( -5) # Exact numeric value with precision 19 (if signed) or 20 (if unsigned) and scale 0 (signed: –2[63] <= n <= 2[63] – 1, unsigned: 0 <= n <= 2[64] – 1)
const SQL_BINARY        = int16( -2) # Binary data of fixed length n.
const SQL_VARBINARY     = int16( -3) # Variable length binary data of maximum length n. The maximum is set by the user.
const SQL_LONGVARBINARY = int16( -4) # Variable length binary data. Maximum length is data source–dependent.
const SQL_TYPE_DATE     = int16( 91) # Year, month, and day fields, conforming to the rules of the Gregorian calendar.
const SQL_TYPE_TIMESTAMP            = int16( 93) # Year, month, day, hour, minute, and second fields, with valid values as defined for the DATE and TIME data types.
const SQL_TYPE_TIME                 = int16( 92) # Hour, minute, and second fields, with valid values for hours of 00 to 23, valid values for minutes of 00 to 59, and valid values for seconds of 00 to 61. Precision p indicates the seconds precision.
#const SQL_INTERVAL_MONTH            = int16(102)
#const SQL_INTERVAL_YEAR             = int16(101)
#const SQL_INTERVAL_YEAR_TO_MONTH    = int16(107)
#const SQL_INTERVAL_DAY              = int16(103)
#const SQL_INTERVAL_HOUR             = int16(104)
#const SQL_INTERVAL_MINUTE           = int16(105)
#const SQL_INTERVAL_SECOND           = int16(106)
#const SQL_INTERVAL_DAY_TO_HOUR      = int16(108)
#const SQL_INTERVAL_DAY_TO_MINUTE    = int16(109)
#const SQL_INTERVAL_DAY_TO_SECOND    = int16(110)
#const SQL_INTERVAL_HOUR_TO_MINUTE   = int16(111)
#const SQL_INTERVAL_HOUR_TO_SECOND   = int16(112)
#const SQL_INTERVAL_MINUTE_TO_SECOND = int16(113)
#const SQL_GUID                      = int16(-11) # Fixed length GUID.

#C Data Types
const SQL_C_CHAR      = int16(  1)
const SQL_C_WCHAR     = int16( -8)
const SQL_C_DOUBLE    = int16(  8)
const SQL_C_SHORT     = int16(  5)
const SQL_C_LONG      = int16(  4)
const SQL_C_FLOAT     = int16(  7)
const SQL_C_BIT       = int16( -7)
const SQL_C_TINYINT   = int16( -6)
const SQL_C_BIGINT    = int16(-27)
const SQL_C_BINARY    = int16( -2)
const SQL_C_TYPE_DATE = int16( 91)
const SQL_C_TYPE_TIMESTAMP            = int16( 93)
const SQL_C_TYPE_TIME                 = int16( 92)
#const SQL_C_INTERVAL_MONTH            = int16(102)
#const SQL_C_INTERVAL_YEAR             = int16(101)
#const SQL_C_INTERVAL_YEAR_TO_MONTH    = int16(107)
#const SQL_C_INTERVAL_DAY              = int16(103)
#const SQL_C_INTERVAL_HOUR             = int16(104)
#const SQL_C_INTERVAL_MINUTE           = int16(105)
#const SQL_C_INTERVAL_SECOND           = int16(106)
#const SQL_C_INTERVAL_DAY_TO_HOUR      = int16(108)
#const SQL_C_INTERVAL_DAY_TO_MINUTE    = int16(109)
#const SQL_C_INTERVAL_DAY_TO_SECOND    = int16(110)
#const SQL_C_INTERVAL_HOUR_TO_MINUTE   = int16(111)
#const SQL_C_INTERVAL_HOUR_TO_SECOND   = int16(112)
#const SQL_C_INTERVAL_MINUTE_TO_SECOND = int16(113)
#const SQL_C_GUID                      = int16(-11)

#Julia structs
immutable SQLDate
	year::Int16
	month::Int16
	day::Int16
end
string(x::SQLDate) = "$(x.year)-$(x.month)-$(x.day)"
immutable SQLTime
	hour::Int16
	minute::Int16
	second::Int16
end
show(io::IO,x::SQLTime) = print(io,"$(x.hour):$(x.minute):$(x.second)")
immutable SQLTimestamp
	year::Int16
	month::Int16
	day::Int16
	hour::Int16
	minute::Int16
	second::Int16
	fraction::Int32 #nanoseconds
end
string(x::SQLTimestamp) = "$(x.year)-$(x.month)-$(x.day) $(x.hour):$(x.minute):$(x.second)"

const SQL2C = [
	SQL_CHAR          => SQL_C_CHAR,
	SQL_VARCHAR       => SQL_C_CHAR,
	SQL_LONGVARCHAR   => SQL_C_CHAR,
	SQL_WCHAR         => SQL_C_WCHAR,
	SQL_WVARCHAR      => SQL_C_WCHAR,
	SQL_WLONGVARCHAR  => SQL_C_WCHAR,
	SQL_DECIMAL       => SQL_C_DOUBLE,
	SQL_NUMERIC       => SQL_C_DOUBLE,
	SQL_SMALLINT      => SQL_C_SHORT,
	SQL_INTEGER       => SQL_C_LONG,
	SQL_REAL          => SQL_C_FLOAT,
	SQL_FLOAT         => SQL_C_DOUBLE,
	SQL_DOUBLE        => SQL_C_DOUBLE,
	SQL_BIT           => SQL_C_BIT,
	SQL_TINYINT       => SQL_C_TINYINT,
	SQL_BIGINT        => SQL_C_BIGINT,
	SQL_BINARY        => SQL_C_BINARY,
	SQL_VARBINARY     => SQL_C_BINARY,
	SQL_LONGVARBINARY => SQL_C_BINARY,
	SQL_TYPE_DATE     => SQL_C_TYPE_DATE,
	SQL_TYPE_TIME 	  => SQL_C_TYPE_TIME,
	SQL_TYPE_TIMESTAMP=> SQL_C_TYPE_TIMESTAMP]

const SQL2Julia = [
	SQL_CHAR          => SQLCHAR,
	SQL_VARCHAR       => SQLVARCHAR,
	SQL_LONGVARCHAR   => SQLVARCHAR,
	SQL_WCHAR         => SQLWCHAR,
	SQL_WVARCHAR      => SQLWCHAR,
	SQL_WLONGVARCHAR  => SQLWCHAR,
	SQL_DECIMAL       => SQLDECIMAL,
	SQL_NUMERIC       => SQLNUMERIC,
	SQL_SMALLINT      => SQLSMALLINT,
	SQL_INTEGER       => SQLINTEGER,
	SQL_REAL          => SQLREAL,
	SQL_FLOAT         => SQLFLOAT,
	SQL_DOUBLE        => SQLDOUBLE,
	SQL_BIT           => Int8,
	SQL_TINYINT       => Int8,
	SQL_BIGINT        => Int64,
	SQL_BINARY        => Uint8,
	SQL_VARBINARY     => Uint8,
	SQL_LONGVARBINARY => Uint8,
	SQL_TYPE_DATE     => SQLDate,
	SQL_TYPE_TIME 	  => SQLTime,
	SQL_TYPE_TIMESTAMP=> SQLTimestamp]

const SQL_TYPES = [
	  1 => "SQL_CHAR",
	 12 => "SQL_VARCHAR",
	 -1 => "SQL_LONGVARCHAR",
	 -8 => "SQL_WCHAR",
	 -9 => "SQL_WVARCHAR",
	-10 => "SQL_WLONGVARCHAR",
	  3 => "SQL_DECIMAL",
	  2 => "SQL_NUMERIC",
	  5 => "SQL_SMALLINT",
	  4 => "SQL_INTEGER",
	  7 => "SQL_REAL",
	  6 => "SQL_FLOAT",
	  8 => "SQL_DOUBLE",
	 -7 => "SQL_BIT",
	 -6 => "SQL_TINYINT",
	 -5 => "SQL_BIGINT",
	 -2 => "SQL_BINARY",
	 -3 => "SQL_VARBINARY",
	 -4 => "SQL_LONGVARBINARY",
	 91 => "SQL_TYPE_DATE",
	 92 => "SQL_TYPE_TIME",
	 93 => "SQL_TYPE_TIMESTAMP",
	102 => "SQL_INTERVAL_MONTH",
	101 => "SQL_INTERVAL_YEAR",
	107 => "SQL_INTERVAL_YEAR_TO_MONTH",
	103 => "SQL_INTERVAL_DAY",
	104 => "SQL_INTERVAL_HOUR",
	105 => "SQL_INTERVAL_MINUTE",
	106 => "SQL_INTERVAL_SECOND",
	108 => "SQL_INTERVAL_DAY_TO_HOUR",
	109 => "SQL_INTERVAL_DAY_TO_MINUTE",
	110 => "SQL_INTERVAL_DAY_TO_SECOND",
	111 => "SQL_INTERVAL_HOUR_TO_MINUTE",
	112 => "SQL_INTERVAL_HOUR_TO_SECOND",
	113 => "SQL_INTERVAL_MINUTE_TO_SECOND",
	-11 => "SQL_GUID"]


text:
function slope(hs,errors)
  x = log10.(hs)
  y = log10.(errors)
  linreg = hcat(fill!(similar(x), 1), x) \ y
  linreg[2]
end

"""
Generates the values of n used for the convergence study within range
[1:min(10^k,n_max)] using a logarithmic distribution of the sample points.
"""
function generate_n_values(k;n_max=120)
  n_values=Int[]
  current=1
  for power=1:k
    for coefficient=1:10
       n=coefficient*current
       if (n>n_max)
        return n_values[2:end]
       end
       append!(n_values,n)
    end
    current=current*10
  end
  return n_values[2:end]
end

"""
   Performs a convergence study.

   Returns a tuple with three entries.
     [1] Values of h.
     [2] Relative errors.
     [3] Slope of the log(h)-log(err) relative error convergence curve.
"""
function convergence_study(f,hs,model_args_series,order,degree,fargs...)
   errors=Float64[]
   for args in model_args_series
     model=CubedSphereDiscreteModel(args...)
     err=f(model,order,degree,fargs...)
     append!(errors,err)
   end
   return hs,errors,slope(hs,errors)
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
function createstiffnessmatrices!(params::Parameters)
  noofnodes = length(params.nodes)
  sb0 = zeros(6, 6)
  for i in 1:noofnodes
    if 1 < i < noofnodes
      ler = params.elements[i-1].er
      lea = params.elements[i-1].ea
      lgj = params.elements[i-1].gj
      l1 = params.nodes[i].z - params.nodes[i-1].z
      l2 = params.nodes[i+1].z - params.nodes[i].z
      params.sa[:, :, i] = stiffnessa!(i, l1, ler, lea, lgj)
      params.sb1[:, :, i] = stiffnessb1!(i, l1, ler, lea, lgj)
      params.sb2[:, :, i] = stiffnessb2!(i, noofnodes, l2, ler, lea, lgj)
      params.sc[:, :, i] = stiffnessc!(i, noofnodes, l2, ler, lea, lgj)
      params.si[:, :, i] = stiffnessi!(i, params.sb1[:, :, i], params.sb2[:, :, i])
    else
      if i == 1
        ler = params.elements[i].er
        lea = params.elements[i].ea
        lgj = params.elements[i].gj
        l2 = params.nodes[i+1].z - params.nodes[i].z
        params.sb2[:, :, i] = stiffnessb2!(i, noofnodes, l2, ler, lea, lgj)
        params.sc[:, :, i] = stiffnessc!(i, noofnodes, l2, ler, lea, lgj)
        params.si[:, :, i] = stiffnessi!(i, sb0, params.sb2[:, :, i])
      else
        ler = params.elements[i-1].er
        lea = params.elements[i-1].ea
        lgj = params.elements[i-1].gj
        l1 = params.nodes[i].z - params.nodes[i-1].z
        params.sa[:, :, i] = stiffnessa!(i, l1, ler, lea, lgj)
        params.sb1[:, :, i] = stiffnessb1!(i, l1, ler, lea, lgj)
        params.si[:, :, i] = stiffnessi!(i, params.sb1[:, :, i], sb0)
      end
    end
  end

  params.sainit[:, :, :] = copy(params.sa)
  params.sb1init[:, :, :] = copy(params.sb1)
  params.sb2init[:, :, :] = copy(params.sb2)
  params.scinit[:, :, :] = copy(params.sc)

  nothing
end

function stiffnessa!(node, l, ler, lea, lgj)
  s = zeros(6, 6)
  if node == 1
    return
  else
    s[1, 1] = -12.0 * ler / l^3
    s[1, 5] =  -6.0 * ler / l^2
    s[2, 2] = -12.0 * ler / l^3
    s[2, 4] =   6.0 * ler / l^2
    s[3, 3] =        -lea / l
    s[4, 2] =  -6.0 * ler / l^2
    s[4, 4] =   2.0 * ler / l
    s[5, 1] =   6.0 * ler / l^2
    s[5, 5] =   2.0 * ler / l
    s[6, 6] =        -lgj / l
  end
  s
end

function stiffnessb1!(node, l, ler, lea, lgj)
  s = zeros(6, 6)
  if node == 1
    return
  else
    s[1, 1] =  12.0 * ler / l^3
    s[1, 5] =  -6.0 * ler / l^2
    s[2, 2] =  12.0 * ler / l^3
    s[2, 4] =   6.0 * ler / l^2
    s[3, 3] =         lea / l
    s[4, 2] =   6.0 * ler / l^2
    s[4, 4] =   4.0 * ler / l
    s[5, 1] =  -6.0 * ler / l^2
    s[5, 5] =   4.0 * ler / l
    s[6, 6] =         lgj / l
  end
  s
end

function stiffnessb2!(node, noofnodes, l, ler, lea, lgj)
  s = zeros(6, 6)
  if node == noofnodes
    return
  else
    s[1, 1] =  12.0 * ler / l^3
    s[1, 5] =   6.0 * ler / l^2
    s[2, 2] =  12.0 * ler / l^3
    s[2, 4] =  -6.0 * ler / l^2
    s[3, 3] =         lea / l
    s[4, 2] =  -6.0 * ler / l^2
    s[4, 4] =   4.0 * ler / l
    s[5, 1] =   6.0 * ler / l^2
    s[5, 5] =   4.0 * ler / l
    s[6, 6] =         lgj / l
  end
  s
end

function stiffnessc!(node, noofnodes, l, ler, lea, lgj)
  s = zeros(6, 6)
  if node == noofnodes
    return
  else
    s[1, 1] = -12.0 * ler / l^3
    s[1, 5] =   6.0 * ler / l^2
    s[2, 2] = -12.0 * ler / l^3
    s[2, 4] =  -6.0 * ler / l^2
    s[3, 3] =        -lea / l
    s[4, 2] =   6.0 * ler / l^2
    s[4, 4] =   2.0 * ler / l
    s[5, 1] =  -6.0 * ler / l^2
    s[5, 5] =   2.0 * ler / l
    s[6, 6] =        -lgj / l
  end
  s
end

function stiffnessi!(node, sb1, sb2)
  s = zeros(6, 6)
  b = sb1 + sb2
  h = b[1, 1] * b[5, 5] - b[1, 5]^2
  ff = b[5, 5] / h
  g = -b[1, 5] / h
  h = b[1, 1] / h
  s[1, 1] = ff
  s[1, 5] = g
  s[2, 2] = ff
  s[2, 4] = -g
  s[3, 3] = node == 1 ? 0.0 : 1.0 / b[3, 3]
  s[4, 2] = -g
  s[4, 4] = h
  s[5, 1] = g
  s[5, 5] = h
  s[6, 6] = node == 1 ? 0.0 : 1.0 / b[6, 6]
  s
end


text:
module LiiBRA

using UnitSystems, Parameters, LinearAlgebra, FFTW
using TSVD, Roots, Statistics, Interpolations
export C_e, Flux, C_se, Phi_s, Phi_e, Phi_se, DRA
export flatten_, R, F, Sim_Model, D_Linear, Construct, tuple_len, interp
export Realise, HPPC, fh!, mag!

include("Functions/Transfer/C_e.jl")
include("Functions/Transfer/C_se.jl")
include("Functions/Transfer/Flux.jl")
include("Functions/Transfer/Phi_s.jl")
include("Functions/Transfer/Phi_e.jl")
include("Functions/Transfer/Phi_se.jl")
include("Methods/DRA.jl")
include("Functions/Sim_Model.jl")

const F,R = faraday(Metric), universal(SI2019) #Faraday Constant / Universal Gas Constant

#---------- Generate Model -----------------#
function Realise(Cell, SList::Array, T::Float64)
    A = B = C = D = tuple()
    for i in SList
        #Arrhenius
        Cell.Const.T = T
        Arr_Factor = (1/Cell.Const.T_ref-1/Cell.Const.T)/R

        #Set Cell Constants
        Cell.Const.SOC = i
        Cell.Const.κ = Cell.Const.κf(Cell.Const.ce0)*exp(Cell.Const.Ea_κ*Arr_Factor)
        Cell.RA.Nfft = Cell.RA.Nfft!(Cell.RA.Fs, Cell.RA.Tlen)
        Cell.RA.f = Cell.RA.f!(Cell.RA.Nfft)
        Cell.RA.s = Cell.RA.s!(Cell.RA.Fs,Cell.RA.Nfft,Cell.RA.f)
        Cell.Neg.β = Cell.Neg.β!(Cell.RA.s)
        Cell.Pos.β = Cell.Pos.β!(Cell.RA.s)

        #Realisation
        Aϕ, Bϕ, Cϕ, Dϕ = DRA(Cell)

        #Flatten output into Tuples
        A = flatten_(A,Aϕ)
        B = flatten_(B,Bϕ)
        C = flatten_(C,Cϕ)
        D = flatten_(D,Dϕ)
    end
return A, B, C, D
end

#---------- Simulate Model -----------------#
function HPPC(Cell,SList::Array,SOC::Float64,λ::Float64,ϕ::Float64,A::Tuple,B::Tuple,C::Tuple,D::Tuple)

    #Set Experiment
    i = Int64(1/Cell.RA.SamplingT) #Sampling Frequency
    Input = [ones(1)*0.; ones(10*i)*λ; ones(40*i)*0.; ones(10*i)*ϕ; ones(40*i)*0.] #1C HPPC Experiment Current Profile
    Tk = ones(size(Input))*Cell.Const.T #Cell Temperature
    t = 0:(1.0/i):((length(Input)-1)/i)

    #Simulate Model
    return Sim_Model(Cell,Input,"Current",Tk,SList,SOC,A,B,C,D,t)
end

#---------- Hankel Formation & SVD -----------------#
function fh!(H,Hlen1,Hlen2,puls,M,Puls_L)
    for lp1 in 1:length(Hlen2), lp2 in 1:length(Hlen1)
            H[Puls_L*(lp2-1)+1:Puls_L*lp2,lp1] .= @view puls[:,Hlen2[lp1]+Hlen1[lp2]]
    end

    #Init = abs.(convert(Vector{float(eltype(H))}, randn(size(H,1))))
    Init = convert(Vector{float(eltype(H))}, ones(size(H,1)))
    U,S,V = tsvd(H, M; initvec=Init, maxiter=1e4)
    println("---")

    for lp1 in 1:length(Hlen2), lp2 in 1:length(Hlen1)
           H[Puls_L*(lp2-1)+1:Puls_L*lp2,lp1] .= @view puls[:,Hlen2[lp1]+Hlen1[lp2]+1]
    end
    return U,S,V'
end

"""
    D_Linear(Cell,ν_neg,ν_pos,σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, κ_eff_Sep)

    Function to linearise D array from input cell type and conductivites.

"""
function D_Linear(Cell, ν_neg, ν_pos, σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, κ_eff_Sep)
    D = Array{Float64}(undef,0,1)
    Dt = Array{Float64}(undef,0,1)
    D_ = Array{Float64}(undef,size(Cell.Transfer.Locs[2],1),1)
    q = Int64(1)
    tfs = Cell.Transfer.tfs
    Elec = Cell.Transfer.Elec
    Locs = Cell.Transfer.Locs
    for i in 1:size(tfs,1)
        if tfs[i] == C_e
            Dt =  zeros(length(Locs[i]))
        elseif tfs[i] == Phi_e
            Dt = Array{Float64}(undef,0,1)
            for pt in Locs[i]
                if pt <= Cell.Neg.L+eps()
                D_[q]  =  @. (Cell.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg*pt/Cell.Neg.L)) - pt*ν_neg*sinh(ν_neg) + Cell.Neg.L*(cosh(ν_neg)-cosh(ν_neg*(Cell.Neg.L-pt)/Cell.Neg.L)))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg)*ν_neg) #Lee. Eqn. 4.22 @ ∞
                elseif pt <= Cell.Neg.L + Cell.Sep.L + eps()
                D_[q] = @. (Cell.Neg.L - pt)/(Cell.Const.CC_A*κ_eff_Sep) + (Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) #Lee. Eqn. 4.23 @ ∞
                else
                D_[q] = @. -Cell.Sep.L/(Cell.Const.CC_A*κ_eff_Sep) + Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg)/(Cell.Const.CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg) + (Cell.Pos.L*(-σ_eff_Pos*cosh(ν_pos) + σ_eff_Pos*cosh((Cell.Const.Ltot-pt)*ν_pos/Cell.Pos.L) +  κ_eff_Pos*(cosh((pt-Cell.Neg.L-Cell.Sep.L)*ν_pos/Cell.Pos.L)-1)) - (pt-Cell.Neg.L-Cell.Sep.L)*κ_eff_Pos*sinh(ν_pos)*ν_pos)/(Cell.Const.CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos*sinh(ν_pos))
                end
                Dt = [Dt; D_[q]]
                q = q+1
            end

        elseif tfs[i] == C_se
            Dt = zeros(length(Locs[i]))
        elseif Elec[i] == "Pos"
            σ_eff = σ_eff_Pos
            κ_eff = κ_eff_Pos
            if tfs[i] == Phi_se
                Dt = @. -1*Cell.Pos.L/(Cell.Const.CC_A*ν_pos*sinh(ν_pos))*((1/κ_eff)*cosh(ν_pos*Locs[i])+(1/σ_eff)*cosh(ν_pos*(Locs[i]-1))) # Contribution to D as G->∞
            elseif tfs[i] == Flux
                Dt = @. -1*ν_pos*(σ_eff*cosh(ν_pos*Locs[i])+κ_eff*cosh(ν_pos*(Locs[i]-1)))/(Cell.Pos.as*F*Cell.Pos.L*Cell.Const.CC_A*(κ_eff+σ_eff)*sinh(ν_pos))
            elseif tfs[i] == Phi_s
                Dt = @. -1*(-Cell.Pos.L*(κ_eff*(cosh(ν_pos)-cosh(Locs[i]-1)*ν_pos))-Cell.Pos.L*(σ_eff*(1-cosh(Locs[i]*ν_pos)+Locs[i]*ν_pos*sinh(ν_pos))))/(Cell.Const.CC_A*σ_eff*(κ_eff+σ_eff)*ν_pos*sinh(ν_pos)) # Contribution to D as G->∞
            end
        elseif Elec[i] == "Neg"
            σ_eff = σ_eff_Neg
            κ_eff = κ_eff_Neg
            if tfs[i] == Phi_se
                Dt = @. Cell.Neg.L/(Cell.Const.CC_A*ν_neg*sinh(ν_neg))*((1/κ_eff)*cosh(ν_neg*Locs[i])+(1/σ_eff)*cosh(ν_neg*(Locs[i]-1))) # Contribution to D as G->∞
            elseif tfs[i] == Flux
                Dt = @. ν_neg*(σ_eff*cosh(ν_neg*Locs[i])+κ_eff*cosh(ν_neg*(Locs[i]-1)))/(Cell.Neg.as*F*Cell.Neg.L*Cell.Const.CC_A*(κ_eff+σ_eff)*sinh(ν_neg))
            elseif tfs[i] == Phi_s
                Dt = @. (-Cell.Neg.L*(κ_eff*(cosh(ν_neg)-cosh(Locs[i]-1)*ν_neg))-Cell.Neg.L*(σ_eff*(1-cosh(Locs[i]*ν_neg)+Locs[i]*ν_neg*sinh(ν_neg))))/(Cell.Const.CC_A*σ_eff*(κ_eff+σ_eff)*ν_neg*sinh(ν_neg)) # Contribution to D as G->∞
            end
        end
        D = [D; Dt]
    end
    return D
end

"""

    flatten_(a::Tuple, b...)

Flattens input Tuple "a" and inserts "b"

"""
function flatten_ end
flatten_() = ()
flatten_(a::Tuple) = Tuple(a)
flatten_(a) = (a,)
flatten_(a::Tuple, b...) = tuple(a..., flatten_(b...)...)
flatten_(a, b...) = tuple(a, flatten_(b...)...)
flatten_tuple(x::Tuple) = flatten_(x...)

"""
    Construct(::String)

    Function to create Construct dictionary.

    Currently supports:

    1. Doyle '94 parameterisation
    2. Chen 2020 parameterisation

"""
function Construct(CellType::String)
    if CellType == "Doyle 94"
        CellType = string(CellType,".jl")
        include(joinpath(dirname(pathof(LiiBRA)), "Data/Doyle_94", CellType))
    elseif CellType == "LG M50"
        CellType = string("LG_M50.jl")
        include(joinpath(dirname(pathof(LiiBRA)), "Data/Chen_2020", CellType))
    end
    return Cell
end

"""
    tuple_len(::NTuple)

    Function to return Tuple length.

"""
tuple_len(::NTuple{N, Any}) where {N} = N #Tuple Size

"""
    interp(MTup::Tuple,SList::Array,SOC)

    Function to interpolate Tuple indices

"""
function interp(MTup::Tuple,SList::Array,SOC)
    T1 = 0
    T2 = 0
    for i in 1:length(SList)-1
        if SList[i] > SOC >= SList[i+1]
            T2 = i
            T1 = i+1
        end
    end
    return M =  @. MTup[T1]+(MTup[T2]-MTup[T1])*(SOC-SList[T1])/(SList[T2]-SList[T1])
end

function mag!(γ::AbstractArray)
    ψ = Array{Float64}(undef,size(γ))
    for j in 1:size(γ,2), i in 1:size(γ,1)
        if real(γ[i,j]) < 0
            ψ[i,j] = -abs(γ[i,j])
        else
            ψ[i,j] = abs(γ[i,j])
        end
    end
    return ψ
end

function mag!(γ::AbstractVector)
    ψ = Vector{Float64}(undef,length(γ))
    for i in 1:length(γ)
        if real(γ[i]) < 0
            ψ[i] = -abs(γ[i])
        else
            ψ[i] = abs(γ[i])
        end
    end
    return ψ
end

end # module

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
@testset "Model permutation functions" begin
    source_pwm = [.7 .1 .1 .1
    .1 .1 .1 .7
    .1 .1 .7 .1]

    source_pwm_2 = [.6 .1 .1 .2
    .2 .1 .1 .6
    .1 .2 .6 .1]

    pwm_to_erode = [.25 .25 .25 .25
                    .97 .01 .01 .01
                    .01 .01 .01 .97
                    .01 .01 .97 .01
                    .25 .25 .25 .25]

    src_length_limits=2:5

    source_priors = assemble_source_priors(3, [source_pwm, source_pwm_2])
    mix_prior=.5

    bg_scores = log.(fill(.25, (12,2)))
    obs=[BioSequences.LongSequence{DNAAlphabet{2}}("ATGATGATGATG")
    BioSequences.LongSequence{DNAAlphabet{2}}("TGATGATGATGA")]
    order_seqs = BioBackgroundModels.get_order_n_seqs(obs, 0)
    coded_seqs = BioBackgroundModels.code_seqs(order_seqs)
    obs=Array(transpose(coded_seqs))
    obsl=[findfirst(iszero,obs[:,o])-1 for o in 1:size(obs)[2]]

    test_model = ICA_PWM_Model("test", source_priors, (falses(0,0),mix_prior), bg_scores, obs, src_length_limits)

    ps_model= permute_source(test_model, Vector{Model_Record}(), obs, obsl, bg_scores, test_model.log_Li, source_priors, iterates=1000,weight_shift_freq=.2,length_change_freq=.2)
    @test ps_model.log_Li > test_model.log_Li
    @test ps_model.sources != test_model.sources
    @test ps_model.mix_matrix == test_model.mix_matrix
    @test ps_model.origin == "PS from test"

    pm_model= permute_mix(test_model, obs, obsl, bg_scores, test_model.log_Li, iterates=1000)
    @test pm_model.log_Li > test_model.log_Li
    @test pm_model.sources == test_model.sources
    @test pm_model.mix_matrix != test_model.mix_matrix
    @test "PM from test" == pm_model.origin

    psfm_model=perm_src_fit_mix(test_model, Vector{Model_Record}(),obs, obsl, bg_scores, test_model.log_Li, source_priors,  iterates=1000)
    @test psfm_model.log_Li > test_model.log_Li
    @test psfm_model.sources != test_model.sources
    @test psfm_model.mix_matrix != test_model.mix_matrix
    @test "PSFM from test" == psfm_model.origin

    fm_model=fit_mix(test_model, Vector{Model_Record}(), obs,obsl,bg_scores)
    @test fm_model.log_Li > test_model.log_Li
    @test fm_model.sources == test_model.sources
    @test fm_model.mix_matrix != test_model.mix_matrix
    @test "FM from test" == fm_model.origin
    @test fit_mix in fm_model.permute_blacklist

    post_fm_psfm=perm_src_fit_mix(fm_model, Vector{Model_Record}(),obs, obsl, bg_scores, test_model.log_Li, source_priors, iterates=1000)
    @test "PSFM from candidate" == post_fm_psfm.origin
    @test fit_mix in post_fm_psfm.permute_blacklist

    rd_model=random_decorrelate(test_model, Vector{Model_Record}(),obs, obsl, bg_scores, test_model.log_Li, source_priors, iterates=1000)
    @test rd_model.log_Li > test_model.log_Li
    @test rd_model.sources != test_model.sources
    @test rd_model.mix_matrix != test_model.mix_matrix
    @test "RD from test" == rd_model.origin

    rs_model=reinit_src(test_model, Vector{Model_Record}(),obs, obsl, bg_scores, test_model.log_Li, source_priors, iterates=1000)
    @test rs_model.log_Li > test_model.log_Li
    @test rs_model.sources != test_model.sources
    @test rs_model.mix_matrix != test_model.mix_matrix
    @test "RS from test" == rs_model.origin

    erosion_sources=[(log.(source_pwm),1),(log.(source_pwm_2),1),(log.(pwm_to_erode),1)]

    eroded_mix=trues(2,3)

    erosion_lh=IPM_likelihood(erosion_sources,obs,obsl, bg_scores, eroded_mix)

    erosion_model=ICA_PWM_Model("erode", "", erosion_sources, test_model.source_length_limits,eroded_mix, erosion_lh)

    eroded_model=erode_model(erosion_model, Vector{Model_Record}(), obs, obsl, bg_scores, erosion_model.log_Li)
    @test eroded_model.log_Li > erosion_model.log_Li
    @test eroded_model.sources != erosion_model.sources
    @test eroded_model.mix_matrix == erosion_model.mix_matrix
    @test "EM from erode" == eroded_model.origin
    @test eroded_model.sources[1]==erosion_model.sources[1]
    @test eroded_model.sources[2]==erosion_model.sources[2]
    @test eroded_model.sources[3]!=erosion_model.sources[3]
    @test eroded_model.sources[3][1]==erosion_model.sources[3][1][2:4,:]

    dbl_eroded=erode_model(eroded_model, Vector{Model_Record}(), obs, obsl, bg_scores, eroded_model.log_Li)
    @test dbl_eroded.permute_blacklist == [erode_model]

    merger_srcs=   [(log.([.1 .7 .1 .1
    .1 .7 .1 .1
    .15 .35 .35 .15]),
    1
    ),
    (log.([.1 .7 .1 .1
        .1 .65 .15 .1
        .1 .7 .1 .1]),
        1
    ),
    (log.([.1 .7 .1 .1
        .1 .7 .1 .1
        .1 .7 .1 .1
        .7 .1 .1 .1]),
        1
    )]

    accurate_srcs=[ (log.([0.90 0.06 0.02 0.02; 0.02 0.02 0.02 0.94; 0.02 0.02 0.94 0.02]), 1),
    (log.([0.94 0.02 0.02 0.02; 0.02 0.07 0.07 0.84; 0.02 0.02 0.94 0.02]), 1),
    (log.([0.95 0.01 0.02 0.02; 0.02 0.02 0.02 0.94; 0.02 0.02 0.94 0.02]), 1)]

    merger_mix = BitMatrix([true false false
    true true false])

    accurate_mix=BitMatrix([true false false
    true true false])

    merger_base=ICA_PWM_Model("merge", "", merger_srcs,src_length_limits, merger_mix, IPM_likelihood(merger_srcs,obs,obsl, bg_scores, merger_mix))

    merger_target=ICA_PWM_Model("target", "", accurate_srcs, src_length_limits, accurate_mix, IPM_likelihood(accurate_srcs,obs,obsl,bg_scores,accurate_mix))

    path=randstring()
    test_record = Model_Record(path, merger_target.log_Li)
    serialize(path, merger_target)

    dm_model=distance_merge(merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li)
    @test dm_model.log_Li > merger_base.log_Li
    @test dm_model.sources != merger_base.sources
    distance_dict=Dict(1=>3,2=>1,3=>1)
    @test all([src==merger_base.sources[n] || src==merger_target.sources[distance_dict[n]] for (n,src) in enumerate(dm_model.sources)])
    @test "DM from merge" == dm_model.origin

    sm_model=similarity_merge(merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li)
    @test sm_model.log_Li > merger_base.log_Li
    @test sm_model.sources != merger_base.sources
    @test all([src==merger_base.sources[n] || src==merger_target.sources[n] for (n,src) in enumerate(sm_model.sources)])
    @test "SM from merge" == sm_model.origin

    shuffled_model=shuffle_sources(merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li)
    @test sum(shuffled_model.sources.==merger_base.sources)==length(shuffled_model.sources)-1
    for (s,src) in enumerate(shuffled_model.sources)
        @test src==merger_base.sources[s] || src==merger_target.sources[s]
        if src == merger_base.sources[s]
            @test shuffled_model.mix_matrix[:,s]==merger_base.mix_matrix[:,s]
        else
            @test shuffled_model.mix_matrix[:,s]==merger_target.mix_matrix[:,s]
        end
    end
    @test "SS from merge" == shuffled_model.origin

    acc_base_mix=falses(2,1);acc_base_mix[1,1]=true
    acc_merge_mix=falses(2,1);acc_merge_mix[2,1]=true

    acc_base=ICA_PWM_Model("accbase", "", [merger_srcs[1]],src_length_limits, acc_base_mix, IPM_likelihood([merger_srcs[1]],obs,obsl, bg_scores, acc_base_mix))
    acc_merge=ICA_PWM_Model("accmerge", "", [accurate_srcs[1]],src_length_limits, acc_merge_mix, IPM_likelihood([accurate_srcs[1]],obs,obsl, bg_scores, acc_merge_mix))

    accpath=randstring()
    acc_record = Model_Record(accpath, acc_merge.log_Li)
    serialize(accpath, acc_merge)

    acc_model=accumulate_mix(acc_base, [acc_record], obs, obsl, bg_scores, acc_merge.log_Li)
    @test acc_model.mix_matrix==trues(2,1)
    @test acc_model.sources == acc_merge.sources
    @test "AM from accbase" == acc_model.origin

    testwk=addprocs(1)[1]
    @everywhere import BioMotifInference

    ddm_model=remotecall_fetch(distance_merge, testwk, merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li, remote=true)
    @test ddm_model.log_Li > merger_base.log_Li
    @test ddm_model.sources != merger_base.sources
    @test all([src==merger_base.sources[n] || src==merger_target.sources[distance_dict[n]] for (n,src) in enumerate(ddm_model.sources)])
    @test "DM from merge" == ddm_model.origin

    dsm_model=remotecall_fetch(similarity_merge, testwk, merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li, remote=true)
    @test dsm_model.log_Li > merger_base.log_Li
    @test dsm_model.sources != merger_base.sources
    @test all([src==merger_base.sources[n] || src==merger_target.sources[n] for (n,src) in enumerate(dsm_model.sources)])
    @test "SM from merge" == dsm_model.origin

    dshuffled_model=remotecall_fetch(shuffle_sources, testwk, merger_base, [test_record], obs, obsl, bg_scores, merger_base.log_Li, remote=true)
    @test sum(dshuffled_model.sources.==merger_base.sources)==length(dshuffled_model.sources)-1
    for (s,src) in enumerate(dshuffled_model.sources)
        @test src==merger_base.sources[s] || src==merger_target.sources[s]
        if src == merger_base.sources[s]
            @test dshuffled_model.mix_matrix[:,s]==merger_base.mix_matrix[:,s]
        else
            @test dshuffled_model.mix_matrix[:,s]==merger_target.mix_matrix[:,s]
        end
    end
    @test "SS from merge" == dshuffled_model.origin

    dacc_model=remotecall_fetch(accumulate_mix, testwk, acc_base, [acc_record], obs, obsl, bg_scores, acc_merge.log_Li)
    @test dacc_model.mix_matrix==trues(2,1)
    @test dacc_model.sources == acc_merge.sources
    @test "AM from accbase" == dacc_model.origin

    info_base=ICA_PWM_Model("info_base", "", [accurate_srcs[1]], src_length_limits, accurate_mix[:,1:1], IPM_likelihood([accurate_srcs[1]],obs,obsl,bg_scores,accurate_mix[:,1:1]))
    info_model=info_fill(info_base,Vector{Model_Record}(),obs,obsl,bg_scores,info_base.log_Li)
    @test info_model.sources[1][1][1,:] == [0., -Inf, -Inf, -Inf]
    @test info_model.log_Li > info_base.log_Li
    @test "IF from info_base" == info_model.origin

    rmprocs(testwk)
    rm(path)
    rm(accpath)
end


text:
__precompile__(true)

module TransformUtils

using LinearAlgebra
using ManifoldsBase
using Manifolds
using StaticArrays
using DocStringExtensions

import Base: convert, promote_rule, *, \, vec
import LinearAlgebra: transpose, adjoint, normalize, normalize!
import ManifoldsBase: vee, vee!, hat, hat!

export
  Quaternion,
  AngleAxis,
  AxisAngle,
  SO3,
  so3,
  SE3,
  Euler,
  skew,
  vee!,
  vee,
  vec,
  setso3!,
  *,
  transpose,
  matrix,
  matrix!,
  inverse,
  compare,
  normalize!,
  normalize,
  q_conj,
  q_conj!,
  convert,
  convert!,
  rotate!,
  rotate,
  wrapRad,
  logmap,
  logmap_SO2,
  rightJacExmap,
  rightJacExmapinv,
  deltaso3vee,
  # Should be good to go
  veeAngleAxis,
  veeQuaternion,
  veeEuler,
  A_invB,
  ominus,
  oplus,
  ⊖,
  ⊕,
  \,
  doinversetransform!,

  # type aliases
  FloatInt,
  VectorFloatInt,

  # basic Taylor exported to try against Pade version in expm
  expmOwn,
  expmOwn1,
  expmOwnT,
  expmOwn2,
  expmOwn3,
  expmOwn4,

  # TODO -- needs refactoring
  R,
  SE2,
  se2vee,
  se2vee!,
  TU,
  veePose3,
  veePose

const TU = TransformUtils
const LA = LinearAlgebra

include("LazyLieManifolds.jl")
include("Legacy.jl")
include("LegacySpecialEuclidean.jl")

end #module


text:
using Test
using IPFitting, NBodyIPs, ProgressMeter, JuLIP
using NBodyIPs: blpolys
using JuLIP: decode_dict
Fit = IPFitting
DB = IPFitting.DB
Data = Fit.Data

function rand_data(sym, N, configtype="rand")
   cubic = (N == 1)
   at = bulk(sym, cubic=cubic) * N
   rattle!(at, 0.1)
   F = (N == 1) ? nothing : rand(3, length(at))
   return Dat(at, configtype; E = rand(), F = F, V = rand(3,3))
end

##
println("Double-Check (de-)dictionisation of basis: ")
basis1 = blpolys(2, ExpTransform(2.0, 3.0), CosCut(5.0, 7.0), 10)
basis2 = blpolys(3, ExpTransform(2.5, 3.0), CosCut2s(2.0, 2.5, 4.0, 5.5), 6)
println(@test decode_dict.( Dict.( basis1 ) ) == basis1)
println(@test decode_dict.( Dict.( basis2 ) ) == basis2)

println("Double-Check (de-)dictionisation of Dat: ")
data1 = [ rand_data(:Ti, 3, "md") for n = 1:10 ]
data2 = [ rand_data(:Ti, 1, "cell") for n = 1:10 ]
println(@test Dat.(Dict.(data1)) == data1)
println(@test Dat.(Dict.(data2)) == data2)

##
println("Create a temporary database.")
tmpdir = mktempdir()
dbpath = joinpath(tmpdir, "temp")
basis = [basis1; basis2]
data = [data1; data2]
db = nothing

try
   global db = DB.LsqDB(dbpath, basis, data)
   println(@test true)
catch
   @info("...something went wrong...")
   println(@test false)
end

println("checking consistency of db")
println(@test DB.dbpath(db) == dbpath)
println(@test DB.kronfile(dbpath) == dbpath * "_kron.h5")
println(@test DB.infofile(dbpath) == dbpath * "_info.json")
println(@test db.basis == basis)
println(@test db.configs == data)

println("re-load the database")
db1 = DB.LsqDB(dbpath)
println(@test DB.dbpath(db1) == dbpath)
println(@test db1.basis == db.basis)
println(@test db1.configs == db.configs)
println(@test db1.Ψ == db.Ψ)

##
println("re-load the database with mmap=true")
db2 = DB.LsqDB(dbpath, mmap=true)
println(@test DB.dbpath(db2) == dbpath)
println(@test db2.basis == db.basis)
println(@test db2.configs == db.configs)
println(@test db2.Ψ == db.Ψ)

##
println("Delete the temporary database")
rm(tmpdir; force=true, recursive=true)


text:
export savefig

const OUTPUT_FORMATS = ("eps", "ps", "pdf", "svg", "jpg", "png")
const OUTPUT_FORMATS_RASTER = ("jpg", "png")

# XXX
# --> make a "path" function
# --> make a "glecom" function

"""
    savefig
"""
function savefig(
            fig::Figure,
            fpath::String = "";
            #
            format::String  = "svg",
            resolution::Int = 200,
            res::Int        = resolution,
            info::Bool      = true
        )::String

    #
    # 1. check name and extension. If an extension is given
    #    it takes precedence over the 'format' keyword.
    #
    odir, fname = splitdir(fpath)
    fname       = ifelse(isempty(fname), fig.id, fname)
    fext        = splitext(fname)[2]
    fext        = ifelse(isempty(fext), format, fext)

    #
    # 2. check format and resolution.
    #
    format = lowercase(strip(fext))
    format in OUTPUT_FORMATS || throw(
        ArgumentError("Unknown output format $format.")
    )
    50 <= res <= 500 || throw(
        ArgumentError("Resolution $res outside of the expected [50, 500] range.")
    )

    #
    # 3. form temporary paths for GLE
    #
    base_path   = GP_ENV["TMP_PATH"] / fig.id
    script_path = base_path * ".gle"
    log_path    = base_path * ".log"

    #
    # 4. form the output path, creating dirs if necessary
    #
    isdir(odir) || mkpath(odir)
    output_dir  = odir
    output_path = output_dir / fname * "." * format

    #
    # 5. assemble the figure proper and write it to script_path
    #
    write_figure(fig, script_path)

    #
    # 6. gle command
    #
    # XXX
    cairo       = ifelse(isdef(fig.transparency), `-cairo`, ``)
    texlabels   = ifelse(isdef(fig.texlabels),    `-tex`,   ``)
    transparent = ifelse(
        isnothing(fig.bgcolor) || Colors.alpha(fig.bgcolor) < 1,
        `-transparent`, ``
    )

    cext  = `-d $format`
    cres  = ifelse(format in OUTPUT_FORMATS_RASTER, `-r $res`, ``)
    copt  = `$cairo $texlabels $transparent`
    cverb = `-vb 0`
    cout  = `-o $output_path $script_path`
    com   = pipeline(
        `gle $cext $cres $copt $cverb $cout`,
        stdout = log_path,
        stderr = log_path
    )

    start = time()

    if !success(com)
        log = read(log_path, String)
        # XXX GP_ENV["CLEAN_TMP"] && cleanup(fig)
        @error """
            GLE error: ...
            $log
            """
    end

    δt     = time() - start
    δt_str = ceil(Int, δt * 1e3)
    info && @info """
        Saved figure $(fig.id) at '$output_path' in $(δt_str)ms.
        """

    return output_path
end


text:
module TestSparse
using LinearAlgebra
using CovarianceFunctions
using SparseArrays
using Test

@testset "sparse" begin
    n = 128
    d = 32
    X = randn(d, n)
    k = CovarianceFunctions.EQ()
    δ = 1e-6
    G = gramian(k, X)
    S = sparse(G, δ)
    @test isapprox(G, S, atol = δ, rtol = δ)
end

end # module


text:
###################
# AbstractQuArray #
###################
    abstract AbstractQuArray{B<:AbstractBasis,T,N}

    # all subtypes of AbstractQuArray have to implement the following functions
    #   - coefftype
    #   - rawcoeffs
    #   - rawbases
    #   - copy
    #   - similar_type

    coeffs(qarr::AbstractQuArray) = rawcoeffs(qarr)

    # works generally as long as the single index form is defined
    bases(qarr::AbstractQuArray, i) = rawbases(qarr, i)
    bases(qarr::AbstractQuArray) = ntuple(ndims(qarr), i->bases(qarr, i))

    ########################
    # Array-like functions #
    ########################
    Base.size(qarr::AbstractQuArray, i...) = size(rawcoeffs(qarr), i...)
    Base.ndims(qarr::AbstractQuArray) = ndims(rawcoeffs(qarr))
    Base.length(qarr::AbstractQuArray) = length(rawcoeffs(qarr))

    Base.getindex(qarr::AbstractQuArray, i...) = getindex(rawcoeffs(qarr), i...)
    Base.setindex!(qarr::AbstractQuArray, i...) = setindex!(rawcoeffs(qarr), i...)

    Base.in(c, qarr::AbstractQuArray) = in(c, rawcoeffs(qarr))

    Base.(:(==))(a::AbstractQuArray, b::AbstractQuArray) = coeffs(a)==coeffs(b) && bases(a)==bases(b)

    typealias AbstractQuVector{B<:AbstractBasis,T} AbstractQuArray{B,T,1}
    typealias AbstractQuMatrix{B<:AbstractBasis,T} AbstractQuArray{B,T,2}

###########
# QuArray #
###########
    type QuArray{B<:AbstractBasis,T,N,A} <: AbstractQuArray{B,T,N}
        coeffs::A
        bases::NTuple{N,B}
        function QuArray(coeffs::AbstractArray{T,N}, bases::NTuple{N,B})
            if checkbases(coeffs, bases)
                new(coeffs, bases)
            else
                error("Coefficient array does not conform to input bases")
            end
        end
    end

    QuArray{B<:AbstractBasis,T,N}(coeffs::AbstractArray{T,N}, bases::NTuple{N,B}) = QuArray{B,T,N,typeof(coeffs)}(coeffs, bases)
    QuArray(coeffs, bases::AbstractBasis...) = QuArray(coeffs, bases)
    QuArray(coeffs) = QuArray(coeffs, map(FiniteBasis, size(coeffs)))

    typealias QuVector{B<:AbstractBasis,T,A} QuArray{B,T,1,A}
    typealias QuMatrix{B<:AbstractBasis,T,A} QuArray{B,T,2,A}

    # returns an appropriate outer constructor for the given (sub)type;
    # the constructor should construct an instance from a coefficient container
    # and a tuple of bases (see  QuArray(coeffs, bases) above)
    similar_type{Q<:QuArray}(::Type{Q}) = QuArray

    ######################
    # Accessor functions #
    ######################
    coefftype{B,T,N,A}(::QuArray{B,T,N,A}) = A

    rawcoeffs(qarr::QuArray) = qarr.coeffs

    rawbases(qarr::QuArray) = qarr.bases
    rawbases(qarr::QuArray, i) = qarr.bases[i]

    Base.copy(qa::QuArray) = QuArray(copy(qa.coeffs), copy(qa.bases))

##############
# CTranspose #
##############
    immutable CTranspose{B<:AbstractBasis,T,N,A} <: AbstractQuArray{B,T,N}
        qarr::QuArray{B,T,N,A}
        CTranspose(qarr::QuVector{B,T,A}) = new(qarr)
        CTranspose(qarr::QuMatrix{B,T,A}) = new(qarr)
        CTranspose(qarr::QuArray) = error("Conjugate transpose is unsupported for QuArrays of dimension $N")
    end

    CTranspose{B<:AbstractBasis,T,N,A}(qa::QuArray{B,T,N,A}) = CTranspose{B,T,N,A}(qa)
    CTranspose{B<:AbstractBasis,T,N}(coeffs::AbstractArray{T,N}, bases::NTuple{N,B}) = CTranspose(QuArray(coeffs, bases))

    typealias DualVector{B,T,A} CTranspose{B,T,1,A}
    typealias DualMatrix{B,T,A} CTranspose{B,T,2,A}

    similar_type{QC<:CTranspose}(::Type{QC}) = CTranspose

    ######################
    # Accessor functions #
    ######################
    coefftype{B,T,N,A}(::CTranspose{B,T,N,A}) = A

    rawcoeffs(ct::CTranspose) = rawcoeffs(ct.qarr)
    coeffs(ct::CTranspose) = rawcoeffs(ct)'

    revind(len, i) = len - (i-1)

    rawbases(ct::CTranspose, i) = rawbases(ct.qarr, i)
    rawbases(ct::CTranspose) =  rawbases(ct.qarr)

    bases(ct::CTranspose, i) = rawbases(ct, revind(ndims(ct), i))

    ########################
    # Array-like functions #
    ########################
    Base.copy(ct::CTranspose) = CTranspose(copy(ct.qarr))

    Base.ndims(ct::CTranspose) = ndims(ct.qarr)
    Base.length(ct::CTranspose) = length(ct.qarr)

    Base.size(ct::CTranspose) = reverse(size(ct.qarr))
    Base.size(ct::CTranspose, i) = size(ct, revind(ndims(ct), i))

    Base.getindex(ct::CTranspose, i, j) = getindex(ct.qarr, j, i)'
    Base.getindex(dv::DualVector, i) = getindex(dv.qarr, i)'

    Base.setindex!(ct::CTranspose, x, i, j) = setindex!(ct.qarr, x', j, i)
    Base.setindex!(dv::DualVector, x, i) = setindex!(dv.qarr, x', i)

    Base.ctranspose(qarr::QuArray) = CTranspose(qarr)
    Base.ctranspose(ct::CTranspose) = ct.qarr

################
# LabelQuArray #
################
    typealias LabelQuArray{B<:LabelBasis,T,N,A} QuArray{B,T,N,A}
    typealias TupleArray{T<:Tuple,N} Array{T,N}

    Base.getindex(larr::LabelQuArray, tups::Union(Tuple,TupleArray)...) = getindex(larr, map(getindex, bases(larr), tups)...)
    Base.setindex!(larr::LabelQuArray, x, tups::Union(Tuple,TupleArray)...) = setindex!(larr, x, map(getindex, bases(larr), tups)...)

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

    Base.promote_rule{B1,B2,T1,T2,N,A1,A2}(::Type{CTranspose{B1,T1,N,A1}}, ::Type{QuArray{B2,T2,N,A2}}) = QuArray{promote_type(B1,B2),promote_type(T1,T2),N,promote_type(A1,A2)}

######################
# Printing Functions #
######################
    Base.summary{B}(qarr::AbstractQuArray{B}) = "$(sizenotation(size(qarr))) $(typerepr(qarr)) in $B"
    Base.summary(larr::LabelQuArray) = "$(sizenotation(size(larr))) $(typerepr(larr)) in labeled bases $(bases(larr))"

    function Base.show(io::IO, qarr::AbstractQuArray)
        println(io, summary(qarr)*":")
        println(io, "...coefficients: $(coefftype(qarr))")
        print(io, repr(coeffs(qarr)))
    end

####################
# Helper Functions #
####################
    typerepr(qa::AbstractQuArray) = "$(typeof(qa).name)"
    typerepr(::QuVector) = "QuVector"
    typerepr(::QuMatrix) = "QuMatrix"
    typerepr(::DualVector) = "DualVector"
    typerepr(::DualMatrix) = "DualMatrix"
    typerepr(::QuArray) = "QuArray"

    sizenotation(@compat(tup::Tuple{Int})) = "$(first(tup))-element"
    sizenotation(@compat(tup::Tuple{Vararg{Int}})) = reduce(*, map(s->"$(s)x", tup))[1:end-1]

    function checkbases{N}(coeffs, bases::NTuple{N,AbstractBasis})
        if ndims(coeffs) == length(bases)
            return reduce(&, [checkcoeffs(coeffs, i, bases[i]) for i=1:N])
        end
        return false
    end

export QuArray, QuVector, QuMatrix, CTranspose, DualVector, DualMatrix,
    rawcoeffs,
    coeffs,
    coefftype,
    rawbases,
    bases


text:
function discreteIntegration(values::Array{Float64,1}, coordinates::AbstractArray{Float64,2}, l::Int, m::Int)

  @polyvar x y z

  sum = 0
  p = ylm(l,m,x,y,z)

  for k in 1:length(values)
    sum += values[k] * p((x,y,z)=>(coordinates[k,1],coordinates[k,2],coordinates[k,3]))
  end
  return sum * ((4*pi)/length(values))
end

function sphericalQuadrature(values::Array{Float64,1}, coordinates::AbstractArray{Float64,2}, L)

  C = [discreteIntegration(values,coordinates,l,m) for l=0:L for m=-l:l]
  return SphericalHarmonicCoefficients(C)
end

# Error propagation
function errorDiscreteIntegration(values::Array{Float64,1}, coordinates::AbstractArray{Float64,2}, l::Int, m::Int)

  @polyvar x y z

  sum = 0
  p = ylm(l,m,x,y,z)

  for k in 1:length(values)
    sum += (values[k] * p((x,y,z)=>(coordinates[k,1],coordinates[k,2],coordinates[k,3])))^2
  end
  return sqrt(sum) * ((4*pi)/length(values))
end

function errorSphericalQuadrature(values::Array{Float64,1}, coordinates::AbstractArray{Float64,2}, L)

  C = [errorDiscreteIntegration(values,coordinates,l,m) for l=0:L for m=-l:l]
  return SphericalHarmonicCoefficients(C)
end


text:
module GlobalMatchingModels
    using Distributions, Parameters, Statistics
    export compute_activation, compute_intensity, compute_similarity, encode, encode!, random_stimulus
    export compute_activations, compute_prob, compute_odds, generate_stimuli
    export MINERVA, REM

    include("MINERVA.jl")
    include("REM.jl")
end


text:
abstract type ActivationFunction end

struct Linear <: ActivationFunction end
(::Linear)(x::T) where {T <: Number} = x

struct ReLU <: ActivationFunction end
(::ReLU)(x::T) where {T <: Number} = max(zero(x), x)

struct Sigmoid <: ActivationFunction end
(::Sigmoid)(x::T) where {T <: Number} = 1 / (1 + exp(-x))

struct Square <: ActivationFunction end
(::Square)(x::T) where {T <: Number} = x^2

struct Tanh <: ActivationFunction end
(::Tanh)(x::T) where {T <: Number} = tanh(x)