text:
# test Petsc matrix functions
facts("\n   ---testing matrix functions---") do
  A = PetscMat(comm)
  PetscMatSetType(A, "mpiaij")

  PetscMatSetSizes(A,sys_size,sys_size, PetscInt(comm_size*sys_size),PetscInt(comm_size*sys_size));
  PetscSetUp(A);

  println("mat_type = ", MatGetType(A))

  B = PetscMat(comm)
  PetscMatSetType(B, "mpiaij")

  PetscMatSetSizes(B,sys_size,sys_size, PetscInt(comm_size*sys_size),PetscInt(comm_size*sys_size));
  PetscSetUp(B);

  x = PetscVec(comm);
  PetscVecSetType(x, VECMPI);
  PetscVecSetSizes(x,sys_size, PetscInt(comm_size*sys_size));

  low, high = PetscVecGetOwnershipRange(x)
  global_indices = Array(low:PetscInt(high - 1))

  for i=1:sys_size
    idxm = [global_indices[i]]   # index
    val = [ rhs[i] ]  # value
    PetscVecSetValues(x, idxm, val, PETSC_INSERT_VALUES)
  end

  PetscVecAssemblyBegin(x)
  PetscVecAssemblyEnd(x)

  xvec_julia = deepcopy(rhs)

  y = PetscVec(comm);
  PetscVecSetType(y, VECMPI);
  PetscVecSetSizes(y,sys_size, PetscInt(comm_size*sys_size));

  for i=1:sys_size
    idxm = [global_indices[i]]   # index
    val = [ rhs[i] ]  # value
    PetscVecSetValues(y, idxm, val, PETSC_INSERT_VALUES)
  end

  PetscVecAssemblyBegin(y)
  PetscVecAssemblyEnd(y)

  y_julia = deepcopy(rhs)

  low, high = PetscMatGetOwnershipRange(A)
  global_indices = Array(low:PetscInt(high - 1))
  println("comm_rank = ", comm_rank, " , global_indices = ", global_indices)

  global_row_ind = zeros(PetscInt, sys_size*sys_size)
  global_col_ind = zeros(PetscInt, sys_size*sys_size)

  # get global row and column indicies, in column major order
  pos = 1
  for i=1:sys_size  # loop over columns
    for j=1:sys_size  # loop over rows
      global_row_ind[pos] = j - 1 + low
      global_col_ind[pos] = i - 1 + low
      pos += 1
    end
  end

  println("global_row_ind = ", global_row_ind)
  println("global_col_ind = ", global_col_ind)
  B_julia = A_julia + PetscScalar(1)

  # test set_values1!
  global_indices1 = global_indices + PetscInt(1)
  global_indices2 = global_indices + PetscInt(1)
  vals = rand(PetscScalar, sys_size, sys_size)

  set_values1!(A, global_indices1, global_indices2, vals)
  PetscMatAssemblyBegin(A)
  PetscMatAssemblyEnd(A)
  vals2 = zeros(vals)
  PetscMatGetValues(A, global_indices, global_indices, vals2)

  @fact vals --> roughly(vals2, atol=1e-13)
  fill!(vals2, 0.0)
  get_values1!(A, global_indices1, global_indices2, vals2)
  @fact vals --> roughly(vals2, atol=1e-13)

  A2s = zeros(3, 3)
  idxm = collect(PetscInt, 1:3)
  idxn = collect(PetscInt, 1:3)
  vals = rand(3,3)
  set_values1!(A2s, idxm, idxn, vals)

  @fact A2s[idxm,  idxn] --> roughly(vals, atol=1e-13)

  set_values1!(A2s, idxm, idxn, vals, PETSC_ADD_VALUES)

  @fact A2s[idxm, idxn] --> roughly(2*vals, atol=1e-13)
  vals2 = zeros(vals)
  get_values1!(A2s, idxm, idxn, vals2)
  @fact vals2 --> roughly(2*vals, atol=1e-13)

  # test the sparse optimized version of +=
  A2ss = sparse(A2s)
  set_values1!(A2ss, idxm, idxn, vals, PETSC_ADD_VALUES)
  fill!(vals2, 0.0)
  get_values1!(A2ss, idxm, idxn, vals2)
  @fact vals2 --> roughly(3*vals, atol=1e-13)

  for i=1:sys_size
    for j = 1:sys_size
      idxm = [ global_indices[i] ] # row index
      idxn = [ global_indices[j] ] # column index
      PetscMatSetValues(A,idxm, idxn, [A_julia[i,j]],PETSC_INSERT_VALUES);
      PetscMatSetValues(B,idxm, idxn, [A_julia[i,j] + PetscScalar(1)],PETSC_INSERT_VALUES);
    end
  end

  B_copy = zeros(PetscScalar, sys_size, sys_size)

  PetscMatAssemblyBegin(A,PETSC_MAT_FINAL_ASSEMBLY);
  PetscMatAssemblyEnd(A,PETSC_MAT_FINAL_ASSEMBLY);

  PetscMatAssemblyBegin(B,PETSC_MAT_FINAL_ASSEMBLY);
  PetscMatAssemblyEnd(B,PETSC_MAT_FINAL_ASSEMBLY);

  # test CreateMatTransposed
  At = MatCreateTranspose(A)

  b2t = PetscVec(comm);
  PetscVecSetType(b2t, VECMPI);
  PetscVecSetSizes(b2t, sys_size, PetscInt(comm_size*sys_size));

  x2 = PetscVec(comm);
  PetscVecSetType(x2, VECMPI);
  PetscVecSetSizes(x2, sys_size, PetscInt(comm_size*sys_size));

  PetscVecSetValues(x2, global_indices, rhs)

  PetscVecAssemblyBegin(x2)
  PetscVecAssemblyEnd(x2)

  PetscMatMult(At, x2, b2t)

  b2t_julia = A_julia.'*rhs

  vals = zeros(PetscScalar, sys_size)
  PetscVecGetValues(b2t, global_indices, vals)

  for i=1:sys_size
    @fact b2t_julia[i] --> roughly(vals[i], atol=1e-12)
  end

  PetscDestroy(At)
  PetscDestroy(b2t)
  PetscDestroy(x2)

  # test MatTranpose
  println("testing out of place transpose")
  At = MatTranspose(A, inplace=false)
  vals = zeros(PetscScalar, sys_size, sys_size)
  valst = zeros(vals)
  PetscMatGetValues(A, global_indices, global_indices, vals)
  PetscMatGetValues(At, global_indices, global_indices, valst)

  for i=1:sys_size
    for j=1:sys_size
      @fact vals[i, j] --> roughly(valst[j, i], atol=1e-12)
    end
  end

  println("testing in place transpose")
  At = MatTranspose(At, inplace=true)
  fill!(vals, 0.0)
  fill!(valst, 0.0)
  PetscMatGetValues(A, global_indices, global_indices, vals)
  PetscMatGetValues(At, global_indices, global_indices, valst)

  println("vals = \n", vals)
  println("valst = \n", valst)

  for i=1:sys_size
    for j=1:sys_size
      @fact vals[i, j] --> roughly(valst[i, j], atol=1e-12)
    end
  end

  PetscDestroy(At)

#=
  for i=1:sys_size
    for j=1:sys_size
      idxm = [global_indices[i] ]  # row index
      idxn = [ global_indices[j] ] # column index
      v = zeros(PetscScalar, 1,1)
      PetscMatGetValues(A, idxm, idxn, v)
#      println("i = ", i, " j = ", j, " v = ", v[1,1], " A[i,j] = ", A_julia[i,j])
      @fact v[1,1] => roughly(A_julia[i,j]) "mismatch at i=$i, j=$j"
      @fact A[i,j] => roughly(A_julia[i,j])
    end
  end

  for i=1:sys_size*sys_size
    A[i] => roughly(A_julia[i])
  end
=#
  @fact size(A) => (sys_size, sys_size)
  @fact size(A, 1) => sys_size
  @fact size(A, 2) => sys_size

  @fact PetscMatGetSize(A) => (comm_size*sys_size, comm_size*sys_size)
  println("Mat local size = ", PetscMatGetLocalSize)
  @fact PetscMatGetLocalSize(A) => (sys_size, sys_size)
  # testing non zero exist code is all we can really do here
  println("Printing A")
  @fact PetscView(A, 0) => 0
  println("Printing B")
  PetscView(B, 0)

  alpha = PetscScalar(2.3)

  PetscMatAXPY(B, alpha, A, DIFFERENT_NONZERO_PATTERN)
  B_julia = alpha*A_julia + B_julia
  PetscMatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size
    for j=1:sys_size
      @fact B_julia[j, i] => roughly(B_copy[i, j])
    end
  end

  println("finished testing MatAXPY")

  PetscMatAYPX(B, alpha, A, DIFFERENT_NONZERO_PATTERN)
  B_julia = alpha*B_julia + A_julia
  PetscMatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size
    for j=1:sys_size
      @fact B_julia[j, i] => roughly(B_copy[i, j])
    end
  end

  println("finished testing MatAYPX")

  PetscMatScale(B, alpha)
  B_julia = alpha*B_julia
  PetscMatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size
    for j=1:sys_size
      @fact B_julia[j, i] => roughly(B_copy[i, j])
    end
  end

  fac = 2.0
  C_julia = rand(3,3)
  C_julia2 = fac*C_julia
  PetscMatScale(C_julia, 2.0)
  @fact C_julia => roughly(C_julia2)

  C_julia = sprand(10, 10, 0.1)
  C_julia2 = fac*C_julia
  PetscMatScale(C_julia, fac)
  @fact C_julia => roughly(C_julia2)

  println("finished testing MatScale")

  PetscMatShift(B, alpha)
  B_julia += alpha*eye(PetscScalar, sys_size)
  PetscMatGetValues(B, global_indices, global_indices, B_copy)
  for i=1:sys_size
    for j=1:sys_size
      @fact B_julia[j, i] => roughly(B_copy[i, j])
    end
  end

  println("finished testing PetscMatShift")

  y_copy = zeros(PetscScalar, sys_size)
  PetscMatMult(B, x, y)
  y_julia = B_julia*xvec_julia

  PetscVecGetValues(y, sys_size, global_indices, y_copy)
  for i=1:sys_size
      @fact y_julia[i] => roughly(y_copy[i])
  end

  PetscMatMultAdd(B, x, y, y)
  y_julia = y_julia + B_julia*xvec_julia
  PetscVecGetValues(y, sys_size, global_indices, y_copy)
  for i=1:sys_size
      @fact y_julia[i] => roughly(y_copy[i])
  end

  PetscMatMultTranspose(A, x, y)
  y_julia = A_julia.'*xvec_julia
  PetscVecGetValues(y, sys_size, global_indices, y_copy)
  for i=1:sys_size
      @fact y_julia[i] => roughly(y_copy[i])
  end

  PetscMatMultHermitianTranspose(A, x, y)
  y_julia = A_julia'*xvec_julia
  PetscVecGetValues(y, sys_size, global_indices, y_copy)
  for i=1:sys_size
      @fact y_julia[i] => roughly(y_copy[i])
  end

  D = PetscMat(comm)
  PetscMatMatMult(A, B, MAT_INITIAL_MATRIX, PetscReal(1.0), D)
  D_julia = A_julia*B_julia
  println("D_julia = ", D_julia)
  D_copy = zeros(PetscScalar, sys_size, sys_size)
  PetscMatGetValues(D, global_indices, global_indices, D_copy)

  for i=1:sys_size
    for j=1:sys_size
      @fact D_copy[j,i] => roughly(D_julia[i,j])
    end
  end

  # test matrix norms
  fnorm = PetscMatNorm(D, NORM_FROBENIUS)
  infnorm = PetscMatNorm(D, NORM_INFINITY)

  @fact fnorm => roughly(sqrt(comm_size)*vecnorm(D_julia), atol=1e-2)
  @fact infnorm => roughly(norm(D_julia, Inf))

  println("testing preallocation")
  nb = PetscInt(100)  # number of times/blocks to insert
  C = PetscMat(comm)
  PetscMatSetType(C, "mpiaij")
#  PetscMatSetFromOptions(C)
  println("nb = ", nb)
  println("comm_size = ", comm_size)
  println("sys_size = ", sys_size)

  PetscMatSetSizes(C, nb*sys_size, nb*sys_size, PetscInt(nb*comm_size*sys_size),PetscInt(nb*comm_size*sys_size));

  # preallocation parameters
  bs = PetscInt(1)
  dnnz = 3*ones(PetscInt, nb*sys_size)  # on diagonal (row + column owned by this process)
  onnz = zeros(PetscInt, nb*sys_size)  # no off diagonal (column not owned by this process)
  dnnzu = Array(PetscInt, 0)  # this is not a symmetric matrix, so unused
  onnzu = Array(PetscInt, 0)  # this is not a symmetric matrix, so unused

  println("dnnz = ", dnnz)
  PetscMatXAIJSetPreallocation(C, bs, dnnz, onnz, dnnzu, onnzu)

