text:
const MACHEP = 1.2e-16
const DWARF = 1.0e-38

function jacFcnDefault!{I <: Integer}(m::I, n::I, x::Vector{Float64}, fvec::Vector{Float64}, fjac::Matrix{Float64}, eps_::Float64, wa::Vector{Float64}, fcn!::Function)
  for j = 1:n
    temp = x[j]
    step = max(eps_ * eps_, eps_ * abs(temp))

    x[j] += step
    fcn!(m, n, x, wa)
    for i = 1:m
      fjac[i, j] = (wa[i] - fvec[i]) / step
    end

    #restore
    x[j] = temp
  end
  return fjac
end

function lm_qrsolv!{I <: Integer}(n::I, r::Matrix{Float64}, ldr::I, pivot::Vector{I}, diagonal::Vector{Float64}, qtb::Vector{Float64}, x::Vector{Float64}, sdiag::Vector{Float64}, W::Vector{Float64})

  # Copy R and Q'* b to preserve input and initialize S
  # In particular, save the diagonal elements of R in x

  # instead of loop...
  # r = triu(r) + triu(r)'
  # x = diag(r)
  # W = qtb

  for j = 1:n
    for i = j:n
      r[(j-1) * ldr + i] = r[(i-1)*ldr + j]
    end

    x[j] = r[(j-1) * ldr + j]
    W[j] = qtb[j]
  end

  # eliminate the diagonal matrix D using a Givens rotation
  for j = 1:n
    if diagonal[pivot[j]] != 0.0
      for k=j:n
        sdiag[k] = 0.0
      end
      sdiag[j] = diagonal[pivot[j]]

      # the transformation to elminate the row of D modify only a single element of Q' * b beyond the first n, which is initially 0
      qtbpj = 0.0
      for k = j:n
        # determine the givens rotation which elminates the appropriate element in the current row of D
        if sdiag[k] == 0.0
          continue
        end

        kk = k + ldr * (k - 1)
        if abs(r[kk]) < abs(sdiag[k])
          _cot = r[kk] / sdiag[k]
          _sin = 1.0 / hypot(1.0, _cot)
          _cos = _sin * _cot
        else
          _tan = sdiag[k] / r[kk]
          _cos = 1.0 / hypot(1.0, _tan)
          _sin = _cos * _tan
        end

        # compute the modified diagonal element of R and the modified element of Q' * b, 0
        r[kk] = _cos * r[kk] + _sin * sdiag[k]
        temp = _cos * W[k] + _sin * qtbpj
        qtbpj = -_sin * W[k] + _cos * qtbpj
        W[k] = temp

        # accumulate the transformation in the row of S
        for i = k+1:n
          temp = _cos * r[(k - 1) * ldr + i] + _sin * sdiag[i]
          sdiag[i] = -_sin * r[(k-1) * ldr + i] + _cos * sdiag[i]
          r[(k-1) * ldr + i] = temp
        end
      end
    end

    sdiag[j] = r[(j -1) * ldr + j]
    r[(j-1) * ldr + j] = x[j]
  end

  # Solve the triangular system for z.  If the system is singular, then obtain a least-squares solution
  nsing = n
  for j=1:n
    if sdiag[j] == 0.0 && nsing ==n
      nsing = j - 1
    end

    if nsing < n
      W[j] = 0.0
    end
  end

  for j = nsing:-1:1
    _sum = 0.0
    for i = j+1:nsing # check this
      _sum += r[(j-1) * ldr + i] * W[i]
    end

    W[j] = (W[j] - _sum) / sdiag[j]
  end

  # permute the components of z back to the components of x
  for j=1:n
    x[pivot[j]] = W[j]
  end
end

function lm_lmpar!{I <: Integer}(n::I, r::Matrix{Float64}, ldr::I, pivot::Vector{I}, diagonal::Vector{Float64}, qtb::Vector{Float64}, delta::Float64, par::Float64,
                  x::Vector{Float64}, sdiag::Vector{Float64}, aux::Vector{Float64}, xdi::Vector{Float64})

  # Compute and store in x the Gauss-Newton direction.  If the jacobian is rank-deficient, obtain a least-squares solution

  # println("here i am")

  nsing = n
  for j = 1:n
    aux[j] = qtb[j]
    if r[(j-1) * ldr + j] == 0.0 && nsing == n
      nsing = j - 1
    end

    if nsing < n
      aux[j] = 0.0
    end
  end

  for j = nsing:-1:1
    aux[j] = aux[j] / r[j + ldr * (j-1)]
    temp = aux[j]
    for i = 1:j - 1  ## CHECK THIS
      aux[i] -= r[(j-1) * ldr + i] * temp
    end
  end

  for j = 1:n
    x[pivot[j]] = aux[j]
  end

  # Initialize the iteration counter, evaluate the function at the origin,
  # and test for acceptance of the Gauss-Newton direction
  xdi = diagonal .* x

  dxnorm = vecnorm(xdi)
  fp = dxnorm - delta
  if fp <= 0.1 * delta
    par = 0.0
    return par
  end

  # if the jacobian is not rank deficient, the newton step provides a lower bound, parl, for the zero of the function.  Otherwise,
  # set this bound to zero

  parl = 0.0
  if nsing >= n
    for j =1:n
      aux[j] = diagonal[pivot[j]] * xdi[pivot[j]] / dxnorm
    end

    for j = 1:n
      _sum = 0.0
      for i = 1:j - 1 # CHECK THIS
        _sum += r[(j-1) * ldr + i] * aux[i]
      end

      aux[j] = (aux[j] - _sum) / r[j + ldr * (j-1)]
    end

    temp = vecnorm(aux)
    parl = fp / delta / temp / temp
  end

  # Calculate the upper bound, paru, for the zero of the function
  for j = 1:n
    _sum = 0.0
    for i = 1:j
      _sum += r[(j-1) * ldr + i] * qtb[i]
    end
    aux[j] = _sum / diagonal[pivot[j]]
  end

  # if the input par lies outside the interval (parl, paru), set par to the closer endpoint

  gnorm = vecnorm(aux)
  paru = gnorm / delta
  if paru == 0.0
    paru = DWARF / min(delta, 0.1)
  end

  # if the input par lies outside of the interval (parl, paru), set par to the closer endpoint

  par = max(par, parl)
  par = min(par, paru)
  if par == 0.0
    par = gnorm / dxnorm
  end

  # iterate
  for iter=1:10
    # evaluate the function at the current value of par
    if par == 0.0
      par = max(DWARF, 0.001 * paru)
    end

    temp = sqrt(par)
    aux = diagonal * temp

    lm_qrsolv!(n, r, ldr, pivot, aux, qtb, x, sdiag, xdi)

    xdi = diagonal .* x
    dxnorm = vecnorm(xdi)
    fp_old = fp
    fp = dxnorm - delta

    if abs(fp) <= 0.1 * delta || (parl == 0.0 && fp <= fp_old && fp_old < 0.0) || iter == 10
      break
    end

    # Compute the Newton correctly
    for j = 1:n
      aux[j] = diagonal[pivot[j]] * xdi[pivot[j]] / dxnorm
    end

    for j=1:n
      aux[j] = aux[j] / sdiag[j]
      for i=j+1:n
        aux[i] -= r[(j-1) * ldr + i] * aux[j]
      end
    end

    temp = vecnorm(aux)
    parc = fp / delta / temp / temp

    # depending on the sign of the function, update parl or paru
    if fp > 0
      parl = max(parl, par)
    else
      paru = min(paru, par)
    end

    par = max(parl, par * parc)
  end
  return par
end

function lmdif2!{I <: Integer}(n::I, m::I, x::Vector{Float64}, mode::I, factor_::Float64, info_::I, epsfcn::Float64, ftol::Float64, xtol::Float64, gtol::Float64,
                maxIter::I, fcn!::Function, jacFcn!::Function = jacFcnDefault!)
  converged = false
  x_converged = false
  g_converged = false
  iterCount = 0
  f_calls = 0
  j_calls = 0
  eps_ = sqrt(max(epsfcn, MACHEP))

  # arrays to create
  fvec = zeros(m)
  diagonal = ones(n)
  qtf = zeros(n)
  fjac = zeros(m, n)
  wa1 = zeros(n)
  wa2 = zeros(n)
  wa3 = zeros(n)
  wf = zeros(m)
  pivot = ones(Int, n)

  # eval the function at the starting point and calculate its norm
  fcn!(m, n, x, fvec)
  f_calls += 1

  fnorm = vecnorm(fvec)

  # initialize the l-m parameter
  par = 0.0
  iterCount += 1

  while ~converged
    # Calculate the jacobian
    jacFcn!(m, n, x, fvec, fjac, eps_, wf, fcn!)
    j_calls += 1

    # println(fjac)

    qr = qrfact!(fjac, Val{true})
    # println(qr)
    # println("=======")
    # println(qr[:Q])
    # println("=======")
    # println(qr[:R])

    q = full(qr[:Q])
    r = qr[:R]
    pivot = qr[:p]
    wa1 = diag(r) # rdiag
    wa2 = r'[:, 1] # first of each column

    # println("fvec: ", fvec)
    #
    # println("next: ", q' * fvec)

    # Form q' * fvec and store in qtf
    #fjac = q
    wf = fjac * fvec
    qtf[1:n] = q' * fvec

    # println("wa1: ", wa1)
    # println("wa2: ", wa2)
    # println("qtf: ", qtf)
    # println("wf: ", wf)
    # println("fjac: ", fjac)
    #
    # error("break")

    # Compute norm of scaled gradient and detect degeneracy
    gnorm = 0.0
    for j = 1:n
      if wa2[pivot[j]] == 0.0
        continue
      end
      _sum = 0.0
      for i = 1:j
        _sum += fjac[i, j] * qtf[i]
      end

      gnorm = max(gnorm, abs(_sum / wa2[pivot[j]] / fnorm))
    end

    # terminate if gnorm <= gtol
    if gnorm <= gtol
      info_ = 4
      converged = true
      break
    end

    if iterCount == 1
      diagonal = copy(wa2)
      wa3 = diagonal .* x
      xnorm = vecnorm(wa3)

      delta = factor_ * xnorm
    else
      diagonal = max(diagonal, wa2)
    end

    inner_success = false

    while ~inner_success
      par = lm_lmpar!(n, fjac, m, pivot, diagonal, qtf, delta, par, wa1, wa2, wf, wa3)

      pnorm = vecnorm(wa3)

      temp2 = par * (pnorm / fnorm)^2

      for j = 1:n
        wa3[j] = 0.0
        for i = 1:j
          wa3[i] -= fjac[i, j] * wa1[pivot[j]]
        end
      end

      temp1 = (vecnorm(wa3) / fnorm)^2

      prered = temp1 + 2.0 * temp2
      dirder = -temp1 + temp2

      # at first call, adjust the initial step bound
      if iterCount == 1 && pnorm < delta
        delta = pnorm
      end

      # Evaluate the function at x + p
      wa2 = x - wa1
      fcn!(m, n, wa2, wf)
      f_calls += 1

      fnorm1 = vecnorm(wf)

      # Evaluate the scaled reduction
      actred = 1.0 - (fnorm1 / fnorm)^2

      # ratio of actual to predicted reduction
      ratio = prered > 0.0 ? actred / prered : 0.0

      # Update the step bound
      if ratio <= 0.25
        if actred >= 0.0
          temp = 0.5
        elseif actred > -99.0
          temp = max(dirder / (2.0 * dirder + actred), 0.1)
        else
          temp = 0.1
        end

        delta = temp * min(delta, pnorm / 0.1)
        par /= temp
      elseif ratio >= 0.75
        delta = 2.0 * pnorm
        par *= 0.5
      elseif par == 0.0
        delta = 2.0 * pnorm
      end

      # on success, update solution and test for convergence
      inner_success = ratio >= 1e-4

      if inner_success
        x = copy(wa2)
        wa2 = diagonal .* x
        fvec = copy(wf)

        xnorm = vecnorm(wa2)
        fnorm = fnorm1

        # Convergence tests
        info_ = 0

        if fnorm <= DWARF
          info_= 10 # success: sum of squares is almost 0
          converged = true
          break
        end

        if abs(actred) <= ftol && prered <= ftol && ratio <= 2.0
          info_ = 1 # success: x is almost stable
        end

        if delta <= xtol * xnorm
          info_ += 2 # success: sum of squares is almost stable
        end

        if info_ != 0
          converged = true
          break
        end

        # Tests for termination and stringent tolerances
        if abs(actred) <= MACHEP && prered <= MACHEP && ratio <= 2.0
          info_ = 6
          converged = true
          break
        end

        if delta <= MACHEP * xnorm
          info_ = 7
          converged = true
          break
        end

        if gnorm <= MACHEP
          info_ = 8
          converged = true
          break
        end
      end # convergence tests
    end # inner loop
    iterCount += 1
  end # outer loop
  return par, info_, f_calls, j_calls
