text:
#
# This File contains the core routines for Slater-Koster type TB models
#
############################################################

"""
Assemble a hamiltonian block for a Slater-Koster type
hamiltonian, with 4 orbitals (s, p).

### Parameters:
* `U` : R / |R|, orientation of the bond, must be a 3-vector
* `hop` : vector of hopping function values, must be a 4-vector
* `mat` : output matrix, must be at least 4 x 4
"""
function _sk4!(U, hop, mat)
   l, m, n = U[1], U[2], U[3]
   # 4 orbitals are s, px, py, pz, these are the mat-indices
   # 4 bond types are : ssσ,spσ,ppσ,ppπ, these are the hop-indices
   mat[1,1] = hop[1]                            # E_ss = V_ssσ
   mat[2,2] = l*l * hop[3] + (1-l*l) * hop[4]   # E_xx = l² V_ppσ + (1-l²) V_ppπ
   mat[3,3] = m*m * hop[3] + (1-m*m) * hop[4]   # E_yy = m² V_ppσ + (1-m²) V_ppπ
   mat[4,4] = n*n * hop[3] + (1-n*n) * hop[4]   # E_zz = n² V_ppσ + (1-n²) V_ppπ
   mat[1,2] = l * hop[2]                        # E_sx = l V_spσ
   mat[1,3] = m * hop[2]                        # E_sy = m V_spσ
   mat[1,4] = n * hop[2]                        # E_sz = n V_spσ
   mat[2,1] = - mat[1,2]                        # E_xs = - E_sx
   mat[3,1] = - mat[1,3]                        # E_ys = - E_sy
   mat[4,1] = - mat[1,4]                        # E_zs = - E_sz
   mat[2,3] = l * m * (hop[3] - hop[4])         # E_xy = l m (V_ppσ - V_ppπ)
   mat[2,4] = l * n * (hop[3] - hop[4])         # E_xz = l n (V_ppσ - V_ppπ)
   mat[3,4] = m * n * (hop[3] - hop[4])         # E_yz = n m (V_ppσ - V_ppπ)
   mat[3,2] =  mat[2,3]                         # E_yx = E_xy
   mat[4,2] =  mat[2,4]                         # E_zx = E_xz
   mat[4,3] =  mat[3,4]                         # E_zy = E_yz
   return mat
end

"""
Assemble a hamiltonian block for a slater-koster type
hamiltonian, with 9 orbitals (s, p, d).

* `U` : orientation of the bond
* `hop` : hopping functions for bond type, 10-vector
* `mat` : output matrix 9 x 9
"""
function _sk9!(U, hop, mat)
   # fill the [1:4, 1:4] block
   _sk4!(U, hop, mat)
   # and then all the rest
   l, m, n = U[1], U[2], U[3]
   # ll, mm, nn = l*l, m*m, n*n    # TODO: for faster evaluation???
   # lm, ln, mn = l*m, l*n, m*n
   # rt3 = √3

   # sd
   mat[1,5] = √3 * l * m * hop[5]
   mat[1,6] = √3 * m * n * hop[5]
   mat[1,7] = √3 * l * n * hop[5]
   mat[1,8] = √3/2 * (l^2 - m^2) * hop[5]
   mat[1,9] = ( n^2 - (l^2 + m^2)/2 ) * hop[5]
   mat[5,1] = mat[1,5]
   mat[6,1] = mat[1,6]
   mat[7,1] = mat[1,7]
   mat[8,1] = mat[1,8]
   mat[9,1] = mat[1,9]

   # pd
   mat[2,5] = √3 * l * l * m * hop[6] + m * (1.0 - 2.0 * l^2) * hop[7]
   mat[2,6] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[2,7] = √3 * l * l * n * hop[6] + n * (1.0 - 2.0 * l^2) * hop[7]
   mat[2,8] = √3/2 * l * (l^2 - m^2) * hop[6] + l * (1.0 - l^2 + m^2) * hop[7]
   mat[2,9] = l * (n^2 - (l^2 + m^2)/2) * hop[6] - √3 * l * n^2 * hop[7]
   mat[5,2] = - mat[2,5]
   mat[6,2] = - mat[2,6]
   mat[7,2] = - mat[2,7]
   mat[8,2] = - mat[2,8]
   mat[9,2] = - mat[2,9]
   mat[3,5] = √3 * l * m * m * hop[6] + l * (1.0 - 2.0 * m^2) * hop[7]
   mat[3,6] = √3 * m * m * n * hop[6] + n * (1.0 - 2.0 * m^2) * hop[7]
   mat[3,7] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[3,8] = √3/2 * m * (l^2 - m^2) * hop[6] - m * (1.0 + l^2 - m^2) * hop[7]
   mat[3,9] = m * (n^2 - (l^2 + m^2)/2) * hop[6] - √3 * m * n^2 * hop[7]
   mat[5,3] = - mat[3,5]
   mat[6,3] = - mat[3,6]
   mat[7,3] = - mat[3,7]
   mat[8,3] = - mat[3,8]
   mat[9,3] = - mat[3,9]
   mat[4,5] = √3 * l * m * n * hop[6] - 2.0 * l * m * n * hop[7]
   mat[4,6] = √3 * m * n * n * hop[6] + m * (1.0 - 2.0 * n^2) * hop[7]
   mat[4,7] = √3 * l * n * n * hop[6] + l * (1.0 - 2.0 * n^2) * hop[7]
   mat[4,8] = √3/2 * n * (l^2 - m^2) * hop[6] - n * (l^2 - m^2) * hop[7]
   mat[4,9] = n * (n^2 - (l^2 + m^2)/2) * hop[6] + √3 * n * (l^2 + m^2) * hop[7]
   mat[5,4] = - mat[4,5]
   mat[6,4] = - mat[4,6]
   mat[7,4] = - mat[4,7]
   mat[8,4] = - mat[4,8]
   mat[9,4] = - mat[4,9]

   # dd
   mat[5,5] = 3.0 * l^2 * m^2 * hop[8] + (l^2 + m^2 - 4.0 * l^2 * m^2) * hop[9] +
               (n^2 + l^2 * m^2) * hop[10]
   mat[6,6] = 3.0 * m^2 * n^2 * hop[8] + (m^2 + n^2 - 4.0 * m^2 * n^2) * hop[9] +
               (l^2 + m^2 * n^2) * hop[10]
   mat[7,7] = 3.0 * l^2 * n^2 * hop[8] + (l^2 + n^2 - 4.0 * l^2 * n^2) * hop[9] +
               (m^2 + l^2 * n^2) * hop[10]
   mat[8,8] = 3.0/4 * (l^2 - m^2)^2 * hop[8] + (l^2 + m^2 - (l^2 - m^2)^2) * hop[9] +
               (n^2 + (l^2 - m^2)^2 /4 ) * hop[10]
   mat[9,9] = (n^2 - (l^2 + m^2) /2)^2 * hop[8] + 3.0 * n^2 * (l^2 + m^2) * hop[9] +
               3.0/4 * (l^2 + m^2)^2 * hop[10]
   mat[5,6] = 3.0 * l * m^2 * n * hop[8] + l * n * (1.0 - 4.0 * m^2) * hop[9] +
               l * n * (m^2 - 1.0) * hop[10]
   mat[5,7] = 3.0 * l^2 * m * n * hop[8] + m * n * (1.0 - 4.0 * l^2) * hop[9] +
               m * n * (l^2 - 1.0) * hop[10]
   mat[5,8] = 3.0/2 * l * m * (l^2 - m^2) * hop[8] + 2.0 * l * m * (m^2 - l^2) * hop[9] +
               1.0/2 * l * m * (l^2 - m^2) * hop[10]
   mat[5,9] = √3 * l * m * (n^2 - (l^2 + m^2)/2) * hop[8] - 2.0*√3 * l * m * n^2 * hop[9] +
               √3/2 * l * m * (1.0 + n^2) * hop[10]
   mat[6,7] = 3.0 * l * m * n^2 * hop[8] + l * m * (1.0 - 4.0 * n^2) * hop[9] +
               l * m * (n^2 - 1.0) * hop[10]
   mat[6,8] = 3.0/2 * m * n * (l^2 - m^2) * hop[8] -
               m * n * (1.0 + 2.0 * (l^2 - m^2)) * hop[9] +
               m * n * (1.0 + (l^2 - m^2) /2) * hop[10]
   mat[6,9] = √3 * m * n * (n^2 - (l^2 + m^2)/2) * hop[8] +
               √3 * m * n * (l^2 + m^2 - n^2) * hop[9] -
               √3/2 * m * n * (l^2 + m^2) * hop[10]
   mat[7,8] = 3.0/2 * l * n * (l^2 - m^2) * hop[8] +
               l * n * (1.0 - 2.0 * (l^2 - m^2)) * hop[9] -
               l * n * (1.0 - (l^2 - m^2) /2) * hop[10]
   mat[7,9] = √3 * l * n * (n^2 - (l^2 + m^2) /2) * hop[8] +
               √3 * l * n * (l^2 + m^2 - n^2) * hop[9] -
               √3/2 * l * n * (l^2 + m^2) * hop[10]
   mat[8,9] = √3/2 * (l^2 - m^2) * (n^2 - (l^2 + m^2) /2) * hop[8] +
               √3 * n^2 * (m^2 - l^2) * hop[9] +
               √3/4 * (1.0 + n^2) * (l^2 - m^2) * hop[10]
   mat[6,5] = mat[5,6]
   mat[7,5] = mat[5,7]
   mat[8,5] = mat[5,8]
   mat[9,5] = mat[5,9]
   mat[7,6] = mat[6,7]
   mat[8,6] = mat[6,8]
   mat[9,6] = mat[6,9]
   mat[8,7] = mat[7,8]
   mat[9,7] = mat[7,9]
   mat[9,8] = mat[8,9]
   return mat
end

# =============== Derivatives of the sk*! functions ================