#  PetscSetUp(C);

  println("A_julia = ", A_julia)
  println("typeof(A_julia) = ", typeof(A_julia))

  A_julia_t = A_julia.'  # transpose because C is row major

  low, high = PetscMatGetOwnershipRange(C)
  idi = zeros(PetscInt, sys_size)  # row indices
  idj = zeros(PetscInt, sys_size)  # column indices

  for i=1:nb

    # get indices
    pos = 1
    for j=1:sys_size # insert block on diagonal
	idi[pos] = low + sys_size*(i - 1) + j - 1
	idj[pos] = low + sys_size*(i - 1) + j - 1
	pos += 1
    end

    # add a random component
    for i=1:sys_size
      for j=1:sys_size
	A_julia_t[i,j] += rand()
      end
    end

#    println("idi = ", idi)
#    println("idj = ", idj)
    PetscMatSetValues(C, idi, idj, A_julia_t, PETSC_INSERT_VALUES)
  end

  PetscMatAssemblyBegin(C, PETSC_MAT_FINAL_ASSEMBLY)
  PetscMatAssemblyEnd(C, PETSC_MAT_FINAL_ASSEMBLY)

  matinfo = PetscMatGetInfo(C, MAT_LOCAL)

  @fact matinfo.mallocs => roughly(0.0)

#  PetscView(C, 0)

  println("finished testing preallocation")

  PetscMatZeroEntries(C)
  cnorm = PetscMatNorm(C, NORM_FROBENIUS)
  @fact cnorm => roughly(0.0, atol=1e-3)  # norm is zero iff matrix is zero

  a = rand(3,3)
  PetscMatZeroEntries(a)
  @fact norm(a) => 0.0

  a = sprand(10, 10, 0.1)
  PetscMatZeroEntries(a)
  @fact sum(a.nzval) => 0.0

  @fact PetscDestroy(A) => 0
  PetscDestroy(B)
  PetscDestroy(C)
  PetscDestroy(D)
  PetscDestroy(y)
  PetscDestroy(x)
end


text:
using Test

## Load in source code
include("../src/Lunettes.jl")
using .Lunettes
@testset "All Tests" begin

@testset "Composition" begin

# set up domain model with structs
@lens struct Leaf
    color::String
end

@lens struct Twig
    leaves::Leaf
end

@lens struct Branch
    twigs::Twig
end
@lens struct Trunk
    branches::Branch
end
@lens struct Tree
    species::String
    trunk::Trunk
end

# build tree
my_tree = Tree("Pohutakawa", Trunk(Branch(Twig(Leaf("Light Green")))))
another_tree = Tree("Pohutakawa", Trunk(Branch(Twig(Leaf("Light Green")))))

@test my_tree === another_tree

# build big Lens
leaf_one = Lens{Tree,:trunk}() ∘
    Lens{Trunk,:branches}() ∘
    Lens{Branch,:twigs}() ∘
    Lens{Twig,:leaves}() ∘
    Lens{Leaf,:color}()

@test getr(leaf_one)(my_tree) == "Light Green"
@test getr(leaf_one, my_tree) == "Light Green"

@test setr(leaf_one, my_tree, "Light Blue") == Tree("Pohutakawa", Trunk(Branch(Twig(Leaf("Light Blue")))))
@test setr(leaf_one, "Light Blue")(my_tree) == Tree("Pohutakawa", Trunk(Branch(Twig(Leaf("Light Blue")))))

end

@testset "Errors" begin

end

@testset "Macro on Parametric types" begin

@lens struct MyTestStruct{T}
   my_field::T
end

@lens struct AnotherTestStruct{T <: Number}
   my_field::T
end

end

end

ext: jl


text:
"""
  (S::NLFunctional)(integrator) -> (z, η, iter, fail_convergence)

Perform functional iteration that is used by implicit methods, where `z` is the
solution, `η` is used to measure the iteration error (see [^HW96]), `iter` is
the number of iteration, and `fail_convergence` reports whether the algorithm
succeed.  It solves

```math
G(z) = dt⋅f(tmp + γ⋅z, p, t+c⋅h)
z = G(z)
```

by iterating

```math
zᵏ⁺¹ = G(zᵏ).
```

[^HW96]: Ernst Hairer and Gerhard Wanner, "Solving Ordinary Differential
Equations II, Springer Series in Computational Mathematics. ISBN
978-3-642-05221-7. Section IV.8.
[doi:10.1007/978-3-642-05221-7](https://doi.org/10.1007/978-3-642-05221-7)
"""
function (S::NLFunctional{false})(integrator)
  nlcache = S.cache
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,tmp,κ,tol,c,γ,min_iter,max_iter = nlcache
  mass_matrix = integrator.f.mass_matrix
  if typeof(integrator.f) <: SplitSDEFunction
    f = integrator.f.f1
  else
    f = integrator.f
  end
  # precalculations
  κtol = κ*tol

  # initial step of functional iteration
  iter = 1
  tstep = t + c*dt
  u = tmp + γ*z
  z₊ = dt*f(u, p, tstep)
  dz = z₊ - z
  ndz = integrator.opts.internalnorm(dz)
  z = z₊

  η = nlcache.ηold
  do_functional = true # TODO: this makes `min_iter` ≥ 2

  # functional iteration
  fail_convergence = false
  while (do_functional || iter < min_iter) && iter < max_iter
    iter += 1
    u = tmp + γ*z
    z₊ = dt*f(u, p, tstep)
    dz = z₊ - z
    ndzprev = ndz
    ndz = integrator.opts.internalnorm(dz)
    θ = ndz/ndzprev
    if θ > 1 || ndz*(θ^(max_iter - iter)/(1-θ)) > κtol
      fail_convergence = true
      break
    end
    η = θ/(1-θ)
    do_functional = (η*ndz > κtol)
    z = z₊
  end

  if (iter >= max_iter && do_functional) || fail_convergence
    integrator.force_stepfail = true
    return (z, η, iter, true)
  end

  return (z, η, iter, false)
end

function (S::NLFunctional{true})(integrator)
  nlcache = S.cache
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,z₊,b,dz,tmp,κ,tol,k,c,γ,min_iter,max_iter = nlcache
  ztmp = b
  mass_matrix = integrator.f.mass_matrix
  if typeof(integrator.f) <: SplitSDEFunction
    f = integrator.f.f1
  else
    f = integrator.f
  end
  # precalculations
  κtol = κ*tol
  # initial step of functional iteration
  iter = 1
  tstep = t + c*dt
  @. u = tmp + γ*z
  f(k, u, p, tstep)
  if mass_matrix == I
    @. z₊ = dt*k
  else
    @. ztmp = dt*k
    mul!(z₊, mass_matrix, ztmp)
  end
  @. dz = z₊ - z
  @. z = z₊
  ndz = integrator.opts.internalnorm(dz)

  η = nlcache.ηold
  do_functional = true # TODO: this makes `min_iter` ≥ 2

  # NLFunctional iteration
  fail_convergence = false
  while (do_functional || iter < min_iter) && iter < max_iter
    iter += 1
    @. u = tmp + γ*z
    f(k, u, p, tstep)
    if mass_matrix == I
      @. z₊ = dt*k
    else
      @. ztmp = dt*k
      mul!(z₊, mass_matrix, ztmp)
    end
    @. dz = z₊ - z
    ndzprev = ndz
    ndz = integrator.opts.internalnorm(dz)
    θ = ndz/ndzprev
    if θ > 1 || ndz*(θ^(max_iter - iter)/(1-θ)) > κtol
      fail_convergence = true
      break
    end
    η = θ/(1-θ)
    do_functional = (η*ndz > κtol)
    @. z = z₊
  end

  if (iter >= max_iter && do_functional) || fail_convergence
    integrator.force_stepfail = true
    return (z, η, iter, true)
  end

  return (z, η, iter, false)
end


text:
using Distributions, InventoryModels
@testset "inventory_system.jl" begin
    item1 = Item(Inventory(1, 11), Supplier(100,1), Supplier(10,1, leadtime = LeadTime(1, 9)), name = "item1")
    item2 = Item(Inventory(1, 12), Supplier(100,1), name = "item2")
    item3 = Item(Inventory(1, 13), Supplier(100,1), name = "item3")
    item4 = Item(Inventory(10,14), Assembly(80,8, item1 => 2, item2 => 4), name = "item4")
    item5 = Item(Inventory(10,15), Assembly(90,2, item2 => 4, item3 =>1), name = "item5")
    item6 = EndProduct(Market(5, Normal, 4, 0, 10, 0.0), Inventory(100,52), Assembly(1280, 1, item4=>1, item5=>2), name = "item6")

    bom = [item3, item2, item5, item4, item6, item1]
    bom_s = InventoryModels.topological_order(bom)
    for (i, el) in enumerate(bom_s)
        for child in InventoryModels.children(el)
            @test child in @view bom_s[i:end]
        end
    end
    is = InventorySystem(52, bom)
    for (i, el) in enumerate(is.bom)
        for child in InventoryModels.children(el)
            @test child in @view is.bom[i:end]
        end
    end
    @test state_size(is) == length(state(is)) == 2+1+1+1+1+1+1+8
    @test state(is) == [0,10,0,10,0,10,0,10,0,52,15,13,14,12,11,9]
    @test action_size(is) == 6*2 + 2

    is2 = sl_sip(1,10,100,0.,0,[20,40,60,40,0,0,0,0], 0.0, horizon = 4)
    @test state(is2) == [0,20,40,60,40,0]
    is2([14,70])
    @test reward(is2) == -(100 + 0 + 50*1)
    @test state(is2) == [0,40,60,40,0,50]
    is2([29,141])
    @test reward(is2) == -10
    is2([58,114])
    @test reward(is2) == -(100 + 114-60)
    is2([28,53])
    @test reward(is2) == -(114-60-40)
    @test state(is2) == [0,0,0,0,0, 14]
    reset!(is2)
    @test state(is2) == [0,20,40,60,40,0]
    @test test_policy(sl_sip(1,10,100,0.25,0,[20,40,60,40], 0.0), [14,70,29,141,58,114,28,53], 10000) ≈ -363 atol = 2
end


text:
## ND LINEAR INTERPOLATION KERNELS

"""
    lerp(t, v0, v1)

Linear interpolation between `f(x0)=v0` and `f(x1)=v1` at the normalized coordinate t
"""
@inline lerp(t, v0, v1) = fma(t, v1, fma(-t, v0, v0))

"""
    bilinear(tx, ty, v00, v10, v01, v11)

Linear interpolation between `f(x0,y0)=v0` and `f(x1,y1)=v1` at the normalized coordinates `tx`, `ty`
"""
@inline function ndlinear(tx, ty, v00, v10, v01, v11)
    return lerp(ty, lerp(tx, v00, v10), lerp(tx, v01, v11))
end

@inline function ndlinear(t::NTuple{2,T1}, v::NTuple{4,T2}) where {T1,T2}
    return ndlinear(t[1], t[2], v[1], v[2], v[3], v[4])
end

bilinear(tx, ty, v00, v10, v01, v11) = ndlinear(tx, ty, v00, v10, v01, v11)

"""
    trilinear(tx, ty, tz, v000, v100, v001, v101, v010, v110, v011, v111)

Linear interpolation between `f(x0,y0,z0)=v0` and `f(x1,y1,z1)=v1` at the normalized coordinates `tx`, `ty`, `tz`
"""
@inline function ndlinear(tx, ty, tz, v000, v100, v001, v101, v010, v110, v011, v111)
    return lerp(
        ty,
        ndlinear(tx, tz, v000, v100, v001, v101),
        ndlinear(tx, tz, v010, v110, v011, v111),
    )
end

@inline function ndlinear(t::NTuple{3,T1}, v::NTuple{8,T2}) where {T1,T2}
    return ndlinear(t[1], t[2], t[3], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8])
end

@inline function trilinear(tx, ty, tz, v000, v100, v001, v101, v010, v110, v011, v111)
    return ndlinear(tx, ty, tz, v000, v100, v001, v101, v010, v110, v011, v111)
end


text:
using Base: Number
using PressureSensitiveMats
using Test