end


text:
"""
    correct_immersed_tendencies!(model, Δt, γⁿ, ζⁿ)

Change the tendency fields to account for the presence of a boundary immersed
within the `model` grid. Does nothing by default.
"""
correct_immersed_tendencies!(model, Δt, γⁿ, ζⁿ) = nothing # fallback function

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 1d24602c-a463-11ec-1b01-4b6d9133a279
using PlutoUI, Dates, SQLiteGraph, CSV, DataFrames, EasyConfig

# ╔═╡ 72412469-67d7-42e3-a7b2-ce05e860e302
md"""
# `reddit_notebook.jl`

- This notebook analyzes the data at [https://snap.stanford.edu/data/soc-RedditHyperlinks.html](https://snap.stanford.edu/data/soc-RedditHyperlinks.html)
"""

# ╔═╡ a78fca5f-fa6a-4929-8861-be95c13492fe
PlutoUI.TableOfContents()

# ╔═╡ a5e78f84-b1b4-4e43-afb3-29cd5ec0e8ff
md"## Load Packages"

# ╔═╡ a51e7694-8209-44e9-bff3-ed2fb607c069
md"## Load Data"

# ╔═╡ 439c9e48-f1b1-4d3f-b1ae-e2b8bbe13c84
url = "https://snap.stanford.edu/data/soc-redditHyperlinks-body.tsv"

# ╔═╡ 678cae6f-d2e8-4d20-85b3-0d5080b07a62
path = joinpath(@__DIR__, "reddit.tsv")

# ╔═╡ b38fb416-419d-40e0-9cb1-d16106969d07
rows = CSV.Rows(path, delim='\t');

# ╔═╡ cf6d196e-9f7e-444b-a88c-78e32bf88fba
first(rows).PROPERTIES

# ╔═╡ be695499-5529-4029-a326-26a9391b2069
function parse_properties(props)
	data = split(props, ',')
	keys = [
		"Number of characters"
		"Number of characters without counting white space"
		"Fraction of alphabetical characters"
		"Fraction of digits"
		"Fraction of uppercase characters"
		"Fraction of white spaces"
		"Fraction of special characters, such as comma, exclamation mark, etc."
		"Number of words"
		"Number of unique works"
		"Number of long words (at least 6 characters)"
		"Average word length"
		"Number of unique stopwords"
		"Fraction of stopwords"
		"Number of sentences"
		"Number of long sentences (at least 10 words)"
		"Average number of characters per sentence"
		"Average number of words per sentence"
		"Automated readability index"
		"Positive sentiment calculated by VADER"
		"Negative sentiment calculated by VADER"
		"Compound sentiment calculated by VADER"
	]
	for i in 1:length(keys)
		keys[i] = replace(keys[i], " " => "_")
	end
	Config(key => parse(Float64, x) for (key,x) in zip(keys, data))
end

# ╔═╡ a291bf9e-48b2-4673-9c0f-1df440ff2db3
df = let
	df = DataFrame(rows)
	df.TIMESTAMP = DateTime.(df.TIMESTAMP, dateformat"yyyy-mm-dd H:M:S")
	df.LINK_SENTIMENT = parse.(Int, df.LINK_SENTIMENT)
	df.PROPERTIES = parse_properties.(df.PROPERTIES)
	df
end

# ╔═╡ d99040c0-eb84-4b98-b4c7-971a3681e0a7
subreddits = sort!(unique!(vcat(df.SOURCE_SUBREDDIT, df.TARGET_SUBREDDIT)))

# ╔═╡ 8e6df98b-22ce-4c1e-bf45-5e154725586a
gb = groupby(df, [:SOURCE_SUBREDDIT, :TARGET_SUBREDDIT]);

# ╔═╡ 0763d60f-82e9-44af-a26b-931e6749d4af
df2 = combine(gb, :PROPERTIES => length => "N_LINKS")

# ╔═╡ 5c222537-f561-4b47-acf5-6731a05d16a9
md"## Make Graph"

# ╔═╡ 30450ece-69b6-45d4-b9ee-8b453cb0f93f
db = let
	db = DB()
	for (i,sub) in enumerate(subreddits)
		insert!(db, Node(i, "Subreddit"; name=sub))
	end
	for row in eachrow(df2)
		i = findfirst(==(row.SOURCE_SUBREDDIT), subreddits)
		j = findfirst(==(row.TARGET_SUBREDDIT), subreddits)
		insert!(db, Edge(i,j,"Has_Link", nlinks=row.N_LINKS))
	end
	db
end

# ╔═╡ 69819a25-1759-4001-a960-44f0dcf88ed2
db[200]

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
EasyConfig = "acab07b0-f158-46d4-8913-50acef6d41fe"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
SQLiteGraph = "de2defd0-d76e-464a-9029-1d71e199ae58"

[compat]
CSV = "~0.10.3"
DataFrames = "~1.3.2"
EasyConfig = "~0.1.10"
PlutoUI = "~0.7.37"
SQLiteGraph = "~0.2.0"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

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

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

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

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

