text:
"""
Handles input coming through Http server requests.
"""
module Input

import HttpCommon, HTTP

export post, files, HttpInput, HttpPostData, HttpFiles, HttpFile

"""
    HttpFile

Represents a file sent over HTTP
"""
mutable struct HttpFile
  name::String
  mime::String
  data::Array{UInt8}
end
HttpFile() = HttpFile("", "", UInt8[])

const HttpPostData  = Dict{String, Union{String, Vector{String}}}
const HttpFiles     = Dict{String,HttpFile}

mutable struct HttpInput
  post::HttpPostData
  files::HttpFiles
end
HttpInput() = HttpInput(HttpPostData(), HttpFiles())

###

mutable struct HttpFormPart
  headers::Dict{String, Dict{String,String}}
  data::Array{UInt8}
end
HttpFormPart() = HttpFormPart(Dict{String,Dict{String,String}}(), UInt8[])

###

function all(request::HTTP.Request) :: HttpInput
  input::HttpInput = HttpInput()
  post_from_request!(request, input)

  input
end

function post(request::HTTP.Request)
  input::HttpInput = all(request)

  input.post
end

function files(request::HTTP.Request)
  input::HttpInput = all(request)

  input.files
end

###

function post_from_request!(request::HTTP.Request, input::HttpInput)
  headers = Dict(request.headers)

  if first(something(findfirst("application/x-www-form-urlencoded", get(headers, "Content-Type", "")), 0:-1)) != 0
    post_url_encoded!(request.body, input.post)
  elseif first(something(findfirst("multipart/form-data", get(headers, "Content-Type", "")), 0:-1)) != 0
    post_multipart!(request, input.post, input.files)
  end

  nothing
end

function post_url_encoded!(http_data::Array{UInt8, 1}, post_data::HttpPostData)
  if occursin("%5B%5D", String(copy(http_data))) || occursin("[]", String(copy(http_data))) # array values []
    for query_part in split(String(http_data), "&")
      qp = split(query_part, "=")
      (size(qp)[1] == 1) && (push!(qp, ""))

      k = HTTP.URIs.unescapeuri(HTTP.URIs.decodeplus(qp[1]))
      v = HTTP.URIs.unescapeuri(HTTP.URIs.decodeplus(qp[2]))
      # collect values like x[] in an array
      if endswith(k, "[]")
        if haskey(post_data, k)
          push!(post_data[k], v)
        else
          post_data[k] = [v]
        end
      else
        post_data[k] = v
      end
    end
  else
    params::Dict{String,String} = HTTP.URIs.queryparams(String(http_data))

    for (key::String, value::String) in params
      post_data[key] = value
    end
  end
end

function post_multipart!(request::HTTP.Request, post_data::HttpPostData, files::HttpFiles) :: Nothing
  headers = Dict(request.headers)
  boundary::String = headers["Content-Type"][(findfirst("boundary=", headers["Content-Type"])[end] + 1):end]

  boundary_length::Int = length(boundary)

  if boundary_length > 0
    form_parts::Array{HttpFormPart} = HttpFormPart[]

    get_multiform_parts!(request.body, form_parts, boundary, boundary_length)

    ### Process form parts

    if length(form_parts) > 0
      for part::HttpFormPart in form_parts
        hasFile::Bool = false
        file::HttpFile = HttpFile()
        fileFieldName::String = ""

        for (field::String, values::Dict{String,String}) in part.headers
          if field == "Content-Disposition" && getkey(values, "form-data", nothing) != nothing

            # Check to see whether this part is a file upload
            # Otherwise, treat as basic POST data

            if getkey(values, "filename", nothing) != nothing
              if length(values["filename"]) > 0
                fileFieldName = values["name"]
                file.name = values["filename"]
                hasFile = true
              end
            elseif getkey(values, "name", nothing) != nothing
              k = values["name"]
              v = String(part.data)

              if endswith(k, "[]")
                if haskey(post_data, k)
                  push!(post_data[k], v)
                else
                  post_data[k] = [v]
                end
              else
                post_data[k] = v
              end
            end
          elseif field == "Content-Type"
            (file.mime, mime) = first(values)
          end
        end # for

        if hasFile
          file.data = part.data

          files[fileFieldName] = file

          fileFieldName = ""
          file = HttpFile()
          hasFile = false
        end # if
      end # for
    end # if
  end

  nothing
end

###

function get_multiform_parts!(http_data::Vector{UInt8}, formParts::Array{HttpFormPart}, boundary, boundaryLength::Int = length(boundary))
  ### Go through each byte of data, parsing it into POST data and files.

  # According to the spec, the boundary chosen by the client must be a unique string
  # i.e. there should be no conflicts with the data within - so it should be safe to just do a basic string search.

  part::HttpFormPart = HttpFormPart()

  headerRaw::Array{UInt8} = UInt8[]

  captureAsData::Bool = false

  crOn::Bool = false
  hadLineEnding::Bool = false
  foundBoundary::Bool = false
  foundFinalBoundary::Bool = false

  bytes::Int = length(http_data)

  byteIndexOffset::Int = 0
  testIndex::Int = 1
  byteTestIndex::Int = 0

  byte::UInt8 = 0x00

  # Skip over the first boundary and CRLF

  byteIndex::Int = boundaryLength + 5

  while !foundFinalBoundary && byteIndex <= bytes
    byte = http_data[byteIndex]

    # Test for boundary.

    if (
      (byte == 0x0d && bytes >= byteIndex+3 && http_data[byteIndex + 1] == 0x0a && Char(http_data[byteIndex + 2]) == '-' && Char(http_data[byteIndex + 3]) == '-')
      || (byte == '-' && bytes >= byteIndex+1 && Char(http_data[byteIndex + 1]) == '-')
      )
      foundBoundary = true
    end

    if byte == 0x0d
      byteIndexOffset = byteIndex + 3
    else
      byteIndexOffset = byteIndex + 1
    end

    byteTestIndex = byteIndexOffset

    testIndex = 1;

    # Find the position of the next char NOT in the boundary
    if foundBoundary
      while testIndex < boundaryLength
        byteTestIndex = byteIndexOffset + testIndex

        if byteTestIndex > bytes || Char(http_data[byteTestIndex]) != boundary[testIndex]
          break
        end

        testIndex = testIndex + 1
      end
    end

    # Check if this boundary is the final one
    if foundBoundary
      if Char(http_data[byteTestIndex + 2]) == '-'
        foundFinalBoundary = true
        byteIndex = byteTestIndex + 5
      else
        byteIndex = byteTestIndex + 3
      end
    end

    ## Otherwise, process data

    if foundBoundary
      captureAsData = false
      crOn = false
      hadLineEnding = false
      foundBoundary = false

      push!(formParts, part)

      part = HttpFormPart()
    else
      if captureAsData
        push!(part.data, byte)
      else
        ## Check for CR

        if byte == 0x0d
          crOn = true
        else
          ## Check for LF and previous CR

          if byte == 0x0a && crOn
            ## Check for CRLFCRLF

            if hadLineEnding
              ## End of headers

              captureAsData = true

              hadLineEnding = false
            else
              ## End of single-line header

              header::String = String(headerRaw)

              if length(header) > 0
                headerParts = split(header, ": "; limit=2)

                valueDecoded = parse_semicolon_fields(String(headerParts[2]));

                if length(valueDecoded) > 0
                  part.headers[headerParts[1]] = valueDecoded
                end
              end

              headerRaw = UInt8[]

              hadLineEnding = true
            end
          else
            if hadLineEnding
              hadLineEnding = false
            end

            push!(headerRaw, byte)
          end

          crOn = false
        end
      end
    end

    byteIndex = byteIndex + 1
  end
end

###

function parse_semicolon_fields(dataString::String)
  dataString = dataString * ";"

  data = Dict{String,String}()

  prevCharacter::Char = 0x00
  inSingleQuotes::Bool = false
  inDoubleQuotes::Bool = false
  ignore::Bool = false
  workingString::String = ""

  dataStringLength::Int = length(dataString)

  dataStringLengthLoop::Int = dataStringLength + 1

  charIndex::Int = 1
  utfIndex::Int = 1 # real index for uft-8

  while charIndex < dataStringLengthLoop
    # character = dataString[charIndex]
    character = dataString[utfIndex]

    if ! inSingleQuotes && character == '"' && prevCharacter != '\\'
      inDoubleQuotes = ! inDoubleQuotes
      ignore = true
    end

    if ! inDoubleQuotes && character == '\'' && prevCharacter != '\\'
      inSingleQuotes = ! inSingleQuotes
      ignore = true
    end

    if charIndex == dataStringLength || (character == ';' && !(inSingleQuotes || inDoubleQuotes))
      workingString = strip(workingString)

      if length(workingString) > 0
        decoded = parse_quoted_params(workingString)

        if decoded != nothing
          (key, value) = decoded

          data[key] = value
        else
          data[workingString] = workingString
        end

        workingString = ""
      end
    elseif ! ignore
      workingString = workingString * string(character)
    end

    prevCharacter = character

    charIndex  = charIndex + 1

    utfIndex = nextind(dataString, utfIndex)   # real index for uft-8

    ignore = false
  end

  return data
end

function parse_quoted_params(data::String)
  tokens = split(data, "="; limit=2)

  if length(tokens) == 2
    return (tokens[1], tokens[2])
  end

  return nothing
end

###

end


text:
include("./imports.jl")

@testset ExtendedTestSet "Basic `Learner` convergence" begin
    learner = testlearner(coeff = 3)
    fit!(learner, repeat([TrainingPhase()], 5))
    @test  learner.model.coeff[1] ≈ 3 atol = 0.1
end

ext: jl


text:
using Oceananigans.Grids: ZDirection, validate_unit_vector

struct Buoyancy{M, G}
    model :: M
    gravity_unit_vector :: G
end

"""
    Buoyancy(; model, gravity_unit_vector=ZDirection())

Uses a given buoyancy `model` to create buoyancy in a model. The optional keyword argument
`gravity_unit_vector` can be used to specify the direction opposite to the gravitational
acceleration (which we take here to mean the "vertical" direction).

Example
=======

```julia
using Oceananigans

grid = RectilinearGrid(size=(1, 8, 8), extent=(1, 1000, 100))
θ = 45 # degrees
g̃ = (0, sind(θ), cosd(θ))

buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃)

model = NonhydrostaticModel(grid=grid, buoyancy=buoyancy, tracers=:b)
```
"""
function Buoyancy(; model, gravity_unit_vector=ZDirection())
    gravity_unit_vector = validate_unit_vector(gravity_unit_vector)
    return Buoyancy(model, gravity_unit_vector)
end

@inline ĝ_x(buoyancy) = @inbounds buoyancy.gravity_unit_vector[1]
@inline ĝ_y(buoyancy) = @inbounds buoyancy.gravity_unit_vector[2]
@inline ĝ_z(buoyancy) = @inbounds buoyancy.gravity_unit_vector[3]