#
# hh[i] = V(r)
# dhh[i] = V'(r)
#
# in other parts of the code hh ≡ bonds ≡ b  and dhh ≡ dbonds ≡ db
#
function _sk4_d!(U, r, hh, dhh, dmat)
   # U = ∇r, so it appears here both as the orientation and as ∇r
   l, m, n = U[1], U[2], U[3]
   ll, lm, ln, mm, mn, nn = l*l, l*m, l*n, m*m, m*n, n*n
   dl = ( (1.0-ll)/r ,     - lm/r ,     - ln/r )
   dm = (     - lm/r , (1.0-mm)/r ,     - mn/r )
   dn = (     - ln/r ,     - mn/r , (1.0-nn)/r )

   for d = 1:3
      dmat[d,1,1] = dhh[1] * U[d]
      dmat[d,2,2] = (ll * dhh[3] + (1. - ll) * dhh[4]) * U[d] + 2*l * hh[3] * dl[d] - 2*l * hh[4] * dl[d]
      dmat[d,3,3] = (mm * dhh[3] + (1. - mm) * dhh[4]) * U[d] + 2*m * hh[3] * dm[d] - 2*m * hh[4] * dm[d]
      dmat[d,4,4] = (nn * dhh[3] + (1. - nn) * dhh[4]) * U[d] + 2*n * hh[3] * dn[d] - 2*n * hh[4] * dn[d]
      dmat[d,1,2] = l * dhh[2] * U[d] + hh[2] * dl[d]
      dmat[d,1,3] = m * dhh[2] * U[d] + hh[2] * dm[d]
      dmat[d,1,4] = n * dhh[2] * U[d] + hh[2] * dn[d]
      dmat[d,2,1] = - dmat[d,1,2]
      dmat[d,3,1] = - dmat[d,1,3]
      dmat[d,4,1] = - dmat[d,1,4]
      dmat[d,2,3] = lm * (dhh[3] - dhh[4]) * U[d] + (dl[d] * m + l * dm[d]) * (hh[3] - hh[4])
      dmat[d,2,4] = ln * (dhh[3] - dhh[4]) * U[d] + (dl[d] * n + l * dn[d]) * (hh[3] - hh[4])
      dmat[d,3,4] = mn * (dhh[3] - dhh[4]) * U[d] + (dm[d] * n + m * dn[d]) * (hh[3] - hh[4])
      dmat[d,3,2] =  dmat[d,2,3]
      dmat[d,4,2] =  dmat[d,2,4]
      dmat[d,4,3] =  dmat[d,3,4]
   end
   return dmat
end