@testset "Utils" begin
    ma1, mv1 = moving_stats(ones(5000), 3)
    ma2, mv2 = moving_stats(ones(5000), fill(3, 5000))
    truth = [
        21   9  22  10  23  11  24  12
        17   5  18   6  19   7  20   8
        13   1  14   2  15   3  16   4
        45  33  46  34  47  35  48  36
        41  29  42  30  43  31  44  32
        37  25  38  26  39  27  40  28
        69  57  70  58  71  59  72  60
        65  53  66  54  67  55  68  56
        61  49  62  50  63  51  64  52]

    @test ma1[3:end] ≈ ones(4998)
    @test mv1[5:end] ≈ zeros(4996)

    @test ma2 ≈ ones(5000)
    @test mv2 ≈ zeros(5000)

    @test mat_shape(ones(72)) isa AbstractVector
    @test mat_shape(ones(5000, 72)) isa AbstractMatrix

    @test sfm(zeros(512)) isa Number
    @test active_sfm(randn(35000, 72), 300) isa AbstractVector

    @test reshape_psm(ones(Int, 400, 72) .* (1:72)')[:, :, rand(1:72)] == truth
end

@testset "Combiners" begin
    multi = fill(313.0, (35000, 72)) .+ randn((35000, 72))

    @test combiner(SNR_MAX(), multi) isa AbstractVector
    @test combiner(PCC(), multi) isa AbstractVector
    @test combiner(PCC2(), multi) isa AbstractVector
    @test combiner(EGC(), multi) isa AbstractVector
    @test apply2seg(s -> combiner(MRC_PSD(10), s), multi, 300) isa AbstractVector
end

@testset "Occupancy Detection" begin
    multi = fill(313.0, (35000, 72)) .+ randn((35000, 72))

    multi_occ = occupancy_detection(multi)
    multi_occ_dist = occupancy_detection(multi, max_dist=true)

    @test multi_occ isa BitVector
    @test multi_occ == zeros(35000)

    @test multi_occ_dist isa BitVector
    @test multi_occ_dist == zeros(35000)
end

@testset "Motion Detection" begin
    single = fill(313.0, 35000) .+ randn(35000)
    multi = fill(313.0, (35000, 72)) .+ randn((35000, 72))

    md_s = move_detect(single, min_samples=10)
    md_m = move_detect(multi, min_samples=10)
    md_solei = move_detect(Solei(), single, min_samples=10)

    @test md_s isa BitVector
    @test md_s == zeros(35000)

    @test  md_m isa BitVector
    @test md_m == zeros(35000)

    @test md_solei isa BitVector
    @test md_solei == zeros(35000)
end


text:
export GaussianMeanVariance

abstract type Gaussian <: SoftFactor end

"""
Description:

    A Gaussian with mean-variance parameterization:

    f(out,m,v) = 𝒩(out|m,v) = (2π)^{-D/2} |v|^{-1/2} exp(-1/2 (out - m)' v^{-1} (out - m))

Interfaces:

    1. out
    2. m (mean)
    3. v (covariance)

Construction:

    GaussianMeanVariance(out, m, v, id=:some_id)
"""
mutable struct GaussianMeanVariance <: Gaussian
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function GaussianMeanVariance(out, m, v; id=generateId(GaussianMeanVariance))
        @ensureVariables(out, m, v)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:m] = self.interfaces[2] = associate!(Interface(self), m)
        self.i[:v] = self.interfaces[3] = associate!(Interface(self), v)

        return self
    end
end

slug(::Type{GaussianMeanVariance}) = "𝒩"

format(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = "$(slug(GaussianMeanVariance))(m=$(format(dist.params[:m])), v=$(format(dist.params[:v])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{GaussianMeanVariance}; m=0.0, v=1.0) = ProbabilityDistribution{Univariate, GaussianMeanVariance}(Dict(:m=>m, :v=>v))
ProbabilityDistribution(::Type{GaussianMeanVariance}; m::Number=0.0, v::Number=1.0) = ProbabilityDistribution{Univariate, GaussianMeanVariance}(Dict(:m=>m, :v=>v))
ProbabilityDistribution(::Type{Multivariate}, ::Type{GaussianMeanVariance}; m=[0.0], v=mat(1.0)) = ProbabilityDistribution{Multivariate, GaussianMeanVariance}(Dict(:m=>m, :v=>v))

dims(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = length(dist.params[:m])

vague(::Type{GaussianMeanVariance}) = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=huge)
vague(::Type{GaussianMeanVariance}, dims::Int64) = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(dims), v=huge*diageye(dims))
vague(::Type{GaussianMeanVariance}, dims::Tuple{Int64}) = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(dims), v=huge*diageye(dims[1]))