@inline ĝ_x(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_y(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_z(::Buoyancy{M, ZDirection}) where M = 1

#####
##### For convenience
#####

@inline required_tracers(bm::Buoyancy) = required_tracers(bm.model)

@inline get_temperature_and_salinity(bm::Buoyancy, C) = get_temperature_and_salinity(bm.model, C)

@inline ∂x_b(i, j, k, grid, b::Buoyancy, C) = ∂x_b(i, j, k, grid, b.model, C)
@inline ∂y_b(i, j, k, grid, b::Buoyancy, C) = ∂y_b(i, j, k, grid, b.model, C)
@inline ∂z_b(i, j, k, grid, b::Buoyancy, C) = ∂z_b(i, j, k, grid, b.model, C)

@inline top_buoyancy_flux(i, j, grid, b::Buoyancy, args...) = top_buoyancy_flux(i, j, grid, b.model, args...)

regularize_buoyancy(b) = b
regularize_buoyancy(b::AbstractBuoyancyModel) = Buoyancy(model=b)

Base.summary(buoyancy::Buoyancy) = string(summary(buoyancy.model), " with -ĝ = ", summary(buoyancy.gravity_unit_vector))

Base.show(io::IO, buoyancy::Buoyancy) = summary(buoyancy)


text:
using ITensors
using Test
using Combinatorics: Combinatorics

digits(::Type{T}, i, j, k) where {T} = T(i * 10^2 + j * 10 + k)

@testset "ITensor $T Contractions" for T in (Float64, ComplexF64)
  mi, mj, mk, ml, mα = 2, 3, 4, 5, 6, 7
  i = Index(mi, "i")
  j = Index(mj, "j")
  k = Index(mk, "k")
  l = Index(ml, "l")
  α = Index(mα, "alpha")
  @testset "Test contract ITensors" begin
    A = randomITensor(T)
    B = randomITensor(T)
    Ai = randomITensor(T, i)
    Bi = randomITensor(T, i)
    Aj = randomITensor(T, j)
    Aij = randomITensor(T, i, j)
    Bij = randomITensor(T, i, j)
    Aik = randomITensor(T, i, k)
    Ajk = randomITensor(T, j, k)
    Ajl = randomITensor(T, j, l)
    Akl = randomITensor(T, k, l)
    Aijk = randomITensor(T, i, j, k)
    Ajkl = randomITensor(T, j, k, l)
    Aikl = randomITensor(T, i, k, l)
    Aklα = randomITensor(T, k, l, α)
    Aijkl = randomITensor(T, i, j, k, l)
    @testset "Test contract ITensor (Scalar*Scalar -> Scalar)" begin
      C = A * B
      @test scalar(C) ≈ scalar(A) * scalar(B)
    end
    @testset "Test contract ITensor (Scalar*Vector -> Vector)" begin
      C = A * Ai
      @test array(C) ≈ scalar(A) * array(Ai)
    end
    @testset "Test contract ITensor (Vector*Scalar -> Vector)" begin
      C = Aj * A
      @test array(C) ≈ scalar(A) * array(Aj)
    end
    @testset "Test contract ITensors (Vectorᵀ*Vector -> Scalar)" begin
      C = Ai * Bi
      CArray = transpose(array(Ai)) * array(Bi)
      @test CArray ≈ scalar(C)
    end
    @testset "Test contract ITensors (Vector*Vectorᵀ -> Matrix)" begin
      C = Ai * Aj
      for ii in 1:dim(i), jj in 1:dim(j)
        @test C[i => ii, j => jj] ≈ Ai[i => ii] * Aj[j => jj]
      end
    end
    @testset "Test contract ITensors (Matrix*Scalar -> Matrix)" begin
      Aij = permute(Aij, i, j)
      C = Aij * A
      @test array(permute(C, i, j)) ≈ scalar(A) * array(Aij)
    end
    @testset "Test contract ITensors (Matrix*Vector -> Vector)" begin
      Aij = permute(Aij, i, j)
      C = Aij * Aj
      CArray = array(permute(Aij, i, j)) * array(Aj)
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrixᵀ*Vector -> Vector)" begin
      Aij = permute(Aij, j, i)
      C = Aij * Aj
      CArray = transpose(array(Aij)) * array(Aj)
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Vector*Matrix -> Vector)" begin
      Aij = permute(Aij, i, j)
      C = Ai * Aij
      CArray = transpose(transpose(array(Ai)) * array(Aij))
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Vector*Matrixᵀ -> Vector)" begin
      Aij = permute(Aij, j, i)
      C = Ai * Aij
      CArray = transpose(transpose(array(Ai)) * transpose(array(Aij)))
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrix*Matrix -> Scalar)" begin
      Aij = permute(Aij, i, j)
      Bij = permute(Bij, i, j)
      C = Aij * Bij
      CArray = LinearAlgebra.tr(array(Aij) * transpose(array(Bij)))
      @test CArray ≈ scalar(C)
    end
    @testset "Test contract ITensors (Matrix*Matrix -> Matrix)" begin
      Aij = permute(Aij, i, j)
      Ajk = permute(Ajk, j, k)
      C = Aij * Ajk
      CArray = array(Aij) * array(Ajk)
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrixᵀ*Matrix -> Matrix)" begin
      Aij = permute(Aij, j, i)
      Ajk = permute(Ajk, j, k)
      C = Aij * Ajk
      CArray = transpose(array(Aij)) * array(Ajk)
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrix*Matrixᵀ -> Matrix)" begin
      Aij = permute(Aij, i, j)
      Ajk = permute(Ajk, k, j)
      C = Aij * Ajk
      CArray = array(Aij) * transpose(array(Ajk))
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrixᵀ*Matrixᵀ -> Matrix)" begin
      Aij = permute(Aij, j, i)
      Ajk = permute(Ajk, k, j)
      C = Aij * Ajk
      CArray = transpose(array(Aij)) * transpose(array(Ajk))
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (Matrix⊗Matrix -> 4-tensor)" begin
      C = Aij * Akl
      for ii in 1:dim(i), jj in 1:dim(j), kk in 1:dim(k), ll in 1:dim(l)
        @test C[i => ii, j => jj, k => kk, l => ll] ≈
              Aij[i => ii, j => jj] * Akl[k => kk, l => ll]
      end
    end
    @testset "Test contract ITensors (3-Tensor*Scalar -> 3-Tensor)" begin
      Aijk = permute(Aijk, i, j, k)
      C = Aijk * A
      @test array(permute(C, i, j, k)) ≈ scalar(A) * array(Aijk) rtol = 1e-12
    end
    @testset "Test contract ITensors (3-Tensor*Vector -> Matrix)" begin
      Aijk = permute(Aijk, i, j, k)
      C = Aijk * Ai
      CArray = reshape(
        reshape(array(permute(Aijk, j, k, i)), dim(j) * dim(k), dim(i)) * array(Ai),
        dim(j),
        dim(k),
      )
      @test CArray ≈ array(permute(C, j, k))
    end
    @testset "Test contract ITensors (Vector*3-Tensor -> Matrix)" begin
      Aijk = permute(Aijk, i, j, k)
      C = Aj * Aijk
      CArray = reshape(
        transpose(array(Aj)) *
        reshape(array(permute(Aijk, j, i, k)), dim(j), dim(i) * dim(k)),
        dim(i),
        dim(k),
      )
      @test CArray ≈ array(permute(C, i, k))
    end
    @testset "Test contract ITensors (3-Tensor*Matrix -> Vector)" begin
      Aijk = permute(Aijk, i, j, k)
      Aik = permute(Aik, i, k)
      C = Aijk * Aik
      CArray =
        reshape(array(permute(Aijk, j, i, k)), dim(j), dim(i) * dim(k)) * vec(array(Aik))
      @test CArray ≈ array(C)
    end
    @testset "Test contract ITensors (3-Tensor*Matrix -> 3-Tensor)" begin
      Aijk = permute(Aijk, i, j, k)
      Ajl = permute(Ajl, j, l)
      C = Aijk * Ajl
      CArray = reshape(
        reshape(array(permute(Aijk, i, k, j)), dim(i) * dim(k), dim(j)) * array(Ajl),
        dim(i),
        dim(k),
        dim(l),
      )
      @test CArray ≈ array(permute(C, i, k, l))
    end
    @testset "Test contract ITensors (Matrix*3-Tensor -> 3-Tensor)" begin
      Aijk = permute(Aijk, i, j, k)
      Akl = permute(Akl, k, l)
      C = Akl * Aijk
      CArray = reshape(
        array(permute(Akl, l, k)) *
        reshape(array(permute(Aijk, k, i, j)), dim(k), dim(i) * dim(j)),
        dim(l),
        dim(i),
        dim(j),
      )
      @test CArray ≈ array(permute(C, l, i, j))
    end
    @testset "Test contract ITensors (3-Tensor*3-Tensor -> 3-Tensor)" begin
      Aijk = permute(Aijk, i, j, k)
      Ajkl = permute(Ajkl, j, k, l)
      C = Aijk * Ajkl
      CArray =
        reshape(array(permute(Aijk, i, j, k)), dim(i), dim(j) * dim(k)) *
        reshape(array(permute(Ajkl, j, k, l)), dim(j) * dim(k), dim(l))
      @test CArray ≈ array(permute(C, i, l))
    end
    @testset "Test contract ITensors (3-Tensor*3-Tensor -> 3-Tensor)" begin
      for inds_ijk in Combinatorics.permutations([i, j, k]),
        inds_jkl in Combinatorics.permutations([j, k, l])

        Aijk = permute(Aijk, inds_ijk...)
        Ajkl = permute(Ajkl, inds_jkl...)
        C = Ajkl * Aijk
        CArray =
          reshape(array(permute(Ajkl, l, j, k)), dim(l), dim(j) * dim(k)) *
          reshape(array(permute(Aijk, j, k, i)), dim(j) * dim(k), dim(i))
        @test CArray ≈ array(permute(C, l, i))
      end
    end
    @testset "Test contract ITensors (4-Tensor*3-Tensor -> 1-Tensor)" begin
      for inds_ijkl in Combinatorics.permutations([i, j, k, l]),
        inds_jkl in Combinatorics.permutations([j, k, l])

        Aijkl = permute(Aijkl, inds_ijkl...)
        Ajkl = permute(Ajkl, inds_jkl...)
        C = Ajkl * Aijkl
        CArray =
          reshape(array(permute(Ajkl, j, k, l)), 1, dim(j) * dim(k) * dim(l)) *
          reshape(array(permute(Aijkl, j, k, l, i)), dim(j) * dim(k) * dim(l), dim(i))
        @test vec(CArray) ≈ array(permute(C, i))
      end
    end
    @testset "Test contract ITensors (4-Tensor*3-Tensor -> 3-Tensor)" begin
      for inds_ijkl in Combinatorics.permutations([i, j, k, l]),
        inds_klα in Combinatorics.permutations([k, l, α])

        Aijkl = permute(Aijkl, inds_ijkl...)
        Aklα = permute(Aklα, inds_klα...)
        C = Aklα * Aijkl
        CArray = reshape(
          reshape(array(permute(Aklα, α, k, l)), dim(α), dim(k) * dim(l)) *
          reshape(array(permute(Aijkl, k, l, i, j)), dim(k) * dim(l), dim(i) * dim(j)),
          dim(α),
          dim(i),
          dim(j),
        )
        @test CArray ≈ array(permute(C, α, i, j))
      end
    end
  end # End contraction testset
end

@testset "Contraction conversions" begin
  @testset "Real scalar * Complex ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    x = rand(Float64)
    A = randomITensor(ComplexF64, i, j)
    B = x * A
    for ii in dim(i), jj in dim(j)
      @test B[i => ii, j => jj] == x * A[i => ii, j => jj]
    end
  end
  @testset "Complex scalar * Real ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    x = rand(ComplexF64)
    A = randomITensor(Float64, i, j)
    B = x * A
    for ii in dim(i), jj in dim(j)
      @test B[i => ii, j => jj] == x * A[i => ii, j => jj]
    end
  end
  @testset "Real ITensor * Complex ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    k = Index(2, "k")
    A = randomITensor(Float64, i, j)
    B = randomITensor(ComplexF64, j, k)
    C = A * B
    @test array(permute(C, i, k)) ≈ array(A) * array(B)
  end
  @testset "Complex ITensor * Real ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    k = Index(2, "k")
    A = randomITensor(ComplexF64, i, j)
    B = randomITensor(Float64, j, k)
    C = A * B
    @test array(permute(C, i, k)) ≈ array(A) * array(B)
  end

  @testset "Outer Product Real ITensor * Complex ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    A = randomITensor(Float64, i)
    B = randomITensor(ComplexF64, j)
    C = A * B
    @test array(permute(C, i, j)) ≈ kron(array(A), transpose(array(B)))
  end

  @testset "Outer Product: Complex ITensor * Real ITensor" begin
    i = Index(2, "i")
    j = Index(2, "j")
    A = randomITensor(ComplexF64, i)
    B = randomITensor(Float64, j)
    C = A * B
    @test array(permute(C, i, j)) ≈ kron(array(A), transpose(array(B)))
  end
end

nothing


text:
#######################################################################################################################################################################################################
#
# Changes to this type
# General
#     2022-Jan-13: migrate abstract temperature dependency type from Photosynthesis.jl
#     2022-Jan-25: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Hierachy of AbstractTemperatureDependency:
- [`Arrhenius`](@ref)
- [`ArrheniusPeak`](@ref)
- [`Q10`](@ref)
"""
abstract type AbstractTemperatureDependency{FT<:AbstractFloat} end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-13: migrate from Photosynthesis.jl, rename to Arrhenius
#     2022-Jan-13: define the struct mutable, use ΔHA directly in the struct, add field T_REF
#     2022-Jan-25: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

An [`Arrhenius`](@ref) type struct using
```math
Y_1 = Y_0 \\cdot \\exp \\left( \\dfrac{H_a}{R T_0} - \\dfrac{H_a}{R T_1} \\right)
```

# Fields

$(TYPEDFIELDS)

"""
mutable struct Arrhenius{FT<:AbstractFloat} <: AbstractTemperatureDependency{FT}
    # parameters that do not change with time
    "Reference temperature `[K]`"
    T_REF::FT
    "Uncorrected vakye at reference temperature"
    VAL_REF::FT
    "Activation energy"
    ΔHA::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-13: migrate from Photosynthesis.jl, rename to ArrheniusPeak
#     2022-Jan-13: define the struct mutable, use ΔHA/ΔHD/ΔSV directly in the struct, add field T_REF/VAL_REF
#     2022-Jan-25: fix documentation
#     2022-Mar-01: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

An [`ArrheniusPeak`](@ref) type struct using
```math
Y_1 = Y_0 \\cdot \\exp \\left( \\dfrac{H_a}{R T_0} - \\dfrac{H_a}{R T_1} \\right)
          \\cdot \\dfrac{ 1 + \\exp \\left( \\dfrac{S_v T_0 - H_d}{R T_0} \\right) }
                        { 1 + \\exp \\left( \\dfrac{S_v T_1 - H_d}{R T_1} \\right) }
```

# Fields

$(TYPEDFIELDS)