function _sk9_d!(U, r, hh, dhh, dmat; dim = 3)
   l, m, n = U[1], U[2], U[3]
   ll, lm, ln, mm, mn, nn = l*l, l*m, l*n, m*m, m*n, n*n
   dl = ( (1.0-ll)/r ,     - lm/r ,     - ln/r )
   dm = (     - lm/r , (1.0-mm)/r ,     - mn/r )
   dn = (     - ln/r ,     - mn/r , (1.0-nn)/r )

   for d = 1 : dim
       # ss
       dmat[d,1,1] = dhh[1] * U[d]
       # sp
       dmat[d,1,2] = l * dhh[2] * U[d] + hh[2] * dl[d]
       dmat[d,1,3] = m * dhh[2] * U[d] + hh[2] * dm[d]
       dmat[d,1,4] = n * dhh[2] * U[d] + hh[2] * dn[d]
       dmat[d,2,1] = - dmat[d,1,2]
       dmat[d,3,1] = - dmat[d,1,3]
       dmat[d,4,1] = - dmat[d,1,4]
       # pp
       dmat[d,2,2] = l*l * dhh[3] * U[d] + (1. -l*l) * dhh[4] * U[d] + 2*l * hh[3] * dl[d] - 2*l * hh[4] * dl[d]
       dmat[d,3,3] = m*m * dhh[3] * U[d] + (1. -m*m) * dhh[4] * U[d] + 2*m * hh[3] * dm[d] - 2*m * hh[4] * dm[d]
       dmat[d,4,4] = n*n * dhh[3] * U[d] + (1. -n*n) * dhh[4] * U[d] + 2*n * hh[3] * dn[d] - 2*n * hh[4] * dn[d]
       dmat[d,2,3] = l * m * (dhh[3] - dhh[4]) * U[d] + (dl[d] * m + l * dm[d]) * (hh[3] - hh[4])
       dmat[d,2,4] = l * n * (dhh[3] - dhh[4]) * U[d] + (dl[d] * n + l * dn[d]) * (hh[3] - hh[4])
       dmat[d,3,4] = m * n * (dhh[3] - dhh[4]) * U[d] + (dm[d] * n + m * dn[d]) * (hh[3] - hh[4])
       dmat[d,3,2] =  dmat[d,2,3]
       dmat[d,4,2] =  dmat[d,2,4]
       dmat[d,4,3] =  dmat[d,3,4]
       # sd
       dmat[d,1,5] = √3 * l * m * dhh[5] * U[d] + √3 * ( dl[d] * m + l * dm[d] ) * hh[5]
       dmat[d,1,6] = √3 * m * n * dhh[5] * U[d] + √3 * ( dm[d] * n + m * dn[d] ) * hh[5]
       dmat[d,1,7] = √3 * l * n * dhh[5] * U[d] + √3 * ( dl[d] * n + l * dn[d] ) * hh[5]
       dmat[d,1,8] = √3/2 * (l^2 - m^2) * dhh[5] * U[d] + √3/2 * ( 2*l * dl[d] - 2*m * dm[d]) * hh[5]
       dmat[d,1,9] = ( n^2 - (l^2 + m^2)/2 ) * dhh[5] * U[d] + ( 2*n * dn[d] - (l * dl[d] + m * dm[d]) ) * hh[5]
       dmat[d,5,1] = dmat[d,1,5]
       dmat[d,6,1] = dmat[d,1,6]
       dmat[d,7,1] = dmat[d,1,7]
       dmat[d,8,1] = dmat[d,1,8]
       dmat[d,9,1] = dmat[d,1,9]
       # pd
       dmat[d,2,5] = √3 * l * l * m * dhh[6] * U[d] + m * (1.0 - 2.0 * l^2) * dhh[7] * U[d] +
                   √3 * (2*l * m * dl[d] + l^2 * dm[d]) * hh[6] +
                   ( dm[d] * (1.0 - 2.0 * l^2) - 4*m*l * dl[d] ) * hh[7]
       dmat[d,2,6] = √3 * l * m * n * dhh[6] * U[d] - 2.0 * l * m * n * dhh[7] * U[d] +
                   √3 * (dl[d] * m * n + l * dm[d] * n + l * m * dn[d]) * hh[6] -
                   2.0 * (dl[d] * m * n + l * dm[d] * n + l * m * dn[d]) * hh[7]
       dmat[d,2,7] = √3 * l * l * n * dhh[6] * U[d] + n * (1.0 - 2.0 * l^2) * dhh[7] * U[d] +
                   √3 * (2*l * dl[d] * n + l*l * dn[d]) * hh[6] +
                   ( dn[d] * (1.0 - 2.0 * l^2) - 4*l*n * dl[d] ) * hh[7]
       dmat[d,2,8] = √3/2 * l * (l^2 - m^2) * dhh[6] * U[d] + l * (1.0 - l^2 + m^2) * dhh[7] * U[d] +
                   √3/2 * (dl[d] * (l^2 - m^2) + l * (2*l * dl[d] - 2*m * dm[d])) * hh[6] +
                   (dl[d] * (1.0 - l^2 + m^2) + l * (-2*l * dl[d] + 2*m * dm[d]) ) * hh[7]
       dmat[d,2,9] = l * (n^2 - (l^2 + m^2)/2) * dhh[6] * U[d] - √3 * l * n^2 * dhh[7] * U[d] +
                   ( dl[d] * (n^2 - (l^2 + m^2)/2) + l * (2*n * dn[d] - l * dl[d] - m * dm[d]) ) * hh[6] -
                   √3 * (dl[d] * n^2 + 2*l*n * dn[d])* hh[7]
       dmat[d,5,2] = - dmat[d,2,5]
       dmat[d,6,2] = - dmat[d,2,6]
       dmat[d,7,2] = - dmat[d,2,7]
       dmat[d,8,2] = - dmat[d,2,8]
       dmat[d,9,2] = - dmat[d,2,9]
       dmat[d,3,5] = √3 * l * m * m * dhh[6] * U[d] + l * (1.0 - 2.0 * m^2) * dhh[7] * U[d] +
                   √3 * (dl[d] * m^2 + 2*l*m * dm[d]) * hh[6] +
                   ( dl[d] * (1.0 - 2.0 * m^2) - 4*l*m * dm[d] ) * hh[7]
       dmat[d,3,6] = √3 * m * m * n * dhh[6] * U[d] + n * (1.0 - 2.0 * m^2) * dhh[7] * U[d] +
                   √3 * (2*m*n * dm[d] + m^2 * dn[d]) * hh[6] +
                   ( dn[d] * (1.0 - 2.0 * m^2) - 4*m*n * dm[d] ) * hh[7]
       dmat[d,3,7] = √3 * l * m * n * dhh[6] * U[d] - 2.0 * l * m * n * dhh[7] * U[d] +
                   ( dl[d] * m * n + l * dm[d] * n + l * m * dn[d] ) * ( √3 * hh[6] - 2.0 * hh[7] )
       dmat[d,3,8] = √3/2 * m * (l^2 - m^2) * dhh[6] * U[d] - m * (1.0 + l^2 - m^2) * dhh[7] * U[d] +
                   √3/2 * ( dm[d] * (l^2 - m^2) + m * (2*l * dl[d] - 2*m * dm[d]) ) * hh[6] -
                   ( dm[d] * (1.0 + l^2 - m^2) + m * (2*l * dl[d] - 2*m * dm[d]) ) * hh[7]
       dmat[d,3,9] = m * (n^2 - (l^2 + m^2)/2) * dhh[6] * U[d] - √3 * m * n^2 * dhh[7] * U[d] +
                   ( dm[d] * (n^2 - (l^2 + m^2)/2) + m * (2*n * dn[d] - l * dl[d] - m * dm[d]) ) * hh[6] -
                   √3 * ( dm[d] * n^2 + 2*m*n * dn[d] ) * hh[7]
       dmat[d,5,3] = - dmat[d,3,5]
       dmat[d,6,3] = - dmat[d,3,6]
       dmat[d,7,3] = - dmat[d,3,7]
       dmat[d,8,3] = - dmat[d,3,8]
       dmat[d,9,3] = - dmat[d,3,9]
       dmat[d,4,5] = √3 * l * m * n * dhh[6] * U[d] - 2.0 * l * m * n * dhh[7] * U[d] +
                   ( dl[d] * m * n + l * dm[d] * n + l * m * dn[d] ) * ( √3 * hh[6] - 2.0 * hh[7] )
       dmat[d,4,6] = √3 * m * n * n * dhh[6] * U[d] + m * (1.0 - 2.0 * n^2) * dhh[7] * U[d] +
                   √3 * ( dm[d] * n^2 + 2*m*n * dn[d] ) * hh[6] +
                   ( dm[d] * (1.0 - 2.0 * n^2) - 4*m*n * dn[d] ) * hh[7]
       dmat[d,4,7] = √3 * l * n * n * dhh[6] * U[d] + l * (1.0 - 2.0 * n^2) * dhh[7] * U[d] +
                   √3 * ( dl[d] * n^2 + 2*l*n * dn[d] ) * hh[6] +
                   ( dl[d] * (1.0 - 2.0 * n^2) - 4*l*n * dn[d] ) * hh[7]
       dmat[d,4,8] = √3/2 * n * (l^2 - m^2) * dhh[6] * U[d] - n * (l^2 - m^2) * dhh[7] * U[d] +
                   √3/2 * ( dn[d] * (l^2 - m^2) + n * (2*l * dl[d] - 2*m * dm[d]) ) * hh[6] -
                   ( dn[d] * (l^2 - m^2) + n * (2*l * dl[d] - 2*m * dm[d]) ) * hh[7]
       dmat[d,4,9] = n * (n^2 - (l^2 + m^2)/2) * dhh[6] * U[d] + √3 * n * (l^2 + m^2) * dhh[7] * U[d] +
                   ( dn[d] * (n^2 - (l^2 + m^2)/2) +  n * (2*n * dn[d] - l * dl[d] - m * dm[d]) ) * hh[6] +
                   √3 * ( dn[d] * (l^2 + m^2) + n * (2*l * dl[d] + 2*m * dm[d]) ) * hh[7]
       dmat[d,5,4] = - dmat[d,4,5]
       dmat[d,6,4] = - dmat[d,4,6]
       dmat[d,7,4] = - dmat[d,4,7]
       dmat[d,8,4] = - dmat[d,4,8]
       dmat[d,9,4] = - dmat[d,4,9]
       # dd
       dmat[d,5,5] = 3.0 * l^2 * m^2 * dhh[8] * U[d] + (l^2 + m^2 - 4.0 * l^2 * m^2) * dhh[9] * U[d] +
             (n^2 + l^2 * m^2) * dhh[10] * U[d] +
             3.0 * ( 2*l * dl[d] * m^2 + l^2 * 2*m * dm[d] ) * hh[8] +
                   (2*l * dl[d] + 2*m * dm[d] - 8*l * dl[d] * m^2 - 8*m * l^2 * dm[d]) * hh[9] +
                   (2*n * dn[d] + 2*l * dl[d] * m^2 + 2*m * l^2 * dm[d]) * hh[10]
       dmat[d,6,6] = 3.0 * m^2 * n^2 * dhh[8] * U[d] + (m^2 + n^2 - 4.0 * m^2 * n^2) * dhh[9] * U[d] +
             (l^2 + m^2 * n^2) * dhh[10] * U[d] +
                   3.0 * ( 2*m * dm[d] * n^2 + m^2 * 2*n * dn[d] ) * hh[8] +
                   (2*m * dm[d] + 2*n * dn[d] - 8*m * dm[d] * n^2 - 8*n * m^2 * dn[d]) * hh[9] +
                   (2*l * dl[d] + 2*m * dm[d] * n^2 + 2*n * m^2 * dn[d]) * hh[10]
       dmat[d,7,7] = 3.0 * l^2 * n^2 * dhh[8] * U[d] + (l^2 + n^2 - 4.0 * l^2 * n^2) * dhh[9] * U[d] +
             (m^2 + l^2 * n^2) * dhh[10] * U[d] +
                   3.0 * ( 2*l * dl[d] * n^2 + l^2 * 2*n * dn[d] ) * hh[8] +
                   (2*l * dl[d] + 2*n * dn[d] - 8*l * dl[d] * n^2 - 8*n * l^2 * dn[d]) * hh[9] +
                   (2*m * dm[d] + 2*l * dl[d] * n^2 + 2*n * l^2 * dn[d]) * hh[10]
       dmat[d,8,8] = 3.0/4 * (l^2 - m^2)^2 * dhh[8] * U[d] + (l^2 + m^2 - (l^2 - m^2)^2) * dhh[9] * U[d] +
             (n^2 + (l^2 - m^2)^2 /4) * dhh[10] * U[d] +
                   3.0 * (l^2 - m^2) * (l * dl[d] - m * dm[d]) * hh[8] +
                   ( 2*l * dl[d] + 2*m * dm[d] - 4.0 * (l^2 - m^2) * (l * dl[d] - m * dm[d]) ) * hh[9] +
                   ( 2*n * dn[d] + (l^2 - m^2) * (l * dl[d] - m * dm[d]) ) * hh[10]
       dmat[d,9,9] = (n^2 - (l^2 + m^2) /2)^2 * dhh[8] * U[d] + 3.0 * n^2 * (l^2 + m^2) * dhh[9] * U[d] +
             3.0/4 * (l^2 + m^2)^2 * dhh[10] * U[d] +
                   2.0 * (n^2 - (l^2 + m^2) /2) * (2*n * dn[d] - (l * dl[d] + m * dm[d])) * hh[8] +
                   3.0 * ( 2*n * dn[d] * (l^2 + m^2) + n^2 * (2*l * dl[d] + 2*m * dm[d]) ) * hh[9] +
                   3.0 * (l^2 + m^2) * (l * dl[d] + m * dm[d]) * hh[10]
       dmat[d,5,6] = 3.0 * l * m^2 * n * dhh[8] * U[d] + l * n * (1.0 - 4.0 * m^2) * dhh[9] * U[d] +
             l * n * (m^2 - 1.0) * dhh[10] * U[d] +
                   3.0 * (dl[d] * m^2 * n + l * 2*m * dm[d] * n + l * m^2 * dn[d]) * hh[8] +
                   ( dl[d] * n * (1.0 - 4.0 * m^2) + l * dn[d] * (1.0 - 4.0 * m^2) - l * n * 8*m * dm[d] ) * hh[9] +
                   ( dl[d] * n * (m^2 - 1.0) + l * dn[d] * (m^2 - 1.0) + l * n * 2*m * dm[d] ) * hh[10]
       dmat[d,5,7] = 3.0 * l^2 * m * n * dhh[8] * U[d] + m * n * (1.0 - 4.0 * l^2) * dhh[9] * U[d] +
             m * n * (l^2 - 1.0) * dhh[10] * U[d] +
                   3.0 * (2*l * dl[d] * m * n + l^2 * dm[d] * n + l^2 * m * dn[d]) * hh[8] +
                   ( dm[d] * n * (1.0 - 4.0 * l^2) + m * dn[d] * (1.0 - 4.0 * l^2) - m * n * 8*l * dl[d] ) * hh[9] +
                   ( dm[d] * n * (l^2 - 1.0) + m * dn[d] * (l^2 - 1.0) + m * n * 2*l * dl[d] ) * hh[10]
       dmat[d,5,8] = 3.0/2 * l * m * (l^2 - m^2) * dhh[8] * U[d] + 2.0 * l * m * (m^2 - l^2) * dhh[9] * U[d] +
             1.0/2 * l * m * (l^2 - m^2) * dhh[10] * U[d] +
                   3.0/2 * ( dl[d] * m * (l^2 - m^2) + l * dm[d] * (l^2 - m^2) + l * m * (2*l * dl[d] - 2*m * dm[d]) ) * hh[8] +
                   2.0 * ( dl[d] * m * (m^2 - l^2) + l * dm[d] * (m^2 - l^2) + l * m * (2*m * dm[d] - 2*l * dl[d]) ) * hh[9] +
                   ( dl[d] * m * (l^2 - m^2)/2 + l * dm[d] * (l^2 - m^2)/2 + l * m * (l*dl[d] - m*dm[d]) ) * hh[10]
       dmat[d,5,9] = √3 * l * m * (n^2 - (l^2 + m^2)/2) * dhh[8] * U[d] - 2.0*√3 * l * m * n^2 * dhh[9] * U[d] +
             √3/2 * l * m * (1.0 + n^2) * dhh[10] * U[d] +
                   √3 * ( dl[d] * m * (n^2 - (l^2 + m^2)/2) + l * dm[d] * (n^2 - (l^2 + m^2)/2) + l * m * (2*n * dn[d] - (l * dl[d] + m * dm[d])) ) * hh[8] -
                   2.0*√3 * ( dl[d] * m * n^2 + l * dm[d] * n^2 + l * m * 2*n * dn[d] ) * hh[9] +
                   √3/2 * ( dl[d] * m * (1.0 + n^2) + l * dm[d] * (1.0 + n^2) + l * m * 2*n * dn[d] ) * hh[10]
       dmat[d,6,7] = 3.0 * l * m * n^2 * dhh[8] * U[d] + l * m * (1.0 - 4.0 * n^2) * dhh[9] * U[d] +
             l * m * (n^2 - 1.0) * dhh[10] * U[d] +
                   3.0 * ( dl[d] * m * n^2 + l * dm[d] * n^2 + l * m * 2*n * dn[d] ) * hh[8] +
                   ( dl[d] * m * (1.0 - 4.0 * n^2) + l * dm[d] * (1.0 - 4.0 * n^2) - l * m * 8*n * dn[d] ) * hh[9] +
                   ( dl[d] * m * (n^2 - 1.0) + l * dm[d] * (n^2 - 1.0) + l * m * 2*n * dn[d] ) * hh[10]
       dmat[d,6,8] = 3.0/2 * m * n * (l^2 - m^2) * dhh[8] * U[d] - m * n * (1.0 + 2.0 * (l^2 - m^2)) * dhh[9] * U[d] +
             m * n * (1.0 + (l^2 - m^2) /2) * dhh[10] * U[d] +
                   3.0/2 * ( dm[d] * n * (l^2 - m^2) + m * dn[d] * (l^2 - m^2) + m * n * (2*l * dl[d] - 2*m * dm[d]) ) * hh[8] -
                   ( dm[d] * n * (1.0 + 2.0 * (l^2 - m^2)) + m * dn[d] * (1.0 + 2.0 * (l^2 - m^2)) + m * n * 4.0 * (l * dl[d] - m * dm[d]) ) * hh[9] +
                   ( dm[d] * n * (1.0 + (l^2 - m^2) /2) + m * dn[d] * (1.0 + (l^2 - m^2) /2) + m * n * (l * dl[d] - m * dm[d]) ) * hh[10]
       dmat[d,6,9] = √3 * m * n * (n^2 - (l^2 + m^2) /2) * dhh[8] * U[d] + √3 * m * n * (l^2 + m^2 - n^2) * dhh[9] * U[d] -
             √3/2 * m * n * (l^2 + m^2) * dhh[10] * U[d] +
                   √3 * ( dm[d] * n * (n^2 - (l^2 + m^2) /2) + m * dn[d] * (n^2 - (l^2 + m^2) /2) +  m * n * (2*n * dn[d] - (l * dl[d] + m * dm[d])) ) * hh[8] +
                   √3 * ( dm[d] * n * (l^2 + m^2 - n^2) + m * dn[d] * (l^2 + m^2 - n^2) + m * n * (2*l * dl[d] + 2*m * dm[d] - 2*n * dn[d]) ) * hh[9] -
                   √3/2 * ( dm[d] * n * (l^2 + m^2) + m * dn[d] * (l^2 + m^2) + m * n * (2*l * dl[d] + 2*m * dm[d]) ) * hh[10]
       dmat[d,7,8] = 3.0/2 * l * n * (l^2 - m^2) * dhh[8] * U[d] + l * n * (1.0 - 2.0 * (l^2 - m^2)) * dhh[9] * U[d] -
             l * n * (1.0 - (l^2 - m^2) /2) * dhh[10] * U[d] +
                   3.0/2 * ( dl[d] * n * (l^2 - m^2) + l * dn[d] * (l^2 - m^2) + l * n * (2*l * dl[d] - 2*m * dm[d]) ) * hh[8] +
                   ( dl[d] * n * (1.0 - 2.0 * (l^2 - m^2)) + l * dn[d] * (1.0 - 2.0 * (l^2 - m^2)) - l * n * 4.0 * (l * dl[d] - m * dm[d]) ) * hh[9] -
                   ( dl[d] * n * (1.0 - (l^2 - m^2) /2) + l * dn[d] * (1.0 - (l^2 - m^2) /2) - l * n * (l * dl[d] - m * dm[d]) ) * hh[10]
       dmat[d,7,9] = √3 * l * n * (n^2 - (l^2 + m^2) /2) * dhh[8] * U[d] + √3 * l * n * (l^2 + m^2 - n^2) * dhh[9] * U[d] -
             √3/2 * l * n * (l^2 + m^2) * dhh[10] * U[d] +
                   √3 * ( dl[d] * n * (n^2 - (l^2 + m^2) /2) + l * dn[d] * (n^2 - (l^2 + m^2) /2) + l * n * (2*n * dn[d] - (l * dl[d] + m * dm[d])) ) * hh[8] +
                   √3 * ( dl[d] * n * (l^2 + m^2 - n^2) + l * dn[d] * (l^2 + m^2 - n^2) + l * n * (2*l * dl[d] + 2*m * dm[d] - 2*n * dn[d]) ) * hh[9] -
                   √3/2 * ( dl[d] * n * (l^2 + m^2) + l * dn[d] * (l^2 + m^2) + l * n * (2*l * dl[d] + 2*m * dm[d]) ) * hh[10]
       dmat[d,8,9] = √3/2 * (l^2 - m^2) * (n^2 - (l^2 + m^2) /2) * dhh[8] * U[d] + √3 * n^2 * (m^2 - l^2) * dhh[9] * U[d] +
             √3/4 * (1.0 + n^2) * (l^2 - m^2) * dhh[10] * U[d] +
                   √3/2 * ( (2*l * dl[d] - 2*m * dm[d]) * (n^2 - (l^2 + m^2) /2) + (l^2 - m^2) * (2*n * dn[d] - (l * dl[d] + m * dm[d])) ) * hh[8] +
                   √3 * ( 2*n * dn[d] * (m^2 - l^2) + n^2 * (2*m * dm[d] - 2*l * dl[d]) ) * hh[9] +
                   √3/4 * ( 2*n * dn[d] * (l^2 - m^2) + (1.0 + n^2) * (2*l * dl[d] - 2*m * dm[d]) ) * hh[10]
       dmat[d,6,5] = dmat[d,5,6]
       dmat[d,7,5] = dmat[d,5,7]
       dmat[d,8,5] = dmat[d,5,8]
       dmat[d,9,5] = dmat[d,5,9]
       dmat[d,7,6] = dmat[d,6,7]
       dmat[d,8,6] = dmat[d,6,8]
       dmat[d,9,6] = dmat[d,6,9]
       dmat[d,8,7] = dmat[d,7,8]
       dmat[d,9,7] = dmat[d,7,9]
       dmat[d,9,8] = dmat[d,8,9]
   end

   return dmat