unsafeMean(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = deepcopy(dist.params[:m]) # unsafe mean

unsafeMode(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = deepcopy(dist.params[:m]) # unsafe mode

unsafeVar(dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}) = dist.params[:v] # unsafe variance
unsafeVar(dist::ProbabilityDistribution{Multivariate, GaussianMeanVariance}) = diag(dist.params[:v])

unsafeCov(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = deepcopy(dist.params[:v]) # unsafe covariance

unsafeMeanCov(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = (deepcopy(dist.params[:m]), deepcopy(dist.params[:v]))

unsafeWeightedMean(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = cholinv(dist.params[:v])*dist.params[:m]

unsafePrecision(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType = cholinv(dist.params[:v])

logPdf(dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}, x) = -0.5*(log(2pi) + log(dist.params[:v]) + (x-dist.params[:m])^2/dist.params[:v])
logPdf(dist::ProbabilityDistribution{Multivariate, GaussianMeanVariance}, x) = -0.5*(dims(dist)*log(2pi) + logdet(dist.params[:v]) + transpose(x-dist.params[:m])*cholinv(dist.params[:v])*(x-dist.params[:m]))

# Converting from m,v to xi,w would require two separate inversions of the covariance matrix;
# this function ensures only a single inversion is performed
function unsafeWeightedMeanPrecision(dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    W = cholinv(dist.params[:v])
    return (W*dist.params[:m], W)
end

isProper(dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}) = (floatmin(Float64) < dist.params[:v] < floatmax(Float64))
isProper(dist::ProbabilityDistribution{Multivariate, GaussianMeanVariance}) = isRoundedPosDef(dist.params[:v])

function ==(t::ProbabilityDistribution{V, GaussianMeanVariance}, u::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    (t === u) && return true
    isApproxEqual(t.params[:m], u.params[:m]) && isApproxEqual(t.params[:v], u.params[:v])
end

function ==(t::ProbabilityDistribution{V, GaussianMeanVariance}, u::ProbabilityDistribution{V, F}) where {V<:VariateType, F<:Gaussian}
    (t === u) && return true
    isApproxEqual(t.params[:m], unsafeMean(u)) && isApproxEqual(t.params[:v], unsafeCov(u))
end

# Average energy functional
function averageEnergy(::Type{GaussianMeanVariance}, marg_out::ProbabilityDistribution{Univariate}, marg_mean::ProbabilityDistribution{Univariate}, marg_var::ProbabilityDistribution{Univariate})
    (m_mean, v_mean) = unsafeMeanCov(marg_mean)
    (m_out, v_out) = unsafeMeanCov(marg_out)

    0.5*log(2*pi) +
    0.5*unsafeLogMean(marg_var) +
    0.5*unsafeInverseMean(marg_var)*(v_out + v_mean + (m_out - m_mean)^2)
end

function averageEnergy(::Type{GaussianMeanVariance}, marg_out::ProbabilityDistribution{Multivariate}, marg_mean::ProbabilityDistribution{Multivariate}, marg_var::ProbabilityDistribution{MatrixVariate})
    (m_mean, v_mean) = unsafeMeanCov(marg_mean)
    (m_out, v_out) = unsafeMeanCov(marg_out)

    0.5*dims(marg_out)*log(2*pi) +
    0.5*unsafeDetLogMean(marg_var) +
    0.5*tr( unsafeInverseMean(marg_var)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)'))
end

function averageEnergy(::Type{GaussianMeanVariance}, marg_out_mean::ProbabilityDistribution{Multivariate, F}, marg_var::ProbabilityDistribution{Univariate}) where F<:Gaussian
    (m, V) = unsafeMeanCov(marg_out_mean)

    0.5*log(2*pi) +
    0.5*unsafeLogMean(marg_var) +
    0.5*unsafeInverseMean(marg_var)*( V[1,1] - V[1,2] - V[2,1] + V[2,2] + (m[1] - m[2])^2 )
end

function averageEnergy(::Type{GaussianMeanVariance}, marg_out_mean::ProbabilityDistribution{Multivariate, F}, marg_var::ProbabilityDistribution{MatrixVariate}) where F<:Gaussian
    (m, V) = unsafeMeanCov(marg_out_mean)
    d = Int64(dims(marg_out_mean)/2)

    0.5*d*log(2*pi) +
    0.5*unsafeDetLogMean(marg_var) +
    0.5*tr( unsafeInverseMean(marg_var)*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ) )
end


text:
# states of request case K
# opened
# closed
function trequest(d, tg, av) # participant iniciates a request
    @unpack T, Av, K, groups = d
    if av["step"]=="0"
        if av["token"]<-200 # suspiciously many requests generated by same person
            msg="Jūsų taškų skaičiaus nepakanka, prašom sudalyvauti talkoje arba pakviesti naują dalyvį."
            sendMessage(tg, chat_id = av["id"], text = msg)
            tbegin(tg, av)
            println("suspiciously many requests generated by same person $(av["id"])")
        else #if !haskey(av,"requestid") ||
            av["path"] = "0"
            av["step"] = "request"
            msg = "Įrašykite klausimą. / Будь ласка, введіть запитання. / Пожалуйста, введите вопрос."
            sendMessage(tg, chat_id = av["id"], text = msg)
        end
        return
    end
    keypath = join(av["path"])
    if av["txt"] == "Parinkti | Виберіть" || av["txt"] == "Atgal | Повернутися"
        tree(tg, av, T)
    elseif av["txt"] == "Pasirinkti" || av["txt"] == "Siųsti | Відправити"
        if av["txt"] == "Siųsti | Відправити"
            keypath = "0Paslaugos | Послуги "
        else
            while isempty(T[keypath]["dav_id"])
                if length(T[keypath]["steps"])<2
                    keypath = "0Paslaugos | Послуги "
                    break
                else
                    keypath = join(T[keypath]["steps"][1:end-1])
                end
            end
        end
        msg = "Lukterėkite, tuoj surasiu talkininką. | Зачекайте, я скоро знайду помічника."
        sendMessage(tg, chat_id = av["id"], text = msg)
        # TODO: redundant keis? Do not understand the question
        keis = K[av["requestid"]]
        subject = join(av["path"][2:end], "/")
        if isempty(subject)
            subject = "bendros"
        end
        keis["subject"] = subject
        keis["dav_id"] = T[keypath]["dav_id"]
        keis["requested_id"] = []
        sendrequest(d, tg, av)
    elseif any(T[keypath]["children"] .== av["txt"])
        av["path"] = vcat(av["path"], av["txt"])
        tree(tg, av, T)
    else #įvedamas užklausimo tekstas
        av["requestid"] = now()
        K[av["requestid"]] = Dict("getter"=>av,"txt"=>deepcopy(av["txt"]),"giver"=>[],"state"=>"opened")
        msg = "Siūskite užklausimą iškart arba pasirinkite iš kurios srities klausimas. | Надішліть запит негайно або виберіть, з якої області запитати."
        k = [["Siųsti | Відправити","Parinkti | Виберіть"]]
        d = Dict(:keyboard => k, :one_time_keyboard => true, :resize_keyboard=>true)
        sendMessage(tg, chat_id = av["id"], text = msg, reply_markup = d)
    end
end

function sendrequest(d::DataBase, tg, av) # the request is broadcasted
    @unpack Av, K = d

    keis = K[av["requestid"]]
    if length(keis["requested_id"]) == length(keis["dav_id"])
        msg = """Pakartotinas užklausimas "$(keis["txt"])" iš "$(keis["subject"])" temos. Gal pažįstat kas gali atsakyti į šį klausimą?"""
        keis["requested_id"] = []
    else
        msg = """Užklausimas "$(keis["txt"])" iš $(keis["subject"]) temos"""
    end
    k = [["Priimti","Atmesti"]]
    kb = Dict(:keyboard => k, :one_time_keyboard => true, :resize_keyboard=>true)
    i = 1
    for id in shuffle(keis["dav_id"])
        if id == av["id"] || any(keis["requested_id"].==id)
            continue
        end
        try
            sendMessage(tg, chat_id = id, text = msg, reply_markup = kb)
        catch
            println("Didn't manage to send to $id, $(Av[id]["first_name"])")
        end
        Av[id]["step"] = "requested"
        Av[id]["getterid"] = av["id"]
        push!(keis["requested_id"], id)
        i=i+1
        if length(i) > 4 #this is one portion of sent requests
            break
        end
    end
end

function dealkeis(d, tg, av)
    @unpack Av, K, groups = d
    if av["txt"] == "Priimti"
        avg = Av[av["getterid"]]
        keis = K[avg["requestid"]]
        if avg["step"] == "request"
            keis["group"] = getgroup(groups)
            keis["group"].second["state"] = "notfree"
            msg = "Aptarkit klausimą prisijungus prie grupės: | Обговоріть проблему, приєднавшись до групи:  $(keis["group"].second["link"])"
            sendMessage(chat_id = av["getterid"], text = msg)
            avg["step"] = "accepted"
            msg = "Atėjus laikui uždarykite klausimą. | Закрийте питання, коли прийде час."
            keyboard = [["Uždaryti klausimą | Закрийте питання","closekeis()"]]
            #sendMessage(chat_id = keis["group"].first, text = msg, reply_markup = tik(keyboard)) # the bottum could be also in a group (at present is too complicated)
            sendMessage(tg, chat_id = avg["id"], text = msg, reply_markup = tik(keyboard))
            if avg["id"] !== 5090964479 # if not owner of the group
                unbanChatMember(tg, chat_id = keis["group"].first, user_id = avg["id"])
            end
        elseif avg["step"] !== "accepted"
            return
        end
        msg = "Aptarkit klausimą prisijungus prie grupės: | Обговоріть проблему, приєднавшись до групи: $(keis["group"].second["link"])"
        sendMessage(chat_id = av["id"], text = msg)
        push!(keis["giver"], av["id"])
        keis["state"] = "accepted"
        if av["id"] !== 5090964479 # owner of the group
            unbanChatMember(chat_id = keis["group"].first,user_id = av["id"])
        end
    elseif av["txt"] == "Atmesti"
        av["getterid"] = nothing
        tbegin(tg, av)
    end
end

function chekeis(d, tg, av)
    @unpack K, Av = d
    n = now()
    for k in keys(K)
        if (k+Minute(30))<n
            if K[k]["state"]=="closed"
                continue
            else
                closekeis(d, tg, av, k)
            end
        end
        if K[k]["state"] == "opened" && (k+Minute(1)) < n
            av = Av[K[k]["getter"]["id"]]
            if (k+Minute(7))<n
                av = closekeis(d, tg, av, k)
            else
                sendrequest(d, tg, av)
            end
        end
    end

    return av
end

function closekeis(d, tg, av, k)
    @unpack K, Av = d
    keis = K[k]
    if keis["state"] == "closed"
        return av
    end
    av = Av[keis["getter"]["id"]]
    av["step"]="0"
    if keis["state"] !== "accepted"
        keis["state"] = "closed"
        sendMessage(tg, chat_id = av["id"], text = "Nepavyko rasti talkininko, bandykite dar kartą. | Не вдалося знайти допомогу, спробуйте ще раз.")
        tbegin(tg, av)
        return av
    end
    keis["state"] = "closed"
    keis["group"].second["state"] = "free"
    chm = getChatMember(tg, chat_id = keis["group"].first, user_id = av["id"])
    if !(chm["status"] == "left") && !(chm["status"] == "creator")
        banChatMember(tg, chat_id = keis["group"].first, user_id = av["id"])
    end
    av["token"] = av["token"]-1
    tfg = 1/(length(keis["giver"])+0.2)
    sendMessage(tg, chat_id = av["id"], text = "Susitikimas uždarytas. | Зустріч закрита.") #(jūsų taškų skaičius $(round(10*av["token"])/10)")
    msg = "Ar talka buvo naudinga? | Чи була допомога корисною?"
    k = [["Taip | Так", "valuableyesno(1)", "Ne | Ні", "valuableyesno(0)"]]
    sendMessage(tg, chat_id = av["id"], text = msg, reply_markup = tik(k))
    for gid in keis["giver"]
        av = Av[gid]
        chm = getChatMember(tg, chat_id = keis["group"].first, user_id = gid)
        if !(chm["status"] == "left") && !(chm["status"] == "creator")
            banChatMember(tg, chat_id = keis["group"].first, user_id = gid)
        end
        Av[gid]["token"] += tfg
        sendMessage(tg, chat_id = gid, text = "Susitikimas uždarytas. | Зустріч закрита.") #(jūsų taškų skaičius $(round(10*Av[gid]["token"])/10))")
        tbegin(tg, av)
        # msg="Galima eiti į pradžią"
        # k=[["Į pradžią", "tbegin()"]]
        # sendMessage(chat_id = gid, text = msg, reply_markup = tik(k))
    end

    return av
end


text:
function settrue(x,b)
    fill!(x,false)
    x[b]=true
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2019 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Function: DataDictionary
# Description: Holds simulation and model parameters as key => value pairs in a Julia Dict()
# Generated on: 2019-10-02T16:15:21.453
#
# Input arguments:
# time_start::Float64 => Simulation start time value (scalar)
# time_stop::Float64 => Simulation stop time value (scalar)
# time_step::Float64 => Simulation time step (scalar)
#
# Output arguments:
# data_dictionary::Dict{String,Any} => Dictionary holding model and simulation parameters as key => value pairs
# ----------------------------------------------------------------------------------- #
function build_data_dictionary(time_span::Tuple{Float64,Float64,Float64}, path_to_biophysical_constants_file::String = "./Default.json", host_type::Symbol = :bacteria)::Dict{String,Any}

	# load the biophysical_constants dictionary
	biophysical_constants_dictionary = build_biophysical_dictionary(path_to_biophysical_constants_file, host_type)

	# stoichiometric_matrix and dilution_matrix -
	stoichiometric_matrix = readdlm("./Network.dat")

	# number of states, and rates -
	(number_of_states,number_of_rates) = size(stoichiometric_matrix)

	# array of species types -
	species_symbol_type_array = [
		:gene	;	# 1	clssrA
		:gene	;	# 2	deGFPssrA
		:gene	;	# 3	sigma_28
		:gene	;	# 4	sigma_70
		:mrna	;	# 5	mRNA_clssrA
		:mrna	;	# 6	mRNA_deGFPssrA
		:mrna	;	# 7	mRNA_sigma_28
		:mrna	;	# 8	mRNA_sigma_70
		:protein	;	# 9	protein_clssrA
		:protein	;	# 10	protein_deGFPssrA
		:protein	;	# 11	protein_sigma_28
		:protein	;	# 12	protein_sigma_70
	]

	# for some whacky reason, we need to add the species_symbol_type_array to the biophysical dictionary -
	biophysical_constants_dictionary["species_symbol_type_array"] = species_symbol_type_array

	# array of gene lengths -
	gene_coding_length_array = [
		1000.0	;	# 1	clssrA
		1000.0	;	# 2	deGFPssrA
		1000.0	;	# 3	sigma_28
		1000.0	;	# 4	sigma_70
	]

	# array of mRNA coding lengths -
	mRNA_coding_length_array = [
		gene_coding_length_array[1]	;	# 5	1	mRNA_clssrA
		gene_coding_length_array[2]	;	# 6	2	mRNA_deGFPssrA
		gene_coding_length_array[3]	;	# 7	3	mRNA_sigma_28
		gene_coding_length_array[4]	;	# 8	4	mRNA_sigma_70
	]

	# array of mRNA coding lengths -
	protein_coding_length_array = [
		round((0.33)*mRNA_coding_length_array[1])	;	# 9	1	protein_clssrA
		round((0.33)*mRNA_coding_length_array[2])	;	# 10	2	protein_deGFPssrA
		round((0.33)*mRNA_coding_length_array[3])	;	# 11	3	protein_sigma_28
		round((0.33)*mRNA_coding_length_array[4])	;	# 12	4	protein_sigma_70
	]

	# array of gene concentrations -
	gene_abundance_array = [
		5.0	;	# (nM) 1	clssrA
		5.0	;	# (nM) 2	deGFPssrA
		5.0	;	# (nM) 3	sigma_28
		0.0	;	# (nM) 4	sigma_70	# we have no sigma 70
	]

	# initial condition array -
	initial_condition_array = [
		gene_abundance_array[1]		;	# 1	clssrA
		gene_abundance_array[2]		;	# 2	deGFPssrA
		gene_abundance_array[3]		;	# 3	sigma_28
		gene_abundance_array[4]		;	# 4	sigma_70
		0.0							;	# 5	mRNA_clssrA
		0.0							;	# 6	mRNA_deGFPssrA
		0.0							;	# 7	mRNA_sigma_28
		0.0							;	# 8	mRNA_sigma_70
		0.0							;	# 9	protein_clssrA
		0.0							;	# 10	protein_deGFPssrA
		0.0							;	# 11	protein_sigma_28
		0.0							;	# 12	protein_sigma_70

		# extra state -
		100.0						;	# 13
	]

	binding_parameter_dictionary = Dict{String,Float64}()
	binding_parameter_dictionary["n_clssrA_sigma_28"] = 1.0
	binding_parameter_dictionary["K_clssrA_sigma_28"] = 0.05
	binding_parameter_dictionary["n_deGFPssrA_sigma_70"] = 1.0
	binding_parameter_dictionary["K_deGFPssrA_sigma_70"] = 0.05
	binding_parameter_dictionary["n_deGFPssrA_clssrA"] = 1.0
	binding_parameter_dictionary["K_deGFPssrA_clssrA"] = 0.05
	binding_parameter_dictionary["n_sigma_28_sigma_70"] = 1.0
	binding_parameter_dictionary["K_sigma_28_sigma_70"] = 0.05
	binding_parameter_dictionary["n_sigma_28_clssrA"] = 1.0
	binding_parameter_dictionary["K_sigma_28_clssrA"] = 0.05

	# Alias the control function parameters -
	control_parameter_dictionary = Dict{String,Float64}()
	control_parameter_dictionary["W_clssrA_RNAP"] = 0.001
	control_parameter_dictionary["W_clssrA_sigma_28"] = 1.0
	control_parameter_dictionary["W_deGFPssrA_RNAP"] = 0.001
	control_parameter_dictionary["W_deGFPssrA_sigma_70"] = 1.0
	control_parameter_dictionary["W_deGFPssrA_clssrA"] = 1.0
	control_parameter_dictionary["W_sigma_28_RNAP"] = 0.001
	control_parameter_dictionary["W_sigma_28_sigma_70"] = 1.0
	control_parameter_dictionary["W_sigma_28_clssrA"] = 1.0
	control_parameter_dictionary["W_sigma_70_RNAP"] = 0.0

	# degradation modifiers -
	degradation_modifier_array = [
		0.0	;	# 1	clssrA
		0.0	;	# 2	deGFPssrA
		0.0	;	# 3	sigma_28
		0.0	;	# 4	sigma_70
		1.0	;	# 5	mRNA_clssrA
		1.0	;	# 6	mRNA_deGFPssrA
		1.0	;	# 7	mRNA_sigma_28
		1.0	;	# 8	mRNA_sigma_70
		1.0	;	# 9	protein_clssrA
		1.0	;	# 10	protein_deGFPssrA
		1.0	;	# 11	protein_sigma_28
		1.0	;	# 12	protein_sigma_70
	]

	# time constant modifiers -
	time_constant_modifier_array = [
		0.0	;	# 1	clssrA
		0.0	;	# 2	deGFPssrA
		0.0	;	# 3	sigma_28
		0.0	;	# 4	sigma_70
		1.0	;	# 5	mRNA_clssrA
		1.0	;	# 6	mRNA_deGFPssrA
		1.0	;	# 7	mRNA_sigma_28
		1.0	;	# 8	mRNA_sigma_70
		1.0	;	# 9	protein_clssrA
		1.0	;	# 10	protein_deGFPssrA
		1.0	;	# 11	protein_sigma_28
		1.0	;	# 12	protein_sigma_70
	]

	# Dilution degrdation matrix -
	dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)

	# Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)

	# Parameter name index array -
	parameter_name_mapping_array = [
		"n_clssrA_sigma_28"	;	# 1
		"K_clssrA_sigma_28"	;	# 2
		"n_deGFPssrA_sigma_70"	;	# 3
		"K_deGFPssrA_sigma_70"	;	# 4
		"n_deGFPssrA_clssrA"	;	# 5
		"K_deGFPssrA_clssrA"	;	# 6
		"n_sigma_28_sigma_70"	;	# 7
		"K_sigma_28_sigma_70"	;	# 8
		"n_sigma_28_clssrA"	;	# 9
		"K_sigma_28_clssrA"	;	# 10
		"W_clssrA_RNAP"	;	# 11
		"W_clssrA_sigma_28"	;	# 12
		"W_deGFPssrA_RNAP"	;	# 13
		"W_deGFPssrA_sigma_70"	;	# 14
		"W_deGFPssrA_clssrA"	;	# 15
		"W_sigma_28_RNAP"	;	# 16
		"W_sigma_28_sigma_70"	;	# 17
		"W_sigma_28_clssrA"	;	# 18
		"W_sigma_70_RNAP"	;	# 19
		"rnapII_concentration"	;	# 20
		"ribosome_concentration"	;	# 21
		"degradation_constant_mRNA"	;	# 22
		"degradation_constant_protein"	;	# 23
		"kcat_transcription"	;	# 24
		"kcat_translation"	;	# 25
		"maximum_specific_growth_rate"	;	# 26
		"saturation_constant_transcription"	;	# 27
		"saturation_constant_translation"	;	# 28
	]

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{String,Any}()
	data_dictionary["number_of_states"] = number_of_states
	data_dictionary["species_symbol_type_array"] = species_symbol_type_array
	data_dictionary["initial_condition_array"] = initial_condition_array
	data_dictionary["gene_coding_length_array"] = gene_coding_length_array
	data_dictionary["mRNA_coding_length_array"] = mRNA_coding_length_array
	data_dictionary["protein_coding_length_array"] = protein_coding_length_array
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix
	data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary
	data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary
	data_dictionary["parameter_name_mapping_array"] = parameter_name_mapping_array
	data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array
	data_dictionary["translation_parameter_array"] = translation_parameter_array
	data_dictionary["degradation_modifier_array"] = degradation_modifier_array
	data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array
	data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary
	data_dictionary["half_life_translation_capacity"] = 8.0

	# extra stuff -
	data_dictionary["R"] = 8.314 			# J mol^-1 K^-1
	data_dictionary["T_K"] = 273.15 + 29.0 	# K
	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
using SafeTestsets

# Start Test Script

@time begin
    @time @safetestset "Utilities" begin
        include("util.jl")
    end

    @time @safetestset "Risk-Adjusted Linearization" begin
        include("risk_adjusted_linearization/constructors.jl")
    end

    @time @safetestset "Numerical algorithms" begin
        include("numerical_algorithms/blanchard_kahn.jl")
        include("numerical_algorithms/compute_psi.jl")
        include("numerical_algorithms/deterministic.jl")
        include("numerical_algorithms/relaxation.jl")
        include("numerical_algorithms/homotopy.jl")
        include("numerical_algorithms/solve.jl")
    end

    @time @safetestset "Examples" begin
        include(joinpath(dirname(@__FILE__), "..", "examples", "rbc_cc", "example_rbc_cc.jl"))
        include(joinpath(dirname(@__FILE__), "..", "examples", "wachter_disaster_risk", "example_wachter.jl"))
        include(joinpath(dirname(@__FILE__), "..", "examples", "crw", "example_crw.jl")) # This example tests case of jump-dependent Σ and Λ
        include(joinpath(dirname(@__FILE__), "..", "examples", "textbook_nk", "example_textbook_nk.jl"))
        include(joinpath(dirname(@__FILE__), "..", "examples", "nk_with_capital", "example_nk_with_capital.jl"))
        include("examples_to_test/rbc_cc.jl")
    end

    @time @safetestset "Sparse Arrays and Jacobians for Objects of Risk-Adjusted Linearizations" begin
        include("sparse_ral_objects/sparse_mu_xi_V_jacobians.jl")
        include("sparse_ral_objects/sparse_array_caches.jl")
    end

    @time @safetestset "Simulations, Impulse Responses, and Simulation-Based Diagnostics" begin
        include("simulation/simulate_rbc_cc.jl")
        include("simulation/simulate_wachter_disaster_risk.jl")
        include("simulation/simulate_crw.jl")
        include("simulation/impulse_responses_rbc_cc.jl")
        include("simulation/impulse_responses_wachter_disaster_risk.jl")
        include("simulation/impulse_responses_crw.jl")
        include("simulation/gausshermite_expectation.jl")
        include("simulation/euler_equation_error.jl")
    end
end


text:
include("./exportData.jl")
using .exportData, DifferentialEquations, DiffEqFlux, PyPlot, DelimitedFiles

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end
function fenkar!(dx, x, p, t)

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi,t0,TI,IA = p[1:16]

	# fenkar dynamics
	dx[1] = H(t-t0;k=1.0).*IA*sin(pi*(t-t0)/TI)^500 - (x[1]*H(uc-x[1])/to + H(x[1]-uc)/tr - x[2]*H(x[1]-uc)*(1.0-x[1])*(x[1]-uc)/td - x[3]*H(x[1]-ucsi;k=xk)/tsi)
	dx[2] = H(uc-x[1])*(1.0-x[2])/(tv1m*H(uv-x[1]) + tv2m*H(x[1]-uv)) - H(x[1]-uc)*x[2]/tvp
	dx[3] = H(uc-x[1])*(1.0-x[3])/twm - H(x[1]-uc)*x[3]/twp

	return nothing
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
# and then p[14:16] = [offset_time, stimulus_period, stimulus_amplitude]
   p = [29.0, 19.6, 1250.0, 3.33, 41.0, 870.0, 0.25, 12.5, 33.3, 10.0, 0.13, 0.04, 0.85,  0.0,   125.0, 0.05]

tspan = (0.0, 3000.0)
u0 = zeros(Float64,3)

#solve and take data
prob = ODEProblem(fenkar!, u0, tspan, p)
sol = solve(prob, Tsit5(); saveat=tspan[begin]:5.0:tspan[end])

# getExpData
ind = 114
t, data, p[15] = getExpData(ind, tInds=1:1000)
p[14] = t[findfirst(data[1,:] .> 0.5)]-p[15]/2

# remake problem
prob = remake(prob, tspan=(0.0, t[end]), p=p)

# use the save_idxs to select the mid-point of the u-variable
sol = solve(prob, Tsit5(); saveat=t, save_idxs=1:1)

function loss(p)
  sol = solve(prob, Tsit5(), p=p, saveat = t, save_idxs=1:1)
  if any((s.retcode != :Success for s in sol))
     loss = Inf
   else
     loss = sum(abs2, sol.-data)
   end
  return loss, sol
end
l0, pred = loss(p)

callback = function (p, l, pred; saving=false)
  @show l
  plt.cla()
  plt.plot(t, data[1,:], "-")
  plt.plot(sol.t, sol[1,:], "--", linewidth=1)
  plt.plot(pred.t, pred[1,:], "--")
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  if saving
  	plt.savefig("./ode/flux/$(ind)_fittings.pdf")
  end
  return false
end

fig = plt.figure(figsize=(4,2))
callback(p, l0, pred)

result = DiffEqFlux.sciml_train(loss, p, cb = callback)
l, pred = loss(result)
callback(result, l, pred; saving=true)
plt.close("all")

open("./ode/flux/$(ind)_params.txt", "w") do io
	writedlm(io, ["tsi" "tv1m" "tv2m" "tvp" "twm" "twp" "td" "to" "tr" "xk" "uc" "uv" "ucsi" "t0" "TI" "IA"])
	writedlm(io, "\n")
	writedlm(io, transpose(p))
	writedlm(io, "\n")
	writedlm(io, transpose(result))
end


text:
module VibronicToolkit

import LinearAlgebra: diag, isdiag
using LinearAlgebra: diagm, Diagonal, eigen, eigvals, I, issymmetric, svd,
                     Symmetric, tr
using Printf: @printf
using Statistics: mean, std

using Distributions: MvNormal, sample
import JSON
using ProgressMeter: Progress, ProgressWrapper
using StatsBase: Weights

# We re-export diag and isdiag from LinearAlgebra.
export
    DegeneracyException,
    SurfaceCouplingException,

    System,
    DenseSystem,
    DiagonalSystem,

    diag,
    isdiag,
    simplify,
    weights,

    Analytical,
    SumOverStates,
    Trotter,
    SamplingFiniteDifference,
    SamplingPrimitiveThermodynamic,

    ground_pes,
    path_mean_std,

    IterativeDecomposition

"""
Solution for vibronic problem.
"""
abstract type Solution end

include("utilities.jl")

include("system.jl")
include("operators.jl")

include("jackknife.jl")

include("analytical.jl")
include("sos.jl")
include("trotter.jl")

include("sampling.jl")
include("sampling_finite_difference.jl")
include("sampling_primitive_thermodynamic.jl")

include("pes.jl")

include("iterative_decomposition.jl")

end


text:
const DEFAULT_ENCODINGS = [
    :domain_wall,
    :one_hot,
]

const DEFAULT_OPTIMIZERS = [
    :descent,
]

const DEFAULT_LEARNING_RATES = [
    0.1,
    0.01,
    0.001,
]

const DEFAULT_DOMAIN_SIZES = [
    3,
    5,
    10,
]

const DEFAULT_TRAINING_SET_SIZES = [
    100,
    1000,
    10000,
    100000,
]

const DEFAULT_ROUNDIND_PRECISIONS = [
    100,
    10,
    5,
    1,
]

# includes
include("all_different.jl")
include("ordered.jl")

const DEFAULT_CONSTRAINTS = [
    :all_different => ALL_DIFFERENT,
    :ordered => ORDERED,
]

const ALL_PARAMETERS = Dict(
    # Search space parameters
    :constraint => DEFAULT_CONSTRAINTS,
    :domain_size => DEFAULT_DOMAIN_SIZES,
    :encoding => DEFAULT_ENCODINGS,
    :training_size_set => DEFAULT_TRAINING_SET_SIZES,

    # Learning parameters
    :learning_rate => DEFAULT_LEARNING_RATES,
    :optimizer => DEFAULT_OPTIMIZERS,
    :rouding_precision => DEFAULT_ROUNDIND_PRECISIONS,
)


text:
using Test, Pkg

const GROUP = get(ENV, "GROUP", "All")

function dev_subpkg(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.develop(PackageSpec(path=subpkg_path))
end

function activate_subpkg_env(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.activate(subpkg_path)
    Pkg.develop(PackageSpec(path=subpkg_path))
    Pkg.instantiate()
end

# Add necessary sub-dependencies
if GROUP == "ArrayInterfaceBlockBandedMatrices"
    dev_subpkg("ArrayInterfaceBandedMatrices")
end

if GROUP == "ArrayInterfaceStaticArrays"
    dev_subpkg("ArrayInterfaceStaticArraysCore")
end

groups = if GROUP == "All"
    ["ArrayInterfaceCore", "ArrayInterface", "ArrayInterfaceBandedMatrices", "ArrayInterfaceBlockBandedMatrices",
     "ArrayInterfaceOffsetArrays", "ArrayInterfaceStaticArrays", "ArrayInterfaceStaticArraysCore"]
else
    [GROUP]
end

@time begin

for g in groups
    if g == "ArrayInterface"

        println("ArrayInterface Tests")

        include("setup.jl")
        include("array_index.jl")
        include("axes.jl")
        include("broadcast.jl")
        include("dimensions.jl")
        include("indexing.jl")
        include("ranges.jl")
        include("size.jl")
        include("stridelayout.jl")
        include("misc.jl")
    else
        dev_subpkg(g)
        subpkg_path = joinpath(dirname(@__DIR__), "lib", g)
        Pkg.test(PackageSpec(name=g, path=subpkg_path))
    end
end

end


text:
using StaticArrays, LinearAlgebra, ForwardDiff, Plots

mutable struct Hopper{T}
    mb::T
    ml::T
    Jb::T
    Jl::T
    r::T
    μ::T
    g::T
end

kinematics(model, q) = @SVector [q[1] + q[3]*sin(q[5]), q[2] - q[3]*cos(q[5])]
# Dimensions
nq = 5 # configuration dim
nu = 2 # control dim
nc = 1 # number of contact points
nf = 2 # number of faces for friction cone pyramid

# Parameters
g = 9.81 # gravity
Δt = 0.1 # time step
μ = 0.5  # coefficient of friction
mb = 10. # body mass
ml = 1.  # leg mass
Jb = 2.5 # body inertia
Jl = 0.25 # leg inertia
r = 0.7 # nominal leg length
model = Hopper(mb,ml,Jb,Jl,r,μ,g)

# from "Dynamically Stable Legged Locomotion p. 144"

function gen_z(func,t,t_flight,t_stance)
    step = 0
    z = zero(t)

    for (i,tt) in enumerate(t)
        if tt > (step+0.5)*(t_flight + t_stance)
            step += 1
        end

        if step == 0
            z[i] = func(tt+0.5*t_flight)
        else
            z[i] = func(tt-((step-1+0.5)*t_flight + step*t_stance))
        end
    end
    return z
end
model.m
function raibert_trajectories(h::Hopper,ẋ,t_step,ρ,tf,T;compress=0.875)
    t_flight = (1-ρ)*t_step
    t_stance = ρ*t_step

    Jeff = h.Jb*h.Jl/(h.Jb+h.Jl)
    @show kh = ((2π/t_step)^2)*Jeff
    ωh = sqrt(kh/Jeff)
    z_min = compress*h.r

	@show kℓ = h.ml * (π / (ρ * t_step))^2.0
    @show dθ0 = ẋ/z_min

    θ_max = dθ0/ωh

    t = range(0,stop=tf,length=T)
    r_param = [1 1; cos(ωh*(t_flight*0.5+t_stance*0.5)) 1]\[h.r;z_min]

    @show z_max = 0.125*h.g*((1-ρ)^2)*t_step^2 + h.r*cos(dθ0/ωh*sin(ωh*t_flight*0.5))

    żlo = 0.5*t_flight*h.g
    zlo = z_max - 0.5/h.g*żlo^2
    function z_traj(i)
        zlo + żlo*i - 0.5*h.g*i^2
    end

    x = range(0,stop=ẋ*tf,length=T)
    z = gen_z(z_traj,t,t_flight,t_stance)
    r_traj = r_param[1]*cos.(ωh.*t) .+ r_param[2]
    θ = θ_max.*sin.(ωh.*t)
    ϕ = -Jl/Jb.*θ

    M_traj = kh * θ
    ṙ_traj = kℓ * (r_traj .- r)

    Q = [[x[t];z[t];r_traj[t]; ϕ[t];θ[t];] for t=1:T]
    U = [[M_traj[t];ṙ_traj[t]] for t=1:T]

    return Q, U, x, z, r_traj, ϕ, θ, M_traj, ṙ_traj, t
end

ẋ = 2
t_step = 0.5
ρ = 0.125
tf = 1.0
T = 1000
h = tf / (T - 1)
Q,U,_x,_z,_r,_θ,_ϕ,_M,_ṙ,_t = raibert_trajectories(model,ẋ,t_step,ρ,tf,T,compress=0.89)

plot(_t,_z,color=:red,width=2.,label="z")
plot!(_t,_r,color=:orange,label="r",width=2.)
plot!(_t,_θ,color=:yellow,label="theta",width=2.)
plot!(_t,_ϕ,color=:cyan,label="phi",width=2.)

plot(_t,_M)
plot!(_t,_ṙ)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, Q, Δt = h)

# Visualization
function visualize!(vis, model::Hopper, q;
		Δt = 0.1, scenario = :vertical)

    r_foot = 0.05
    r_leg = 0.5 * r_foot

    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))

    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

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

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

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

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

        z_shift = [0.0; 0.0; r_foot]

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

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

	# settransform!(vis["/Cameras/default"],
	# 	compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end


text:
## Adapted from sundialsTB/kinsol/examples_ser/mkinTest_nds.m

## %mkinTest_dns - KINSOL example problem (serial, dense)
## %   Simple test problem for the Dense linear solver in KINSOL
## %   This example solves the system
## %       y(1)^2 + y(2)^2 = 1
## %       y(2) = y(1)^2
## %

## % Radu Serban <radu@llnl.gov>

import Sundials

## function to be optimized
function sysfn(y_in, fy_in, a_in)
    y = Sundials.asarray(y_in)
    fy = Sundials.asarray(fy_in)
    fy[1] = y[1]^2 + y[2]^2 - 1.0
    fy[2] = y[2] - y[1]^2
    return int32(0)   # indicates normal return
end

## Initialize problem
neq = 2
kmem = Sundials.KINCreate()
flag = Sundials.KINSetFuncNormTol(kmem, 1.0e-5)
flag = Sundials.KINSetScaledStepTol(kmem, 1.0e-4)
flag = Sundials.KINSetMaxSetupCalls(kmem, 1)
y = ones(neq)
flag = Sundials.KINInit(kmem, sysfn, y)
flag = Sundials.KINDense(kmem, neq)
## Solve problem
scale = ones(neq)
strategy = 1   # KIN_LINESEARCH
flag = Sundials.KINSol(kmem,
                       y,
                       strategy,
                       scale,
                       scale)

println("Solution: ", y)
residual = ones(2)
sysfn(y, residual, [1,2])
println("Residual: ", residual)

## Sundials.KINFree(kmem)   # segfaults, probably because Julia arrays and NVectors share memory


text:
struct ThermalStandardUCOutages <: PSI.AbstractStandardUnitCommitment end
struct ThermalDispatchOutages <: PSI.AbstractThermalDispatchFormulation end

function PSI.time_constraints!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::PSI.DeviceModel{T, ThermalStandardUCOutages},
    ::Type{S},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    parameters = PSI.model_has_parameters(optimization_container)
    resolution = PSI.model_resolution(optimization_container)
    initial_conditions_on =
        PSI.get_initial_conditions(optimization_container, PSI.ICKey(PSI.TimeDurationON, T))
    initial_conditions_off = PSI.get_initial_conditions(
        optimization_container,
        PSI.ICKey(PSI.TimeDurationOFF, T),
    )
    ini_conds, time_params =
        PSI._get_data_for_tdc(initial_conditions_on, initial_conditions_off, resolution)
    forecast_label = "outage"
    constraint_infos = Vector{DeviceDurationConstraintInfo}()
    for (ix, ic) in enumerate(ini_conds[:, 1])
        name = PSI.device_name(ic)
        info = DeviceDurationConstraintInfo(
            name,
            time_params[ix],
            Tuple(ini_conds[ix, :]),
            1.0,
            PSI.get_time_series(optimization_container, ic.device, forecast_label),
        )
        push!(constraint_infos, info)
    end

    if !(isempty(ini_conds))
        if parameters
            device_duration_parameters_outage!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(PSI.DURATION, T),
                (
                    PSI.make_variable_name(PSI.OnVariable, T),
                    PSI.make_variable_name(PSI.StartVariable, T),
                    PSI.make_variable_name(PSI.StopVariable, T),
                ),
                PSI.UpdateRef{T}(OUTAGE, forecast_label),
            )
        else
            device_duration_look_ahead_outage!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(PSI.DURATION, T),
                (
                    PSI.make_variable_name(PSI.OnVariable, T),
                    PSI.make_variable_name(PSI.StartVariable, T),
                    PSI.make_variable_name(PSI.StopVariable, T),
                ),
            )
        end
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end
    return