"""
mutable struct ArrheniusPeak{FT<:AbstractFloat} <: AbstractTemperatureDependency{FT}
    # parameters that do not change with time
    "Reference temperature `[K]`"
    T_REF::FT
    "Uncorrected vakye at reference temperature"
    VAL_REF::FT
    "Activation energy"
    ΔHA::FT
    "Deactivation energy"
    ΔHD::FT
    "Entropy factor"
    ΔSV::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-13: migrate from Photosynthesis.jl, rename to Q10
#     2022-Jan-14: make structure mutable
#     2022-Jan-25: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

An [`Q10`](@ref) type struct using
```math
Y_1 = Y_0 \\cdot Q_{10} ^ \\dfrac{T_1 - T_0}{10}
```

# Fields

$(TYPEDFIELDS)

"""
mutable struct Q10{FT<:AbstractFloat} <: AbstractTemperatureDependency{FT}
    # parameters that do not change with time
    "Reference temperature `[K]`"
    T_REF::FT
    "Uncorrected vakye at reference temperature"
    VAL_REF::FT
    "Power of Q10 correction"
    Q_10::FT
end

#######################################################################################################################################################################################################
#
# Changes to the constructors
# General
#     2022-Jan-14: migrate from Photosynthesis.jl
#     2022-Feb-11: add temperature dependent Jmax and Vcmax TD from CLM
# Sources
#     Lavigne and Ryan (1997) Growth and maintenance respiration rates of aspen, blackspruce and jack pine stems at northern and southern BOREAS sites
#     Bernacchi et al. (2001) Improved temperature response functions for models of Rubisco‐limited photosynthesis
#     Boyd et al. (2001) Temperature responses of C4 photosynthesis: biochemical analysis of Rubisco, phosphoenolpyruvate carboxylase, and carbonic anhydrase in Setaria viridis
#     Leuning (2002) Temperature dependence of two parameters in a photosynthesis model
#     Kattge et al. (2007) Temperature acclimation in a biochemical model of photosynthesis: a reanalysis of data from 36 species
#     Sperry et al. (2019) The impact of rising CO2 and acclimation on the response of US forests to global warming
#     Johnson et al. (2021) The limiting factors and regulatory processes that control the environmental responses of C3, C3–C4 intermediate, and C4 photosynthesis
#     CLM5 Documentation. Chapter 9 Page 106
#
#######################################################################################################################################################################################################
KcTDBernacchi(FT)          = Arrhenius{FT}(T_25(), 41.0264925, 79430.0);
KcTDCLM(FT)                = Arrhenius{FT}(T_25(), 40.49     , 79430.0);
KoTDBernacchi(FT)          = Arrhenius{FT}(T_25(), 28208.88  , 36380.0);
KoTDCLM(FT)                = Arrhenius{FT}(T_25(), 27840.0   , 36380.0);
KpepTDCLM(FT)              = Arrhenius{FT}(T_25(), 8.0       , 36000.0);
KpepTDBoyd(FT)             = Arrhenius{FT}(T_25(), 16.0      , 36300.0);
KqTDJohnson(FT)            = Arrhenius{FT}(T_25(), 300       , 37000.0);
RespirationTDBernacchi(FT) = Arrhenius{FT}(T_25(), NaN       , 46390.0);
VcmaxTDBernacchi(FT)       = Arrhenius{FT}(T_25(), NaN       , 65330.0);
VomaxTDBernacchi(FT)       = Arrhenius{FT}(T_25(), NaN       , 60110.0);
ΓStarTDBernacchi(FT)       = Arrhenius{FT}(T_25(), 4.33164375, 37830.0);
ΓStarTDCLM(FT)             = Arrhenius{FT}(T_25(), 4.275     , 37830.0);

JmaxTDBernacchi(FT)                = ArrheniusPeak{FT}(T_25(), NaN , 57500.0, 439000.0, 1400.0);
JmaxTDCLM(FT, t::Number = T_25())  = ArrheniusPeak{FT}(T_25(), NaN , 50000.0, 200000.0, 659.70 - 0.75 * (t - T_0()) );
JmaxTDLeuning(FT)                  = ArrheniusPeak{FT}(T_25(), NaN , 50300.0, 152044.0, 495.0 );
RespirationTDCLM(FT)               = ArrheniusPeak{FT}(T_25(), NaN , 46390.0, 150650.0, 490.0 );
VcmaxTDCLM(FT, t::Number = T_25()) = ArrheniusPeak{FT}(T_25(), NaN , 72000.0, 200000.0, 668.39 - 1.07 * (t - T_0()) );
VcmaxTDLeuning(FT)                 = ArrheniusPeak{FT}(T_25(), NaN , 73637.0, 149252.0, 486.0 );
VpmaxTDBoyd(FT)                    = ArrheniusPeak{FT}(T_25(), NaN , 94800.0, 73300.0 , 250.0 );
ΗCTDJohnson(FT)                    = ArrheniusPeak{FT}(T_25(), 1.0 , 0.0    , 220000.0, 710.0 );
ΗLTDJohnson(FT)                    = ArrheniusPeak{FT}(T_25(), 0.75, 0.0    , 220000.0, 710.0 );

Q10TDAngiosperm(FT) = Q10{FT}(T_25(), 0.0140/8760, 1.4);
Q10TDGymnosperm(FT) = Q10{FT}(T_25(), 0.0425/8760, 1.7);


text:
using Oceananigans.BoundaryConditions: Flux, Value, Gradient, flip, BoundaryCondition, ContinuousBoundaryFunction
using Oceananigans.BoundaryConditions: getbc, regularize_boundary_condition, LeftBoundary, RightBoundary
using Oceananigans.BoundaryConditions: DefaultBoundaryCondition
using Oceananigans.TurbulenceClosures: AbstractScalarDiffusivity, h_diffusivity, z_diffusivity
using Oceananigans.Operators: index_left, index_right, Δx, Δy, Δz, div

import Oceananigans.BoundaryConditions: regularize_immersed_boundary_condition, bc_str

struct ImmersedBoundaryCondition{W, E, S, N, B, T}
    west :: W
    east :: E
    south :: S
    north :: N
    bottom :: B
    top :: T
end

const IBC = ImmersedBoundaryCondition

bc_str(::IBC) = "ImmersedBoundaryCondition"

Base.summary(ibc::IBC) =
    string(bc_str(ibc), " with ",
           "west=", bc_str(ibc.west), ", ",
           "east=", bc_str(ibc.east), ", ",
           "south=", bc_str(ibc.south), ", ",
           "north=", bc_str(ibc.north), ", ",
           "bottom=", bc_str(ibc.bottom), ", ",
           "top=", bc_str(ibc.top))

Base.show(io::IO, ibc::IBC) =
    print(io, "ImmersedBoundaryCondition:", '\n',
              "├── west: ", summary(ibc.west), '\n',
              "├── east: ", summary(ibc.east), '\n',
              "├── south: ", summary(ibc.south), '\n',
              "├── north: ", summary(ibc.north), '\n',
              "├── bottom: ", summary(ibc.bottom), '\n',
              "└── top: ", summary(ibc.top))

"""
    ImmersedBoundaryCondition(; interfaces...)

Return an ImmersedBoundaryCondition with conditions on individual
cell `interfaces ∈ (west, east, south, north, bottom, top)`
between the fluid and immersed boundary.
"""
function ImmersedBoundaryCondition(; west = nothing,
                                     east = nothing,
                                     south = nothing,
                                     north = nothing,
                                     bottom = nothing,
                                     top = nothing)

    return ImmersedBoundaryCondition(west, east, south, north, bottom, top)
end

#####
##### Nothing and FluxBoundaryCondition.
#####
##### Very Important Note: For FluxBoundaryCondition,
##### we assume fluxes are directed along the "inward-facing normal".
##### For example, east_immersed_flux = - user_flux.
##### With this convention, positive fluxes _increase_ boundary-adjacent
##### cell values, and negative fluxes _decrease_ them.
#####

for side in (:west, :south, :bottom)
    side_ib_flux = Symbol(side, :_ib_flux)
    @eval begin
        @inline $side_ib_flux(i, j, k, ibg, ::Nothing, args...) = zero(eltype(ibg))
        @inline $side_ib_flux(i, j, k, ibg, bc::FBC, loc, c, closure, K, id, args...) = + getbc(bc, i, j, k, ibg, args...)
    end
end

for side in (:east, :north, :top)
    side_ib_flux = Symbol(side, :_ib_flux)
    @eval begin
        @inline $side_ib_flux(i, j, k, ibg, ::Nothing, args...) = zero(eltype(ibg))
        @inline $side_ib_flux(i, j, k, ibg, bc::FBC, loc, c, closure, K, id, args...) = - getbc(bc, i, j, k, ibg, args...)
    end
end

#####
##### Value boundary condition. It's harder!
#####

# Harder in some ways... ValueBoundaryCondition...
const VBC = BoundaryCondition{Value}
const GBC = BoundaryCondition{Gradient}
const VBCorGBC = Union{VBC, GBC}
const ASD = AbstractScalarDiffusivity

# "Gradient" utility for Value or Gradient boundary conditions
@inline right_gradient(i, j, k, ibg, κ, Δ, bc::GBC, c, clock, fields) = getbc(bc, i, j, k, ibg, clock, fields)
@inline left_gradient(i, j, k, ibg, κ, Δ, bc::GBC, c, clock, fields)  = getbc(bc, i, j, k, ibg, clock, fields)

@inline function right_gradient(i, j, k, ibg, κ, Δ, bc::VBC, c, clock, fields)
    cᵇ = getbc(bc, i, j, k, ibg, clock, fields)
    cⁱʲᵏ = @inbounds c[i, j, k]
    return 2 * (cᵇ - cⁱʲᵏ) / Δ
end

@inline function left_gradient(i, j, k, ibg, κ, Δ, bc::VBC, c, clock, fields)
    cᵇ = getbc(bc, i, j, k, ibg, clock, fields)
    cⁱʲᵏ = @inbounds c[i, j, k]
    return 2 * (cⁱʲᵏ - cᵇ) / Δ
end

# Metric and index gymnastics for the 6 facets of the cube

@inline function _west_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δx(index_left(i, LX), j, k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock)
    ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

@inline function _east_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δx(index_right(i, LX), j, k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock)
    ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

@inline function _south_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δy(i, index_left(j, LY), k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock)
    ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

@inline function _north_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δy(i, index_right(j, LY), k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock)
    ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

@inline function _bottom_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δz(i, j, index_left(k, LZ), ibg, LX, LY, LZ)
    κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock)
    ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

@inline function _top_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δz(i, j, index_right(k, LZ), ibg, LX, LY, LZ)
    κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock)
    ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

sides = [:west, :east, :south, :north, :bottom, :top]

for side in sides
    flux = Symbol(side, "_ib_flux")
    _flux = Symbol("_", flux)

    @eval begin
        @inline $flux(i, j, k, ibg, bc::VBCorGBC, args...) = $_flux(i, j, k, ibg, bc::VBCorGBC, args...)
        @inline $_flux(i, j, k, ibg, bc::VBCorGBC, args...) = zero(ibg) # fallback for non-ASD closures

        @inline $flux(i, j, k, ibg, bc::VBCorGBC, loc, c, closures::Tuple{<:Any}, Ks, id, clock, fields) =
            $_flux(i, j, k, ibg, bc, loc, c, closures[1], Ks[1], id, clock, fields)

        @inline $flux(i, j, k, ibg, bc::VBCorGBC, loc, c, closures::Tuple{<:Any, <:Any}, Ks, id, clock, fields) =
            $_flux(i, j, k, ibg, bc, loc, c, closures[1], Ks[1], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[2], Ks[2], id, clock, fields)

        @inline $flux(i, j, k, ibg, bc::VBCorGBC, loc, c, closures::Tuple{<:Any, <:Any, <:Any}, Ks, id, clock, fields) =
            $_flux(i, j, k, ibg, bc, loc, c, closures[1], Ks[1], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[2], Ks[2], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[3], Ks[3], id, clock, fields)

        @inline $flux(i, j, k, ibg, bc::VBCorGBC, loc, c, closures::Tuple{<:Any, <:Any, <:Any, <:Any}, Ks, id, clock, fields) =
            $_flux(i, j, k, ibg, bc, loc, c, closures[1], Ks[1], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[2], Ks[2], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[3], Ks[3], id, clock, fields) +
            $_flux(i, j, k, ibg, bc, loc, c, closures[4], Ks[4], id, clock, fields)

        @inline $flux(i, j, k, ibg, bc::VBCorGBC, loc, c, closures::Tuple, Ks, id, clock, fields) =
            $_flux(i, j, k, ibg, bc, loc, c, closures[1], Ks[1], id, clock, fields) +
             $flux(i, j, k, ibg, bc, loc, c, closures[2:end], Ks[2:end], id, clock, fields)
    end
end

#####
##### Immersed flux divergence
#####