end


text:
using Base
using Base.Test

@testset "Random Integer" begin
    src::Vector{Int64} = [62,21,65,27,78,23,97,53,38,67]
    count::Int64 = 10
    seed::UInt32 = 17
    r::UnitRange{Int64} = 1:100

    # - Generate a new array - #
    srand(seed) # - Seed the GLOBAL_RNG
    x::Vector{Int64} = rand(r, count)
    @test x == src

    # - In-place generation - #
    srand(seed) # - Seed the GLOBAL_RNG
    x = ones(x) # - Reset
    rand!(x, r)
    @test x == src
end

@testset "Random Float64" begin
    src::Vector{Float64} = [
        0.08110913680564935, 0.5205151824383689, 0.3607870878129784,
        0.5461467030095473, 0.9463823998310394, 0.004177010042717644,
        0.5210377372245407, 0.6744743681550878, 0.8823899341292194,
        0.3139038622430719]

    count::Int64 = 10
    seed::UInt32 = 17

    # - new array - #
    srand(seed)
    x::Vector{Float64} = rand(Float64, count)
    @test x == src

    # - In-place generation - #
    srand(seed)
    rand!(x)
    @test x == src
end

ext: jl


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

abstract type AbstractCoordinateVector{T, N, S} <: StaticArrays.FieldVector{N, T} end

struct CartesianVector{ T <: RealQuantity } <: AbstractCoordinateVector{T, 3, :cartesian}
    x::T
    y::T
    z::T
end

function (+)(p::CartesianPoint{T}, v::CartesianVector{T})::CartesianPoint{T} where {T <: Real}
    return CartesianPoint{T}( p.x + v.x, p.y + v.y, p.z + v.z )
end
function (+)(u::CartesianVector{T}, v::CartesianVector{T})::CartesianVector{T} where {T <: Real}
    return CartesianVector{T}( u.x + v.x, u.y + v.y, u.z + v.z )
end
function broadcast(f::Function, v::Vector{CartesianPoint{T}}, p::CartesianVector{T}) where {T}
    f(v ,[p for i in eachindex(v)])
end

function (-)(p::CartesianPoint{T}, v::CartesianVector{T})::CartesianPoint{T} where {T <: Real}
    return CartesianPoint{T}( p.x - v.x, p.y - v.y, p.z - v.z )
end
function (-)(u::CartesianVector{T}, v::CartesianVector{T})::CartesianVector{T} where {T <: Real}
    return CartesianVector{T}( u.x - v.x, u.y - v.y, u.z - v.z )
end
function (-)(p1::CartesianPoint{T}, p2::CartesianPoint{T})::CartesianVector{T} where {T <: Real}
    return CartesianVector{T}( p1.x - p2.x, p1.y - p2.y, p1.z - p2.z )
end

struct CylindricalVector{ T <: RealQuantity } <: AbstractCoordinateVector{T, 3, :cylindrical}
    r::T
    φ::T # in radian
    z::T
end

function (+)(p::CylindricalPoint{T}, v::CylindricalVector{T})::CylindricalPoint{T} where {T <: Real}
    error("Not yet defined")
end
function (-)(p::CylindricalPoint{T}, v::CylindricalVector{T})::CylindricalPoint{T} where {T <: Real}
    error("Not yet defined")
end
function (-)(p1::CylindricalPoint{T}, p2::CylindricalPoint{T})::CylindricalPoint{T} where {T <: Real}
    error("Not yet defined")
end


text:
# AutoRegressive Normal Model.
#

using LinearAlgebra
using PaddedViews

export ARNormal1D
export ARNormal
export meanfilter

#######################################################################
# Regressors{X}D is an iterable to access dynamically the context of a
# given sample.
#

struct Regressors1D{T <: AbstractFloat, K}
    x::PaddedView{<: AbstractFloat}
end

function Base.getindex(reg::Regressors1D{T, K}, n::Int) where {T <: AbstractFloat, K}
    r = reg.x[n-K:n]
    r[end, end] = 1 # "bias" regressor
    return r
end

#######################################################################
# 1D AutoRegressive Normal.

struct ARNormal1D{K} <: Model
    # IIR Filter
    h::ConjugateParameter{<:Normal}

    # Noise precision
    λ::ConjugateParameter{<:Gamma}

    function ARNormal1D(h, λ) where K
        μ₀, _ = stdparam(h.prior)
        model = new{length(μ₀) - 1}(h, λ)
        h.stats = (X) -> stats1D_h(model, X)
        λ.stats = (X) -> stats1D_λ(model, X)
        return model
    end
end

function ARNormal1D(μ₀::Vector{T}, Σ₀::Matrix{T}, μ::Vector{T},
                    Σ::Matrix{T}, a₀::Real, b₀::Real, a::Real, b::Real) where T <: AbstractFloat
    ARNormal1D(
        ConjugateParameter(Normal(μ₀, Σ₀), Normal(μ, Σ)),
        ConjugateParameter(Gamma([T(a₀)], [T(b₀)]), Gamma([T(a)], [T(b)])),
    )
end

function ARNormal1D(T; order::Integer, initσ::Real = 0.0,
                    pseudocounts::Real = 1.0)
    # Parameters of the prior over the filter `h`
    μ₀ = zeros(T, order+1)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, order+1, order+1)

    # Initial parameters of the posterior over the filter `h`
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    # We use `pseudocounts` as the default parameterization of the
    # Gamma prior/posterior
    ARNormal1D(μ₀, Σ₀, μ, Σ, pseudocounts, pseudocounts, pseudocounts,
               pseudocounts)
end

meanfilter(model::ARNormal1D{K}) where K = reverse(mean(model.h.posterior)[1:K])