end

function outage_constraints!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::PSI.DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel, D <: Union{ThermalStandardUCOutages, ThermalDispatchOutages}}
    parameters = PSI.model_has_parameters(optimization_container)
    resolution = PSI.model_resolution(optimization_container)
    initial_conditions =
        PSI.get_initial_conditions(optimization_container, PSI.ICKey(OutageStatus, T))
    forecast_label = "outage"
    constraint_infos = Vector{DeviceOutageConstraintInfo}()
    for (ix, ic) in enumerate(initial_conditions)
        name = PSI.device_name(ic)
        info = DeviceOutageConstraintInfo(
            name,
            ic,
            1.0,
            PSI.get_time_series(optimization_container, ic.device, forecast_label),
        )
        push!(constraint_infos, info)
    end

    if !(isempty(initial_conditions))
        if parameters
            device_outage_ub_parameter!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ON, T),
                PSI.UpdateRef{T}(OUTAGE, forecast_label),
            )
        else
            device_outage_ub!(
                optimization_container,
                constraint_infos,
                PSI.make_constraint_name(OUTAGE, T),
                PSI.make_variable_name(PSI.ON, T),
            )
        end
    else
        @warn "Data doesn't contain generators with initial condition for outage status, consider adjusting your formulation"
    end

    return
end