@inline function immersed_flux_divergence(i, j, k, ibg::GFIBG, bc, loc, c, closure, K, id, clock, fields)
    # Fetch fluxes across immersed boundary
    q̃ᵂ =   west_ib_flux(i, j, k, ibg, bc.west,   loc, c, closure, K, id, clock, fields)
    q̃ᴱ =   east_ib_flux(i, j, k, ibg, bc.east,   loc, c, closure, K, id, clock, fields)
    q̃ˢ =  south_ib_flux(i, j, k, ibg, bc.south,  loc, c, closure, K, id, clock, fields)
    q̃ᴺ =  north_ib_flux(i, j, k, ibg, bc.north,  loc, c, closure, K, id, clock, fields)
    q̃ᴮ = bottom_ib_flux(i, j, k, ibg, bc.bottom, loc, c, closure, K, id, clock, fields)
    q̃ᵀ =    top_ib_flux(i, j, k, ibg, bc.top,    loc, c, closure, K, id, clock, fields)

    iᵂ, jˢ, kᴮ = index_right.((i, j, k), loc)
    iᴱ, jᴺ, kᵀ = index_left.((i, j, k), loc)
    LX, LY, LZ = loc

    # Impose i) immersed fluxes if we're on an immersed boundary or ii) zero otherwise.
    qᵂ = conditional_flux(iᵂ, j, k, ibg, flip(LX), LY, LZ, q̃ᵂ, zero(eltype(ibg)))
    qᴱ = conditional_flux(iᴱ, j, k, ibg, flip(LX), LY, LZ, q̃ᴱ, zero(eltype(ibg)))
    qˢ = conditional_flux(i, jˢ, k, ibg, LX, flip(LY), LZ, q̃ˢ, zero(eltype(ibg)))
    qᴺ = conditional_flux(i, jᴺ, k, ibg, LX, flip(LY), LZ, q̃ᴺ, zero(eltype(ibg)))
    qᴮ = conditional_flux(i, j, kᴮ, ibg, LX, LY, flip(LZ), q̃ᴮ, zero(eltype(ibg)))
    qᵀ = conditional_flux(i, j, kᵀ, ibg, LX, LY, flip(LZ), q̃ᵀ, zero(eltype(ibg)))

    return div(i, j, k, ibg, loc, qᵂ, qᴱ, qˢ, qᴺ, qᴮ, qᵀ)
end

@inline immersed_∂ⱼ_τ₁ⱼ(i, j, k, ibg::GFIBG, U, u_bc::IBC, closure, K, clock, fields) =
    immersed_flux_divergence(i, j, k, ibg, u_bc, (f, c, c), U.u, closure, K, nothing, clock, fields)

@inline immersed_∂ⱼ_τ₂ⱼ(i, j, k, ibg::GFIBG, U, v_bc::IBC, closure, K, clock, fields) =
    immersed_flux_divergence(i, j, k, ibg, v_bc, (c, f, c), U.v, closure, K, nothing, clock, fields)

@inline immersed_∂ⱼ_τ₃ⱼ(i, j, k, ibg::GFIBG, U, w_bc::IBC, closure, K, clock, fields) =
    immersed_flux_divergence(i, j, k, ibg, w_bc, (c, c, f), U.w, closure, K, nothing, clock, fields)

@inline immersed_∇_dot_qᶜ(i, j, k, ibg::GFIBG, C, c_bc::IBC, closure, K, id, clock, fields) =
    immersed_flux_divergence(i, j, k, ibg, c_bc, (c, c, c), C, closure, K, id, clock, fields)

#####
##### Boundary condition "regularization"
#####

const ZFBC = BoundaryCondition{Flux, Nothing}
regularize_immersed_boundary_condition(ibc::ZFBC, ibg::GFIBG, args...) = ibc # keep it

# Compiler hint
@inline immersed_flux_divergence(i, j, k, ibg::GFIBG, bc::ZFBC, loc, c, closure, K, id, clock, fields) = zero(ibg)

regularize_immersed_boundary_condition(default::DefaultBoundaryCondition, ibg::GFIBG, loc, field_name, args...) =
    regularize_immersed_boundary_condition(default.boundary_condition, ibg, loc, field_name, args...)

# Convert certain non-immersed boundary conditions to immersed boundary conditions
function regularize_immersed_boundary_condition(ibc::Union{VBC, GBC, FBC}, ibg::GFIBG, loc, field_name, args...)
    ibc = ImmersedBoundaryCondition(Tuple(ibc for i=1:6)...)
    regularize_immersed_boundary_condition(ibc, ibg, loc, field_name, args...)
end

"""
    regularize_immersed_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
                                           topo, loc, dim, I, prognostic_field_names) where C
"""
function regularize_immersed_boundary_condition(bc::IBC, grid, loc, field_name, prognostic_field_names)

    topo = topology(grid)

    west   = loc[1] === Face ? nothing : regularize_boundary_condition(bc.west,   topo, loc, 1, LeftBoundary,  prognostic_field_names)
    east   = loc[1] === Face ? nothing : regularize_boundary_condition(bc.east,   topo, loc, 1, RightBoundary, prognostic_field_names)
    south  = loc[2] === Face ? nothing : regularize_boundary_condition(bc.south,  topo, loc, 2, LeftBoundary,  prognostic_field_names)
    north  = loc[2] === Face ? nothing : regularize_boundary_condition(bc.north,  topo, loc, 2, RightBoundary, prognostic_field_names)
    bottom = loc[3] === Face ? nothing : regularize_boundary_condition(bc.bottom, topo, loc, 3, LeftBoundary,  prognostic_field_names)
    top    = loc[3] === Face ? nothing : regularize_boundary_condition(bc.top,    topo, loc, 3, RightBoundary, prognostic_field_names)

    return ImmersedBoundaryCondition(; west, east, south, north, bottom, top)
end

#####
##### Alternative implementation for immersed flux divergence
#####

#=
# Another idea...
# loc is the field location
# These are evaluated on both sides of a cell (eg left and right)
function immersed_flux_x(i, j, k, ibg, bc, loc, c, closure, K, id, args...)
    qᵂ = west_ib_flux(i, j, k, ibg, bc.west, loc, c, closure, K, id, args...)
    qᴱ = east_ib_flux(i, j, k, ibg, bc.east, loc, c, closure, K, id, args...)

    LX, LY, LZ = loc
    west_boundary = immersed_peripheral_node(flip(LX), LY, LZ, i, j, k, ibg) & !inactive_node(LX, LY, LZ, i, j, k, ibg)
    east_boundary = immersed_peripheral_node(flip(LX), LY, LZ, i, j, k, ibg) & !inactive_node(LX, LY, LZ, i-1, j, k, ibg)

    return ifelse(west_boundary, qᵂ, zero(ibg)) + ifelse(east_boundary, qᴱ, zero(ibg))
end

function immersed_flux_y(i, j, k, ibg, bc, loc, c, closure, K, id, args...)
    qˢ = south_ib_flux(i, j, k, ibg, bc.south, loc, c, closure, K, id, args...)
    qᴺ = north_ib_flux(i, j, k, ibg, bc.north, loc, c, closure, K, id, args...)

    LX, LY, LZ = loc
    south_boundary = immersed_peripheral_node(LX, flip(LY), LZ, i, j, k, ibg) & !inactive_node(LX, LY, LZ, i, j, k, ibg)
    north_boundary = immersed_peripheral_node(LX, flip(LY), LZ, i, j, k, ibg) & !inactive_node(LX, LY, LZ, i, j-1, k, ibg)

    return ifelse(south_boundary, qˢ, zero(ibg)) + ifelse(north_boundary, qᴺ, zero(ibg))
end

function immersed_flux_z(i, j, k, ibg, bc, loc, c, closure, K, id, args...)
    qᴮ = bottom_ib_flux(i, j, k, ibg, bc.bottom, loc, c, closure, K, id, args...)
    qᵀ =    top_ib_flux(i, j, k, ibg, bc.top,    loc, c, closure, K, id, args...)

    LX, LY, LZ = loc
    bottom_boundary = immersed_peripheral_node(LX, LY, flip(LZ), i, j, k, ibg) & !inactive_node(LX, LY, LZ, i, j, k, ibg)
    top_boundary    = immersed_peripheral_node(LX, LY, flip(LZ), i, j, k, ibg) & !inactive_node(LX, LY, LZ, i, j, k-1, ibg)

    return ifelse(bottom_boundary, qᴮ, zero(ibg)) + ifelse(top_boundary, qᵀ, zero(ibg))
end

@inline immersed_∂ⱼ_τ₁ⱼ(i, j, k, ibg::GFIBG, U, u_bc::IBC, closure, K, args...) =
    return 1/Vᶠᶜᶜ(i, j, k, ibg) * (δxᶠᵃᵃ(i, j, k, ibg, Ax_qᶜᶜᶜ, immersed_flux_x, u_bc, (f, c, c), U.u, closure, K, nothing, args...) +
                                   δyᵃᶜᵃ(i, j, k, ibg, Ax_qᶜᶜᶜ, immersed_flux_y, u_bc, (f, c, c), U.u, closure, K, nothing, args...) +
                                   δzᵃᵃᶜ(i, j, k, ibg, Az_qᶠᶜᶠ, immersed_flux_z, u_bc, (f, c, c), U.u, closure, K, nothing, args...))

@inline immersed_∂ⱼ_τ₂ⱼ(i, j, k, ibg::GFIBG, U, v_bc::IBC, closure, K, args...) =
    return 1/Vᶜᶠᶜ(i, j, k, ibg) * (δxᶜᵃᵃ(i, j, k, ibg, Ax_qᶠᶠᶜ, immersed_flux_x, v_bc, (c, c, c), U.v, closure, K, nothing, args...) +
                                   δyᵃᶠᵃ(i, j, k, ibg, Ax_qᶜᶜᶜ, immersed_flux_y, v_bc, (c, f, c), U.v, closure, K, nothing, args...) +
                                   δzᵃᵃᶜ(i, j, k, ibg, Az_qᶜᶠᶠ, immersed_flux_z, v_bc, (c, f, c), U.v, closure, K, nothing, args...))

    immersed_flux_divergence(i, j, k, ibg, v_bc, (c, f, c), U.v, closure, K, nothing, args...)

@inline immersed_∂ⱼ_τ₃ⱼ(i, j, k, ibg::GFIBG, U, w_bc::IBC, closure, K, args...) =
    return 1/Vᶜᶜᶠ(i, j, k, ibg) * (δxᶜᵃᵃ(i, j, k, ibg, Ax_qᶠᶜᶠ, immersed_flux_x, w_bc, (c, c, f), U.w, closure, K, nothing, args...) +
                                   δyᵃᶜᵃ(i, j, k, ibg, Ax_qᶜᶠᶠ, immersed_flux_y, w_bc, (c, c, f), U.w, closure, K, nothing, args...) +
                                   δzᵃᵃᶠ(i, j, k, ibg, Az_qᶜᶜᶜ, immersed_flux_z, w_bc, (c, c, f), U.w, closure, K, nothing, args...))

@inline immersed_∇_dot_qᶜ(i, j, k, ibg::GFIBG, C, c_bc::IBC, closure, K, id, args...) =
    return 1/Vᶜᶜᶜ(i, j, k, ibg) * (δxᶜᵃᵃ(i, j, k, ibg, Ax_qᶠᶜᶜ, immersed_flux_x, c_bc, (c, c, c), U.u, closure, K, id, args...) +
                                   δyᵃᶜᵃ(i, j, k, ibg, Ax_qᶜᶠᶜ, immersed_flux_y, c_bc, (c, c, c), U.u, closure, K, id, args...) +
                                   δzᵃᵃᶜ(i, j, k, ibg, Az_qᶜᶜᶠ, immersed_flux_z, c_bc, (c, c, c), U.u, closure, K, id, args...))

=#


text:
# Generic operators (e.g. advection and diffusion, explicit/implicit)
# using building blocks:
# - difference operator
# - gradient
# - divergence
# - laplacian

# Central and upwind types
# struct Central  end; central  = Central()
# struct Upwind   end; upwind   = Upwind()
# struct Downwind end; downwind = Downwind()

############### 2d difference operators ##################
# Difference operators in x- and y-direction:
# Central-x
@inline δfδx_2dᶜ(f, i, j, Δx) = (f[i + 1,j] - f[i - 1,j]) / (2Δx)
# Upwind-x
@inline δfδx_2dᵘ(f, i, j, Δx) = (f[i,j] - f[i - 1,j]) / Δx
# Downwind-x
@inline δfδx_2dᵈ(f, i, j, Δx) = (f[i + 1,j] - f[i,j]) / Δx
# Central-y
@inline δfδy_2dᶜ(f, i, j, Δy) = (f[i,j + 1] - f[i,j - 1]) / (2Δy)
# Upwind-y
@inline δfδy_2dᵘ(f, i, j, Δy) = (f[i,j] - f[i,j - 1]) / Δy
# Downwind-y
@inline δfδy_2dᵈ(f, i, j, Δy) = (f[i,j + 1] - f[i,j]) / Δy