function (model::ARNormal1D{K})(x::Vector{T}) where {T <: AbstractFloat, K}
    px = PaddedView(0., x, (1-K:length(x), ))
    regressors = Regressors1D{T, K}(px)

    # Expectation of the natural parameters of the model.
    E_λ, E_lnλ = gradlognorm(model.λ.posterior)
    E_Th = gradlognorm(model.h.posterior)

    N = length(px) - K
    llh = zeros(T, N)
    for n in 1:N
        r = vec(regressors[n])
        stats = vcat(r * px[n], -.5 * vec(r * r'))
        llh[n] = -.5 * (log(2π) - E_lnλ + E_λ * px[n]^2) + E_λ * dot(E_Th, stats)
    end
    llh
end

function stats1D_h(model::ARNormal1D{K}, x::Vector{T}) where {T <: AbstractFloat, K}
    px = PaddedView(0., x, (1-K:length(x), ))
    regressors = Regressors1D{T, K}(px)
    E_λ, E_lnλ = gradlognorm(model.λ.posterior)
    retval = zeros(T, (K + 1) * (K + 2), length(x))
    N = length(px) - K
    for n in 1:N
        r = vec(regressors[n])
        retval[:, n] = vcat(E_λ * r * px[n], -.5 * E_λ * vec(r * r'))
    end
    retval
end

function stats1D_λ(model::ARNormal1D{K}, x::Vector{T}) where {T <: AbstractFloat, K}
    px = PaddedView(0., x, (1-K:length(x), ))
    regressors = Regressors1D{T, K}(px)
    E_Th = gradlognorm(model.h.posterior)
    retval = zeros(T, 2, length(x))
    N = length(px) - K
    for n in 1:N
        r = vec(regressors[n])
        stats = vcat(r * px[n], -.5 * vec(r * r'))
        retval[:, n] = [-.5 * px[n]^2 + dot(E_Th, stats), .5]
    end
    retval
end

#######################################################################
# ARNormal: concatenation of D independent ARNormal1D

struct ARNormal{D, K} <: Model
    filters::Vector{ARNormal1D}

    function ARNormal(filters::Vector{ARNormal1D{K}}) where K
        model = new{length(filters), K}(filters)

        # We override the stats
        for (d, filter) in enumerate(filters)
            filter.h.stats = data -> stats_h(model, d, data)
            filter.λ.stats = data -> stats_λ(model, d, data)
        end
        model
    end
end

function ARNormal(μ₀::Vector{T}, Σ₀::Matrix{T}, μ::Vector{T},
                  Σ::Matrix{T}, a₀::Real, b₀::Real, a::Real, b::Real,
                  dim::Integer) where T <: AbstractFloat
    ARNormal([ARNormal1D(μ₀, Σ₀, deepcopy(μ), deepcopy(Σ), a₀, b₀,
                         deepcopy(a), deepcopy(b)) for i in 1:dim])
end

function ARNormal(T = Float64; dim::Integer, order::Integer, initσ::Real = 0.0,
                    pseudocounts::Real = 1.0)
    # Parameters of the prior over the filter `h`
    μ₀ = zeros(T, order+1)
    Σ₀ = T(1/pseudocounts) * Matrix{T}(I, order+1, order+1)

    # Initial parameters of the posterior over the filter `h`
    μ = μ₀ .+ T(initσ) .* randn(T)
    Σ = deepcopy(Σ₀)

    # We use `pseudocounts` as the default parameterization of the
    # Gamma prior/posterior
    ARNormal(μ₀, Σ₀, μ, Σ, pseudocounts, pseudocounts, pseudocounts,
             pseudocounts, dim)
end

getconjugateparams(model::ARNormal) = vcat([m.h for m in model.filters],
                                           [m.λ for m in model.filters])

function (model::ARNormal)(X::Matrix{T}) where T <: AbstractFloat
    retval = zeros(T, size(X, 2))
    for (d, m) in enumerate(model.filters)
        retval .+= m(X[d, :])
    end
    retval
end

function stats_h(model::ARNormal, d::Integer, X::Matrix{T}) where T <: AbstractFloat
    stats1D_h(model.filters[d], X[d, :])
end

function stats_λ(model::ARNormal, d::Integer, X::Matrix{T}) where T <: AbstractFloat
    stats1D_λ(model.filters[d], X[d, :])
end

max_stars_repo_name: fnlandini/BayesianModels

max_issues_repo_name: fnlandini/BayesianModels

max_forks_repo_name: fnlandini/BayesianModels


text:
function parser(julia,jupyter)
    f = open("test\\Test_1.jl","r")
    line = 0
    file = open("writing files\\Test_1.ipynb","w")
    write(file,"{ \"cells\":[ ")
    while !eof(f)
        line+=1
        s = readline(f)
        if (s[1]=='#')
            md = open("templates\\markdown_cell.json","r")
            println("Its a comment!!")
            println("Line Number: $line")
            while !eof(md)
                mdr = readline(md)
                write(file,"$mdr")
            end
            write(file,",")
            close(md)
        else
            cc = open("templates\\code_cell.json","r")
            println("Its a code line")
            println("Line Number: $line")
            while !eof(cc)
                ccr = readline(cc)
                write(file,"$ccr")
            end
            write(file,",")
            close(cc)
        end

    end
    write(file,"],")
    template = open("templates\\metadata.json","r")
    while !eof(template)
        metaR = readline(template)
        write(file,"$metaR")
    end
    write(file,"}")
    close(file)
end


text:
"""
    rewire_edges(g, p)

Rewire each edge of `g` with probability `p`.

When `g` is a cyclic k-lattice, this essentially produces Watts-Strogatz networks.
"""
function rewire_edges(g::T, p::Real) where T <: AbstractGraph
    h = typeof(g)(nv(g))
    for e in edges(g)
        (u,v) = src(e), dst(e)
        # this is somewhat inefficient, but at least it's correct
        rewirings = append!(
            [(u,w) for w in setdiff(vertices(g), outneighbors(h, u)) if w != v],
            [(w,v) for w in setdiff(vertices(g), inneighbors(h, v)) if w != u]
            )
        if !isempty(rewirings) && rand() < p
            (u, v) = rand(rewirings)
        end
        add_edge(h, u, v)
    end
    return h
end

rewire_edges(p::Real) = g -> rewire_edges(g,p)

"""
    drop_edges!(g, p)

Remove each edge in `g` with probability `p`.
"""
function drop_edges!(g, p::Real)
    for e in collect(edges(g))
        if rand() < p
            rem_edge(g,p)
        end
    end
    return g
end

"""
    drop_edges(g, p)

Remove each edge in `g` with probability `p`.

This function creates a copy of `g`. See [drop_edges!](@ref) for the in-place version.
"""
drop_edges(g, p::Real) = drop_edges!(copy(g), p)

drop_edges!(p::Real) = g -> drop_edges!(g,p)
drop_edges(p::Real) = g -> drop_edges!(g,p)


text:
module Solutions

using HDF5
using OffsetArrays
using SharedArrays

using GeometricBase
using GeometricBase.Config
using GeometricBase.Utils
using GeometricEquations

using GeometricBase: fromarray

export DEFAULT_NSAVE, DEFAULT_NWRITE

const DEFAULT_NSAVE = 1
const DEFAULT_NWRITE = 0

export DataSeries, TimeSeries, Solution, DeterministicSolution
export get_solution, get_solution!, set_solution!

include("solutions/solution.jl")

export AtomicSolution,
    AtomicSolutionODE, AtomicSolutionPODE,
    AtomicSolutionDAE, AtomicSolutionPDAE
export update!, cut_periodic_solution!

include("solutions/atomic_solution.jl")
include("solutions/atomic_solution_ode.jl")
include("solutions/atomic_solution_pode.jl")
include("solutions/atomic_solution_dae.jl")
include("solutions/atomic_solution_pdae.jl")

export SolutionIO, SolutionHDF5, hdf5

include("solutions/io.jl")
include("solutions/io_hdf5.jl")

export SolutionODE, SSolutionODE, PSolutionODE, SolutionPODE, SSolutionPODE, PSolutionPODE
export SolutionDAE, SSolutionDAE, PSolutionDAE, SolutionPDAE, SSolutionPDAE, PSolutionPDAE
export get_initial_conditions, get_initial_conditions!, set_initial_conditions!

include("solutions/solution_ode.jl")
include("solutions/solution_pode.jl")
include("solutions/solution_dae.jl")
include("solutions/solution_pdae.jl")
include("solutions/solutions_hdf5.jl")

include("solutions/atomic_solution_constructors.jl")

end


text:
#=
    There are six general cases for processing a/b (when b excludes 0)
    and six special cases (where the denominator includes 0 as a bound)
    and two exceptional cases (where the denominator is 0 or straddles 0)
                  b <= 0       b >= 0
       a <=  0    LteLte       LteGte
       a <0< a    ZerLte       ZerGte
       a >=  0    GteLte       GteGte
=#

function (/){S<:Sculpt,W<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::Flex{W,Q,C})
    z = zero(C)
    if ((b.lo < z) & (b.hi > z))
      divisorContainsZero(a,b)
   elseif a.hi <= z
       if     b.hi <= z
           divLteLte(a,b)
       else # b.lo.fp >= z
           divLteGte(a,b)
       end
   elseif a.lo >= z
       if     b.hi <= z
           divGteLte(a,b)
       else # b.lo >= 0
           divGteGte(a,b)
       end
   else   # a straddles 0
       if     b.hi <= z
           divZerLte(a,b)
       else # b.lo >= 0
           divZerGte(a,b)
       end
   end
end

for (fn,loa,lob,hia,hib) in [ (:divLteLte, :(a.hi), :(b.lo), :(a.lo), :(b.hi)),
                              (:divLteGte, :(a.lo), :(b.lo), :(a.hi), :(b.hi)),
                              (:divGteLte, :(a.hi), :(b.hi), :(a.lo), :(b.lo)),
                              (:divGteGte, :(a.lo), :(b.hi), :(a.hi), :(b.lo)),
                              (:divZerLte, :(a.hi), :(b.hi), :(a.lo), :(b.hi)),
                              (:divZerGte, :(a.lo), :(b.lo), :(a.hi), :(b.lo)),
                            ]
  @eval begin
    function ($fn){S<:Sculpt,W<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::Flex{W,Q,C})
        aLoIsOpen, aHiIsOpen = boundries(S)
        bLoIsOpen, bHiIsOpen = boundries(W)
        cType = boundries( (aLoIsOpen|bLoIsOpen), (aHiIsOpen|bHiIsOpen) )

        lo = (/)(($loa), ($lob), RoundDown)
        hi = (/)(($hia), ($hib), RoundUp)

        Flex{cType,Q,C}(lo, hi)
    end
  end
end

function divisorContainsZero{S<:Sculpt,W<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::Flex{W,Q,C})
    throw(ErrorException("Divisor contains zero: $(a) / $(b)"))
end

(/){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::C) = (/)(a, Flex{S,Q,C}(b))
(/){S<:Sculpt,Q<:Qualia,C<:Clay}(a::C, b::Flex{S,Q,C}) = (/)(Flex{S,Q,C}(a), b)

(/){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::Flex{S,Q,C}, b::T) = (/)(a, Flex{S,Q,C}(convert(C,b)))
(/){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::T, b::Flex{S,Q,C}) = (/)(Flex{S,Q,C}(convert(C,a)), b)

(/){S<:Sculpt,T<:Sculpt,Q<:Qualia,R<:Qualia,C<:Clay,D<:Clay}(a::Flex{S,Q,C}, b::Flex{T,R,D}) = (/)(promote(a,b)...)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec "${JULIA:-julia}" "$@" ${BASH_SOURCE[0]}
=#

using Libdl
using PackageCompiler
using Pkg

Pkg.activate(@__DIR__)
Pkg.develop(PackageSpec(name="PyCall", path=dirname(@__DIR__)))
Pkg.build("PyCall")
Pkg.activate()

sysimage_path = joinpath(@__DIR__, "sys.$(Libdl.dlext)")
create_sysimage(
    [:PyCall],
    sysimage_path = sysimage_path,
    project = @__DIR__(),
    precompile_execution_file = joinpath(@__DIR__, "precompile.jl"),
)

write(joinpath(@__DIR__, "_julia_path"), Base.julia_cmd().exec[1])

@info "System image: $sysimage_path"


text:
include("Hilbert.jl")
include("OffHilbert.jl")
include("LowRankIntegralOperator.jl")
include("HierarchicalDSF.jl")
include("HierarchicalOperator.jl")
include("hierarchicalsolve.jl")


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# ╔═╡ 219d899e-8e54-437e-b20b-840e923ed0c7
begin

	# load -
	using PlutoUI
	using Plots
	using LinearAlgebra
	using CSV
	using DataFrames

	# setup paths -
	_PATH_TO_ROOT = pwd()
	_PATH_TO_DATA = joinpath(_PATH_TO_ROOT, "data")

	# show -
	nothing
end

# ╔═╡ 83a9dd26-8ab0-11ec-05a7-97fba528df76
md"""
### Multiple Substrates and Inhibitor Models

In this lecture, we continue our discussion of enzyme kinetics. Today we will:

* Discuss competitive and non-competitive inhibitor models
* Introduce phenomenological and mechanistic multi-substrate enzyme kinetics models
* Incorporating different levels of information into the bounds constraints
"""

# ╔═╡ b722071c-4256-4fea-84c2-3a209822a1c1
md"""
### Competitive, non- and uncompetitive inhibitors
Inhibitors are critical research and therapeutic tools; they are often a go-to strategy for manipulating biological function (and are widely used in the allosteric regulation context).

* [Clinical application: Literature search on small-molecule checkpoint inhibitors](https://pubmed.ncbi.nlm.nih.gov/?term=small-molecule%20checkpoint%20inhibitors)

Let's revisit the previous feedback inhibition examples, getting more granular about the type of inhibition (and inhibitor) involved. There are three types of inhibitors:

* __Non-competitive inhibition__: A non-competitive inhibitor reduces the enzyme's activity and binds equally well to the enzyme whether or not it has already bound the substrate. This type of inhibition reduces the maximum rate of a chemical reaction without changing the apparent binding affinity of the catalyst for the substrate. No chemistry is possible for the E:I or E:S:I complexes. [Derivation for non-competitive inhibition here can be found here.](https://www.youtube.com/watch?v=pq5TtzVoZdM)
* __Competitive inhibition__: A competitive inhibitor prevents the binding of the substrate to the enzyme. In competitive inhibition, the maximum reaction velocity $V_{max}$ is unchanged. At the same time, the apparent affinity of the substrate to the binding site is decreased (decreasing affinity leads to increased $K_{m}$). [Derivation for competitive inhibition can be found here.](https://www.youtube.com/watch?v=3-mj0Cx0JbA)
* __Uncompetitive inhibition__: An uncompetitive inhibitor binds only to the enzyme-substrate complex. Uncompetitive inhibition typically occurs in reactions with two or more substrates or products. Both $V_{max}$ and $K_{m}$ decreased by the inhibitor. [Derivation for uncompetitive inhibition can be found here.](https://www.youtube.com/watch?v=EyAln4sDxmg)
"""

# ╔═╡ 681d0325-c997-49bd-ad5a-da2e62a4fade
md"""
__Example__: Non-competitive inhibition.
The form for an enzyme catalyzed reaction of a single substrate $S$ to product $P$ by enzyme $E$ in the presense of a non-compeitive inhibitor $I$ is given by:

$$v=V_{max}^{app}\left(\frac{S}{K_{m}+S}\right)$$

where the _appearant_ maximum reaction velocity is given by:

$$V_{max}^{app} = \frac{k_{cat}E}{1+I/K_{I}}$$

where $K_{I}$ denotes an inhibition constant (units: Inhibitor concentration), and $k_{cat}$ denotes the turnover number for the enzyme of interest (units: 1/time).
"""

# ╔═╡ 3e0457b7-9d6a-488a-9ff6-06cb90bde738
let

	# setup some parameters -
	kcat = 13.7 # units s^-1
	Eₒ = 2.0 	# units: *mol/L
	Km = 5.0 	# units: *mol/L
	KI = 20.0 	# units: *mol/L

	# set the substrate and inhibitor concentration, and initialize some space -
	I_array = range(0.0,step = 10,stop = 50.0) |> collect
	S_array = range(0.0,step = 0.01,stop = 100.0) |> collect
	Nₛ = length(S_array);
	Nᵢ = length(I_array)

	# initialize -
	v_array = Array{Float64, 2}(undef, Nₛ, Nᵢ)

	# main loop -
	for (i,S) ∈ enumerate(S_array)
		for (j,Iₒ) ∈ enumerate(I_array)

			# compute the apparent Vmax -
			Vmax_app = (kcat*Eₒ)/(1 + Iₒ/KI)
			sat_term = (S/(Km+S))
			v_array[i,j] = Vmax_app*sat_term
		end
	end

	# plots -
	for (j,Iₒ) ∈ enumerate(I_array)
		if (j == 1)
			plot(S_array,v_array[:,j], c=:blue,lw=2, label="I = $(Iₒ) mmol/L")
		else
			plot!(S_array,v_array[:,j], c=:red,lw=2, label="I = $(Iₒ) mmol/L")
		end
	end

	xlabel!("Substrate S (mM)", fontsize=18)
	ylabel!("v (*mol/s)", fontsize=18)
end

# ╔═╡ 8375496a-c9fb-43a3-84b9-138edf08b92c
md"""
__Example__: Competitive inhibition
The form for an enzyme-catalyzed reaction of a single substrate $S$ to a product $P$ by enzyme $E$ in the presence of a competitive inhibitor $I$ is given by:

$$v=V_{max}\left(\frac{S}{K_{m}^{app}+S}\right)$$

where the _appearant_ saturation constant is given by:

$$K_{m}^{app} = K_{m}\left(1+I/K_{I}\right)$$

where $K_{I}$ denotes an inhibition constant (units: Inhibitor concentration).

"""

# ╔═╡ 80202c43-6dae-4e53-8161-516fbe36a094
let

	# setup some parameters -
	kcat = 13.7 # units s^-1
	Eₒ = 2.0 	# units: *mol/L
	Km = 5.0 	# units: *mol/L
	KI = 20.0 	# units: *mol/L
	n = 10.0

	# set the substrate and inhibitor concentration, and initialize some space -
	I_array = range(0.0,step = 10,stop = 50.0) |> collect
	S_array = range(0.0,step = 0.01,stop = 100.0) |> collect
	Nₛ = length(S_array);
	Nᵢ = length(I_array)

	# initialize -
	v_array = Array{Float64, 2}(undef, Nₛ, Nᵢ)

	# main loop -
	for (i,S) ∈ enumerate(S_array)
		for (j,Iₒ) ∈ enumerate(I_array)

			# compute the apparent Vmax -
			Vmax = (kcat*Eₒ)
			Km_app = Km*(1+Iₒ/KI)
			sat_term = (S^n)/(Km_app^n+S^n)
			v_array[i,j] = Vmax*sat_term
		end
	end

	# plots -
	for (j,Iₒ) ∈ enumerate(I_array)
		if (j == 1)
			plot(S_array,v_array[:,j], c=:blue,lw=2, label="I = $(Iₒ) mmol/L",
				legend=:bottomright)
		else
			plot!(S_array,v_array[:,j], c=:red,lw=2, label="I = $(Iₒ) mmol/L")
		end
	end

	xlabel!("Substrate S (mM)", fontsize=18)
	ylabel!("Rate (*mol/L-time)", fontsize=18)
end

# ╔═╡ b92c2a56-018f-4e27-acd0-bf2488545cba
md"""
##### Discrete state models of inhbitors

###### Theory
Suppose we model the rate $v_{j}$ as the product of a kinetic limit (a simple model of the rate) and a correction term that accounts for the missing regulation:

$$v_{j} = r_{j}\theta\left(...\right)_{j}$$

where $v_{j}$ denotes the overall rate (units: $\mu$M/time), $r_{j}$ denotes the kinetic limit i.e., the maximum rate of conversion (units: $\mu$M/time) and
$0\leq \theta\left(...\right)_{j}\leq 1$ (units: dimensionless) is a control function that describes the influence of effector molecules.

Suppose an enzyme $E$ can exits in one of $s=1,2\dots,\mathcal{S}$ possible microstates, where each microstate $s$ has some pseudo energy $\epsilon_{s}$. Some microstates will lead to activity (the ability to carry out the chemical reactions, while others will not). For each microstate $s$, let's assign a pseudo energy $\epsilon_{s}$, where by definition $\epsilon_{1}=0$; we assume the base state has the lowest energy. Next, suppose the probability that enzyme $E$ is in microstate $s$ follows a [Boltzmann distribution](https://en.wikipedia.org/wiki/Boltzmann_distribution) which says:

$$p_{i} = \frac{1}{Z} \times f_{i}\exp\left(-\beta\epsilon_{i}\right)\qquad{i=1,2,\dots,\mathcal{S}}$$

where $p_{i}$ denotes the probability that enzyme $E$ is in microstate $i=1,2,\dots,\mathcal{S}$, $f_{i}$ denotes a state-specific factor $f_{i}\in\left[0,1\right]$, $\beta$ denotes the [thermodynamic beta](https://en.wikipedia.org/wiki/Thermodynamic_beta) and $Z$ denotes a normalization factor (called the [Partiton function](https://en.wikipedia.org/wiki/Partition_function_(statistical_mechanics)) in the statistical physics community). We can find $Z$ using the summation law of discrete probolity e.g.,  $\sum_{s}p_{s} = 1$ which gives:

$$Z = \sum_{s=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)$$

which gives:

$$p_{i} = \frac{f_{i}\exp\left(-\beta\epsilon_{i}\right)}{\displaystyle \sum_{s=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)}\qquad{i=1,2,\dots,\mathcal{S}}$$.

Finally, we relate the probability that enzyme $E$ is in microstate $s$ back to the $\theta$ control function by computing the overall probability that the desired event happens, e.g., enzyme $E$ catalyzes the reaction of interests. We know if $\Omega = \left\{1,2,\dots,\mathcal{S}\right\}$, then we can define the subset $\mathcal{A}\subseteq\Omega$ in which the desired event happens. Given $\mathcal{A}$, the $\theta$ function becomes:

$$\theta=\sum_{s\in{\mathcal{A}}}p_{s}$$
"""

# ╔═╡ ad1609c2-b63e-410b-b119-2c4b327cd87a
md"""
###### A DSM model for non-competitive inhibition

Enzyme $E$, inhibited by a non-competitive inhibitor $I$, can exist in one of four possible microstates:

* __s = 1__: No substrate $S$ is bound. No chemical reaction is possible.
* __s = 2__: Substrate $S$ is bound to enzyme $E$, but inhibitor $I$ is not bound. A chemical reaction is possible.
* __s = 3__: Both the substrate $S$ and inhibitor $I$ are bound to enzyme $E$. No chemical reaction is possible
* __s = 4__: Inhibitor $I$ is bound to enzyme $E$. No chemical reaction is possible

There are four possible states $\Omega=\left\{1,2,3,4\right\}$, while the subset of states resulting in a chemical reaction $\mathcal{A}\subseteq\Omega$ is given by $\mathcal{A}=\left\{2\right\}$. Thus, the $\theta$ funtion is given by:

$$\theta = p_{2}$$

or:

$$\theta = \frac{f_{2}\exp\left(-\beta\epsilon_{2}\right)}
{\displaystyle \sum_{i=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)}$$

"""

# ╔═╡ 4acfa066-8422-462d-966f-d025abf8d5b2
let

	# setup constants -
	kcat = 13.7 			# units: sec^-1
	Eₒ = 2.0 				# units: *mol/L
	R = 8.314 				# units: J/mol-K
	T = (273.15 + 25.0) 	# units: K
	β = (1/(R*T))			# units: mol/J
	K₂ = 25.0 				# units: mmol/L
	K₃ = 10.0 				# units: mmol/L
	K₄ = 20.0 				# units: mmol/L
	n₂ = 1.0 				# units: dimensionless
	n₃ = n₂					# units: dimensionless
	n₄ = 1.0 	 			# units: dimensionless

	# setup ϵ array -
	ϵ_array = [
		0.0 		; # 1 s₁ units: J/mol
		-5000.0 	; # 2 s₂ units: J/mol
		-1000.0 	; # 3 s₃ units: J/mol
		-1000.0 	; # 4 s₄ units: J/mol
	];

	# compute the W_array -
	W_array = exp.(-β*ϵ_array)

	# set the inhibitor concentration, and initialize some space -
	I_array = range(0.0,step = 10,stop = 50.0) |> collect
	S_array = range(0.0,step = 0.01,stop = 100.0) |> collect
	Nₛ = length(S_array);
	Nᵢ = length(I_array)

	# initialize -
	θ_array = Array{Float64, 2}(undef, Nₛ, Nᵢ)

	# main loop -
	for (i,S) ∈ enumerate(S_array)
		for (j,Iₒ) ∈ enumerate(I_array)

			# compute fᵢ -
			f₄ = ((Iₒ/K₄)^(n₄))/(1 + (Iₒ/K₄)^(n₄))
			f₃ = (Iₒ/K₄)*(((S/K₂)^(n₃))/(1 + (S/K₂)^(n₃)))
			f₂ = (((S/K₂)^(n₂))/(1 + (S/K₂)^(n₂)))
			f₁ = 1

			# populate the f array -
			f_array = [f₁,f₂,f₃,f₄];

			# compute θ -
			D = dot(f_array,W_array)
			N = f₂*W_array[2]
			θ_array[i,j] = kcat*Eₒ*(N/D)
		end
	end

	# plots -
	for (j,Iₒ) ∈ enumerate(I_array)
		if (j == 1)
			plot(S_array,θ_array[:,j], c=:blue,lw=2, label="I = $(Iₒ) mmol/L",
				legend=:bottomright)
		else
			plot!(S_array,θ_array[:,j], c=:red,lw=2, label="I = $(Iₒ) mmol/L")
		end
	end

	xlabel!("Substrate S (mM)", fontsize=18)
	ylabel!("Rate (*mol/L-time)", fontsize=18)
end

# ╔═╡ 7ffdd795-395e-449f-9b4c-cc88208f56f3
md"""
###### A DSM model for competitive inhibition

Enzyme $E$, inhibited by a competitive inhibitor $I$, can exist in one of three possible microstates:

* __s = 1__: No substrate $S$ is bound. No chemical reaction is possible.
* __s = 2__: Substrate $S$ is bound to enzyme $E$, but inhibitor $I$ is not bound. A chemical reaction is possible.
* __s = 3__: Inhibitor $I$ is bound to enzyme $E$. No chemical reaction is possible

There are three possible states $\Omega=\left\{1,2,3\right\}$, while the subset of states resulting in reaction $\mathcal{A}\subseteq\Omega$ is given by $\mathcal{A}=\left\{2\right\}$. Thus, the $\theta$ funtion is given by:

$$\theta = p_{2}$$

or:

$$\theta = \frac{f_{2}\exp\left(-\beta\epsilon_{2}\right)}
{\displaystyle \sum_{i=1}^{\mathcal{S}}f_{i}\exp\left(-\beta\epsilon_{i}\right)}$$
"""

# ╔═╡ 0dbe24c4-224b-4fd8-9f9c-37f6223bd030
let

	# setup constants -
	kcat = 13.7 			# units: sec^-1
	Eₒ = 2.0 				# units: *mol/L
	R = 8.314 				# units: J/mol-K
	T = (273.15 + 25.0) 	# units: K
	β = (1/(R*T))			# units: mol/J
	K₂ = 10.0 				# units: mmol/L
	K₃ = 5.0 				# units: mmol/L
	n₂ = 10.0 				# units: dimensionless
	n₃ = 1.0				# units: dimensionless

	# setup ϵ array -
	ϵ_array = [
		0.0 		; # 1 s₁ units: J/mol
		-4000.0 	; # 2 s₂ units: J/mol
		-100.0 		; # 3 s₃ units: J/mol
	];

	# compute the W_array -
	W_array = exp.(-β*ϵ_array)

	# set the inhibitor concentration, and initialize some space -
	I_array = range(0.0,step = 10,stop = 50.0) |> collect
	S_array = range(0.0,step = 0.01,stop = 100.0) |> collect
	Nₛ = length(S_array);
	Nᵢ = length(I_array)

	# initialize -
	θ_array = Array{Float64, 2}(undef, Nₛ, Nᵢ)

	# main loop -
	for (i,S) ∈ enumerate(S_array)
		for (j,Iₒ) ∈ enumerate(I_array)

			# compute fᵢ -
			f₃ = ((Iₒ/K₃)^(n₃))/(1 + (Iₒ/K₃)^(n₃))
			f₂ = ((S^(n₂))/((K₂*(1+Iₒ/K₃))^(n₂) + S^(n₂)))
			f₁ = 1

			# populate the f array -
			f_array = [f₁,f₂,f₃];

			# compute θ -
			D = dot(f_array,W_array)
			N = f₂*W_array[2]
			θ_array[i,j] = kcat*Eₒ*(N/D)
		end
	end

	# plots -
	for (j,Iₒ) ∈ enumerate(I_array)
		if (j == 1)
			plot(S_array,θ_array[:,j], c=:blue,lw=2, label="I = $(Iₒ) mmol/L",
				legend=:bottomright)
		else
			plot!(S_array,θ_array[:,j], c=:red,lw=2, label="I = $(Iₒ) mmol/L")
		end
	end

	xlabel!("Substrate S (mM)", fontsize=18)
	ylabel!("Rate (*mol/L-time)", fontsize=18)
end

# ╔═╡ 967500aa-664a-4bca-81b6-1ee95bec595a
md"""
### Kinetics of multiple substrate reactions

##### Power-law kinetics and biochemical systems theory (BST)
Power-law kinetics are a flexible tool to describe multiple substrate kinetics. Suppose reaction $v_{i}$ depends upon $j=1,2,\dots,\mathcal{F}$ factors. These factors can be concentration e.g., substrates or products well as other type of data e.g., categorical data:

$$v_{i} = \alpha_{i}\prod_{j=1}^{\mathcal{F}}X_{j}^{f_{ij}}\qquad{i=1,2,\dots,\mathcal{R}}$$

where $\alpha_{j}$ denotes the rate constant for reaction $j$, $X_{j}$ denotes the abundance of
factor $j$ and $f_{ij}$ denotes the kinetic order of factor $j$ in reaction $j$. Power-law kinetics are a prominent feature of Biochemical Systems Theory (BST). Biochemical systems theory has been developed since the 1960s by Michael Savageau, Eberhard Voit, and others for the systems analysis of biochemical processes:

* [Atkinson MR, Savageau MA, Myers JT, Ninfa AJ. Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell. 2003 May 30;113(5):597-607. doi: 10.1016/s0092-8674(03)00346-5. PMID: 12787501.](https://pubmed.ncbi.nlm.nih.gov/12787501/)
* [Alvarez-Vasquez F, Sims KJ, Cowart LA, Okamoto Y, Voit EO, Hannun YA. Simulation and validation of modeled sphingolipid metabolism in Saccharomyces cerevisiae. Nature. 2005 Jan 27;433(7024):425-30. doi: 10.1038/nature03232. PMID: 15674294.](https://pubmed.ncbi.nlm.nih.gov/15674294/)
* [Goel G, Chou IC, Voit EO. Biological systems modeling and analysis: a biomolecular technique of the twenty-first century. J Biomol Tech. 2006;17(4):252-269.](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC2291792/?report=classic)

##### General multisubstrate kinetics
Suppose the irreversible rate $v_{i}$ is dependent upon susbtrates $S_{j},j=1,2,\dots,\mathcal{S}$, then the multiple saturation kinetic form is given by:

$$v_{i} = V_{max,i}\left[\frac{\prod_{j}\frac{S_{j}}{K_{j}}}{\prod_{j}\left(1+\frac{S_{j}}{K_{j}}\right) - 1}\right]\qquad{i=1,2,\dots,\mathcal{R}}$$

where $V_{max,i}$ denote the maximum reaction rate (units: concentration/time), $S_{j}$ denotes the
substrate concentration (units: concentration) and $K_{j}$ denotes the saturation constant for substrate $j$.

* [Liebermeister W, Klipp E. Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theor Biol Med Model. 2006;3:41. Published 2006 Dec 15. doi:10.1186/1742-4682-3-41](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC1781438/)

##### Ping–pong, sequential- and random-order mechanistic models
Multi substrate enzyme kinetics models can be broadly broken down into classes:

* __Ping-pong mechanisms__: In a [ping-pong mechanism](https://chem.libretexts.org/Courses/University_of_Arkansas_Little_Rock/CHEM_4320_5320%3A_Biochemistry_1/05%3A_Michaelis-Menten_Enzyme_Kinetics/5.6%3A_Multi-Substrate_Sequential_Mechanisms), also called a double-displacement reaction, the enzyme is transformed into an intermediate state where the first substrate to product reaction occurs, followed by a second reaction. Enzymes with ping–pong mechanisms include [oxidoreductases](https://en.wikipedia.org/wiki/Oxidoreductase), [transferases](https://en.wikipedia.org/wiki/Transferase), and [serine proteases](https://en.wikipedia.org/wiki/Serine_protease) such as [trypsin](https://en.wikipedia.org/wiki/Trypsin), [chymotrypsin](https://en.wikipedia.org/wiki/Chymotrypsin) and several enzymes of the blood clotting cascade.
* __Random and sequential order mechanisms__: In a [random or sequential order model](https://chem.libretexts.org/Courses/University_of_Arkansas_Little_Rock/CHEM_4320_5320%3A_Biochemistry_1/05%3A_Michaelis-Menten_Enzyme_Kinetics/5.6%3A_Multi-Substrate_Sequential_Mechanisms) both substrates bind to the enzyme before the chemistry can occur. Products can be released in random or sequential order.
"""

# ╔═╡ 3f4ae7c0-f425-49c9-9df7-7424eaca7aa2
md"""
### Kinetics and flux balance analysis calculations

Kinetic rate equations (and the associated regulation of enzyme activity, enzyme state, and overall biochemical state) are included in the flux balance analysis calculation through the flux bounds equation. There are different levels of information included in the bounds. These levels are partitioned (broadly) into metabolite and metabolite-free groups.

##### General bounds model
Flux bounds constrain the values that each reaction in a metabolic network can take. A general model for these bounds is given by:

$$-\delta_{j}\left[{V_{max,j}^{\circ}}\left(\frac{e_{j}}{e^{\circ}}\right)\theta_{j}\left(\dots\right){f_{j}\left(\dots\right)}\right]\leq{v_{j}}\leq\left[{V_{max,j}^{\circ}}\left(\frac{e_{j}}{e^{\circ}}\right)\theta_{j}\left(\dots\right){f_{j}\left(\dots\right)}\right]$$

where $V_{max,j}^{\circ}$ denotes the maximum reaction velocity (units: flux) computed at some characteristic enzyme abundance (units: concentration), the ratio $e/e^{\circ}$ is a correction for enzyme abundance (units: dimensioness), $\theta_{j}\left(\dots\right)\in\left[0,1\right]$ is the fraction of maximial enzyme activity (a function or measurement producing units: dimensionless), and $f_{j}\left(\dots\right)$ is a function describing the substrate dependence of the reaction rate $j$ (units: dimensionless). Both $\theta_{j}\left(\dots\right)$ and $f_{j}\left(\dots\right)$ could have associated parameters, e.g., saturation or binding constants, etc.  Finally, the quanity $\delta_{j}\in\left\{0,1\right\}$ is a _binary_ variable:
* If reaction $j$ is __reversible__ $\delta_{j}=1$ or,
* If reaction $j$ is __irreversible__ $\delta_{j}=0$

##### Metabolite free bounds models
In metabolite-free models, we, as the simulation engineer, do not have access to any information about the concentrations of the substrates that appear in the various form of the kinetics. In these cases, $f_{j}\left(\dots\right)=1,\forall{j}$ and $\theta_{j}=1,\forall{j}$ which gives bounds of the form:

$$-\delta_{j}\left[{V_{max,j}^{\circ}}\left(\frac{e_{j}}{e^{\circ}}\right)\right]\leq{v_{j}}\leq\left[{V_{max,j}^{\circ}}\left(\frac{e_{j}}{e^{\circ}}\right)\right]$$

Finally, if we do not have access to enzyme abundance data, then $\left(e/e^{\circ}\right)=1$ and $V_{max,j}^{\circ}$ is calculated using a _characteristic_ enzyme abundance.

##### Metabolite-dependent bounds models
In this class of bounds, we assume that we have some (or full) information about the concentrations of the metabolites. Assuming general multi-substrate kinetics, the bounds would take the form (shown for irreversible):

$$0\leq{v_{j}}\leq{V_{max,i}^{\circ}\left(\frac{e_{j}}{e^{\circ}}\right)\theta_{j}\left(\cdots\right)\left[\frac{\prod_{k}\frac{S_{j}}{K_{j}}}{\prod_{k}\left(1+\frac{S_{j}}{K_{j}}\right) - 1}\right]}$$

where the products are carried out over the reactants of the reaction $j$. Note: we could have chosen a different form for the kinetics, e.g., a power-law formulation, and the choice of enzyme activity model is up to the simulation designer.

##### Sources of data for bounds models

* [Chang A, Jeske L, Ulbrich S, Hofmann J, Koblitz J, Schomburg I, Neumann-Schaal M, Jahn D, Schomburg D. BRENDA, the ELIXIR core data resource in 2021: new developments and updates. Nucleic Acids Res. 2021 Jan 8;49(D1): D498-D508. doi: 10.1093/nar/gkaa1025. PMID: 33211880; PMCID: PMC7779020.](https://pubmed.ncbi.nlm.nih.gov/33211880/)
* [Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers--the database of key numbers in molecular and cell biology. Nucleic Acids Res. 2010 Jan;38(Database issue): D750-3. doi: 10.1093/nar/gkp889. Epub 2009 Oct 23. PMID: 19854939; PMCID: PMC2808940.](https://pubmed.ncbi.nlm.nih.gov/19854939/)
* [Park JO, Rubin SA, Xu YF, Amador-Noguez D, Fan J, Shlomi T, Rabinowitz JD. Metabolite concentrations, fluxes and free energies imply efficient enzyme usage. Nat Chem Biol. 2016 Jul;12(7):482-9. doi: 10.1038/nchembio.2077. Epub 2016 May 2. PMID: 27159581; PMCID: PMC4912430.](https://pubmed.ncbi.nlm.nih.gov/27159581/)

"""

# ╔═╡ f7a88d98-7149-485f-ade1-b1d39f28daaa
md"""
### Summary and Conclusions

In this lecture we:

* Discussed competitive and non-competitive inhibitor models
* Introduced phenomenological and mechanistic multi-substrate enzyme kinetics models
* Incorporated different levels of information into the bounds constraints
"""

# ╔═╡ 5391cecc-bbd5-44e4-9a92-2184b4252f20
md"""
### Next Time

In-class problem set in which we estimate the flux through the Urea cycle with different types of bounds constraints.
"""

# ╔═╡ 518e00f4-16a2-4ca1-a365-d9ef17437d97
TableOfContents(title="📚 Table of Contents", indent=true, depth=5, aside=true)

# ╔═╡ c375b63a-c383-4476-84fc-4b3078de8c64
html"""
<script>
	// initialize -
	var section = 0;
	var subsection = 0;
	var subsubsection = 0;
	var headers = document.querySelectorAll('h3, h5, h6');

	// main loop -
	for (var i=0; i < headers.length; i++) {

		var header = headers[i];
	    var text = header.innerText;
	    var original = header.getAttribute("text-original");
	    if (original === null) {

			// Save original header text
	        header.setAttribute("text-original", text);
	    } else {

			// Replace with original text before adding section number
	        text = header.getAttribute("text-original");
	    }

	    var numbering = "";
	    switch (header.tagName) {
	        case 'H3':
	            section += 1;
	            numbering = section + ".";
	            subsection = 0;
	            break;
	        case 'H5':
	            subsection += 1;
	            numbering = section + "." + subsection;
	            break;
			case 'H6':
	            subsubsection += 1;
	            numbering = section + "." + subsection + "." + subsubsection;
	            break;
	    }
		// update the header text
		header.innerText = numbering + " " + text;
	};
</script>"""

# ╔═╡ 9780c35f-f044-410c-a16d-a38378c6feed
html"""
<style>
main {
    max-width: 860px;
    width: 70%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
CSV = "~0.10.2"
DataFrames = "~1.3.2"
Plots = "~1.25.9"
PlutoUI = "~0.7.34"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "9519274b50500b8029973d241d32cfbf0b127d97"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.2"

[[deps.Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f9982ef575e19b0e5c7a98c6e75ee496c0f73a93"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.12.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ae13fcbc7ab8f16b0856729b050ef0c446aa3492"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.4.4+0"

[[deps.FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[deps.FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[deps.FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "51d2dfe8e590fbd74e7a842cf6d13d8a2f45dc01"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.6+0"

[[deps.GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "RelocatableFolders", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "4a740db447aae0fbeb3ee730de1afbb14ac798a1"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.63.1"

[[deps.GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "aa22e1ee9e722f1da183eb33370df4c1aeb6c2cd"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.63.1+0"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[deps.Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[deps.Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+2"

[[deps.Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[deps.Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.2+0"

[[deps.LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[deps.LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a8f4f279b6fa3c3c4f1adadd78a621b13a506bce"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.9"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[deps.Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[deps.Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[deps.Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[deps.Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[deps.Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[deps.Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.7"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+1"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "648107615c15d4e09f7eca16307bc821c1f718d8"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.13+0"

[[deps.Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[deps.PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "6f1b25e8ea06279b5689263cc538f51331d7ca17"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.1.3"

[[deps.Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"]
git-tree-sha1 = "1d0a11654dbde41dc437d6733b68ce4b28fbe866"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.25.9"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "8979e9802b4ac3d58c503a20f2824ad67f9074dd"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.34"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "37c1631cb3cc36a535105e6d5557864c82cd8c2b"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.5.0"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.RelocatableFolders]]
deps = ["SHA", "Scratch"]
git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca"
uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
version = "0.1.3"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "a635a9333989a094bddc9f940c04c549cd66afcf"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.3.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "d88665adc9bcf45903013af0982e2fd05ae3d0a6"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "118e8411d506d583fbbcf4f3a0e3c5a9e83370b8"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.15"

[[deps.StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "d21f2c564b21a202f4677c0fba5b5ee431058544"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.4"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[deps.Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

[[deps.Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[deps.Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "66d72dc6fcc86352f01676e8f0f698562e60510f"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.23.0+0"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[deps.XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[deps.Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[deps.Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[deps.Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[deps.Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[deps.Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[deps.Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[deps.Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[deps.Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[deps.Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[deps.Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[deps.Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[deps.Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[deps.Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[deps.Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[deps.Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[deps.Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[deps.Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[deps.Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.2+0"

[[deps.libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[deps.libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[deps.libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+1"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[deps.x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[deps.x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[deps.xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─83a9dd26-8ab0-11ec-05a7-97fba528df76
# ╟─b722071c-4256-4fea-84c2-3a209822a1c1
# ╟─681d0325-c997-49bd-ad5a-da2e62a4fade
# ╟─3e0457b7-9d6a-488a-9ff6-06cb90bde738
# ╟─8375496a-c9fb-43a3-84b9-138edf08b92c
# ╠═80202c43-6dae-4e53-8161-516fbe36a094
# ╟─b92c2a56-018f-4e27-acd0-bf2488545cba
# ╟─ad1609c2-b63e-410b-b119-2c4b327cd87a
# ╠═4acfa066-8422-462d-966f-d025abf8d5b2
# ╟─7ffdd795-395e-449f-9b4c-cc88208f56f3
# ╠═0dbe24c4-224b-4fd8-9f9c-37f6223bd030
# ╟─967500aa-664a-4bca-81b6-1ee95bec595a
# ╟─3f4ae7c0-f425-49c9-9df7-7424eaca7aa2
# ╟─f7a88d98-7149-485f-ade1-b1d39f28daaa
# ╟─5391cecc-bbd5-44e4-9a92-2184b4252f20
# ╠═518e00f4-16a2-4ca1-a365-d9ef17437d97
# ╠═219d899e-8e54-437e-b20b-840e923ed0c7
# ╟─c375b63a-c383-4476-84fc-4b3078de8c64
# ╟─9780c35f-f044-410c-a16d-a38378c6feed
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Libdl

mods = ["MODULES.f90"];

routines = ["CUMULATE.f90", "DRIVE.f90", "FSM.f90", "INITIALIZE.f90", "OUTPUT.f90", "PHYSICS.f90", "QSAT.f90", "SET_PARAMETERS.f90", "SNOW.f90", "SOIL.f90", "SURF_EBAL.f90", "SURF_EXCH.f90", "SURF_PROPS.f90", "TRIDIAG.f90"];

if Sys.iswindows()
    flags = ["-m$(Sys.WORD_SIZE)", "-shared", "-O3"]
else
    flags = ["-m$(Sys.WORD_SIZE)", "-shared", "-O3", "-fPIC"]
end

run(`gfortran $flags $mods $routines -o FSM.$(Libdl.dlext)`)

rm.(filter!(f -> endswith(f, "mod"), readdir()))