function PSI.initial_conditions!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    formulation::D,
) where {T <: PSY.ThermalGen, D <: Union{ThermalStandardUCOutages, ThermalDispatchOutages}}
    PSI.status_initial_condition!(optimization_container, devices, formulation)
    PSI.output_initial_condition!(optimization_container, devices, formulation)
    PSI.duration_initial_condition!(optimization_container, devices, formulation)
    outage_status_initial_condition!(optimization_container, devices, formulation)
    return
end

function outage_status_initial_condition!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.ThermalGen, D <: Union{ThermalStandardUCOutages, ThermalDispatchOutages}}
    PSI._make_initial_conditions!(
        optimization_container,
        devices,
        ThermalStandardUCOutages(),
        nothing,
        PSI.ICKey(OutageStatus, T),
        _make_initial_condition_outage_status,
        _get_outage_initial_value,
    )

    return
end


text:
@testset "jd2cal" begin
    y, m, d, fd = jd2cal(2400000.5, 50123.9999)
    @test (y, m, d) == (1996, 2, 10)
    @test fd ≈ 0.9999 atol=1e-7
    @test_throws ERFAException jd2cal(typemax(Cdouble), 0.0)
end

@testset "jdcalf" begin
    y, m, d, fd = jdcalf(4, 2400000.5, 50123.9999)
    @test (y, m, d, fd) == (1996, 2, 10, 9999)
    @test_throws ERFAException jdcalf(4, 2e9, 0.0)
    @test_logs (:warn,) jdcalf(12, 2400000.5, 50123.9999)
end


text:
using GraphTheoryMips
#using Base.Test
using FactCheck
using JuMP
using Cbc
using LightGraphs

function make_simple_undirected_graph{T<:Integer}(n::T, edgelist::Vector{Tuple{T,T}})
    g = Graph(n)
    for (s,d) in edgelist
        add_edge!(g, Edge(s,d))
    end
    return g
end

function PetersenGraph()
    e = [
        (1, 2), (1, 5), (1, 6),
        (2, 3), (2, 7),
        (3, 4), (3, 8),
        (4, 5), (4, 9),
        (5, 10),
        (6, 8), (6, 9),
        (7, 9), (7, 10),
        (8, 10)
    ]
    return make_simple_undirected_graph(10,e)
end

facts("Testing graph mips") do
    petersen = PetersenGraph()
    k_3_3 = CompleteBipartiteGraph(3,3)

    @fact length(bin_decisions(max_independent_set(petersen))) --> 4
    @fact length(bin_decisions(min_dominating_set(petersen)))  --> 3
    @fact length(bin_decisions(min_vertex_cover(petersen)))    --> 6
    @fact distances(petersen, 1)[3]                            --> 2
    @fact length(bin_decisions(bipartite_set(k_3_3)))          --> 3
end

# write your own tests here
#@test 1 == 1


text:
using ITensors, ITensorNetworkAD
using ITensorNetworkAD.ITensorNetworks: PEPS, Models
using ITensorNetworkAD.Optimizations: gd_error_tracker

function test_peps_gd()
  Ny, Nx = 3, 4
  num_sweeps = 20
  cutoff = 1e-15
  maxdim = 1000
  sites = siteinds("S=1/2", Ny, Nx)
  peps = PEPS(sites; linkdims=2)
  randn!(peps)
  H_line = Models.lineham(Models.Model("tfim"), sites; h=1.0)
  return gd_error_tracker(
    peps, H_line; stepsize=0.005, num_sweeps=num_sweeps, cutoff=cutoff, maxdim=maxdim
  )
end

test_peps_gd()


text:
# check that all interace functions are provided

# --- LazySet ---

# support vector
@test check_method_implementation(LazySet, σ,
                                  Function[S -> (Vector{Float64}, S{Float64})])
# dimension
check_method_implementation(LazySet, dim, Function[S -> (S{Float64},)])

# --- AbstractPolyhedron ---

# constraints list
@test check_method_implementation(AbstractPolyhedron, constraints_list,
                                  Function[S -> (S{Float64},)])

# --- AbstractPolytope ---

# vertices list
@test check_method_implementation(AbstractPolytope, vertices_list,
                                  Function[S -> (S{Float64},)])

# --- AbstractCentrallySymmetric ---

# center
@test check_method_implementation(AbstractCentrallySymmetric, center,
                                  Function[S -> (S{Float64},)])

# --- AbstractCentrallySymmetricPolytope ---
# (copies methods from AbstractCentrallySymmetric)

# center (from AbstractCentrallySymmetric)
@test check_method_implementation(AbstractCentrallySymmetricPolytope, center,
                                  Function[S -> (S{Float64},)])

# --- AbstractZonotope ---

# genmat
@test check_method_implementation(AbstractZonotope, genmat,
                                  Function[S -> (S{Float64},)])
# generators
@test check_method_implementation(AbstractZonotope, generators,
                                  Function[S -> (S{Float64},)])

# --- AbstractHyperrectangle ---

# radius_hyperrectangle (2x)
@test check_method_implementation(AbstractHyperrectangle, radius_hyperrectangle,
                                  Function[S -> (S{Float64},)])
@test check_method_implementation(AbstractHyperrectangle, radius_hyperrectangle,
                                  Function[S -> (S{Float64}, Int)])
# isflat
@test check_method_implementation(AbstractHyperrectangle, isflat,
                                  Function[S -> (S{Float64},)])

# --- AbstractSingleton ---

# element (2x)
@test check_method_implementation(AbstractSingleton, element,
                                  Function[S -> (S{Float64},)])
@test check_method_implementation(AbstractSingleton, element,
                                  Function[S -> (S{Float64}, Int)])

# --- AbstractPolygon ---

# tohrep
@test check_method_implementation(AbstractPolygon, tohrep,
                                  Function[S -> (S{Float64},)])
# tovrep
@test check_method_implementation(AbstractPolygon, tovrep,
                                  Function[S -> (S{Float64},)])

# --- AbstractHPolygon ---

# no functions yet


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule openPMD_api_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("openPMD_api")
JLLWrappers.@generate_main_file("openPMD_api", UUID("0cbd5639-4ed4-5caa-9af9-db1a1b49a12d"))
end  # module openPMD_api_jll


text:
# ======================================================================================
# ======================================================================================
# ======================================================================================
function Copy(primary::Dict, secondary::Dict, attribute::String)
      if !haskey(primary, attribute) && haskey(secondary, attribute)
          primary[attribute] = secondary[attribute]
      end
end
# ======================================================================================
function CopyLoop(primary, secondary, attributelist)
    for attribute in attributelist
        if !haskey(primary, attribute) && haskey(secondary, attribute)
            primary[attribute] = secondary[attribute]
        end
    end
end


text:
export FewShotDataLoader

using Serialization
import Base: size
import StatsBase: sample, Weights

struct FewShotDataLoader
    labels::AbstractArray{L,1} where {L}
    samples::AbstractArray
    label2Idices::Dict
    dims_feature::Integer

    function FewShotDataLoader(path::String)
        data = deserialize(path)
        @assert typeof(data) <: Dict "$path should containing a Dict, get $(typeof(data))"
        labels = data["labels"]
        samples = data["data"]
        @assert ndims(samples) >= 2 "ndims of samples should be at least 2, get $(ndims(samples))"
        N = size(labels)[end]
        label2Idices = Dict{eltype(labels),Array{Integer,1}}()
        for i = 1:N
            label = labels[i]
            if !haskey(label2Idices, label)
                label2Idices[label] = []
            end
            push!(label2Idices[label], i)
        end
        dims_feature = reduce((a, b) -> a * b, size(samples)[1:end-1])
        new(labels, samples, label2Idices, dims_feature)
    end
end