############### 3d difference operators ##################
# Difference operators in x- and y-direction:
# Central-x
@inline δfδx_3dᶜ(f, i, j, k, Δx)  = (f[i + 1,j,k] - f[i - 1,j,k]) / (2Δx)
# Upwind-x
@inline δfδx_3dᵘ(f, i, j, k, Δx)   = (f[i,j,k] - f[i - 1,j,k]) / Δx
# Downwind-x
@inline δfδx_3dᵈ(f, i, j, k, Δx) = (f[i + 1,j,k] - f[i,j,k]) / Δx
# Central-y)
@inline δfδy_3dᶜ(f, i, j, k, Δy)  = (f[i,j + 1,k] - f[i,j - 1,k]) / (2Δy)
# Upwind-y
@inline δfδy_3dᵘ(f, i, j, k, Δy)   = (f[i,j,k] - f[i,j - 1,k]) / Δy
# Downwind-y
@inline δfδy_3dᵈ(f, i, j, k, Δy) = (f[i,j + 1,k] - f[i,j,k]) / Δy
# Central-z
@inline δfδz_3dᶜ(f, i, j, k, Δz)  = (f[i,j,k + 1] - f[i,j,k - 1]) / (2Δz)
# Upwind-z
@inline δfδz_3dᵘ(f, i, j, k, Δz)   = (f[i,j,k] - f[i,j,k - 1]) / Δz
# Downwind-z
@inline δfδz_3dᵈ(f, i, j, k, Δz) = (f[i,j,k + 1] - f[i,j,k]) / Δz

############### 2d gradient operators ##################
# central gradient
@inline ∇2dᶜ(f, i, j, (Δx, Δy)) = (δfδx_2dᶜ(f, i, j, Δx), δfδy_2dᶜ(f, i, j, Δy))
# upwind gradient
@inline ∇2dᵘ(f, i, j, (Δx, Δy)) = (δfδx_2dᵘ(f, i, j, Δx), δfδy_2dᵘ(f, i, j, Δy))
# downwind gradient
@inline ∇2dᵈ(f, i, j, (Δx, Δy)) = (δfδx_2dᵈ(f, i, j, Δx), δfδy_2dᵈ(f, i, j, Δy))

############### 3d gradient operators ##################
# central gradient
@inline ∇3dᶜ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᶜ(f, i, j, k, Δx), δfδy_3dᶜ(f, i, j, k, Δy), δfδz_3dᶜ(f, i, j, k, Δz))
# upwind gradient
@inline ∇3dᵘ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᵘ(f, i, j, k, Δx), δfδy_3dᵘ(f, i, j, k, Δy), δfδz_3dᵘ(f, i, j, k, Δz))
# downwind gradient
@inline ∇3dᵈ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᵈ(f, i, j, k, Δx), δfδy_3dᵈ(f, i, j, k, Δy), δfδz_3dᵈ(f, i, j, k, Δz))

############### 2d divergence operators ##################
# central divergence
@inline div2dᶜ(f, i, j, (Δx, Δy)) = δfδx_2dᶜ(f, i, j, Δx) + δfδy_2dᶜ(f, i, j, Δy)
# upwind divergence
@inline div2dᵘ(f, i, j, (Δx, Δy)) = δfδx_2dᵘ(f, i, j, Δx) + δfδy_2dᵘ(f, i, j, Δy)
# downwind divergence
@inline div2dᵈ(f, i, j, (Δx, Δy)) = δfδx_2dᵈ(f, i, j, Δx) + δfδy_2dᵈ(f, i, j, Δy)

############### 3d divergence operators ##################
# central divergence
@inline div3dᶜ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᶜ(f, i, j, k, Δx) + δfδy_3dᶜ(f, i, j, k, Δy) + δfδz_3dᶜ(f, i, j, k, Δz))
# upwind divergence
@inline div3dᵘ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᵘ(f, i, j, k, Δx) + δfδy_3dᵘ(f, i, j, k, Δy) + δfδz_3dᵘ(f, i, j, k, Δz))
# downwind divergence
@inline div3dᵈ(f, i, j, k, (Δx, Δy, Δz)) = (δfδx_3dᵈ(f, i, j, k, Δx) + δfδy_3dᵈ(f, i, j, k, Δy) + δfδz_3dᵈ(f, i, j, k, Δz))

# Laplacian = div(gradient(f))
# central Laplacian
# @inline Δˣˣ(f,i,j,central)  = div(∇(f,i,j,central ),i,j,central)
# @inline Δˣˣ(f,i,j,upwind)   = div(∇(f,i,j,upwind  ),i,j,upwind)
# @inline Δˣˣ(f,i,j,downwind) = div(∇(f,i,j,downwind),i,j,downwind)


text:
module OpenCV

import Base.size

include("cv_cxx.jl")

include("cv_wrap.jl")

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    abstract type AbstractFunctionConversionBridge{F, S} <: AbstractBridge end

Abstract type to support writing bridges in which the function changes but the
set does not.

By convention, the transformed function is stored in the `.constraint` field.
"""
abstract type AbstractFunctionConversionBridge{F,S} <: AbstractBridge end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    bridge::AbstractFunctionConversionBridge,
)
    if !invariant_under_function_conversion(attr)
        throw(
            MOI.UnsupportedAttribute(
                attr,
                "Bridge of type `$(typeof(bridge))` does not support getting " *
                "the attribute `$attr` because " *
                "`MOIB.Constraint.invariant_under_function_conversion($attr)` " *
                "returns `false`.",
            ),
        )
    end
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    ::Type{<:AbstractFunctionConversionBridge{F,S}},
) where {F,S}
    return invariant_under_function_conversion(attr) &&
           MOI.supports(model, attr, MOI.ConstraintIndex{F,S})
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.AbstractConstraintAttribute,
    bridge::AbstractFunctionConversionBridge,
    value,
)
    if !invariant_under_function_conversion(attr)
        throw(
            MOI.UnsupportedAttribute(
                attr,
                "Bridge of type `$(typeof(bridge))` does not support setting " *
                "the attribute `$attr` because " *
                "`MOIB.Constraint.invariant_under_function_conversion($attr)` " *
                "returns `false`.",
            ),
        )
    end
    MOI.set(model, attr, bridge.constraint, value)
    return
end

"""
    invariant_under_function_conversion(attr::MOI.AbstractConstraintAttribute)

Returns whether the value of the attribute does not change if the constraint
`F`-in-`S` is transformed into a constraint `G`-in-`S` where `F` and `G` are
equivalent representations of the same function.

If it returns true, then subtypes of
[`Constraint.AbstractFunctionConversionBridge`](@ref) such as
[`Constraint.ScalarFunctionizeBridge`](@ref) and
[`Constraint.VectorFunctionizeBridge`](@ref) will automatically support
[`MOI.get`](@ref) and [`MOI.set`](@ref) for `attr`.
"""
invariant_under_function_conversion(::MOI.AbstractConstraintAttribute) = false

function invariant_under_function_conversion(
    ::Union{
        MOI.ConstraintSet,
        MOI.ConstraintBasisStatus,
        MOI.ConstraintPrimal,
        MOI.ConstraintPrimalStart,
        MOI.ConstraintDual,
        MOI.ConstraintDualStart,
    },
)
    return true
end


text:
"""
    Exponential(θ)

The *Exponential distribution* with scale parameter `θ` has probability density function

```math
f(x; \\theta) = \\frac{1}{\\theta} e^{-\\frac{x}{\\theta}}, \\quad x > 0
```

```julia
Exponential()      # Exponential distribution with unit scale, i.e. Exponential(1)
Exponential(b)     # Exponential distribution with scale b

params(d)          # Get the parameters, i.e. (b,)
scale(d)           # Get the scale parameter, i.e. b
rate(d)            # Get the rate parameter, i.e. 1 / b
```

External links