[[deps.BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

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

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

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

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

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

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

[[deps.DBInterface]]
git-tree-sha1 = "9b0dc525a052b9269ccc5f7f04d5b3639c65bca5"
uuid = "a10d1c49-ce27-4219-8d33-6db1a4562965"
version = "2.5.0"

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

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

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

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

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

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

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

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

[[deps.EasyConfig]]
deps = ["JSON3", "OrderedCollections", "StructTypes"]
git-tree-sha1 = "c070b3c48a8ba3c6e6507997f0a7f5ebf85c3600"
uuid = "acab07b0-f158-46d4-8913-50acef6d41fe"
version = "0.1.10"

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "8c1f668b24d999fb47baf80436194fdccec65ad2"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.4"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

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

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

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

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

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

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

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

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

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

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

[[deps.SQLite]]
deps = ["BinaryProvider", "DBInterface", "Dates", "Libdl", "Random", "SQLite_jll", "Serialization", "Tables", "Test", "WeakRefStrings"]
git-tree-sha1 = "8e14d9b200b975e93a0ae0e5d17dea1c262690ee"
uuid = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
version = "1.4.0"

[[deps.SQLiteGraph]]
deps = ["DBInterface", "DataFrames", "EasyConfig", "JSON3", "SQLite"]
git-tree-sha1 = "004f75e30ea2b98d6ab487a758816ddae8c96cd1"
uuid = "de2defd0-d76e-464a-9029-1d71e199ae58"
version = "0.2.0"

[[deps.SQLite_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "f79c1c58951ea4f5bb63bb96b99bf7f440a3f774"
uuid = "76ed43ae-9a5d-5a62-8c75-30186b810ce8"
version = "3.38.0+0"

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

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

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

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

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

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

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

[[deps.StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─72412469-67d7-42e3-a7b2-ce05e860e302
# ╟─a78fca5f-fa6a-4929-8861-be95c13492fe
# ╟─a5e78f84-b1b4-4e43-afb3-29cd5ec0e8ff
# ╠═1d24602c-a463-11ec-1b01-4b6d9133a279
# ╟─a51e7694-8209-44e9-bff3-ed2fb607c069
# ╠═439c9e48-f1b1-4d3f-b1ae-e2b8bbe13c84
# ╠═678cae6f-d2e8-4d20-85b3-0d5080b07a62
# ╠═b38fb416-419d-40e0-9cb1-d16106969d07
# ╠═cf6d196e-9f7e-444b-a88c-78e32bf88fba
# ╠═be695499-5529-4029-a326-26a9391b2069
# ╠═a291bf9e-48b2-4673-9c0f-1df440ff2db3
# ╠═d99040c0-eb84-4b98-b4c7-971a3681e0a7
# ╠═8e6df98b-22ce-4c1e-bf45-5e154725586a
# ╠═0763d60f-82e9-44af-a26b-931e6749d4af
# ╟─5c222537-f561-4b47-acf5-6731a05d16a9
# ╠═30450ece-69b6-45d4-b9ee-8b453cb0f93f
# ╠═69819a25-1759-4001-a960-44f0dcf88ed2
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
## Adjoints Delta
@adjoint function evaluate(s::Delta, x::AbstractVector, y::AbstractVector)
  evaluate(s, x, y), Δ -> begin
    (nothing, nothing, nothing)
  end
end

@adjoint function Distances.pairwise(d::Delta, X::AbstractMatrix, Y::AbstractMatrix; dims=2)
  D = Distances.pairwise(d, X, Y; dims = dims)
  if dims == 1
      return D, Δ -> (nothing, nothing, nothing)
  else
      return D, Δ -> (nothing, nothing, nothing)
  end
end

@adjoint function Distances.pairwise(d::Delta, X::AbstractMatrix; dims=2)
  D = Distances.pairwise(d, X; dims = dims)
  if dims == 1
      return D, Δ -> (nothing, nothing)
  else
      return D, Δ -> (nothing, nothing)
  end
end

## Adjoints DotProduct
@adjoint function evaluate(s::DotProduct, x::AbstractVector, y::AbstractVector)
  dot(x, y), Δ -> begin
    (nothing, Δ .* y, Δ .* x)
  end
end

@adjoint function Distances.pairwise(d::DotProduct, X::AbstractMatrix, Y::AbstractMatrix; dims=2)
  D = Distances.pairwise(d, X, Y; dims = dims)
  if dims == 1
      return D, Δ -> (nothing, Δ * Y, (X' * Δ)')
  else
      return D, Δ -> (nothing, (Δ * Y')', X * Δ)
  end
end

@adjoint function Distances.pairwise(d::DotProduct, X::AbstractMatrix; dims=2)
  D = Distances.pairwise(d, X; dims = dims)
  if dims == 1
      return D, Δ -> (nothing, 2 * Δ * X)
  else
      return D, Δ -> (nothing, 2 * X * Δ)
  end
end

## Adjoints Sinus
@adjoint function evaluate(s::Sinus, x::AbstractVector, y::AbstractVector)
  d = (x - y)
  sind = sinpi.(d)
  val = sum(abs2, sind ./ s.r)
  gradx = 2π .* cospi.(d) .* sind ./ (s.r .^ 2)
  val, Δ -> begin
    ((r = -2Δ .* abs2.(sind) ./ s.r,), Δ * gradx, - Δ * gradx)
  end
end

@adjoint function ColVecs(X::AbstractMatrix)
    back(Δ::NamedTuple) = (Δ.X,)
    back(Δ::AbstractMatrix) = (Δ,)
    function back(Δ::AbstractVector{<:AbstractVector{<:Real}})
        throw(error("In slow method"))
    end
    return ColVecs(X), back
end

@adjoint function RowVecs(X::AbstractMatrix)
    back(Δ::NamedTuple) = (Δ.X,)
    back(Δ::AbstractMatrix) = (Δ,)
    function back(Δ::AbstractVector{<:AbstractVector{<:Real}})
        throw(error("In slow method"))
    end
    return RowVecs(X), back
end

@adjoint function Base.map(t::Transform, X::ColVecs)
    pullback(_map, t, X)
end

@adjoint function Base.map(t::Transform, X::RowVecs)
    pullback(_map, t, X)
end


text:
stress_scale = 1.0e5
strain_scale = 1.0

tid= 200

include("nnutil.jl")

# H0 = constant(H1/stress_scale)
testtype = "NeuralNetwork2D"

force_scales = [5.0]

force_scale = 5.0

H0 = [1335174.0968380707  326448.3267263398 0.0
       326448.3267263398 1326879.2022994285 0.0
            0.0                0.0     526955.763626241]/stress_scale

H0inv = inv(H0)

n_data = [200]
porder = 2
# density 4.5*(1 - 0.25) + 3.2*0.25
fiber_fraction = 0.25
#todo
#fiber_fraction = 1.0
prop = Dict("name"=> testtype, "rho"=> 4.5*(1 - fiber_fraction) + 3.2*fiber_fraction, "nn"=>nn)

T = 0.1
NT = 200

# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5

porder = 2

nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny,porder; force_scale=force_scale)

ndofs=2
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]

        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]

        else
            error("polynomial order error, porder= ", porder)

        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords, elnodes, prop, 3))
    end
end

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

nstrain = 3
ngps_per_elem = length(domain.elements[1].weights)
neles = domain.neles

obs_state = zeros(Float64, NT+1, domain.neqs)
state = zeros(Float64, NT+1, domain.neqs)
strain = zeros(Float64,NT+1, neles * ngps_per_elem, nstrain)
stress = zeros(Float64,NT+1, neles * ngps_per_elem, nstrain)

config = [9, 20, 20, 20, 20, 4]
@load  "Data/theta_ite10_bfgs.jld2" last_theta

ForwardNewmarkSolver(globdat, domain, last_theta, T, NT, strain_scale, stress_scale,
                     obs_state, state, strain, stress, -1.0, 0.0,  1e-4,
                     1e-6, 10)

ts = LinRange(0, T, NT+1)

state_all = zeros(Float64, NT+1, 2*(nx*porder + 1)*(ny*porder + 1))
for i = 1:NT+1
    state_all[i, :] = domain.state
    state_all[i, domain.eq_to_dof] = state[i,:]
end

close("all")
ux = [reshape(state_all[i,1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[1,end] for i = 1:NT+1]
plot(ts, ux)
savefig("Debug/order$porder/test_ux_$(tid)_$force_scale.png")

close("all")
uy = [reshape(state_all[i,(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[1,end] for i = 1:NT+1]
plot(ts, uy)
savefig("Debug/order$porder/test_uy_$(tid)_$force_scale.png")


text:
using OhMyREPL
import REPL

include(Pkg.dir("VT100","test","TerminalRegressionTests.jl"))

TerminalRegressionTests.automated_test(
    joinpath(@__DIR__,"flicker/simple.multiout"),
    # Test a bunch of intersting cases:
    # - Typing out keywords
    # - Inserting/deleting braces
    # - Newlines
    # and try to make sure nothing flickers
    [c for c in "function foo(a = b())\nend\n\x4"],
    aggressive_yield=true) do emuterm
    repl = REPL.LineEditREPL(emuterm, false)
    repl.specialdisplay = val->display(REPL.REPLDisplay(repl), val)
    repl.interface = REPL.setup_interface(repl)
    OhMyREPL.Prompt.insert_keybindings(repl)
    REPL.run_repl(repl)
end


text:
module Dual

using LinearAlgebra
using LogExpFunctions
import OptimalTransport.add_singleton
import OptimalTransport.dot_vecwise

export ot_entropic_semidual
export ot_entropic_semidual_grad
export getprimal_ot_entropic_semidual
export ot_entropic_dual
export ot_entropic_dual_grad
export getprimal_ot_entropic_dual

function logKexp(logK::AbstractMatrix, x::AbstractVector)
    return logsumexp(logK .+ add_singleton(x, Val(1)); dims=2)
end

function logKexp(logK::AbstractMatrix, x::AbstractMatrix; dims)
    return logsumexp(logK + x; dims=dims)
end

"""
    ot_entropic_semidual(μ, v, eps, K)

Computes the semidual (in the second argument) of the entropic optimal transport loss, with source marginal `μ`, regularization parameter `ε`, and Gibbs kernel `K`.

That is,
```math
    \\operatorname{OT}_{\\varepsilon}(\\mu, \\nu) = \\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle + \\varepsilon \\Omega(\\gamma)
```
with ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{ij} \\log \\gamma_{ij}``, then the semidual in the second argument ν is [^Z21]
```math
\\begin{aligned}
    \\operatorname{OT}_{\\varepsilon}^*(\\mu, v) &= \\sup_{\\nu} \\langle v, \\nu \\rangle - \\operatorname{OT}_{\\varepsilon}(\\mu, \\nu) \\ \\
    &= -\\varepsilon \\left\\langle \\mu, \\log\\left( \\dfrac{\\mu}{K e^{v/\\varepsilon}} \\right) - 1\\right\\rangle.
\\end{aligned}
```
Notably, the semidual is computationally advantageous for solving variational problems since it is a smooth and unconstrained function of `v` since it admits a closed form gradient. See [^CP16] for a detailed discussion of dual methods for variational problems in optimal transport.

[^CP16]: Cuturi, Marco, and Gabriel Peyré. "A smoothed dual approach for variational Wasserstein problems." SIAM Journal on Imaging Sciences 9.1 (2016): 320-343.
[^Z21]: Zhang, Stephen Y. “A Unified Framework for Non-Negative Matrix and Tensor Factorisations with a Smoothed Wasserstein Loss.” ArXiv: Machine Learning, 2021.
"""
function ot_entropic_semidual(μ, v, eps, K)
    return eps *
           (-dot_vecwise(xlogx.(μ) .- μ, one.(μ)) + dot_vecwise(μ, log.(K * exp.(v / eps))))
end

"""
    ot_entropic_semidual_grad(μ, v, eps, K)

Computes the gradient with respect to `v` of the semidual of the entropic optimal transport loss. That is,
```math
\\nabla_v \\operatorname{OT}^*_{\\varepsilon}(\\mu, v) = K^\\top \\left( \\dfrac{\\mu}{K e^{v/\\varepsilon}} \\right) \\odot e^{v/\\varepsilon}.
```

See also: [`ot_entropic_semidual`](@ref)
"""
function ot_entropic_semidual_grad(μ, v, eps, K)
    return K' * (μ ./ (K * exp.(v / eps))) .* exp.(v / eps)
end

"""
    getprimal_ot_entropic_semidual(μ, v, eps, K)

Computes the the primal variable `ν` corresponding to the dual variable `v` at optimality. That is,
```math
\\nu^\\star = e^{v^\\star/\\varepsilon} \\odot K^\\top \\dfrac{\\mu}{K e^{v^\\star/\\varepsilon}}.
```

See also: [`ot_entropic_semidual`](@ref)
"""
function getprimal_ot_entropic_semidual(μ, v, eps, K)
    return exp.(v / eps) .* (K' * (μ ./ (K * exp.(v / eps))))
end

"""
    ot_entropic_dual(u, v, eps, K)

Computes the dual in both arguments of entropic optimal transport loss, where `u` and `v` are the dual variables associated with the source and target marginals respectively.

That is,
```math
    \\begin{aligned}
    \\operatorname{OT}_{\\varepsilon}^*(u, v) &= \\sup_{\\mu, \\nu} \\langle u, \\mu \\rangle + \\langle v, \\nu \\rangle - \\operatorname{OT}_\\varepsilon(\\mu, \\nu) \\ \\
    &= \\varepsilon \\log \\langle e^{u/\\varepsilon}, K e^{v/\\varepsilon} \\rangle.
    \\end{aligned}
```
"""
function ot_entropic_dual(u, v, eps, K)
    # (μ, ν) → min_{γ ∈ Π(μ, ν)} ε H(γ | K)
    # has Legendre transform
    # (u, v) → ε log < exp(u/ε), K exp(v/ε) >
    return eps * log.(dot_vecwise(exp.(u / eps), K * exp.(v / eps)))
end

"""
    ot_entropic_dual_grad(u, v, eps, K)

Computes the gradient with respect to `u` and `v` of the dual of the entropic optimal transport loss. That is,
```math
\\begin{aligned}
\\nabla_u \\operatorname{OT}^*_{\\varepsilon}(u, v) &= \\dfrac{e^{u/\\varepsilon} \\odot K e^{v/\\varepsilon}}{\\langle e^{u/\\varepsilon}, K e^{v/\\varepsilon} \\rangle} \\ \\
\\nabla_v \\operatorname{OT}^*_{\\varepsilon}(u, v) &= \\dfrac{e^{v/\\varepsilon} \\odot K^\\top e^{u/\\varepsilon}}{\\langle e^{v/\\varepsilon}, K^\\top e^{u/\\varepsilon} \\rangle}.
\\end{aligned}
```

See also: [`ot_entropic_dual`](@ref)
"""
function ot_entropic_dual_grad(u, v, eps, K)
    U = exp.(u / eps)
    V = exp.(v / eps)
    grad_u = (U .* (K * V)) ./ dot_vecwise(U, K * V)'
    grad_v = (V .* (K' * U)) ./ dot_vecwise(V, K' * U)'
    return grad_u, grad_v
end

"""
    getprimal_ot_entropic_dual(u, v, eps, K)

Computes the the primal variable `γ` corresponding to the dual variable `u, v` at optimality. That is,
```math
    \\gamma = \\operatorname{softmax}(\\mathrm{diag}(e^{u/\\varepsilon}) K \\mathrm{diag}(e^{v/\\varepsilon}))
```

See also: [`ot_entropic_dual`](@ref)
"""
function getprimal_ot_entropic_dual(u, v, eps, K)
    γ = K .* add_singleton(exp.(-u / eps), Val(2)) .* add_singleton(exp.(-v / eps), Val(1))
    return γ
end

end


text:
const CACHE_MODE = Symbol(uppercase(@load_preference("cache_mode", "DEFAULT")))
const CACHE_DIR = @load_preference(
    "cache_dir", joinpath(DEPOT_PATH[1], "datadeps", "JuliaConSchedule")
)
const TIMEOUT = parse(Float64, @load_preference("timeout", "5.0"))
const TERMINAL_LINKS = parse(Bool, @load_preference("terminal_links", "false"))
const JULIACON_TIMEZONE = tz"UTC"

function set_cachemode(mode::Symbol)
    @assert mode in (:DEFAULT, :NEVER, :ALWAYS)
    @set_preferences!("cache_mode" => uppercase(string(mode)))
    @info("New cache mode set; restart your Julia session for this change to take effect!")
end
function set_cachedir(path::String)
    @set_preferences!("cache_dir" => string(path))
    @info("New cache dir set; restart your Julia session for this change to take effect!")
end
function set_timeout(val::Union{Float64,Int})
    @assert val > 0
    @set_preferences!("timeout" => string(val))
    @info("New timeout set; restart your Julia session for this change to take effect!")
end
function set_terminallinks(on::Bool)
    @set_preferences!("terminal_links" => string(on))
    @info(
        "New terminal links preference set; restart your Julia session for this change to take effect!"
    )
end

const PRETALX_JSON_URL = "https://pretalx.com/juliacon2021/schedule/export/schedule.json"
const DATA_ARCHIVE_JSON_URL = "https://raw.githubusercontent.com/JuliaCon/JuliaConDataArchive/master/juliacon2021_schedule/schedule.json"
const jcon = Ref{JuliaConSchedule}()

default_json_url() = DATA_ARCHIVE_JSON_URL
default_now() = TimeZones.now(localzone())

"""
    debugmode(on::Bool=true)

Simulates that we are live / in the middle of JuliaCon.
"""
function debugmode(on::Bool=true)
    if on
        # @eval JuliaCon default_now() = ZonedDateTime(Dates.DateTime("2021-07-20T18:00:00.000"), tz"MET") # JuliaCon2021 workshops
        # @eval JuliaCon default_now() = ZonedDateTime(Dates.DateTime("2021-07-30T22:00:00.000"), tz"MET") # JuliaCon2021 posters
        @eval JuliaCon function default_now()
            return ZonedDateTime(Dates.DateTime("2021-07-30T21:30:00.000"), tz"MET")
        end # JuliaCon2021
    else
        @eval JuliaCon default_now() = TimeZones.now(localzone())
    end
    return nothing
end

"""
    set_json_source(src::Symbol)

Anticipated input: `:pretalx`, `:github` (JuliaConDataArchive)
"""
function set_json_source(src::Symbol)
    if src == :pretalx
        @eval JuliaCon default_json_url() = PRETALX_JSON_URL
    else
        @eval JuliaCon default_json_url() = DATA_ARCHIVE_JSON_URL
    end
    return nothing
end


text:
abstract type AbstractEventHandler end

mutable struct AllToAll{T} <: AbstractEventHandler
    next_collision_time::T
    partner1::Int
    partner2::Int
    collision_type::Symbol
end

AllToAll{T}() where {T} = AllToAll{T}(0, -1, -1, :none)

function AllToAll(fluid::HardSphereFluid{N,T}, flow_type) where {N,T}
	event_handler = AllToAll{T}()
	find_collision!(event_handler, fluid, flow_type)

	return event_handler
end

function find_collision(::AllToAll, balls, box, flow::AbstractFlowDynamics)

	partner1 = -1
	partner2 = -1
	min_collision_time = Inf
	collision_type = :none

	for i in 1:length(balls), j in 1:length(box.walls)

		t = collision_time(balls[i], box.walls[j], flow)

		if t < min_collision_time
			partner1 = i
			partner2 = j
			collision_type = :wall_collision
			min_collision_time = t
		end

	end

	for i in 1:length(balls), j in i+1:length(balls)

		t = collision_time(balls[i], balls[j], flow)

		if t < min_collision_time
			partner1 = i
			partner2 = j
			collision_type = :disc_collision
			min_collision_time = t
		end

	end

	return partner1, partner2, collision_type, min_collision_time
end

find_collision(event_handler, fluid::HardSphereFluid, flow_type) = find_collision(event_handler, fluid.balls, fluid.box, flow_type)

function find_collision!(event_handler::AllToAll, fluid::HardSphereFluid, flow_type)
	partner1, partner2, collision_type, min_collision_time = find_collision(event_handler, fluid, flow_type)

	event_handler.partner1 = partner1
	event_handler.partner2 = partner2
	event_handler.collision_type = collision_type
	event_handler.next_collision_time += min_collision_time
end

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
# Degree-optimal polynomials
export Degopt,
    grow!, get_degopt_crefs, scale!, normalize!, square!, get_degopt_coeffs;

struct Degopt{T}
    x::Vector{Tuple{Vector{T},Vector{T}}}
    y::Vector{T}
end

"""
    degopt=Degopt(x,y)
    degopt=Degopt(HA,HB,y)
    degopt=Degopt(graph)

    struct Degopt{T}
        x::Vector{Tuple{Vector{T},Vector{T}}}
        y::Vector{T}
    end

Creates an object representing a degree-optimal polynomial, i.e., the
coefficients in an evaluation scheme which maximizes the degree for a fixed
number of multiplications. The object represents the coefficients in

    A0=I
    A1=A
    A2=(x *A0+x *A1)(x *A0+x *A1)
    A4=(x *A0+x *A1+x*A2)(x *A0+x *A1+x *A2)
    A8=(x *A0+x *A1+x*A2)(x *A0+x *A1+x *A2)
     ..

and

    Out=y*A0+y*A1+y*A4+y*A8+y*B4...

The `x`-values are given in the argument `x`, which is a
`Vector{Tuple{Vector{T},Vector{T}}}`, containing the elements of each sum. The
`y`-vector contains the elements to form the output. In the constructor with
`HA`, `HB` and `y` the elements of `x` are stored as matrices.

The coefficients in graphs generated by [`graph_degopt`](@ref) can be recovered
by using `Degopt(graph)`.
"""
function Degopt(x, y)
    T = eltype(y)
    return Degopt{T}(x, y)
end

function Degopt(graph)
    # Number of multiplications
    k = count(values(graph.operations) .== :mult)
    # Check so it's not a complicated graph
    if (count(values(graph.operations) .== :ldiv) > 0)
        error(
            "The graph must originate from a Degopt form " *
            "and cannot have divisions",
        )
    end
    (x_cref, y_cref) = get_degopt_crefs(k)
    T = eltype(graph)
    xv = Vector{Tuple{Vector{T},Vector{T}}}(undef, k)
    y = Vector{T}(undef, k + 2)

    for (i, x) in enumerate(x_cref)
        xv[i] = (get_coeffs(graph, x[1]), get_coeffs(graph, x[2]))
    end
    y = get_coeffs(graph, y_cref)
    return Degopt(xv, y)
end

function Degopt(HA::AbstractMatrix, HB::AbstractMatrix, y::AbstractVector)
    m = size(HA, 1)
    T = eltype(HA)
    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, m)
    for k = 1:m
        x[k] = (HA[k, 1:k+1], HB[k, 1:k+1])
    end
    return Degopt(x, y)
end

"""
    grow!(degopt::Degopt)

Increases the [`Degopt`](@ref) by one multiplication without modifying the
function values.
"""
function grow!(degopt::Degopt{T}) where {T}
    k = size(degopt.x[end][1], 1)
    push!(degopt.x, (zeros(T, k + 1), zeros(T, k + 1)))
    push!(degopt.y, zero(T))
    return nothing
end

"""
    (x,z)=get_degopt_crefs(k)
    (x,z)=get_degopt_crefs(graph)

Returns linear combination references (crefs) related
to [`graph_degopt`](@ref). Specifically
`x` is a `Vector{Tuple{Vector{Tuple{Symbol,Int}},Vector{Tuple{Symbol,Int}}}}`
such that `x[2][1]` corresponds to the coefficients of the left hand side of the
multiplication

    B2=(α_2_1 *I + α_2_2 *A)(β_2_1 *I + β_2_2 *A)

i.e., the crefs corresponding to `[α_2_1, α_2_2]`. See [`graph_degopt`](@ref).
Hence, `get_coeffs(graph,x[2][1])` returns the corresponding numerical values of
the coefficients.
"""
function get_degopt_crefs(k)
    TT = Tuple{Symbol,Int}
    x = Vector{Tuple{Vector{TT},Vector{TT}}}(undef, k)
    for i = 1:k
        x[i] = (Vector{TT}(undef, i + 1), Vector{TT}(undef, i + 1))
    end

    for s = 2:k+1
        for (b, base) in enumerate(["Ba$s", "Bb$s"])
            if s == 2
                x[s-1][b][1] = (Symbol(base), 1)
            end
            for i = 2:s-1
                if i == 2
                    x[s-1][b][1] = (Symbol("$(base)_2"), 1)
                end
                x[s-1][b][i] = (Symbol("$(base)_$(i)"), 2)
            end
            x[s-1][b][s] = (Symbol(base), 2)
        end
    end

    zlength = k + 2
    z = Vector{TT}(undef, zlength)
    z[1] = (:T2k2, 1)
    z[2] = (:T2k2, 2)
    for i = 3:zlength-1
        z[i] = (Symbol("T2k$(i)"), 2)
    end
    z[k+2] = (Symbol("T2k$(k+3)"), 2)

    return (x, z)
end

function get_degopt_crefs(graph::Compgraph)
    return get_degopt_crefs(count(values(graph.operations) .== :mult))
end

"""
    scale!(degopt::Degopt,α)

Effectively change a [`Degopt`](@ref) such that the input is scaled by `α`. If
`p` is the original function, `p(α x)` will be the new function.
"""
function scale!(degopt::Degopt, α)
    for x in degopt.x
        x[1][2] *= α
        x[2][2] *= α
    end
    degopt.y[2] *= α
    return nothing
end

"""
    square!(degopt::Degopt)

Effectively square a [`Degopt`](@ref) in the sense that the output is square. If
`p` is the original function, `p(x)^2` will be the new function.
"""
function square!(degopt::Degopt)
    y0 = degopt.y
    push!(degopt.x, (deepcopy(y0), deepcopy(y0)))
    degopt.y[:] = zero(y0)
    push!(degopt.y, 1)
    return nothing
end

function row1_normalize!(xx, a, b, c, d)
    q3 = xx[3]
    xx[1] += q3 * a * c
    xx[2] += q3 * (a * d + b * c)
    xx[3] = q3 * b * d
    return nothing
end

import LinearAlgebra.normalize!; # for overloading
"""
    normalize!(degopt::Degopt,tp=:row1)

Normalizes the [`Degopt`](@ref) coefficients, in the way specified by `tp`. If
the `rp==:row1` the `degopt` will be transformed to an equivalent
[`Degopt`](@ref) with first row equal to `(0 1) (0 1)`.
"""
function normalize!(degopt::Degopt, tp = :row1)
    if (tp == :row1)
        a = degopt.x[1][1][1]
        b = degopt.x[1][1][2]
        c = degopt.x[1][2][1]
        d = degopt.x[1][2][2]

        for (i, x) in enumerate(degopt.x)
            if (i > 1)
                for k = 1:2
                    row1_normalize!(x[k], a, b, c, d)
                end
            end
        end
        row1_normalize!(degopt.y, a, b, c, d)

        degopt.x[1][1][1] = 0
        degopt.x[1][1][2] = 1
        degopt.x[1][2][1] = 0
        degopt.x[1][2][2] = 1

        return degopt
    else
        error("Unknown normalization")
    end
end

"""
    get_degopt_coeffs(degopt) -> (HA,HB,y)
    get_degopt_coeffs(graph) -> (HA,HB,y)

Returns the coefficients of the degree-optimal polynomial using two matrices and
one vectors as described in [`Degopt`](@ref).
"""
function get_degopt_coeffs(degopt::Degopt)
    p = size(degopt.x, 1)
    T = eltype(degopt.y)
    HA = zeros(T, p, p + 1)
    HB = zeros(T, p, p + 1)
    for i = 1:p
        HA[i, 1:(i+1)] = degopt.x[i][1]
        HB[i, 1:(i+1)] = degopt.x[i][2]
    end
    y = degopt.y
    return (HA, HB, y)
end
get_degopt_coeffs(graph::Compgraph) = get_degopt_coeffs(Degopt(graph));


text:
# Test by creating a network, writing it, then reading it back in and make sure that all weights and biases match
import NeuralVerification: Layer, ReLU, Id

@testset "Write + Read nnet test" begin

    # 3 --> 3 --> 2 --> 5
    l1 = Layer([3.0 2.0 1.0; 5.0 6.0 7.0; 8.0 9.0 10.0], [0.8; 1.0; 1.2], ReLU())
    l2 = Layer([1.5 2.5 3.5; 4.5 6.5 7.5], [-1.0; -3.0], ReLU())
    # The .nnet file doesn't store an activation at each layer,
    l3 = Layer([10.0 -1.0; -2.0 3.0; 4.0 5.0; 10.0 7.0; -3.5 -4.5], [0.0; -1.0; 0.0; 10.0; -10.0], Id())

    # Write out the network
    network_file = string(tempname(), ".nnet")
    println("File: ", network_file)
    nnet = Network([l1, l2, l3])
    write_nnet(network_file, nnet; header_text="test_header")

    # Read back in the network
    new_nnet = read_nnet(network_file)

    # Test that all weights, biases, and activations are the same
    @test new_nnet.layers[1].weights == l1.weights;
    @test new_nnet.layers[1].bias == l1.bias;
    @test new_nnet.layers[1].activation == l1.activation;

    @test new_nnet.layers[2].weights == l2.weights;
    @test new_nnet.layers[2].bias == l2.bias;
    @test new_nnet.layers[2].activation == l2.activation;

    @test new_nnet.layers[3].weights == l3.weights;
    @test new_nnet.layers[3].bias == l3.bias;
    @test new_nnet.layers[3].activation == l3.activation;
end


text:
using FastUnitfulOperations
using Test

using SortingAlgorithms
using Unitful

@testset "FastUnitfulOperations.jl" begin
    @testset "Product" begin
        a = rand(200, 3)
        b = rand(3, 24)
        au = a * u"m"
        bu = b * u"s"

        @test uproduct(au, bu) == a * b * u"m*s"
        @test uproduct(a, bu) == a * b * u"s"
        @test uproduct(au, b) == a * b * u"m"
        @test uproduct(a, b) == a * b

        @test au *ᵘ bu == a * b * u"m*s"
    end

    @testset "Sorting" begin
        a = rand(1000)
        au = a * u"m"

        @test usort(au) |> issorted
        @test usort(au; alg=RadixSort) |> issorted
        @test usort(a; alg=QuickSort) |> issorted

        auc = copy(au)
        usort!(auc)
        @test auc |> issorted

        auc = copy(au)
        usort!(auc; alg=RadixSort)
        @test auc |> issorted

        ac = copy(a)
        usort!(ac; alg=QuickSort)
        @test ac |> issorted

        @test au[usortperm(au)] |> issorted
        @test au[usortperm(au; alg=RadixSort)] |> issorted
        @test a[usortperm(a; alg=QuickSort)] |> issorted

        ix = Vector{Int}(undef, length(au))
        usortperm!(ix, au)
        @test au[ix] |> issorted

        ix = collect(1:length(au))
        usortperm!(ix, au; alg=RadixSort, initialized=true)
        @test au[ix] |> issorted

        ix = collect(1:length(a))
        usortperm!(ix, a; alg=QuickSort, initialized=true)
        @test a[ix] |> issorted
    end
end


text:
isinteractive() && include("common.jl")

####
#### test directions mechanism
####

@testset "directions" begin
    directions = Directions(0b110101)
    is_forwards = Vector{Bool}()
    for i in 1:6
        is_forward, directions = next_direction(directions)
        push!(is_forwards, is_forward)
    end
    @test collect(is_forwards) == [true, false, true, false, true, true]
    @test rand(RNG, Directions).flags isa UInt32
end

####
#### dummy trajectory for unit testing
####

"""
Trajectory type that is easy to inspect and reason about, for unit testing.

The field `visited` keeps track of visited nodes, can be reset with `empty!`.
"""
struct DummyTrajectory{L,T,D}
    "Positions that trigger turning."
    turning::T
    "Positions that trigger/mimic divergence."
    divergent::D
    "Log density."
    ℓ::L
    "Visited positions."
    visited::Vector{Int}
end

function DummyTrajectory(ℓ; turning = Set(Int[]), divergent = Set(Int[]))
    DummyTrajectory(turning, divergent, ℓ, Int[])
end

Base.empty!(trajectory) = empty(trajectory.visited)

move(::DummyTrajectory, z, is_forward) = z + (is_forward ? 1 : -1)

function is_turning(::DummyTrajectory, τ)
    turn_flag, positions = τ
    @test length(positions) > 1 # not called on a leaf
    turn_flag
end

function combine_turn_statistics(::DummyTrajectory, τ₁, τ₂)
    turn_flag1, positions1 = τ₁
    turn_flag2, positions2 = τ₂
    @test last(positions1) + 1 == first(positions2) # adjacency and order
    (turn_flag1 && turn_flag2, first(positions1):last(positions2))
end

function combine_visited_statistics(::DummyTrajectory, v₁, v₂)
    a1, s1 = v₁
    a2, s2 = v₂
    (a1 + a2, s1 + s2)
end

# copyied from StatsFuns.jl
log1mexp(x::Real) = x < log(0.5) ? log1p(-exp(x)) : log(-expm1(x))

function combine_proposals(_, ::DummyTrajectory, zeta1, zeta2, logprob2, is_forward)
    lp2 = logprob2 > 0 ? 0.0 : logprob2
    lp1 = logprob2 > 0 ? oftype(lp2, -Inf) : log1mexp(lp2)
    if !is_forward
        # exchange so that we can test for adjacency, and join as a UnitRange
        zeta1, zeta2 = zeta2, zeta1
        lp1, lp2 = lp2, lp1
    end
    z1, p1 = zeta1
    z2, p2 = zeta2
    @test last(z1) + 1 == first(z2) # adjacency and order
    (first(z1):last(z2), vcat(p1 .+ lp1, p2 .+ lp2))
end

function calculate_logprob2(::DummyTrajectory, is_doubling, ω₁, ω₂, ω)
    biased_progressive_logprob2(is_doubling, ω₁, ω₂, ω)
end

function leaf(trajectory::DummyTrajectory, z, is_initial)
    @unpack turning, divergent, ℓ, visited = trajectory
    d = z ∈ divergent
    is_initial && @argcheck !d                              # don't start with divergent
    Δ = ℓ(z)
    v = is_initial ? (0.0, 0) : (min(exp(Δ), 1), 1)
    !is_initial && push!(visited, z)                        # save position
    if d
        nothing, v
    else
        (((z:z, [0.0]),          # ζ = nodes, log prob. within tree
          Δ,                     # ω = Δ for leaf
          (z ∈ turning, z:z)),   # τ
         v)
    end
end

"A log density for testing."
testℓ(z) = -abs2(z - 3) * 0.1

"Total acceptance rate of `ℓ` over `z`"
testA(ℓ, z) = sum(min.(exp.(ℓ.(z)), 1))

"sum of acceptance rates for trajectory."
testA(trajectory::DummyTrajectory) = testA(trajectory.ℓ, trajectory.visited)

@testset "dummy adjacent tree full" begin
    trajectory = DummyTrajectory(testℓ)
    (ζ, ω, τ, z′, i′), v = adjacent_tree(nothing, trajectory, 0, 0, 2, true)
    @test first(ζ) == 1:4
    @test sum(exp, last(ζ)) ≈ 1
    @test trajectory.visited == 1:4
    @test !is_turning(trajectory, τ)
    @test v[1] ≈ testA(trajectory)
    @test v[2] == 4
    @test z′ = i′ == 4
end

@testset "dummy adjacent tree turning" begin
    trajectory = DummyTrajectory(testℓ; turning = 5:7)
    t, v = adjacent_tree(nothing, trajectory, 0, 0, 3, true)
    @test trajectory.visited == 1:6 # [5,6] is turning
    @test t == InvalidTree(5, 6)
    @test v[1] == testA(trajectory)
    @test v[2] == 6
end

@testset "dummy adjacent tree divergent" begin
    trajectory = DummyTrajectory(testℓ; divergent = 5:7)
    t, v = adjacent_tree(nothing, trajectory, 0, 0, 3, true)
    @test trajectory.visited == 1:5 # 5 is divergent
    @test t == InvalidTree(5)
    @test v[1] ≈ testA(testℓ, 1:5)
    @test v[2] == 5
end

@testset "dummy adjacent tree full backward" begin
    trajectory = DummyTrajectory(testℓ)
    (ζ, ω, τ, z′, i′), v = adjacent_tree(nothing, trajectory, 0, 0, 3, false)
    @test first(ζ) == -8:-1
    @test sum(exp, last(ζ)) ≈ 1
    @test trajectory.visited == -(1:8)
    @test !is_turning(trajectory, τ)
    @test v[1] ≈ testA(testℓ, -(1:8))
    @test v[2] == 8
    @test z′ == i′ == -8
end

@testset "dummy sampled tree" begin
    trajectory = DummyTrajectory(testℓ)
    ζ, v, termination, depth = sample_trajectory(nothing, trajectory, 0, 3, Directions(0b101))
    @test trajectory.visited == [1, -1, -2, 2, 3, 4, 5]
    @test first(ζ) == -2:5
    @test sum(exp, last(ζ)) ≈ 1
    @test termination == DynamicHMC.REACHED_MAX_DEPTH
    @test v[1] ≈ testA(trajectory)
    @test v[2] == 7             # initial node does not participate in acceptance rate
end

####
#### Detailed balance tests
####

"An accumulator for log probabilities associated with a position on a trajectory."
empty_accumulator() = Dict{Int,Float64}()

"Add log probabilities `πs` at positions `zs` into `accumulator`."
function add_log_probabilities!(accumulator, zs, πs)
    for (z, π) in zip(zs, πs)
        accumulator[z] = haskey(accumulator, z) ? logaddexp(accumulator[z], π) : π
    end
    accumulator
end

"Normalize an accumulator by depth."
function normalize_accumulator(accumulator, depth)
    D = log(0.5) * depth
    Dict((k => v + D) for (k, v) in pairs(accumulator))
end

"""
An accumulator with the probability of visiting nodes for all trees with `depth`, strarting
from `z`, on `trajectory`.
"""
function visited_log_probabilities(trajectory, z, depth)
    accumulator = empty_accumulator()
    for flags in 0:(2^depth - 1)
        ζ = first(sample_trajectory(nothing, trajectory, z, depth, Directions(UInt32(flags))))
        add_log_probabilities!(accumulator, ζ...)
    end
    normalize_accumulator(accumulator, depth)
end

"""
The probability of visiting node `z′` for all trees with `depth`, strarting from `z`, on
`trajectory`.
"""
function transition_log_probability(trajectory, z, z′, depth)
    p = -Inf
    for flags in 0:(2^depth - 1)
        zs, πs = first(sample_trajectory(nothing, trajectory, z, depth,
                                         Directions(UInt32(flags))))
        ix = findfirst(isequal(z′), zs)
        if ix ≠ nothing
            p = logaddexp(p, πs[ix])
        end
    end
    p + depth * log(0.5)
end

@testset "transition calculations consistency check" begin
    trajectory = DummyTrajectory(testℓ)
    depth = 5
    z = 9
    for (z′, π) in pairs(visited_log_probabilities(trajectory, z, depth))
        @test π ≈ transition_log_probability(trajectory, z, z′, depth)
    end
end

"""
$(SIGNATURES)

For all transitions from `z`, test the detailed balance condition, ie

``ℙ(z) ℙ(j ∣ z) = ℙ(j) ℙ(z ∣ j)``

where ``ℙ(z) = exp(ℓ(z))`` and the transition probabilities ``ℙ(⋅∣⋅)`` are
calculated using `visited_log_probabilities` and `transition_log_probability`.

(We use logs for more accurate calculation.)
"""
function test_detailed_balance(trajectory, z, depth; atol = √eps())
    @unpack ℓ = trajectory
    ℓz = ℓ(z)
    for (z′, π) in pairs(visited_log_probabilities(trajectory, z, depth))
        π′ = transition_log_probability(trajectory, z′, z, depth)
        @test (π + ℓz) ≈ (π′ + ℓ(z′)) atol = atol
    end
end

@testset "detailed balance" begin
    for max_depth in 1:5
        test_detailed_balance(DummyTrajectory(testℓ), 0, max_depth)
    end
    for max_depth in 1:5
        test_detailed_balance(DummyTrajectory(testℓ; turning = 1:2), 3, max_depth)
    end
    for max_depth in 1:6
        test_detailed_balance(DummyTrajectory(testℓ; divergent = 10:11), 3, max_depth)
    end
    for max_depth in 1:6
        test_detailed_balance(DummyTrajectory(testℓ; divergent = 10:12, turning = -3:-2), 3,
                              max_depth)
    end
end


text:
println("Dg) Test of the PhotoIonizationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
wa = Atomic.Computation("xx",  Nuclear.Model(26.),
                        initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                        finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                        process = JAC.PhotoIonFluor,
                        processSettings=PhotoIonizationFluores.Settings()  )

wb = perform(wa)

## [E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], false, false, false, true, true, Tuple{Int64,Int64,Int64}[(1,1,0), (1,2,0), (1,2,0)], 0., 0., 0., "Coulomb"


text:
##########################################
########   GaGe Low-level API    #########
##########################################

function CsInitialize()
    ccall((:CsInitialize, :CsSsm), Int32, ())
end

function CsGetSystem( ; u32BoardType=UInt32(0), u32Channels=UInt32(0),
                        u32SampleBits=UInt32(0), i16Index=Int16(0))

    phSystem = Ref{CSHANDLE}(0)

    ccall((:CsGetSystem, :CsSsm), Int32,
          (Ref{CSHANDLE}, UInt32, UInt32, UInt32, Int16),
          phSystem, u32BoardType, u32Channels, u32SampleBits, i16Index)

    return phSystem[]
end

function CsFreeSystem(arg1::CSHANDLE)
    ccall((:CsFreeSystem, :CsSsm), Int32, (CSHANDLE,), arg1)
end

function CsGet(hSystem::CSHANDLE, nIndex::Integer, nConfig::Integer)

    nIndex == CS_ACQUISITION ? pData = CSACQUISITIONCONFIG() :
    nIndex == CS_CHANNEL     ? pData = CSCHANNELCONFIG() :
                               pData = CSTRIGGERCONFIG()

    ccall((:CsGet, :CsSsm), Int32,
          (CSHANDLE, Int32, Int32, Ref{CSACQUISITIONCONFIG}),
          hSystem, nIndex, nConfig, pData)

    return pData
end

function CsSet(hSystem::CSHANDLE, nIndex::Integer)
    nIndex == CS_ACQUISITION ? pData = CSACQUISITIONCONFIG() :
    nIndex == CS_CHANNEL     ? pData = CSCHANNELCONFIG() :
                               pData = CSTRIGGERCONFIG()

    ccall((:CsSet, :CsSsm), Int32,
          (CSHANDLE, Int32, Ref{CSACQUISITIONCONFIG}), hSystem, nIndex, pData)

    return pData
end

function CsGetSystemInfo(hSystem::CSHANDLE)
    SystemInfo = CSSYSTEMINFO()
    ccall((:CsGetSystemInfo, :CsSsm), Int32, (CSHANDLE, PCSSYSTEMINFO),
                                              hSystem, pSystemInfo)
end

function CsGetSystemCaps(hSystem::CSHANDLE, CapsId::UInt32, pBuffer, BufferSize)
    ccall((:CsGetSystemCaps, :CsSsm), Int32, (CSHANDLE, UInt32, Ptr{Cvoid}, Ptr{UInt32}), hSystem, CapsId, pBuffer, BufferSize)
end

function CsDo(hSystem::CSHANDLE, i16Operation::Integer)
    ccall((:CsDo, :CsSsm), Int32, (CSHANDLE, Int16), hSystem, i16Operation)
end

function CsTransfer(hSystem::CSHANDLE)
    pInData::PIN_PARAMS_TRANSFERDATA
    outData::POUT_PARAMS_TRANSFERDATA

    ccall((:CsTransfer, :CsSsm), Int32, (CSHANDLE, PIN_PARAMS_TRANSFERDATA,
                        POUT_PARAMS_TRANSFERDATA), hSystem, pInData, outData)

end

function CsTransferEx(hSystem::CSHANDLE)
    pInData::PIN_PARAMS_TRANSFERDATA_EX
    outData::POUT_PARAMS_TRANSFERDATA_EX

    ccall((:CsTransferEx, :CsSsm), Int32,
            (CSHANDLE, PIN_PARAMS_TRANSFERDATA_EX, POUT_PARAMS_TRANSFERDATA_EX),
            hSystem, pInData, outData)

end

function CsGetEventHandle(hSystem::CSHANDLE, u32EventType::UInt32, phEvent)
    ccall((:CsGetEventHandle, :CsSsm), Int32, (CSHANDLE, UInt32, Ptr{Cint}),
        hSystem, u32EventType, phEvent)
end

function CsGetStatus(hSystem::CSHANDLE)
    ccall((:CsGetStatus, :CsSsm), Int32, (CSHANDLE,), hSystem)
end

function CsGetErrorString(i32ErrorCode::Integer)
    errstring = Vector{UInt8}(undef,255)
    ccall((:CsGetErrorStringA, :CsSsm), Int32, (Int32, Ref{UInt8}, Cint),
                                        i32ErrorCode, errstring, 255)

    return errstring |> pointer |> unsafe_string
end

function CsTransferAS(hSystem::CSHANDLE)
    pInData::PIN_PARAMS_TRANSFERDATA
    pOutParams::POUT_PARAMS_TRANSFERDATA
    pToken
    ccall((:CsTransferAS, :CsSsm), Int32,
    (CSHANDLE, PIN_PARAMS_TRANSFERDATA, POUT_PARAMS_TRANSFERDATA, Ptr{Int32}),
    hSystem, pInData, pOutParams, pToken)

end

function CsGetTransferASResult(hSystem::CSHANDLE)
    nChannelIndex::Int32
    pi64Results
    ccall((:CsGetTransferASResult, :CsSsm),Int32, (CSHANDLE, Int32, Ptr{int64}),
                                            hSystem, nChannelIndex, pi64Results)

end

function CsExpertCall(hSystem::CSHANDLE, pFunctionParams::Cint)
    ccall((:CsExpertCall, :CsSsm), Int32, (CSHANDLE, Cint), hSystem, pFunctionParams)
end

function CsConvertToSigHeader(pHeader::PCSDISKFILEHEADER,
                              pSigStruct::PCSSIGSTRUCT, szComment, szName)
    ccall((:CsConvertToSigHeader, :CsSsm), Int32,
                        (PCSDISKFILEHEADER, PCSSIGSTRUCT, Ptr{Cint}, Ptr{Cint}),
                        pHeader, pSigStruct, szComment, szName)
end

function CsConvertFromSigHeader(pHeader::PCSDISKFILEHEADER,
                                pSigStruct::PCSSIGSTRUCT, szComment, szName)
    ccall((:CsConvertFromSigHeader, :CsSsm), Int32,
        (PCSDISKFILEHEADER, PCSSIGSTRUCT, Ptr{Cint}, Ptr{Cint}),
        pHeader, pSigStruct, szComment, szName)
end

function CsRetrieveChannelFromRawBuffer(pMulrecRawDataBuffer::Cint,
        i64RawDataBufferSize::int64, u32SegmentId::UInt32, ChannelIndex::uInt16,
        i64Start::int64, i64Length::int64, pNormalizedDataBuffer::Cint,
        i64TrigTimeStamp, i64ActualStart, i64ActualLength)
    ccall((:CsRetrieveChannelFromRawBuffer, :CsSsm), Int32,
(Cint, int64, UInt32, uInt16, int64, int64, Cint, Ptr{int64}, Ptr{int64}, Ptr{int64}),
    pMulrecRawDataBuffer, i64RawDataBufferSize, u32SegmentId, ChannelIndex,
    i64Start, i64Length, pNormalizedDataBuffer, i64TrigTimeStamp, i64ActualStart,
    i64ActualLength)
end

function SetAcquisitionParameters(hSystem::CSHANDLE, pCsAcquisitionCfg)
    ccall((:SetAcquisitionParameters, :CsSsm), Int32,
        (CSHANDLE, Ref{CSACQUISITIONCONFIG}), hSystem, pCsAcquisitionCfg)
end

function SetChannelParameters(hSystem::CSHANDLE, pCsChannelCfg)
    ccall((:SetChannelParameters, :CsSsm), Int32,
            (CSHANDLE, Ptr{CSCHANNELCONFIG}), hSystem, pCsChannelCfg)
end

function SetTriggerParameters(hSystem::CSHANDLE, pCsTriggerCfg)
    ccall((:SetTriggerParameters, :CsSsm), Int32,
                    (CSHANDLE, Ptr{CSTRIGGERCONFIG}), hSystem, pCsTriggerCfg)
end

##########################################
########     Custom Wrappers     #########
##########################################
function CsGetAcq(hSystem::CSHANDLE, nConfig::Integer)
    pData = CSACQUISITIONCONFIG()
    err = ccall((:CsGet, :CsSsm), Int32,
          (CSHANDLE, Int32, Int32, Ref{CSACQUISITIONCONFIG}),
          hSystem, CS_ACQUISITION, nConfig, pData)
    if err < 0
        return error(CsGetErrorString(err))
    else
        return pData
    end
end

function CsGetTrig(hSystem::CSHANDLE, nConfig::Integer)
    pData = CSTRIGGERCONFIG()
    err = ccall((:CsGet, :CsSsm), Int32,
          (CSHANDLE, Int32, Int32, Ref{CSTRIGGERCONFIG}),
          hSystem, CS_TRIGGER, nConfig, pData)
    if err < 0
        return error(CsGetErrorString(err))
    else
        return pData
    end
end

function CsGetChan(hSystem::CSHANDLE, nConfig::Integer)
    pData = CSCHANNELCONFIG()
    err = ccall((:CsGet, :CsSsm), Int32,
          (CSHANDLE, Int32, Int32, Ref{CSCHANNELCONFIG}),
          hSystem, CS_CHANNEL, nConfig, pData)
    if err < 0
        return error(CsGetErrorString(err))
    else
        return pData
    end
end


text:
module TestForeach

using Folds
using Referenceables
using Test
using Transducers: PreferParallel
include("utils.jl")

raw_testdata_executors = """
PreferParallel()
SequentialEx()
ThreadedEx()
ThreadedEx(basesize = 1)
"""

testdata_executors = map(split(raw_testdata_executors, "\n", keepempty = false)) do x
    x => Base.include_string((@__MODULE__), x)
end

@testset "increment" begin
    @testset "default" begin
        xs = zeros(3)
        Folds.foreach(referenceable(xs)) do x
            x[] += 1
        end
        @test xs == ones(3)
    end
    @testset "$label" for (label, ex) in testdata_executors
        xs = zeros(3)
        Folds.foreach(referenceable(xs), ex) do x
            x[] += 1
        end
        @test xs == ones(3)
    end
end

@testset "DistributedEx should throw" begin
    xs = ones(3)
    err = @test_error Folds.foreach(identity, referenceable(xs), DistributedEx())
    @test "used with a referenceable" ∈ᵉʳʳᵒʳ err
end

end  # module


text:
# function fullGreen(para, site::AbstractVector, orbital::AbstractVector, extT::AbstractVector = collect(1:length(orbital)), subdiagram = false; name = Symbol("Gn$(length(site))"), resetuid = false, even = true)
#     # @assert para.diagType == GreenNDiag
#     @assert length(extT) == length(orbital) == length(site)
#     if even
#         @assert length(extT) % 2 == 0
#     end

#     resetuid && uidreset()

#     gn = []
#     uniqueR = Set(site)
#     permutation = [] # keep track of the permutation after the site index rearrangement
#     for r in uniqueR
#         ind = findall(x -> x == r, site)
#         if even && (length(ind) % 2 == 1)
#             continue
#         end
#         t = extT[ind]
#         o = orbital[ind]
#         bareGId = BareGreenNId(para, o, t, r)
#         push!(gn, Diagram(bareGId, name = Symbol("gn$(length(t))")))
#         append!(permutation, ind)
#     end

#     if isempty(gn)
#         return nothing
#     else
#         return Diagram(GreenNId(para, orbital, extT, site), Prod(), gn, name = name, factor = parity(permutation))
#     end
# end

function fullGreen(para, hop::Vector{BareHoppingId}, subdiagram = false; name = Symbol("Gn$(length(hop)*2)"), resetuid = false, even = true)
    # @assert para.diagType == GreenNDiag
    # @assert length(extT) == length(orbital) == length(site)
    # if even
    #     @assert length(extT) % 2 == 0
    # end
    extT, orbital, site, _creation = [], [], [], []
    for h in hop
        append!(extT, h.extT)
        append!(site, h.site)
        append!(_creation, [true, false])
        append!(orbital, h.orbital)
    end
    # println("calculate: ", hop, " . site: ", site)

    resetuid && uidreset()

    gn = []
    uniqueR = Set(site)
    permutation = [] # keep track of the permutation after the site index rearrangement
    for r in uniqueR
        ind = findall(x -> x == r, site)
        if even && (length(ind) % 2 == 1)
            return nothing
        end
        t = extT[ind]
        o = orbital[ind]
        c = _creation[ind]
        bareGId = BareGreenNId(para, orbital = o, t = t, r = r, creation = c)
        push!(gn, Diagram(bareGId, name = Symbol("gn$(length(t))")))
        append!(permutation, ind)
    end

    if isempty(gn)
        return nothing
    else
        for h in hop
            push!(gn, Diagram(h, name = :hop))
        end
        # println(permutation)
        return Diagram(GreenNId(para, orbital = orbital, t = extT, r = site, creation = _creation), Prod(), gn, name = name, factor = parity(permutation))
    end
end


text:
#! format: off

using SIIPExamples # Only needed for the tutorial, comment if you want to run
import DisplayAs # Only needed for the tutorial
using PowerSimulationsDynamics
using PowerSystems
using Logging
using Sundials
using Plots
gr()
PSD = PowerSimulationsDynamics

file_dir = joinpath(
    dirname(dirname(pathof(SIIPExamples))),
    "script",
    "4_PowerSimulationsDynamics_examples",
    "Data",
)

sys = System(joinpath(file_dir, "14bus.raw"), joinpath(file_dir, "dyn_data.dyr"))

sim = PSD.Simulation(
    ResidualModel, #Type of model used
    sys,         #system
    file_dir,       #path for the simulation output
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "BUS 02-BUS 04-i_4");
    console_level = Logging.Info,
)

show_states_initial_value(sim)

PSD.execute!(sim, IDA(); abstol = 1e-8)

result = read_results(sim)
p = plot()
for b in get_components(Bus, sys)
    voltage_series = get_voltage_magnitude_series(result, get_number(b))
    plot!(
        p,
        voltage_series;
        xlabel = "Time",
        ylabel = "Voltage Magnitude [pu]",
        label = "Bus - $(get_name(b))",
    )
end
img = DisplayAs.PNG(p) # This line is only needed because of literate use display(p) when running locally

p2 = plot()
for g in get_components(ThermalStandard, sys)
    state_series = get_state_series(result, (get_name(g), :ω))
    plot!(
        p2,
        state_series;
        xlabel = "Time",
        ylabel = "Speed [pu]",
        label = "$(get_name(g)) - ω",
    )
end
img = DisplayAs.PNG(p2) # This line is only needed because of literate use display(p2) when running locally

res = small_signal_analysis(sim; reset_simulation = true)

scatter(res.eigenvalues; legend = false)

sys = System(joinpath(file_dir, "14bus.raw"), joinpath(file_dir, "dyn_data.dyr"))

thermal_gen = get_component(ThermalStandard, sys, "generator-6-1")
remove_component!(sys, get_dynamic_injector(thermal_gen))
remove_component!(sys, thermal_gen)

storage = GenericBattery(
    name = "Battery",
    bus = get_component(Bus, sys, "BUS 06"),
    available = true,
    prime_mover = PrimeMovers.BA,
    active_power = 0.6,
    reactive_power = 0.16,
    rating = 1.1,
    base_power = 25.0,
    initial_energy = 50.0,
    state_of_charge_limits = (min = 5.0, max = 100.0),
    input_active_power_limits = (min = 0.0, max = 1.0),
    output_active_power_limits = (min = 0.0, max = 1.0),
    reactive_power_limits = (min = -1.0, max = 1.0),
    efficiency = (in = 0.80, out = 0.90),
)

add_component!(sys, storage)

res = solve_powerflow(sys)
res["bus_results"]

inverter = DynamicInverter(
    name = get_name(storage),
    ω_ref = 1.0, # ω_ref,
    converter = AverageConverter(rated_voltage = 138.0, rated_current = 100.0),
    outer_control = OuterControl(
        VirtualInertia(Ta = 2.0, kd = 400.0, kω = 20.0),
        ReactivePowerDroop(kq = 0.2, ωf = 1000.0),
    ),
    inner_control = CurrentControl(
        kpv = 0.59,     #Voltage controller proportional gain
        kiv = 736.0,    #Voltage controller integral gain
        kffv = 0.0,     #Binary variable enabling the voltage feed-forward in output of current controllers
        rv = 0.0,       #Virtual resistance in pu
        lv = 0.2,       #Virtual inductance in pu
        kpc = 1.27,     #Current controller proportional gain
        kic = 14.3,     #Current controller integral gain
        kffi = 0.0,     #Binary variable enabling the current feed-forward in output of current controllers
        ωad = 50.0,     #Active damping low pass filter cut-off frequency
        kad = 0.2,
    ),
    dc_source = FixedDCSource(voltage = 600.0),
    freq_estimator = KauraPLL(
        ω_lp = 500.0, #Cut-off frequency for LowPass filter of PLL filter.
        kp_pll = 0.084,  #PLL proportional gain
        ki_pll = 4.69,   #PLL integral gain
    ),
    filter = LCLFilter(lf = 0.08, rf = 0.003, cf = 0.074, lg = 0.2, rg = 0.01),
)
add_component!(sys, inverter, storage)

sys

sim = PSD.Simulation(
    ResidualModel, #Type of model used
    sys,         #system
    file_dir,       #path for the simulation output
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "BUS 02-BUS 04-i_4");
    console_level = Logging.Info,
)

res = small_signal_analysis(sim)

scatter(res.eigenvalues)

PSD.execute!(sim, IDA(); abstol = 1e-8)

result = read_results(sim)
p = plot()
for b in get_components(Bus, sys)
    voltage_series = get_voltage_magnitude_series(result, get_number(b))
    plot!(
        p,
        voltage_series;
        xlabel = "Time",
        ylabel = "Voltage Magnitude [pu]",
        label = "Bus - $(get_name(b))",
    )
end
img = DisplayAs.PNG(p) # This line is only needed because of literate use display(p) when running locally

p2 = plot()
for g in get_components(ThermalStandard, sys)
    state_series = get_state_series(result, (get_name(g), :ω))
    plot!(
        p2,
        state_series;
        xlabel = "Time",
        ylabel = "Speed [pu]",
        label = "$(get_name(g)) - ω",
    )
end
state_series = get_state_series(result, ("Battery", :ω_oc))
plot!(p2, state_series; xlabel = "Time", ylabel = "Speed [pu]", label = "Battery - ω")
img = DisplayAs.PNG(p2) # This line is only needed because of literate use display(p2) when running locally

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
export iauPvmpv
"""
Subtract one pv-vector from another.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a       double[2][3]      first pv-vector
   b       double[2][3]      second pv-vector

Returned:
   amb     double[2][3]      a - b

Note:
   It is permissible to re-use the same array for any of the
   arguments.

Called:
   iauPmp       p-vector minus p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvmpv(double a[2][3], double b[2][3], double amb[2][3])
function iauPvmpv(a::AbstractMatrix{<:Real}, b::AbstractMatrix{<:Real})
   # Preallocate return values
   amb = zeros(Float64, 3, 2)

   ccall((:iauPvmpv, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Ref{Cdouble}),
          convert(Array{Float64, 2}, a'),
          convert(Array{Float64, 2}, b'),
          amb)

   return SMatrix{2,3}(amb')
end


text:
# Model: Y. S. Cheng, C. S. Wang (Motion of Particles in Bends of Circular Pipes, 1980)
# Based on: bendTubePenetrationCheng1980.m

dir = "./src"
push!(LOAD_PATH,dir)

module tubeBend

import physicalConstants: k, g

struct Parameters
#CarrierFluidProfile
     T     :: Float64                                        # Carrier fluid temperature (T) SI units: K
     Q     :: Float64                                        # Carrier fluid flow rate (Q) SI units: m3 s-1
     mu    :: Float64                                        # Carrier fluid dynamic viscosity (mu) SI units: N s m-2
end
struct Properties
#Fluid
     ro_d  :: Float64                                        # Droplet density (ro_d) SI units: kg m-3
     mu_d  :: Float64                                        # Droplet dynamic viscosity (mu_d) SI units: N s m-2
     d_d   :: Float64                                        # Droplet diameter (d_d) SI units: m
#Tube
     Rb    :: Float64                                        # Radius of the bend (Rb) SI units: m
     theta :: Float64                                        # Bend angle (theta) SI units: radians
     d     :: Float64                                        # Tube inner diameter (d) SI units: m
end

function AerosolPenetration(d :: Properties,
                           Rb :: Properties,
                            Q :: Parameters) :: Float64
     Rt = d/2                                                # Tube inner radius (d) SI units: m
     Ro = Rb/Rt                                              # Curvature ratio (Ro) SI units: dimensionless
     U = 4*Q/(pi*d^2)                                        # Fluid flow mean velocity (U) SI units: m s-1
end

function Cunninghamslipcorrection(mu :: Parameters,
                                   T :: Parameters,
                                 d_d :: Properties) :: Float64
     lambda = (mu/p)*sqrt((pi*k*T)/(2*MM))                   # Model for the mean free path of the carrier fluid (lambda)
     C = 1 + (lambda/dd)*(2.34+1.05*exp(-.39*d_d/lambda))    #Cunningham slip correction (C)
end

function

#Debug test
# println(k)

end


text:
tree_hash_fmt(d) = bytes2hex(Pkg.GitTools.tree_hash(d))

EMPTY_TREE_HASH = "none (package dev'ed)"

function project_and_tree_hash(pkg::Base.PkgId)
    for env in Base.load_path()
        project_file = Base.env_project_file(env)
        project_file isa String || continue
        mth = explicit_manifest_uuid_path_tree_hash(project_file,pkg)
        mth === nothing || return env, mth
    end
    return nothing, nothing
end

# Analogous to explicit_manifest_uuid_path from base/loading.jl,
# but returning the tree_hash in the manifest instead
function explicit_manifest_uuid_path_tree_hash(project_file::String, pkg::Base.PkgId)::Union{Nothing,String}
    manifest_file = Base.project_file_manifest_path(project_file)
    manifest_file === nothing && return nothing # no manifest, skip env

    d = Base.get_deps(Base.parsed_toml(manifest_file))
    entries = get(d, pkg.name, nothing)::Union{Nothing, Vector{Any}}
    entries === nothing && return nothing # TODO: allow name to mismatch?
    for entry in entries
        entry = entry::Dict{String, Any}
        uuid = get(entry, "uuid", nothing)::Union{Nothing, String}
        uuid === nothing && continue
        if Base.UUID(uuid) === pkg.uuid
            return get(entry, "git-tree-sha1", EMPTY_TREE_HASH)::Union{Nothing, String}
        end
    end
    return nothing
end


text:
using Printf
inputF = ARGS[1]

function getP_mutations(infile)
    f = open( infile, "r")
    n = countlines(f)
    seekstart(f)

    germMutation_dict = Dict{String, String}()
    somaticMutation_dict = Dict{String, String}()

    ## println(string("CHR\tPOS\tFREQ"))
    for i = 1:n
        line = readline(f)
        #if !( startswith( line, "序号" ) )
        if !(occursin("序号", line))
            s = split(line, "\t")
            #freq = parse(Float64, s[1])

            patho_HRR = split(s[9], ";")

            vus_HRR = split(s[10], ";")
            for v in patho_HRR
                if ( startswith( v, "g" ) )
                    geneInfo = split(v, ":")
                    if length(geneInfo) == 2
                        println(@sprintf("%s\t%s\t胚系\t良性", geneInfo[1][2:end], geneInfo[2]))
                    else
                        println(@sprintf("%s\t非点突变\t胚系\t良性", geneInfo[1][2:end] ))
                    end
                elseif (v != "" && v != "/" )
                    println(@sprintf("%s\t非点突变\t体系\t致病", v ))
                end

            end
            for v in vus_HRR
                if ( startswith( v, "g" ) )
                    geneInfo = split(v, ":")
                    if length(geneInfo) == 2
                        println(@sprintf("%s\t%s\t胚系\t意义不明", geneInfo[1][2:end], geneInfo[2]))
                    else
                        println(@sprintf("%s\t非点突变\t胚系\t意义不明", geneInfo[1][2:end] ))
                    end
                elseif (v != "" && v != "/" )
                    println(@sprintf("%s\t非点突变\t体系\t意义不明", v ))
                end
            end

        end
    end
    return 1
end

convertF = "mutation_a1.txt"

## special characters "#{}()[]<>|&*?~;" should now be quoted in commands
##

##source is built-in command

##try like this:

##   run(`zsh -c "source $cmd"`)

##   run(`bash -c source $cmd`)

##windows:

##   run(`cmd /c 'echo xxx'`)

run(`powershell /c Get-Content -Path $inputF "|" Out-File -FilePath .\\$convertF -Encoding utf8`)

#run(`powershell /c Get-Content -Path D:\\docs_x\\$inputF`)

#a = getP_mutations(inputF)
a = getP_mutations(convertF)


text:
@testset "OpenCL.CmdQueue" begin
    @testset "OpenCL.CmdQueue constructor" begin
        has_warned = false
        @test_throws MethodError cl.CmdQueue(nothing, nothing)
        for platform in cl.platforms()
            for device in cl.devices(platform)
                ctx = cl.Context(device)
                @test cl.CmdQueue(ctx) != nothing
                @test cl.CmdQueue(ctx, device) != nothing
                @test cl.CmdQueue(ctx, :profile) != nothing
                try
                    cl.CmdQueue(ctx, device, :out_of_order)
                    cl.CmdQueue(ctx, device, (:profile, :out_of_order))
                catch err
                    if !has_warned
                        warn("Platform $(device[:platform][:name]) does not seem to " *
                             "suport out of order queues: \n$err")
                        has_warned = true
                    end
                end
                @test_throws ArgumentError cl.CmdQueue(ctx, device, :unrecognized_flag)
                for flag in [:profile, :out_of_order]
                    @test_throws ArgumentError cl.CmdQueue(ctx, (flag, :unrecognized_flag))
                    @test_throws ArgumentError cl.CmdQueue(ctx, device, (:unrecognized_flag, flag))
                    @test_throws ArgumentError cl.CmdQueue(ctx, (flag, flag))
                    @test_throws ArgumentError cl.CmdQueue(ctx, device, (flag, flag))
                end
            end
        end
    end

    @testset "OpenCL.CmdQueue info" begin
        for platform in cl.platforms()
            for device in cl.devices(platform)
                ctx = cl.Context(device)
                q1 = cl.CmdQueue(ctx)
                q2 = cl.CmdQueue(ctx, device)
                for q in (q1, q2)
                    @test q[:context] == ctx
                    @test q[:device] == device
                    @test q[:reference_count] > 0
                    @test typeof(q[:properties]) == cl.CL_command_queue_properties
                end
            end
        end
    end
end


text:
mutable struct CachedFactor{T<:Unsigned} <: AbstractFactor{T}
    cache::Union{Vector{T}, Nothing}
    basefactor::AbstractFactor{T}
    lockobj::Threads.TatasLock
end

Base.length(var::CachedFactor) = length(var.basefactor)

getname(var::CachedFactor) = getname(var.basefactor)

getlevels(var::CachedFactor) = getlevels(var.basefactor)

function CachedFactor(factor::AbstractFactor{T}) where {T<:Unsigned}
    CachedFactor{T}(nothing, factor, Threads.TatasLock())
end

function cache(basefactor::AbstractFactor{T}) where {T<:Unsigned}
    CachedFactor{T}(nothing, basefactor, Threads.TatasLock())
end

function slice(factor::CachedFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    basefactor = factor.basefactor
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if isa(basefactor, CachedFactor) || isa(basefactor, Factor)
        slice(basefactor, fromobs, toobs, slicelength)
    else
        lockobj = factor.lockobj
        lock(lockobj)
        try
            if factor.cache === nothing
                slices = slice(basefactor, 1, length(basefactor), SLICELENGTH)
                cachedata = Vector{T}(undef, length(basefactor))
                fold(0, slices) do offset, slice
                    n = length(slice)
                    view(cachedata, (1 + offset):(n + offset)) .= slice
                    offset += n
                end
                factor.cache = cachedata
            end
        finally
            unlock(lockobj)
        end
        slice(factor.cache, fromobs, toobs, slicelength)
    end
end

function Base.map(factor::CachedFactor{T}, dataframe::AbstractDataFrame; permute::Bool = false) where {T<:Unsigned}
    CachedFactor(map(factor.basefactor, dataframe; permute = permute))
end

function isordinal(factor::CachedFactor{T}) where {T<:Unsigned}
    isordinal(factor.basefactor)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# This file was generated, do not modify it. # hide
using Gnuplot
ecycl_x(r,k,θ) = r*(k .+ 1).*cos.(θ) .- r*cos.((k .+ 1) .* θ)
ecycl_y(r,k,θ) = r*(k .+ 1).*sin.(θ) .- r*sin.((k .+ 1) .* θ)
θ = LinRange(0,6.2π,1000)
@gp(ecycl_x(1,1,θ), ecycl_y(1,1,θ), "w l lw 2 t '1'",
    "set size square")
for k in 2:0.5:5.
    @gp(:-, ecycl_x(2k,k,θ), ecycl_y(2k,k,θ), "w l lw 2 t '$(k)' ",
        "set key outside title 'k, r=2k' box opaque",
        xlabel = "x(θ) = r(k+1)cos(θ) -rcos((k+1)θ)",
        ylabel = "y(θ) = r(k+1)cos(θ) -rcos((k+1)θ) ",
        title = "Epicycloid",:-)
end
@gp
save(term="pngcairo font 'Consolas, 12' size 600,600", output="plt2_ex6dot1.png") # hide


text:
import Term.box: ROUNDED, get_row, get_title_row, fit
import Term
println("get row")
@time get_row(ROUNDED, [3], :top)
@time get_row(ROUNDED, 3, :top)

println("get title row")
@time get_title_row(:top, Term.box.ROUNDED, nothing, width=12, style="red")

for justify in (:left, :center, :right)
    @time get_title_row(:top, Term.box.ROUNDED, "test", width=12,
                                justify=justify, title_style="blue", style="red")
end


text:
"""
    Univariate optimization and root finding examples with a few different packages
"""

using Plots
#=
 You should load all packages at beginning of file!
 I just load the optimizers later to make it clear which package is being demonstrated.
=#

# Function definition
# f(x) = x[1]^2 - 0.5
f(x) = x[1]^3 - x[1]
xlim_low = -1.5
xlim_high = 1.5
plot1 = plot(f, xlim_low, xlim_high, label = "f(x)")
plot1 = plot!(x -> 0, xlim_low, xlim_high, label = "")
plot(plot1)

# Optimization
initial_value = -0.01
scatter!([initial_value],
         [0.],
         label = "Initial Value",
         markershape = :x)

## Optim
using Optim

optim_result = Optim.optimize(f, [initial_value], BFGS())
print("Optim minimum: ")
if optim_result.ls_success
    println(optim_result.minimizer)
    scatter!([optim_result.minimizer], [optim_result.minimum], label = "Optim minimum")
else
    println("Failed to converge")
end

# Root finding

## Roots
using Roots
roots_zero = find_zero(f, initial_value)
println("Roots 'find_zero' zero: $(roots_zero)")
scatter!([roots_zero],[f(roots_zero)], label = "Roots 'find_zero' zero")

## Alternative syntax
roots_fzero_zero = fzero(f, initial_value)
println("Roots 'fzero' zero: $(roots_fzero_zero)")
scatter!([roots_fzero_zero],[f(roots_zero)], label = "Roots 'fzero' zero")

## Alternative syntax
fx = ZeroProblem(f, initial_value)
roots_solve_zero = solve(fx)
println("Roots 'solve' zero: $(roots_solve_zero)")
scatter!([roots_solve_zero],[f(roots_zero)], label = "Roots 'solve' zero")

## NLsolve
using NLsolve
f!(F, x) = F[1] = f(x)
nls_zero = nlsolve(f!, [initial_value])
println("NLsolve zero: $(nls_zero.zero)")
scatter!([nls_zero.zero], [f(nls_zero.zero)], label = "NLsolve zero")

## LeastSquaresOptim
### Nonlinear least squares
### minₓ F(x)ᵀF(x)
using LeastSquaresOptim

ls_optim_result = LeastSquaresOptim.optimize(f, [initial_value], Dogleg())
print("LeastSquaresOptim minimum: ")
if ls_optim_result.converged
    println(ls_optim_result.minimizer)
    scatter!([ls_optim_result.minimizer], [f(ls_optim_result.minimizer)], label = "LeastSquaresOptim minimum")
else
    println("Failed to converge")
end


text:
# functions for mating

function mate_genotypes(gmap,sire,dam)

end

# chelen: cM
function number_of_crosses(chrlen)
   # to Morgan
   lambda = chrlen/100.0
   return rand(Poisson(lambda))
end

# chelen: cM
function location_of_crossover(chrlen,pos,ncross)
   nloci = length(pos)
   loccross = zeros(Int,ncross)
   ntrials = 0
   while(true)
      ntrials = ntrials +1
      poscross = sort(rand(Uniform(0,chrlen),ncross))
      loccross .= 0
      k = 1
      for i=1:nloci
         if poscross[k]<=pos[i]
            loccross[k] = i
            if k==ncross
              break
            end
            k = k + 1
         end
      end
      if sum(loccross .> 0)==ncross
         break
      end
      if ntrials>10
         # too many trials: no crossover
         @warn "too many trials to find crossovers"
         return (0,Vector{Tuple{Int,Int}}(),Vector{Float64}())
      end
   end
   # which chromosomes
   gameteid = Vector{Tuple{Int,Int}}(undef,ncross)
   for i=1:ncross
      whichp = ifelse(rand()<0.5, 1, 2)
      whichm = ifelse(rand()<0.5, 3, 4)
      gameteid[i] = (whichp,whichm)
   end
   return adjusted_location_of_crossover(gameteid,loccross)
end

# removing double cross etc
function adjusted_location_of_crossover(gameteid,loccross)
   # effective location
   ncross = length(loccross)
   effective = zeros(Int,ncross)
   n = zeros(Int,4)
   prev = zeros(Int,4)
   which = zeros(Int,ncross)
   k = 0
   # remove double cross
   adj_ncross = 0
   for i=1:ncross
      # (1,3)->1, (1,4)->2, (2,3)->3, (2,4)->4
      which[i] = ifelse(gameteid[i][1]==1,gameteid[i][1]+gameteid[i][2]-3,gameteid[i][1]+gameteid[i][2]-2)
      n[which[i]] = n[which[i]] + 1
      if n[which[i]] == 2
         n[which[i]] = 0
         which[prev[which[i]]] = 0
         prev[which[i]] = 0
         which[i] = 0
         adj_ncross = adj_ncross - 1
      else
         prev[which[i]] = i
         adj_ncross = adj_ncross + 1
      end
   end

   # adjusted crossover
   if adj_ncross < ncross
      new_gameteid = Array{Tuple{Int,Int}}(undef,adj_ncross)
      new_loccross = zeros(Int, adj_ncross)
      k = 0
      for i=1:ncross
         if which[i]>0
            k = k +1
            new_gameteid[k] = gameteid[i]
            new_loccross[k] = loccross[i]
         end
      end
      return adj_ncross, new_gameteid, new_loccross
   else
      return ncross, gameteid, loccross
   end
end

# sample a gamete from an individual's paternal and maternal haplotypes
#   gamete1: copy of gp ----------
#   gamete2: copy of gp ----------
#   gamete3: copy of gm ==========
#   gamete4: copy of gm ==========
# crossover occurring between (1,3), (1,4), (2,3), or (2,4)
function generate_gamete(chrlen,pos,gp,gm)
   nloci = length(pos)
   ncross = number_of_crosses(pos[end])
   chosen = rand(1:4)
   if ncross==0 || nloci<ncross+1
      # no crossover; either gp or gm
      if chosen<=2
         return copy(gp)
      else
         return copy(gm)
      end
   else
      adj_ncross,gameteid,loccross = location_of_crossover(chrlen,pos,ncross)
      return generate_gamete_with_crossover(gameteid,loccross,chosen,gp,gm)
   end
end

function generate_gamete_with_crossover(gameteid,loccross,chosen,gp,gm)
   # make a copy
   gp_copy1 = copy(gp)
   gp_copy2 = copy(gp)
   gm_copy1 = copy(gm)
   gm_copy2 = copy(gm)
   tmp = copy(gp)
   tmp .= 0
   # pointers
   gam = [gp_copy1, gp_copy2, gm_copy1, gm_copy2]
   # simulate crossover
   ncross = length(loccross)
   for i=1:ncross
      # which chromosome
      (whichp,whichm) = gameteid[i]
      # swap genotypes between gam[whichp][loc:nloci] and gam[whichm][loc:nloci]
      loc = loccross[i]
      tmp[loc:end] = gam[whichp][loc:end]
      gam[whichp][loc:end] = gam[whichm][loc:end]
      gam[whichm][loc:end] = tmp[loc:end]
   end
   # choose the gamete
   return gam[chosen]
end

"""
    progeny = mating(map,sire,dam)

Given parents `sire` and `dam` as `QMSimIndividualGenome` and a map structure `QMSimMap`, this function produces one progeny as `QMSimIndividualGenome`.
"""
function mating(gmap::QMSimMap,sire::QMSimIndividualGenome,dam::QMSimIndividualGenome)
   chromosome_set = generate_chromosome_set(gmap)
   for i in 1:gmap.nchr
      chrlen = gmap.chr[i].pos[end]
      chromosome_set[i].gp .= generate_gamete(chrlen, gmap.chr[i].pos, sire.chr[i].gp, sire.chr[i].gm)
      chromosome_set[i].gm .= generate_gamete(chrlen, gmap.chr[i].pos, dam.chr[i].gp, dam.chr[i].gm)
   end
   tbv = get_true_breeding_value(gmap,chromosome_set)
   return QMSimIndividualGenome(tbv,chromosome_set)
end


text:
using OrderedCollections
using Test
using Random, Serialization

@test isempty(detect_ambiguities(Base, Core, OrderedCollections))

tests = [
         "ordered_dict",
         "ordered_set",
        ]

if length(ARGS) > 0
    tests = ARGS
end

@testset "OrderedCollections" begin

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end

end # @testset


text:
"""
converts a primitive into a drawable for `window`
"""
function Drawable(window::WT, primitive::PT) where {WT, PT}
    error("Can't draw $PT with window type: $WT")
end

"""
Creates a rasterizer clojure for a certain window type
"""
function rasterizer{WT}(
        window::WT,
        vertexarray::AbstractArray,
        uniforms::Tuple,
        vertexshader::Function,
        fragmentshader::Function;
        geometryshader = nothing,
        max_primitives = 4,
        primitive_in = :points,
        primitive_out = :triangle_strip,
    )
    error("No rasterizer available for window type: $WT")
end

"Enables debugging information for a composable"
@field Debugging = false

"Name of something, usually a String or Symbol"
@field Name = gensym(:Visualize)

"Defines if an object is visible or not"
@field Visible = true

"Defines if an object is focused or not"
@field Focused = false

@field Rotation = Quaternion(1, 0, 0, 0)
@field Scale = Vec3f0(1)
@field Position = Vec3f0(0)
@field Translation = Vec3f0(0)

abstract type Transformation <: ReactiveComposable end

@composed type Transform <: Transformation
    Scale
    Rotation
    Translation
end

@field XAxis
@field YAxis
@field ZAxis
@field Primitive
"""
The Area of a viewable. Can be a Rectangle or a tuple with two elements.
"""
@field Area = IRect(0, 0, 500, 500)

"""
Resolution is the smaller brother of Area. It can be NTuple{2, <: Number} or a Vec{2, <: Number}.
"""
@field Resolution

function default(::Type{Resolution}, x::Partial)
    area = get(x, Area) # APIDOC this gets the area if in `x`, or creates a default area
    widths(area)
end

function default(::Type{Area}, x::Partial)
    haskey(x, Resolution) && return convertfor(Area, x, get(x, Resolution))
    default(Area)
end
function convertfor(::Type{Area}, x::Partial, value)
    isa(value, Tuple) && length(value) == 2 && return IRect(0, 0, value...)
    isa(value, SimpleRectangle) && return IRect(origin(value), widths(value))
    throw(UsageError(Area, value))
end

@field Spacetransform = identity

@field Projection = eye(Mat4f0)
@field View = eye(Mat4f0)
@field ProjectionView = eye(Mat4f0)

@field BoundingBox = AABB{Float32}()

@composed type Shared
    Transform
    BoundingBox
end

@field Color = Vec4f0(1f0, 1f0, 1f0, 1f0)

"""
Thickness e.g. of a line or stroke.
"""
@field Thickness = 2f0
@field Vertices = error("Needs default")