"""
Fields
======

`support_samples`: array of support samples with size (`feature_dim`, `num_support`, `num_tasks`)
`support_labels`: array of labels of support samples with size (`num_support`, `num_tasks`)
`query_samples`: array of query samples with size (`feature_dim`, `num_query`, `num_tasks`)
`query_labels`: array of labels of query samples with size (`num_query`, `num_tasks`)
"""
struct MetaDataSample{T,L}
    support_n_ways::Integer
    support_k_shots::Integer
    query_n_ways::Integer
    query_k_shots::Integer
    support_samples::AbstractArray{T,3}
    support_labels::AbstractArray{L,2}
    query_samples::AbstractArray{T,3}
    query_labels::AbstractArray{L,2}
end

function MetaDataSample(;
    support_n_ways,
    support_k_shots,
    query_n_ways,
    query_k_shots,
    support_samples,
    support_labels,
    query_samples,
    query_labels,
)
    @assert support_n_ways > 0 "support_n_ways must be positive: $support_n_ways"
    @assert support_k_shots > 0 "support_k_shots must be positive: $support_k_shots"
    @assert query_n_ways > 0 "query_n_ways must be positive: $query_n_ways"
    @assert query_k_shots > 0 "query_k_shots must be positive: $query_k_shots"
    @assert ndims(support_samples) >= 2 "support_samples must be at least 2-dim"
    @assert ndims(query_samples) >= 2 "query_samples must be at least 2-dim"
    @assert size(support_samples) == size(query_samples) "size of support samples and query samples should be the same"
    n_tasks = size(support_samples)[end]
    if ndims(support_samples) > 2
        support_samples =
            reshape(support_samples, :, support_n_ways * support_k_shots, n_tasks)
        query_samples = reshape(query_samples, :, query_n_ways * query_k_shots, n_tasks)
    end
    for i = 1:n_tasks
        support_label = support_labels[:, i]
        query_label = query_labels[:, i]
        @assert isempty(intersect(support_label, query_label)) "support_labels and query_labels is not disjoint at task $i"
    end
    MetaDataSample(
        support_n_ways,
        support_k_shots,
        query_n_ways,
        query_k_shots,
        support_samples,
        support_labels,
        query_samples,
        query_labels,
    )
end

function size(meta_data_sample::MetaDataSample)
    return size(meta_data_sample.support_samples)[end]
end

function _sample(
    dloader::FewShotDataLoader;
    support_n_ways = 2,
    support_k_shots = 5,
    query_n_ways = support_n_ways,
    query_k_shots = support_k_shots,
)
    # find exclusive train/test samples
    uniq_labels = unique(dloader.labels)
    idx_map = Dict([(k, i) for (i, k) in enumerate(uniq_labels)])
    support_target_labels = sample(uniq_labels, support_n_ways, replace = false)
    idxs = [idx_map[label] for label in support_target_labels]
    weight_vs = ones(size(uniq_labels, 1))
    weight_vs[idxs] .= 0
    query_target_labels =
        sample(uniq_labels, Weights(weight_vs), query_n_ways, replace = false)
    support_candidate_idxs =
        vcat([dloader.label2Idices[label] for label in support_target_labels]...)
    query_candidate_idxs =
        vcat([dloader.label2Idices[label] for label in query_target_labels]...)
    train_idxs =
        sample(support_candidate_idxs, support_n_ways * support_k_shots, replace = false)
    test_idxs = sample(query_candidate_idxs, query_n_ways * query_k_shots, replace = false)
    support_samples = reshape(
        dloader.samples[repeat([:], ndims(dloader.samples) - 1)..., train_idxs],
        dloader.dims_feature,
        support_n_ways * support_k_shots,
        1,
    )
    support_labels = Utils.add_dim(dloader.labels[train_idxs,])
    query_samples = reshape(
        dloader.samples[repeat([:], ndims(dloader.samples) - 1)..., test_idxs],
        dloader.dims_feature,
        query_n_ways * query_k_shots,
        1,
    )
    query_labels = Utils.add_dim(dloader.labels[test_idxs])
    return support_samples, support_labels, query_samples, query_labels
end

function sample(
    dloader::FewShotDataLoader;
    support_n_ways = 2,
    support_k_shots = 5,
    query_n_ways = support_n_ways,
    query_k_shots = support_k_shots,
)
    support_samples, support_labels, query_samples, query_labels = _sample(
        dloader,
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_n_ways = query_n_ways,
        query_k_shots = query_k_shots,
    )
    MetaDataSample(
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_n_ways = query_n_ways,
        query_k_shots = query_k_shots,
        support_samples = support_samples,
        support_labels = support_labels,
        query_samples = query_samples,
        query_labels = query_labels,
    )
end

function sample(
    dloader::FewShotDataLoader,
    n_tasks::Integer;
    support_n_ways = 2,
    support_k_shots = 5,
    query_n_ways = support_n_ways,
    query_k_shots = support_k_shots,
)
    support_sets = []
    support_labels = []
    query_sets = []
    query_labels = []
    for i = 1:n_tasks
        support, support_label, query, query_label = _sample(
            dloader,
            support_n_ways = support_n_ways,
            support_k_shots = support_k_shots,
            query_n_ways = query_n_ways,
            query_k_shots = query_k_shots,
        )
        push!(support_sets, support)
        push!(support_labels, support_label)
        push!(query_sets, query)
        push!(query_labels, query_label)
    end
    MetaDataSample(
        support_n_ways = support_n_ways,
        support_k_shots = support_k_shots,
        query_n_ways = query_n_ways,
        query_k_shots = query_k_shots,
        support_samples = cat(support_sets..., dims = Val(3)),
        support_labels = cat(support_labels..., dims = Val(2)),
        query_samples = cat(query_sets..., dims = Val(3)),
        query_labels = cat(query_labels..., dims = Val(2)),
    )
end