* [Exponential distribution on Wikipedia](http://en.wikipedia.org/wiki/Exponential_distribution)

"""
struct Exponential{T<:Real} <: ContinuousUnivariateDistribution
    θ::T        # note: scale not rate
    Exponential{T}(θ::T) where {T} = new{T}(θ)
end

function Exponential(θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Exponential, θ > zero(θ))
    return Exponential{T}(θ)
end

Exponential(θ::Integer) = Exponential(float(θ))
Exponential() = Exponential(1.0, check_args=false)

@distr_support Exponential 0.0 Inf

### Conversions
convert(::Type{Exponential{T}}, θ::S) where {T <: Real, S <: Real} = Exponential(T(θ))
convert(::Type{Exponential{T}}, d::Exponential{S}) where {T <: Real, S <: Real} = Exponential(T(d.θ), check_args=false)

#### Parameters

scale(d::Exponential) = d.θ
rate(d::Exponential) = inv(d.θ)

params(d::Exponential) = (d.θ,)
partype(::Exponential{T}) where {T<:Real} = T

#### Statistics

mean(d::Exponential) = d.θ
median(d::Exponential) = logtwo * d.θ
mode(::Exponential{T}) where {T<:Real} = zero(T)

var(d::Exponential) = d.θ^2
skewness(::Exponential{T}) where {T} = T(2)
kurtosis(::Exponential{T}) where {T} = T(6)

entropy(d::Exponential{T}) where {T} = one(T) + log(d.θ)

#### Evaluation

zval(d::Exponential, x::Real) = x / d.θ
xval(d::Exponential, z::Real) = z * d.θ

pdf(d::Exponential, x::Real) = (λ = rate(d); x < 0 ? zero(λ) : λ * exp(-λ * x))
function logpdf(d::Exponential{T}, x::Real) where T<:Real
    λ = rate(d)
    x < 0 ? -T(Inf) : log(λ) - λ * x
end

cdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -expm1(-zval(d, x)) : zero(T)
ccdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? exp(-zval(d, x)) : zero(T)
logcdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? log1mexp(-zval(d, x)) : -T(Inf)
logccdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -zval(d, x) : zero(T)

quantile(d::Exponential, p::Real) = -xval(d, log1p(-p))
cquantile(d::Exponential, p::Real) = -xval(d, log(p))
invlogcdf(d::Exponential, lp::Real) = -xval(d, log1mexp(lp))
invlogccdf(d::Exponential, lp::Real) = -xval(d, lp)

gradlogpdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -rate(d) : zero(T)

mgf(d::Exponential, t::Real) = 1/(1 - t * scale(d))
cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))

#### Sampling
rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))

#### Fit model

struct ExponentialStats <: SufficientStats
    sx::AbstractFloat   # (weighted) sum of x
    sw::AbstractFloat   # sum of sample weights

    ExponentialStats(sx::Real, sw::Real) = new(sx, sw)
end

suffstats(::Type{<:Exponential}, x::AbstractArray{T}) where {T<:Real} = ExponentialStats(sum(x), length(x))
suffstats(::Type{<:Exponential}, x::AbstractArray{T}, w::AbstractArray) where {T<:Real} = ExponentialStats(dot(x, w), sum(w))

fit_mle(::Type{<:Exponential}, ss::ExponentialStats) = Exponential(ss.sx / ss.sw)


text:
function ctx_init(context, flags=0)
  ccall((:grn_ctx_init, libgroonga),
        Ptr{Void}, (Ptr{Void}, Int,), context, flags)
end

function ctx_fin(context, flags=0)
  ccall((:grn_ctx_fin, libgroonga),
        Ptr{Void}, (Ptr{Void}, Int,), context, flags)
end

function ctx_open(flags=0)
  ccall((:grn_ctx_open, libgroonga),
        Ptr{Void}, (Int,), flags)
end

function ctx_use(context, database)
  ccall((:grn_ctx_use, libgroonga),
        Ptr{Void}, (Ptr{Void}, Ptr{Void},), context, database)
end

function ctx_db(context)
  ccall((:grn_ctx_db, libgroonga),
        Ptr{Void}, (Ptr{Void},), context)
end

function ctx_get(context, command)
  ccall((:grn_ctx_get, libgroonga),
        Int,
        (Ptr{Void}, Ptr{UInt8}, Int),
        context,
        bytestring(command),
        length(command))
end

function ctx_send(context, command)
  ccall((:grn_ctx_send, libgroonga),
        Int,
        (Ptr{Void}, Ptr{UInt8}, Int, Int),
        context,
        bytestring(command),
        length(command),
        0)
end

function ctx_recv(context)
  if Groonga.ctx_db(context) == C_NULL
    error("Database is not opened.")
  end

  str = Array(Ptr{UInt8})
  str_len = Array(Int32)
  flags = Array(Int)
  ccall((:grn_ctx_recv, libgroonga),
        Void,
        (Ptr{Void}, Ptr{Ptr{UInt8}}, Ptr{Int32}, Ptr{Int}),
        context,
        str,
        str_len,
        flags)
  if str[] == C_NULL
    return nothing
  end
  bytestring(str[], str_len[])
end


text:
module LightGaffes

function generate_scenario()
    subject = ["grandma", "I", "mom"]
    verb = ["tripped", "insulted", "spilled coffee on", "picked my nose in front of"]
    object = ["my teacher", "the cashier", "my kid", "the president"]

    print(subject[Int(ceil(rand()*length(subject)))], " ",
        verb[Int(ceil(rand()*length(verb)))], " ",
        object[Int(ceil(rand()*length(subject)))])

end

end # module


text:
"""
    ConstituencyTree{T}

Constituency parse tree of a natural language sentence.
"""
struct ConstituencyTree{T}
    node::T
    children
end

const Tree = ConstituencyTree

Tree{T}(label) where T = Tree{T}(label, [])
Tree() = Tree{Any}(nothing, [])
Tree(node) = Tree(node, [])
Tree(node::String, child::String) = Tree(node, [child])

AbstractTrees.children(tree::Tree) = tree.children
AbstractTrees.printnode(io::IO, tree::Tree) =
    print(io, ifelse(tree.node === nothing, "", tree.node))

height(tree) = isleaf(tree) ? 1 : 1 + maximum(height.(children(tree)))

isleaf(args...)    = true
isleaf(tree::Tree) = isempty(tree.children)

isterminal(x) = false
isterminal(tree::Tree) = length(tree.children) == 1 && isempty(children(tree.children[1]))

label(tree::Tree{Nothing}) = ""
label(tree::Tree) = tree.node
label(x::String) = x

function production(tree; nonterminal=identity, terminal=identity)
    lhs, rhs = label(tree), label.(children(tree))
    if isterminal(tree)
        return nonterminal(lhs), terminal.(rhs)
    else
        return nonterminal(lhs), nonterminal.(rhs)
    end
end

"""
    productions(tree; search=PreOrderDFS, nonterminal=identity, terminal=identity)

Return a vector of (lhs, rhs) productions from a constituency parse tree.

# Arguments

- `tree`: the tree to search

# Keywords

- `nonterminal`: a function to call on nonterminal symbols.
- `terminal`: a function to call on terminal symbols.

# Returns

- a `Vector` of (lhs, rhs) tuples

"""
function productions(tree; search=PreOrderDFS, ks...)
    ps = [production(p; ks...) for p in search(tree)]
    return [p for p in ps if p[2] != ()]
end

function Base.show(io::IO, tree::ConstituencyTree)
    print(io, typeof(tree))
    pprint(io, tree; multiline=false)
end

Base.:(==)(t1::Tree, t2::Tree) = label(t1) == label(t2) && children(t1) == children(t2)

"""
    pprint([io::IO,] tree; indent=2, multiline=true)

Print a constituency parse tree in bracketed format.

# Arguments

- `io`: `IO` stream to write the tree to
- `tree`: the tree to search

# Keywords

- `multiline=true`: whether to include newlines in string representation
- `indent=2`: how many characters of whitespace to use in indentation (ignored if `multiline` is `false`)

# Returns

- bracketed `String` representation of the parse tree
"""
pprint(tree; kws...) = pprint(stdout, tree; kws...)
pprint(io::IO, tree; kws...) = print(io, brackets(tree; kws...))

function brackets(tree; depth=0, indent=2, multiline=true)
    s = "(" * label(tree)
    if !isleaf(tree)
        for child in children(tree)
            if !isleaf(child)
                if multiline
                    s *= "\n"
                    s *= repeat(" ", (depth + 1) * indent)
                else
                    if !isempty(label(tree))
                        s *= " "
                    end
                end
                kws = (depth=depth + 1, indent=indent, multiline=multiline)
                s *= brackets(child; kws...)
            else
                if !isempty(label(child))
                    s *= " " * label(child)
                end
            end
        end
    end
    return s * ")"
end

# iteration

abstract type ConstituencyTreeIterator end

Base.IteratorSize(::Type{<:ConstituencyTreeIterator}) = Base.SizeUnknown()

struct POSIterator{T} <: ConstituencyTreeIterator
    iterator::T
end
Base.iterate(pos::POSIterator, state...) = iterate(pos.iterator, state...)

"""
    POS(tree)

Iterator for part-of-speech tagged words.

```jldoctest
julia> POS(tree"(S (NP (DT the) (N cat)) (VP (V ate)))") |> collect
3-element Array{Any,1}:
 ("DT", "the")
 ("N", "cat")
 ("V", "ate")
```

# Arguments

- `tree`: the tree to search

# Returns

- `POSIterator`- a lazy iterator over (POS, token) pairs
"""
POS(tree::Tree) =
    POSIterator((label(node), node.children[1]) for node in PreOrderDFS(tree)
                if length(children(node)) == 1 && isleaf(node.children[1]))

struct WordsIterator{T} <: ConstituencyTreeIterator
    iterator::T
end
Base.iterate(w::WordsIterator, state...) = iterate(w.iterator, state...)

"""
    Words(tree)

Iterator for words in a sentence.

```jldoctest
julia> Words(tree"(S (NP (DT the) (N cat)) (VP (V ate)))") |> collect
3-element Array{Any,1}:
 "the"
 "cat"
 "ate"
```

# Arguments

- `tree`: the tree to search

# Returns

- `WordsIterator`: an iterator over tokens
"""
Words(tree::Tree) = WordsIterator(Leaves(tree))


text:
# Prepare for plotting in test environment

# Default is to run offline tests, producing PNG files in temp. dir.,
# with default plotting package (Plots w/pyplot backend).

# To test other packages/backends, set them up first and define `psplotter`;
# e.g. do
# `using Plots; gr(); psplotter = :Plots`
# before running the test suite.

# To get figures on screen instead of PNG files, do
# `displayplots=true`
# before running the test suite.

if isdefined(Main,:psplotter)
    psplotter = Main.psplotter
    # warning: assume plotting package is properly initialized
else
    psplotter = Symbol(get(ENV,"PSPLOTTER","Plots"))
    if psplotter == :Plots
        using Plots
#        pyplot() # seems most reliable for png output
    elseif psplotter == :PyPlot
        using PyPlot
    elseif psplotter == :Makie
        using GLMakie
    else
        error("invalid ENV[\"PSPLOTTER\"]")
    end
end
@info("Plotting package is $psplotter")
setpsplotter(psplotter)
headless = isdefined(Main,:displayplots) ? !Main.displayplots : true
gs = setgs(headless=headless, savefigs=headless, fig_id=1)


text:
# Estimate polynomial linear regression model with a half-T prior.

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Import the dataset.

data = CSV.read(joinpath("..", "..", "data", "Howell1.csv"), DataFrame)

# Use only adults and standardize

df = filter(row -> row[:age] >= 18, data);
df[!, :weight] = convert(Vector{Float64}, df[:, :weight]);
df[!, :weight_s] = (df[:, :weight] .- mean(df[:, :weight])) / std(df[:, :weight]);
df[!, :weight_s2] = df[:, :weight_s] .^ 2;

# Define a structure to hold the data: observables, covariates,
# and the degrees of freedom for the prior.

"""
Linear regression model ``y ∼ Xβ + ϵ``, where ``ϵ ∼ N(0, σ²)`` IID.
Flat prior for `β`, half-T for `σ`.
"""
Base.@kwdef mutable struct LinearRegressionModel{Ty <: AbstractVector, Tx <: AbstractMatrix,
Tv <: Real}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
    "Degrees of freedom for prior."
    v::Tv
end

# Write a function to return a properly dimensioned transformation.

function make_transformation(model::LinearRegressionModel)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

N = size(df, 1)
x = hcat(ones(N), hcat(df[:, :weight_s], df[:, :weight_s2]));
model = LinearRegressionModel(;y = df[:, :height], x=x, v=1.0)

# Pack parameters *as a single argument*.

function (model::LinearRegressionModel)(θ)
    @unpack y, x, v = model   # extract data
    @unpack β, σ = θ                 # extract parameters
    loglikelihood(Normal(0, σ), y .- x*β) + logpdf(TDist(v), σ)
end

# Evaluate at model function at some initial valuues

println()
model((β = [1.0, 2.0, 3.0], σ = 1.0)) |> display
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

p = as_particles(posterior)
display(p)

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
           Mean         SD       Naive SE       MCSE      ESS
    a 154.609019750 0.36158389 0.0057171433 0.0071845548 1000
   b1   5.838431778 0.27920926 0.0044146860 0.0048693502 1000
   b2  -0.009985954 0.22897191 0.0036203637 0.0047224478 1000
sigma   5.110136300 0.19096315 0.0030193925 0.0030728192 1000

Quantiles:
          2.5%        25.0%        50.0%       75.0%        97.5%
    a 153.92392500 154.3567500 154.60700000 154.8502500 155.32100000
   b1   5.27846200   5.6493250   5.83991000   6.0276275   6.39728200
   b2  -0.45954687  -0.1668285  -0.01382935   0.1423620   0.43600905
sigma   4.76114350   4.9816850   5.10326000   5.2300450   5.51500975
";

# end of m4.5d.jl


text:
# Implements some well known distillation/filtering protocols
#
# It will create a QC state of the form
# rho_hatA, hatB, F_A, F_B where hatA, and hatB contain the outputs,
# and FA, FB are the flag registers
#

export filtering, filteringMakeChoi, measureScheme, measureSchemeMakeChoi;

"""
``(rhoQC,P,F,psucc) = filtering(rhoAB,eps)`

Implements the filtering protocol on input state *rho* with filtering parameter *eps*. This filtering protocol is designed
to work well on states of the form p * EPR + (1-p) |10><10|.

Input:
- *rhoAB* 2x2 state to be filtered
- *eps* filtering parameter, smaller eps typically means higher fidelity but lower probability of success

Output:
- *rhoQC* filtered state with flags: Q=hatA, hatB contains the output state, C=F_A,F_B contains the flags (1 success, 0 failure)
- *P* projector specifying the global success condition on the flags; for this scheme it is given by |11><11|
- *F* fidelity achieved for the given input state
- *psucc* probability of success
"""

function filtering(rho::AbstractMatrix, eps::Number; makeChoi::Bool=false)

	(d1,d2) = size(rho);
	@assert d1 == 4 "Filtering takes a 2x2 input state."
	@assert eps <= 1 "Filtering parameter eps needs to be less than 1."
	@assert eps >= 0 "Filtering parameter eps must be positive."

	# Measurement operator for Alice and Bob
	MokAlice = sqrt(eps) * eVec(2,1)*eVec(2,1)' + eVec(2,2)*eVec(2,2)';
	MfailAlice = sqrt(1-eps) * eVec(2,1) * eVec(2,1)';
	MokBob = sqrt(eps) * eVec(2,2)*eVec(2,2)' + eVec(2,1)*eVec(2,1)';
	MfailBob= sqrt(1-eps) * eVec(2,2) * eVec(2,2)';

	# How success and failure are indicated on the local flag registers
	# Here the same for both Alice and Bob
	# success = |1><1|, failure = |0><0|
	fok = eVec(2,2)*eVec(2,2)';
	ffail = eVec(2,1) * eVec(2,1)';

	# Overall success only when both succeed
	P = kron(fok,fok);

	# Compute the state rhoQC
	(rhoQC, F, psucc) = measureScheme(rho, 2, MokAlice, MfailAlice, MokBob, MfailBob, fok, ffail, fok, ffail, P);

	return(rhoQC, P, F, psucc);
end

"""
``(choiA, choiB) = filteringMakeChoi(eps)`

Computes the choi states for Alice and Bob in the filtering protocol
with filtering parameter *eps*. This filtering protocol is designed
to work well on states of the form p * EPR + (1-p) |10><10|.

Input:
- *eps* filtering parameter, smaller eps typically means higher fidelity but lower probability of success

Output:
- *choi* state for Alice hatA,FA,A
- *choi* state for Bob hatB,FB,B
"""

function filteringMakeChoi(eps::Number)

	@assert eps <= 1 "Filtering parameter eps needs to be less than 1."
	@assert eps >= 0 "Filtering parameter eps must be positive."

	# Measurement operator for Alice and Bob
	MokAlice = sqrt(eps) * eVec(2,1)*eVec(2,1)' + eVec(2,2)*eVec(2,2)';
	MfailAlice = sqrt(1-eps) * eVec(2,1) * eVec(2,1)';
	MokBob = sqrt(eps) * eVec(2,2)*eVec(2,2)' + eVec(2,1)*eVec(2,1)';
	MfailBob= sqrt(1-eps) * eVec(2,2) * eVec(2,2)';

	# How success and failure are indicated on the local flag registers
	# Here the same for both Alice and Bob
	# success = |1><1|, failure = |0><0|
	fok = eVec(2,2)*eVec(2,2)';
	ffail = eVec(2,1) * eVec(2,1)';

	choiA = measureSchemeMakeChoi(MokAlice,MfailAlice,fok, ffail);
	choiB = measureSchemeMakeChoi(MokBob,MfailBob,fok, ffail);

	return(choiA, choiB);
end

""" `(rhoQC, F, psucc) = measureScheme(rhoAB, k, MokAlice, MfailAlice, MokBob, MfailBob, fokAlice, ffailAlice, fokBob, ffailBob, P)`

Implements a distillation scheme in which Alice and Bob perform the indicated measurement and then output a success or failure flag individually. Global success is determined by the projector P which may compare the to local flag registers.

Input:
- *rhoAB* input state to be distilled
- *k* local dimension of the max entangled state to be produced
- *MokAlice* Kraus operator corresponding to Alice success
- *MfailAlice* Kraus operator corresponding to Alice failure
- *MokBob* Kraus operator corresponding to Bob success
- *MfailBob* Kraus operator corresponding to Bob failure
- *fokAlice* flag Alice outputs for success
- *ffailAlice* flag Alice outputs for failure
- *fokBob* flag Bob outputs for success
- *ffailBob* flag Bob outputs for failure
- *P* projector determining for which joint flags they declare global success

Output:
- *rhoQC* output state where Q=hatA, hatB is the output state and C=F_A,F_B the flag register
- *F* fidelity achieved
- *psucc* success probability achieved

"""

function measureScheme(rho::AbstractMatrix, k::Number, MokAlice::AbstractMatrix, MfailAlice::AbstractMatrix, MokBob::AbstractMatrix, MfailBob::AbstractMatrix, fokAlice::AbstractMatrix, ffailAlice::AbstractMatrix, fokBob::AbstractMatrix, ffailBob::AbstractMatrix, P::AbstractMatrix)

	@assert isQuantumState(rho) "Input is not a quantum state."
	@assert k >= 2 "Output dimension is at least 2."

	# Desired output state
	want = maxEnt(k);

	# Construct the QC state of the output
	rhoQC = kron(kron(MokAlice,MokBob) * rho * kron(MokAlice,MokBob)', kron(fokAlice, fokBob));
	rhoQC += kron(kron(MokAlice,MfailBob) * rho * kron(MokAlice,MfailBob)', kron(fokAlice, ffailBob));
	rhoQC += kron(kron(MfailAlice,MokBob) * rho * kron(MfailAlice,MokBob)', kron(ffailAlice, fokBob));
	rhoQC += kron(kron(MfailAlice,MfailBob) * rho * kron(MfailAlice,MfailBob)', kron(ffailAlice, ffailBob));

	# Compute prob. of success and the post-measurement state conditioned on success, includes flag register
	rhoSuccUN = kron(eye(k^2),P) * rhoQC * kron(eye(k^2),P)';
	psucc = trace(rhoSuccUN);
	rhoSucc  = rhoSuccUN/psucc;

	# Tracing out the flag register
	(d1,d2) = size(P);
	rhoSuccNoFlag = partialtrace(rhoSucc,[2], [k^2, d1]);

	# Compute the fidelity with max entangled state
	F = entFidelity(rhoSuccNoFlag);

	return(rhoQC, F, psucc);
end

""" `choiState = measureSchemeMakeChoi(MokAlice, MfailAlice, fokAlice, ffailAlice)`

Computes the choi state of a measuring distillation/filtering protocol of one party: Outputs choi_QA where A is a copy of the input and Q is the output of the distillation map A->hatA, fA

Input:
- *MokAlice* Kraus operator corresponding to Alice success
- *MfailAlice* Kraus operator corresponding to Alice failure
- *fokAlice* flag Alice outputs for success
- *ffailAlice* flag Alice outputs for failure

"""

function measureSchemeMakeChoi(MokAlice::AbstractMatrix, MfailAlice::AbstractMatrix, fokAlice::AbstractMatrix, ffailAlice::AbstractMatrix)

	@assert isQuantumState(fokAlice) "Success tag is not a valid quantum state."
	@assert isQuantumState(ffailAlice) "Failure tag is not a valid quantum state."

	# Compute input dimensions
	(outA, inA) = size(MokAlice);
	(fA, fA2) = size(fokAlice);

	# We'll produce the choi state starting with the maximally entangled one.
	rho = maxEnt(inA);

	# Construct the QC state of the output
	id = eye(inA);
	choi = kron(kron(MokAlice,id) * rho * kron(MokAlice,id)', fokAlice);
	choi += kron(kron(MfailAlice,id) * rho * kron(MfailAlice,id)', ffailAlice);

	# Reorder to hatA, F_A, A
	choi = permutesystems(choi, [1, 3, 2], [outA, inA, fA]);

	return choi;
end


text:
using Unitful
using DataStructures

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t*u"s"
result["r"]    = [[0.4 * cos(t[i]),
                   0.5 * sin(t[i]),
                   0.3 * cos(t[i])] for i in eachindex(t)]*u"m"

plot(result, ["r", "r[2]", "r[2:3]"], heading="Plot vector signals")

println("\n... result info:")
printResultInfo(result)


text:
@testset "Test BC handlers" begin

    using jFEMTools
    using SparseArrays
    using LinearAlgebra

    mesh = unitSquareMesh(TriangleCell, (2,2));
    dim = 2
    element = PoissonVirtualElement(dim,1);
    u = TrialFunction(VEMFunctionSpace(mesh,element))
    dh = DofHandler(mesh, [u]);

    # ### Boundary conditions
    dbc = Dirichlet(dh,u,"boundary",x -> 0);
    @test dbc.prescribed_dofs == [1,2,3,5,6,7,8,9]
    @test dbc.values == zeros(8)
    K = sparse(I(9)); b = ones(9);
    apply!(K,b,dbc);
    @test b == [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    @test K == sparse(I(9))
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
module TestArrayMethods

using Timestamps
using Test

# clstamp = convert(TimestampArray, cl)

# facts("array methods work") do

#     context("aggregating methods") do
#         @fact sum(clstamp).value       --> roughly(23095.24, atol=.01)
#         @fact prod(clstamp[1:2]).value --> 11473.85
#     end

#     context("statistical methods") do
#         @fact mean(clstamp).value    --> roughly(46.190479999999994)
#         @fact var(clstamp).value     --> roughly(1496.90588754469)
#         @fact std(clstamp).value     --> roughly(38.68986802180501)
#         @fact maximum(clstamp).value --> 144.19
#         @fact minimum(clstamp).value --> 14.0
#     end
# end

end # module


text:
"""
    SVESMtrain(svr::LIBSVM.AbstractSVR, esn::AbstractLeakyESN [, y_target::AbstractArray{Float64}])

Train the ESN using SVM regression, as described in [1].

[1] Shi, Zhiwei, and Min Han. "Support vector echo-state machine for chaotic time-series prediction." IEEE Transactions on Neural Networks 18.2 (2007): 359-372.

"""
function SVESMtrain(svr::LIBSVM.AbstractSVR,
    esn::AbstractLeakyESN;
    y_target::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)

    if size(y_target, 1) == 1
        fitted_svr = LIBSVM.fit!(svr, states_new', vec(y_target))
    else
        fitted_svr = []
        for i=1:size(y_target, 1)
            push!(fitted_svr, LIBSVM.fit!(svr, states_new', y_target[i,:]))
        end
    end
    return fitted_svr
end

"""
    SVESM_direct_predict(esn::AbstractLeakyESN, test_in::AbstractArray{Float64}, fitted_svr::LIBSVM.AbstractSVR)

Given the input data return the corresponding predicted output, as described in [1].

[1] Shi, Zhiwei, and Min Han. "Support vector echo-state machine for chaotic time-series prediction." IEEE Transactions on Neural Networks 18.2 (2007): 359-372.

"""
function SVESM_direct_predict(esn::AbstractLeakyESN,
    test_in::AbstractArray{Float64},
    fitted_svr::LIBSVM.AbstractSVR)

    x = esn.states[:, end]
    prediction_states = zeros(Float64, size(esn.states, 1), size(test_in, 2))

    if esn.extended_states == false
        for i=1:size(test_in, 2)
            x_new = nla(esn.nla_type, x)
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x_new, test_in[:,i])
            prediction_states[:, i] = x
        end
    else
        for i=1:size(test_in, 2)
            x_new = nla(esn.nla_type, x)
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x_new[1:esn.res_size], test_in[:,i]), test_in[:, i])
            prediction_states[:,i] = x
        end
    end
    output = LIBSVM.predict(fitted_svr, prediction_states')
    return output
end

#predict if one-variable timeseries is provided
function SVESMpredict(esn::AbstractLeakyESN,
    predict_len::Int,
    fitted_svr::LIBSVM.AbstractSVR)

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]
    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = LIBSVM.predict(fitted_svr, x_new)
            output[:, i] = out
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = LIBSVM.predict(fitted_svr, x_new)
            output[:, i] = out
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output
end

#predict for multidimensional timeseries
"""
    SVESMpredict(esn::AbstractLeakyESN, predict_len::Int, fitted_svr::AbstractArray{Any})

Return the prediction for a given length of the constructed ESN struct using SVMs.
"""
function SVESMpredict(esn::AbstractLeakyESN,
    predict_len::Int,
    fitted_svr::AbstractArray{Any})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = []
            for i=1:size(fitted_svr, 1)
                push!(out, LIBSVM.predict(fitted_svr[i], x_new)[1])
            end
            output[:, i] = out
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = Array{Float64}(undef, esn.out_size)
            for i=1:size(fitted_svr, 1)
                out[i] = LIBSVM.predict(fitted_svr[i], x_new)[1]
            end
            output[:, i] = out
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output
end

"""
    SVESMpredict_h_steps(esn::AbstractLeakyESN, predict_len::Int, h_steps::Int, test_data::AbstractArray{Float64}, fitted_svr::LIBSVM.AbstractSVR)

Return the prediction for h steps ahead of the constructed ESN struct using SVMs.
"""
function SVESMpredict_h_steps(esn::AbstractLeakyESN,
    predict_len::Int,
    h_steps::Int,
    test_data::AbstractArray{Float64},
    fitted_svr::LIBSVM.AbstractSVR)

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]
    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = LIBSVM.predict(fitted_svr, x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, test_data[:,i])
            else
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
            end
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)
            out = LIBSVM.predict(fitted_svr, x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = vcat(ReservoirComputing.leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], test_data[:,i]), test_data[:,i])
            else
                x = vcat(ReservoirComputing.leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
            end
        end
    end
    return output
end


text:
"""
    dominate(p, q)

Returns `1` if `p` is dominated by `q`, `-1` if otherwise, and `0` if dominance cannot be determined.
"""
function dominate(p::T, q::T) where {T <: AbstractArray}
    ret = 0
    for (i,j) in zip(p,q)
        if i < j
            ret == -1 && return 0
            ret = 1
        elseif j < i
            ret == 1 && return 0
            ret = -1
        end
    end
    return ret
end

"""
dominations(P::AbstractVector)

Returns a domination matrix of all elements in the input collection `P`.
"""
function dominations(P::AbstractVector{T}) where {T <: AbstractArray}
    l = length(P)
    D = zeros(Int8, l, l)
    for i in 1:l
        for j in (i+1):l
            D[i,j] = dominate(P[i],P[j])
            D[j,i] = -D[i,j]
        end
    end
    D
end

"""
    nondominatedsort!(R, F)

Calculate fronts for fitness values `F`, and store ranks of the individuals into `R`.
"""
function nondominatedsort!(R, P)
    n = size(P,2)
    @assert length(R) == n "Ranks must be defined for the whole population"

    Sₚ = Dict(i=>Set() for i in 1:n)
    C = zeros(Int, n)

    # construct first front
    F =[Int[]]
    for i in 1:n
        for j in i+1:n
            r = dominate(view(P,:,i), view(P,:,j)) #M[i,j]
            if r == 1
                push!(Sₚ[i], j)
                C[j] += 1
            elseif r == -1
                push!(Sₚ[j], i)
                C[i] += 1
            end
        end
        if C[i] == 0
            R[i] = 1
            push!(F[1], i)
        end
    end

    # construct rest of the fronts
    while !isempty(last(F))
        Q = Int[]
        for i in last(F)
            for j in Sₚ[i]
                C[j] -= 1
                if C[j] == 0
                    push!(Q, j)
                    R[j] = length(F) + 1
                end
            end
        end
        push!(F, Q)
    end
    isempty(last(F)) && pop!(F)

    F #, R #, Sₚ
end

"""
    crowding_distance!((C, F, fronts)

Calculate crowding distance for individuals and save the results into `C`
given the fitness values `F` and collection of `fronts`.
"""
function crowding_distance!(C::AbstractVector, F::AbstractMatrix{T}, fronts) where {T}
    for f in fronts
        cf = @view C[f]
        if length(cf) <= 2
            cf .= typemax(T)
        else
            # sort front by each objective value
            SF = F[:, f]
            d = size(SF,1)
            IX = zeros(Int, size(SF))
            IIX = zeros(Int, size(SF))
            for i in 1:d
                irow, iirow, row = view(IX,i,:), view(IIX,i,:), view(SF,i,:)
                sortperm!(irow, row)
                sortperm!(iirow, irow)
                permute!(row, irow)
            end
            nrm = SF[:,end] - SF[:,1]
            dst = (hcat(SF, fill(typemax(T), d)) - hcat(fill(typemin(T), d), SF)) ./ nrm
            dst[isnan.(dst)] .= zero(T)
            ss = sum(mapslices(v->diag(dst[:,v]) + diag(dst[:,v.+1]), IIX, dims=1), dims=1)
            cf .= vec(ss)/d
        end
    end
    C
end


text:
# Test for tree-based analysis tools found in `AnalysisTools.jl`.
using Test
using IncrementalInference

@testset "Number of non-zero calculation for frontals." begin
    # Alternative way of calculating number of upper triangular matrix elements.
    nnzFrontalsRecursive(dim) = dim==1 ? 1 : dim + nnzFrontalsRecursive(dim-1)
    # Both must agree for any integer number dimension.
    for dim in 1:100
        @test nnzFrontalsRecursive(dim) == nnzFrontals(dim)
    end
end

@testset "Number of non-zero calculation for full cliques." begin
    fg = generateGraph_Kaess()
    vo = [:l1, :l2, :x1, :x2, :x3]
    tree = buildTreeReset!(fg, vo)
    # Must agree with hand-calculated values, iSAM2 paper.
    @test nnzClique(tree.cliques[1]) == 3
    @test nnzClique(tree.cliques[2]) == 5
    @test nnzClique(tree.cliques[3]) == 2
end

@testset "Number of non-zero calculation for full trees." begin
    fg = generateGraph_Kaess()
    vo = [:l1, :l2, :x1, :x2, :x3]
    tree = buildTreeReset!(fg, vo)
    # Must agree with hand-calculated values, iSAM2 paper.
    @test nnzTree(tree) == 10
end

@testset "Test drawTree" begin
    fg = generateGraph_Kaess(graphinit=false)
    vo = [:l1, :l2, :x1, :x2, :x3]
    tree = buildTreeReset!(fg, vo)

    drawTree(tree, show=false)
end

##


text:
"""
    ahp_RI(n)
Return AHP Random Index for a given n.
# Arguments:
 - `n::Int64`: Index value

 # Description
 The function returns the random index for a given integer between 3 and 15.

# Output
- `::Float64`: Random Index (RI) value.

# Examples
```julia-repl
julia> ahp_RI(3)
0.58
```

# References

Saaty, Thomas L. "Decision making with the analytic hierarchy process." International journal of services sciences 1.1 (2008): 83-98.
```
"""
function ahp_RI(n::Int64)::Float64
    # First index is n = 3
    # RI[3] = 0.58
    # RI[4] = 0.90
    ris = [0.58, 0.90, 1.12, 1.24, 1.32,
          1.41, 1.45, 1.49, 1.51, 1.53,
          1.56, 1.57, 1.59]
    if n < 3
        return 0
    elseif n <= 15
        return ris[n - 2]
    else
        return ris |> last
    end

end

"""
    ahp_consistency(comparisonMatrix)
    Test if a comparison matrix whether consistent or not.

 # Arguments:
 - `comparisonMatrix::DataFrame`: Comparison matrix for AHP.

 # Description
 AHP is based on subjective comparison between criteria. The success of the method highly depends on
 consistency of these comparison.

# Output
- `::AHPConsistencyResult`: AhpConsistencyResult object that holds multiple outputs including the boolean that indicates the matrix is whether consistent or not.

# Examples
```julia-repl
julia> K = [
               1 7 (1 / 5) (1 / 8) (1 / 2) (1 / 3) (1 / 5) 1;
               (1 / 7) 1 (1 / 8) (1 / 9) (1 / 4) (1 / 5) (1 / 9) (1 / 8);
               5 8 1 (1 / 3) 4 2 1 1;
               8 9 3 1 7 5 3 3;
               2 4 (1 / 4) (1 / 7) 1 (1 / 2) (1 / 5) (1 / 5);
               3 5 (1 / 2) (1 / 5) 2 1 (1 / 3) (1 / 3);
               5 9 1 (1 / 3) 5 3 1 1;
               1 8 1 (1 / 3) 5 3 1 1
           ];

julia> K
8×8 Array{Float64,2}:
 1.0       7.0  0.2    0.125     0.5   0.333333  0.2       1.0
 0.142857  1.0  0.125  0.111111  0.25  0.2       0.111111  0.125
 5.0       8.0  1.0    0.333333  4.0   2.0       1.0       1.0
 8.0       9.0  3.0    1.0       7.0   5.0       3.0       3.0
 2.0       4.0  0.25   0.142857  1.0   0.5       0.2       0.2
 3.0       5.0  0.5    0.2       2.0   1.0       0.333333  0.333333
 5.0       9.0  1.0    0.333333  5.0   3.0       1.0       1.0
 1.0       8.0  1.0    0.333333  5.0   3.0       1.0       1.0

julia> dmat = makeDecisionMatrix(K);
julia> ahp_consistency(dmat).isConsistent
true
```

# References
Saaty, Thomas L. "Decision making with the analytic hierarchy process." International journal of services sciences 1.1 (2008): 83-98.
"""
function ahp_consistency(comparisonMatrix::DataFrame)::AHPConsistencyResult

    n, m = size(comparisonMatrix)

    csums = colsums(comparisonMatrix)
    normalizedComparisonMatrix = zeros(Float64, n, n)

    for i in 1:n
        for j in 1:n
            normalizedComparisonMatrix[i, j] = comparisonMatrix[i, j] / csums[j]
        end
    end

    priority_vector = rowmeans(normalizedComparisonMatrix)
    #mcomparisonMatrix = convert(Array{Float64,2}, comparisonMatrix)
    mcomparisonMatrix = Matrix{Float64}(comparisonMatrix)

    consistency_vector = mcomparisonMatrix * priority_vector

    pc_matrix = consistency_vector ./ priority_vector

    lambda_max = sum(pc_matrix) / n

    CI  = (lambda_max - n) / (n - 1)
    ri  = ahp_RI(n)
    CR  = CI / ri

    isConsistent = (CR < 0.1)

    result = AHPConsistencyResult(
        comparisonMatrix,
        makeDecisionMatrix(normalizedComparisonMatrix),
        consistency_vector,
        priority_vector,
        pc_matrix,
        lambda_max,
        CI,
        ri,
        CR,
        isConsistent
    )

    return result

end

"""
    AHP(comparisonMatrixList, criteriaComparisonMatrix)

Apply AHP (Analytical Hierarchy Process) for a given list of comparison matrices and criteria comparison matrix.

# Arguments:
 - `comparisonMatrixList::Array{DataFrame,1}`: Array of comparison matrices for all of the criteria.
 - `criteriaComparisonMatrix::DataFrame`: Criteria comparison matrix for AHP (Comparison of columns).

 # Description
 AHP is based on subjective comparison between criteria. The success of the method highly depends on
 consistency of these comparison. The method test the consistency first. At the next step, weights are
 calculated. The ordering of rows is determined by these weights.

# Output
- `::AHPResult`: AhpResult object that holds multiple outputs including calculated weights and scores.

# Examples
```julia-repl
julia> K
8×8 Array{Float64,2}:
 1.0       7.0  0.2    0.125     0.5   0.333333  0.2       1.0
 0.142857  1.0  0.125  0.111111  0.25  0.2       0.111111  0.125
 5.0       8.0  1.0    0.333333  4.0   2.0       1.0       1.0
 8.0       9.0  3.0    1.0       7.0   5.0       3.0       3.0
 2.0       4.0  0.25   0.142857  1.0   0.5       0.2       0.2
 3.0       5.0  0.5    0.2       2.0   1.0       0.333333  0.333333
 5.0       9.0  1.0    0.333333  5.0   3.0       1.0       1.0
 1.0       8.0  1.0    0.333333  5.0   3.0       1.0       1.0

julia> A1
4×4 Array{Float64,2}:
 1.0       3.0  0.2       2.0
 0.333333  1.0  0.142857  0.333333
 5.0       7.0  1.0       4.0
 0.5       3.0  0.25      1.0

julia> A2
4×4 Array{Float64,2}:
 1.0   0.5       4.0       5.0
 2.0   1.0       6.0       7.0
 0.25  0.166667  1.0       3.0
 0.2   0.142857  0.333333  1.0

julia> A3
4×4 Array{Float64,2}:
 1.0       0.5  0.166667  3.0
 2.0       1.0  0.25      5.0
 6.0       4.0  1.0       9.0
 0.333333  0.2  0.111111  1.0

julia> A4
4×4 Array{Float64,2}:
 1.0       7.0  0.25      2.0
 0.142857  1.0  0.111111  0.2
 4.0       9.0  1.0       5.0
 0.5       5.0  0.2       1.0

julia> A5
4×4 Array{Float64,2}:
 1.0       6.0  2.0   3.0
 0.166667  1.0  0.25  0.333333
 0.5       4.0  1.0   2.0
 0.333333  3.0  0.5   1.0

julia> A6
4×4 Array{Float64,2}:
 1.0  0.25  0.5  0.142857
 4.0  1.0   2.0  0.333333
 2.0  0.5   1.0  0.2
 7.0  3.0   5.0  1.0

julia> A7
4×4 Array{Float64,2}:
 1.0       3.0   7.0  1.0
 0.333333  1.0   4.0  0.333333
 0.142857  0.25  1.0  0.142857
 1.0       3.0   7.0  1.0

julia> A8
4×4 Array{Float64,2}:
 1.0    2.0       5.0       8.0
 0.5    1.0       3.0       6.0
 0.2    0.333333  1.0       3.0
 0.125  0.166667  0.333333  1.0

julia> km = makeDecisionMatrix(K);
julia> as = map(makeDecisionMatrix, [A1, A2, A3, A4, A5, A6, A7, A8]);
julia> result = ahp(as, km);
julia> result.bestIndex
3
julia> result.scores
4-element Array{Float64,1}:
 0.2801050163111839
 0.14822726478768022
 0.3813036392434616
 0.19036407965767424
```

# References
Saaty, Thomas L. "Decision making with the analytic hierarchy process." International journal of services sciences 1.1 (2008): 83-98.
"""
function ahp(comparisonMatrixList::Array{DataFrame,1}, criteriaComparisonMatrix::DataFrame)::AHPResult

    result_list = map(ahp_consistency, comparisonMatrixList)

    n = length(result_list)

    ncriteria, _ = size(criteriaComparisonMatrix)

    ncandidates, _ = size(comparisonMatrixList[1])

    decision_matrix = zeros(Float64, ncandidates, ncriteria)

    @inbounds for i in 1:n
        decision_matrix[:,i] = result_list[i].priority
    end

    criteria_consistency = ahp_consistency(criteriaComparisonMatrix)

    weights = criteria_consistency.priority

    decision_matrix_df = makeDecisionMatrix(decision_matrix)
    ordering_result =  decision_matrix * weights

    bestIndex = sortperm(ordering_result) |> last

    result = AHPResult(
        comparisonMatrixList,
        criteriaComparisonMatrix,
        criteria_consistency,
        decision_matrix_df,
        ordering_result,
        weights,
        bestIndex
    )

    return result
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
export LayerMaterial, Conductor, Dielectric, Stackup
export copper, fr408

"""
    LayerMaterial

Subtypes of `LayerMaterial` hold the material properties for the layers of the
PCB.
"""
abstract type LayerMaterial end
name(::LayerMaterial) = ""
ρ_20(::LayerMaterial) = NaN
tc(::LayerMaterial) = NaN
ϵ(::LayerMaterial) = NaN

"""
    Conductor(name, ρ, tc)

**Fields**
- `name`  -- name of the material
- `ρ_20`  -- conductivity at 20C
- `tc`    -- temperature coefficient

see constant `copper` as example.
"""
struct Conductor <: LayerMaterial
  name  :: String
  ρ_20  :: Float64
  tc    :: Float64
end
conductivity(c::Conductor, temperature) = conductivity(c.ρ_20,c.tc,temperature)
const copper = Conductor("copper",1.68e-8, 0.003862) # (Ωm, 1/K)

"""
    Dielectric(name, ϵ)

**Fields**
- `name`    -- name of the material
- `ϵ`       -- permittivity
"""
struct Dielectric <: LayerMaterial
  name :: String
  ϵ    :: Float64
end
const ϵ0 = 8.854187817e-12
const fr408 = Dielectric("fr408",3.67*ϵ0)

"""
    Stackup(material,thickness)

Specify the PCB stackup.

example:
for double sided FR408 PCB with 1oz copper 63mil thick.
`Stackup([copper,fr408,copper],[0.48e-3,1.6e-3 - (2*0.48e-3),0.48e-3])`

copper and fr408 are exported constants.
"""
struct Stackup
  material  :: Tuple
  thickness :: Tuple
  function Stackup(s::Tuple,t::Tuple)
    if length(s) != length(t)
      throw(ArgumentError("length of material and thickness must match"))
    end
    if length(s) < 3
      throw(ArgumentError("stackup must have at least three layers"))
    end
    if mod(length(s),2) == 0
      throw(ArgumentError("stackup must have odd number of layers"))
    end
    for i in eachindex(s)
      if mod(i,2) == 0
        if typeof(s[i]) != Dielectric
          throw(ArgumentError("even layers must be Dielectric"))
        end
      else
        if typeof(s[i]) != Conductor
          throw(ArgumentError("odd layers must be Conductor"))
        end
      end
    end
    new(s,t)
  end
end
function Stackup(s::Vector{T}, t::Vector) where T<:LayerMaterial
  Stackup(ntuple(x->s[x],length(s)), ntuple(x->t[x],length(t)))
end


text:
# First Half

function isvalid1(line)
    m = match(r"^(\d+)-(\d+) ([a-z]{1}): (.*)$", line)
    lwst, high = parse(Int, m[1]), parse(Int, m[2])
    letter = m[3]
    pw = m[4]
    eachmatch(Regex(letter), pw) |> collect |> length ∈ lwst:high
end

puzzleinput = readlines("inputs/day02.txt")
count(isvalid1, puzzleinput)

# Second Half

function isvalid2(line)
    m = match(r"^(\d+)-(\d+) ([a-z]{1}): (.*)$", line)
    p1, p2 = parse(Int, m[1]), parse(Int, m[2])
    letter = m[3][1]
    pw = m[4]
    (pw[p1] == letter) ⊻ (pw[p2] == letter)
end

puzzleinput = readlines("inputs/day02.txt")
count(isvalid2, puzzleinput)


text:
"""
    NexUpload(nexSerial, name; pid=pid, cid=cid)

A Nextion control to download the tft file.
"""
struct NexUpload <: AbstractNexObject
    _nid::NexID

    function NexUpload(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid)
    end
end