text:
"""
    TDist(ν)

The *Students T distribution* with `ν` degrees of freedom has probability density function

```math
f(x; d) = \\frac{1}{\\sqrt{d} B(1/2, d/2)}
\\left( 1 + \\frac{x^2}{d} \\right)^{-\\frac{d + 1}{2}}
```

```julia
TDist(d)      # t-distribution with d degrees of freedom

params(d)     # Get the parameters, i.e. (d,)
dof(d)        # Get the degrees of freedom, i.e. d
```

External links

[Student's T distribution on Wikipedia](https://en.wikipedia.org/wiki/Student%27s_t-distribution)

"""
struct TDist{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    TDist{T}(ν::T) where {T <: Real} = new{T}(ν)
end

function TDist(ν::T) where {T <: Real}
    @check_args(TDist, ν > zero(ν))
    return TDist{T}(ν)
end

TDist(ν::T, ::NoArgCheck) where {T<:Real} = TDist{T}(ν)
TDist(ν::Integer) = TDist(float(ν))

@distr_support TDist -Inf Inf

#### Conversions
convert(::Type{TDist{T}}, ν::Real) where {T<:Real} = TDist(T(ν))
convert(::Type{TDist{T}}, d::TDist{S}) where {T<:Real, S<:Real} = TDist(T(d.ν), NoArgCheck())

#### Parameters

dof(d::TDist) = d.ν
params(d::TDist) = (d.ν,)
@inline partype(d::TDist{T}) where {T<:Real} = T

#### Statistics

mean(d::TDist{T}) where {T<:Real} = d.ν > 1 ? zero(T) : T(NaN)
median(d::TDist{T}) where {T<:Real} = zero(T)
mode(d::TDist{T}) where {T<:Real} = zero(T)

function var(d::TDist{T}) where T<:Real
    ν = d.ν
    isinf(ν) && return one(T)
    ν > 2 ? ν / (ν - 2) :
    ν > 1 ? T(Inf) : T(NaN)
end

skewness(d::TDist{T}) where {T<:Real} = d.ν > 3 ? zero(T) : T(NaN)

function kurtosis(d::TDist{T}) where T<:Real
    ν = d.ν
    ν > 4 ? 6 / (ν - 4) :
    ν > 2 ? T(Inf) : T(NaN)
end

function entropy(d::TDist{T}) where T <: Real
    isinf(d.ν) && return entropy( Normal(zero(T), one(T)) )
    h = d.ν/2
    h1 = h + 1//2
    h1 * (digamma(h1) - digamma(h)) + log(d.ν)/2 + lbeta(h, 1//2)
end

#### Evaluation & Sampling

@_delegate_statsfuns TDist tdist ν

rand(rng::AbstractRNG, d::TDist) = randn(rng) / ( isinf(d.ν) ? 1 : sqrt(rand(rng, Chisq(d.ν))/d.ν) )

function cf(d::TDist{T}, t::Real) where T <: Real
    isinf(d.ν) && return cf(Normal(zero(T), one(T)), t)
    t == 0 && return complex(1)
    h = d.ν/2
    q = d.ν/4
    complex(2(q*t^2)^q * besselk(h, sqrt(d.ν) * abs(t)) / gamma(h))
end

gradlogpdf(d::TDist{T}, x::Real) where {T<:Real} = isinf(d.ν) ? gradlogpdf(Normal(zero(T), one(T)), x) : -((d.ν + 1) * x) / (x^2 + d.ν)


text:
#### PlotFuncs

# Provides a set of plot functions that operate on StateVec.

using Requires

@init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
  using .Plots
  export plot_state, plot_states
end

export nice_string

k_start = 1
k_stop_min = 1000
markershapes = Symbol[:utriangle, :circle, :rect, :star5,
                      :diamond, :hexagon, :cross, :xcross,
                      :dtriangle, :rtriangle, :ltriangle,
                      :pentagon, :heptagon, :octagon,
                      :star4, :star6, :star7, :star8,
                      :vline, :hline, :+, :x]
markershapes_net = Symbol[:hline, :+, :x]
markersize = [6, 4, 2]

function nice_string(name)
  friendly_name = string(name)
  friendly_name = replace(friendly_name, "θ" => "theta")
  friendly_name = replace(friendly_name, "ρ" => "rho")
  friendly_name = replace(friendly_name, "α" => "alpha")
  friendly_name = replace(friendly_name, "∇" => "grad")
  friendly_name = replace(friendly_name, "εδ" => "entr-detr")
  return friendly_name
end

@init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin

"""
    plot_state(sv::StateVec,
                    grid::Grid,
                    name_s::Symbol = nothing,
                    directory::AbstractString,
                    i = 0,
                    include_ghost = false,
                    xlims::Union{Nothing, Tuple{R, R}} = nothing,
                    ylims::Union{Nothing, Tuple{R, R}} = nothing
                    ) where R

Save the plot variable along the z-direction in `StateVec` given the
grid, `grid`, variable name `name_s`, directory `directory`,
sub-domain `i`, and a `Bool`, `include_ghost`, indicating
to include include or exclude the ghost points.
"""
function plot_state(sv::StateVec,
                    grid::Grid,
                    directory::AbstractString,
                    name_s::Symbol;
                    i = 0,
                    include_ghost = true,
                    filename = nothing,
                    i_Δt = 1,
                    xlims::Union{Nothing, Tuple{R, R}} = nothing,
                    ylims::Union{Nothing, Tuple{R, R}} = nothing
                    ) where R
  domain_range = include_ghost ? over_elems(grid) : over_elems_real(grid)
  k_stop_local = min(length(domain_range), k_stop_min)
  domain_range = domain_range[k_start:k_stop_local]
  x = [grid.zc[k] for k in domain_range]
  name = nice_string(name_s)
  if i==0
    plot()
    for i in over_sub_domains(sv, name_s)
      name_i = nice_string(var_string(sv, name_s, i))
      y = [sv[name_s, k, i] for k in domain_range]
      plot!(y, x, markershapes = markershapes[i], label = name_i)
    end
    plot!(title = name * " vs z", xlabel = name, ylabel = "z")
  else
    name = nice_string(var_string(sv, name_s, i))
    y = [sv[name_s, k, i] for k in domain_range]
    plot(y, x, markershapes = markershapes[i], label = name)
    plot!(title = name * " vs z", xlabel = name, ylabel = "z")
  end
  if xlims != nothing; plot!(xlims = xlims); end
  if ylims != nothing; plot!(ylims = xlims); end
  filename == nothing && (filename = name)
  mkpath(directory)
  savefig(joinpath(directory, filename))
end

"""
    plot_states(sv::StateVec,
                grid::Grid,
                name_ids,
                directory::AbstractString,
                include_ghost = false,
                xlims::Union{Nothing, Tuple{R, R}} = nothing,
                ylims::Union{Nothing, Tuple{R, R}} = nothing
                ) where R

Save the plot variable along the z-direction in `StateVec` given the
grid, `grid`, variable name `name_s`, directory `directory`,
sub-domain `i`, and a `Bool`, `include_ghost`, indicating to
include include or exclude the ghost points.
"""
function plot_states(sv::StateVec,
                     grid::Grid,
                     directory::AbstractString,
                     name_ids;
                     include_ghost = true,
                     filename = nothing,
                     i_Δt = 1,
                     sources = false,
                     xlims::Union{Nothing, Tuple{R, R}} = nothing,
                     ylims::Union{Nothing, Tuple{R, R}} = nothing
                     ) where R
  if sources
    domain_range = over_elems_real(grid)
  else
    domain_range = include_ghost ? over_elems(grid) : over_elems_real(grid)
  end
  k_stop_local = min(length(domain_range), k_stop_min)
  domain_range = domain_range[k_start:k_stop_local]

  x = [grid.zc[k] for k in domain_range]
  plot()
  gm, en, ud, sd, al = allcombinations(DomainIdx(sv))
  if sources
    @inbounds for name_id in name_ids
      @inbounds for i in sd
        source_term_names = unique([s.name for k in domain_range for s in sv[name_id, k, i]])
        data = [Dict(s.name => s.value for s in sv[name_id, k, i]) for k in domain_range]
        @inbounds for s in source_term_names
          y = [x[s] for x in data]
          if any([abs(v) > eps(typeof(v)) for v in y])
            plot!(y, x, label = nice_string(s)*var_suffix(sv, name_id, i), markershapes = markershapes[i], markersize = markersize[i])
          end
        end
        y = [sum([x[s] for s in source_term_names]) for x in data]
        if any([abs(v) > eps(typeof(v)) for v in y])
          plot!(y, x, label = "net_"*nice_string(var_string(sv, name_id, i)), markershapes = markershapes_net[i], markersize = markersize[i])
        end
      end
    end
  else
    @inbounds for name_id in name_ids
      @inbounds for i in over_sub_domains(sv, name_id)
        y = [sv[name_id, k, i] for k in domain_range]
        plot!(y, x, label = nice_string(s)*var_suffix(sv, name_id, i), markershapes = markershapes[i], markersize = markersize[i])
      end
    end
  end
  filename == nothing && (filename = nice_string(name_ids[1]))
  plot!(title = nice_string(filename)*" vs z", xlabel = nice_string(filename), ylabel = "z")
  if xlims != nothing; plot!(xlims = xlims); end
  if ylims != nothing; plot!(ylims = xlims); end
  mkpath(directory)
  savefig(joinpath(directory, filename))
end

end # @require

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
export
ruleSPDirichletOutNP,
ruleVBDirichletOut,
ruleVBDirichletIn1

ruleSPDirichletOutNP(msg_out::Nothing, msg_a::Message{PointMass, V}) where V<:VariateType = Message(V, Dirichlet, a=deepcopy(msg_a.dist.params[:m]))

ruleVBDirichletOut(marg_out::Any, marg_a::ProbabilityDistribution{V}) where V<:VariateType = Message(V, Dirichlet, a=unsafeMean(marg_a))

function ruleVBDirichletIn1(marg_out::ProbabilityDistribution{Multivariate}, marg_a::Any)
    Message(Multivariate, Function, log_pdf = (a)->sum((a.-1).*unsafeLogMean(marg_out)) - sum(loggamma.(a)) + loggamma(sum(a)))
end


text:
###
# Creating an archive
#
# To create an archive in C:
#     1. Ask `archive_write_new` for an archive writer object.
#     2. Set any global properties. In particular, you should set
#        the compression and format to use.
#     3. Call `archive_write_open` to open the file (most people
#         will use `archive_write_open_file` or `archive_write_open_fd`,
#         which provide convenient canned I/O callbacks for you).
#     4. For each entry:
#         * construct an appropriate `struct archive_entry` structure
#         * `archive_write_header` to write the header
#         * `archive_write_data` to write the entry data
#     5. `archive_write_close` to close the output
#     6. `archive_write_free` to cleanup the writer and release resources
#
# To create an archive in Julia:
#     1. Create a LibArchive.Writer
#     2. Set any global properties. In particular, you should set
#        the compression and format to use.
#     3. Write each entries. LibArchive.jl will call libarchive
#        open functions if needed.
#     4. Call LibArchive.free to end processing.

@compat abstract type WriterData end

type Writer{T<:WriterData} <: Archive
    data::T
    ptr::Ptr{Void}
    opened::Bool
    @static if isdefined(Base, :UnionAll)
        function Writer{T}(data::T) where
            T
            ptr = ccall((:archive_write_new, libarchive), Ptr{Void}, ())
            ptr == C_NULL && throw(OutOfMemoryError())
            obj = new(data, ptr, false)
            finalizer(obj, free)
            obj
        end
    else
        function Writer(data::T)
            ptr = ccall((:archive_write_new, libarchive), Ptr{Void}, ())
            ptr == C_NULL && throw(OutOfMemoryError())
            obj = new(data, ptr, false)
            finalizer(obj, free)
            obj
        end
    end
end

Writer{T<:WriterData}(data::T) = Writer{T}(data)

function free(archive::Writer)
    ptr = archive.ptr
    ptr == C_NULL && return
    ccall((:archive_write_free, libarchive), Cint, (Ptr{Void},), ptr)
    archive.ptr = C_NULL
    nothing
end

Base.close(archive::Writer) =
    @_la_call(archive_write_close, (Ptr{Void},), archive)

function Base.cconvert(::Type{Ptr{Void}}, archive::Writer)
    archive.ptr == C_NULL && error("archive already freed")
    archive
end
Base.unsafe_convert(::Type{Ptr{Void}}, archive::Writer) = archive.ptr

function ensure_open(archive::Writer)
    archive.opened && return
    archive.opened = true
    try
        do_open(archive)
    catch
        archive.opened = false
        rethrow()
    end
    nothing
end

set_bytes_per_block(archive::Writer, bytes_per_block) =
    @_la_call(archive_write_set_bytes_per_block, (Ptr{Void}, Cint),
              archive, bytes_per_block)
get_bytes_per_block(archive::Writer) =
    ccall((:archive_write_get_bytes_per_block, libarchive),
          Cint, (Ptr{Void},), archive)

set_bytes_in_last_block(archive::Writer, bytes_in_last_block) =
    @_la_call(archive_write_set_bytes_in_last_block, (Ptr{Void}, Cint),
              archive, bytes_in_last_block)
get_bytes_in_last_block(archive::Writer) =
    ccall((:archive_write_get_bytes_in_last_block, libarchive),
          Cint, (Ptr{Void},), archive)

"""
The dev/ino of a file that won't be archived. This is used
to avoid recursively adding an archive to itself.
"""
set_skip_file(archive::Writer, dev, ino) =
    @_la_call(archive_write_set_skip_file, (Ptr{Void}, Int64, Int64),
              archive, dev, ino)

###
# Open filename

immutable WriteFileName{T} <: WriterData
    name::T
end

Writer(fname::AbstractString) = Writer(WriteFileName(fname))
Writer() = Writer(WriteFileName(nothing))

function do_open{T}(archive::Writer{WriteFileName{T}})
    data = archive.data
    @_la_call(archive_write_open_filename, (Ptr{Void}, Cstring),
              archive, data.name)
end

function do_open(archive::Writer{WriteFileName{Void}})
    @_la_call(archive_write_open_filename, (Ptr{Void}, Ptr{Void}),
              archive, C_NULL)
end

immutable WriteFD <: WriterData
    fd::Cint
end

Writer{T<:Integer}(fd::T) = Writer(WriteFD(Cint(fd)))

function do_open(archive::Writer{WriteFD})
    data = archive.data
    @_la_call(archive_write_open_fd, (Ptr{Void}, Cint), archive, data.fd)
end

###
# Open memory

immutable WriteMemory{T,TO} <: WriterData
    ptr::T
    obj::TO # Reference for GC
    size::Csize_t
    used::typeof(Ref{Csize_t}())
end

Writer{T<:Ptr,TO}(ptr::T, size, obj::TO=nothing) =
    Writer(WriteMemory{T,TO}(ptr, obj, size, Ref(Csize_t(0))))
Writer(ary::Vector, size=sizeof(ary)) = Writer(pointer(ary), size, ary)

Base.unsafe_convert(::Type{Ptr{Void}}, mem_data::WriteMemory) = mem_data.ptr

function do_open{T<:WriteMemory}(archive::Writer{T})
    data = archive.data
    @_la_call(archive_write_open_memory,
              (Ptr{Void}, Ptr{Void}, Csize_t, Ptr{Csize_t}),
              archive, data, data.size, data.used)
end

get_used{T<:WriteMemory}(archive::Writer{T}) = archive.data.used[]

###
# Generic writer

immutable GenericWriteData{T} <: WriterData
    data::T
end

Writer{T}(data::T) = Writer(GenericWriteData{T}(data))

writer_open(archive::Writer, data) = nothing
function writer_writebytes end
writer_close(archive::Writer, data) = nothing

function writer_open_callback(c_archive, archive)
    try
        clear_error(archive)
        writer_open(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

function writer_write_callback(c_archive, archive, buff::Ptr{Void},
                               length::Csize_t)
    try
        clear_error(archive)
        ary = unsafe_wrap(Array, Ptr{UInt8}(buff), length)
        return Cssize_t(writer_writebytes(archive, archive.data.data, ary))
    catch ex
        set_exception(archive, ex)
        return Cssize_t(-1)
    end
end

function writer_close_callback(c_archive, archive)
    try
        clear_error(archive)
        writer_close(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

writer_writebytes(archive, io::IO, ary) = write(io, ary)

function do_open{T<:GenericWriteData}(archive::Writer{T})
    @_la_call(archive_write_open,
              (Ptr{Void}, Any, Ptr{Void}, Ptr{Void}, Ptr{Void}),
              archive, archive,
              to_open_callback(writer_open_callback, Writer{T}),
              to_write_callback(writer_write_callback, Writer{T}),
              to_close_callback(writer_close_callback, Writer{T}))
end

function write_header(archive::Writer, entry::Entry)
    ensure_open(archive)
    @_la_call(archive_write_header, (Ptr{Void}, Ptr{Void}), archive, entry)
end

@inline function unsafe_archive_write(archive::Writer, buf, len)
    ensure_open(archive)
    nb = ccall((:archive_write_data, libarchive),
               Cssize_t, (Ptr{Void}, Ptr{Void}, Csize_t), archive, buf, len)
    nb < 0 && _la_error(Cint(nb), archive)
    nb % Csize_t
end

unsafe_write(archive::Writer, buf::Ptr{UInt8}, len::UInt) =
    unsafe_archive_write(archive, buf, len)
Base.write(archive::Writer, c::UInt8) =
    unsafe_archive_write(archive, Ref(c), 1)

if !isdefined(Base, :unsafe_write)
    Base.write(archive::Writer, data::Array) =
        unsafe_archive_write(archive, data, sizeof(data))
end

function finish_entry(archive::Writer)
    ensure_open(archive)
    @_la_call(archive_write_finish_entry, (Ptr{Void},), archive)
end

"""
Marks the archive as FATAL so that a subsequent free() operation
won't try to close() cleanly.  Provides a fast abort capability
when the client discovers that things have gone wrong.
"""
write_fail(archive::Writer) =
    @_la_call(archive_write_fail, (Ptr{Void},), archive)

Writer(f::Function) = archive_guard(f, Writer())
Writer(f::Function, args...; kws...) =
    archive_guard(f, Writer(args...; kws...))

@deprecate file_writer(fname::AbstractString) Writer(fname)
@deprecate file_writer() Writer()
@deprecate file_writer{T<:Integer}(fd::T) Writer(fd)
@deprecate mem_writer(data::Vector) Writer(data)
@deprecate mem_writer(obj) Writer(pointer(obj), sizeof(obj), obj)
@deprecate gen_writer(data) Writer(data)
@deprecate write_data(archive::Writer, data) write(archive, data)

# /*
#  * Set write options.
#  */
# /* Apply option to the format only. */
# int archive_write_set_format_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option to the filter only. */
# int archive_write_set_filter_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option to both the format and the filter. */
# int archive_write_set_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option string to both the format and the filter. */
# int archive_write_set_options(struct archive *_a,
# 			    const char *opts);

# /* This interface is currently only available for archive_write_disk handles.  */
# ssize_t	 archive_write_data_block(struct archive *,
# 				    const void *, size_t, int64_t);


text:
"""
Utilities to find and plot the biggest changing parameters
"""

"""
    find parameter indices of the biggest changing parametesr in the curve
"""
function biggest_movers(mdc::AbstractCurveSolution, num::Integer; rev=false)
    diff = trajectory(mdc)[:,end] - trajectory(mdc)[:,1]
    ids = sortperm(diff, by=abs, rev=!rev)
    ids = ids[1:num]
end


text:
export Beta

import Distributions: Beta, params
import SpecialFunctions: digamma

vague(::Type{ <: Beta }) = Beta(1.0, 1.0)

prod_analytical_rule(::Type{ <: Beta }, ::Type{ <: Beta }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::Beta, right::Beta)
    left_a, left_b   = params(left)
    right_a, right_b = params(right)
    T                = promote_type(eltype(left), eltype(right))
    return Beta(left_a + right_a - one(T), left_b + right_b - one(T))
end

function logmean(dist::Beta)
    a, b = params(dist)
    return digamma(a) - digamma(a + b)
end

function mirroredlogmean(dist::Beta)
    a, b = params(dist)
    return digamma(b) - digamma(a + b)
end