text:
using KrylovKit: OrthonormalBasis, householder, rows, cols, hschur!, schur2eigvals, schur2eigvecs, permuteschur!

@testset "Orthonormalize with algorithm $alg" for alg in (cgs, mgs, cgs2, mgs2, cgsr, mgsr)
    @testset for S in (Float32, Float64, ComplexF32, ComplexF64)
        b = OrthonormalBasis{Vector{S}}()
        A = randn(S,(n,n))
        v, r, x = orthonormalize(A[:,1], b, alg)
        @test r ≈ norm(A[:,1])
        @test norm(v) ≈ 1
        @test length(x) == 0
        push!(b, v)
        for i = 2:n
            v, r, x = orthonormalize(A[:,i], b, alg)
            @test norm([r,norm(x)]) ≈ norm(A[:,i])
            @test norm(v) ≈ 1
            @test length(x) == i-1
            push!(b, v)
        end
        U = hcat(b...)
        @test U'*U ≈ I
        v = randn(S, n)
        @test U*v ≈ b*v
    end
end

@testset "Givens and Householder" begin
    @testset for S in (Float32, Float64, ComplexF32, ComplexF64)
        U, = svd!(randn(S,(n,n)))
        b = OrthonormalBasis(map(copy,cols(U)))
        v = randn(S,(n,))
        g, r = givens(v, 1, n)
        @test rmul!(U,g) ≈ hcat(rmul!(b,g)...)
        h, r = householder(v, axes(v,1), 3)
        @test r ≈ norm(v)
        v2 = lmul!(h, copy(v))
        v3 = zero(v2); v3[3] = r;
        @test v2 ≈ v3
        @test lmul!(h, one(U)) ≈ rmul!(one(U), h)
        @test lmul!(h, one(U))' ≈ lmul!(h', one(U))
        @test rmul!(U, h) ≈ hcat(rmul!(b,h)...)
    end
end

@testset "Rows and cols iterator" begin
    @testset for S in (Float32, Float64, ComplexF32, ComplexF64)
        A = randn(S, (n,n))
        rowiter = rows(A)
        @test typeof(first(rowiter)) == eltype(rowiter)
        @test all(t->t[1]==t[2], zip(rowiter, [A[i,:] for i=1:n]))
        coliter = cols(A)
        @test typeof(first(coliter)) == eltype(coliter)
        @test all(t->t[1]==t[2], zip(coliter, [A[:,i] for i=1:n]))
    end
end

@testset "Dense Schur factorisation and associated methods" begin
    @testset for S in (Float32, Float64, ComplexF32, ComplexF64)
        H = hessenberg(rand(S,n,n)).H

        # schur factorisation of Hessenberg matrix
        T,U,w = hschur!(copy(H))
        @test H*U ≈ U*T
        @test schur2eigvals(T) ≈ w

        # full eigenvector computation
        V = schur2eigvecs(T)
        @test T*V ≈ V*Diagonal(w)

        # selected eigenvector computation
        p = randperm(n)
        select = p[1:(n>>1)]
        V2 = schur2eigvecs(T, select)
        @test T*V2 ≈ V2*Diagonal(w[select])

        # permuting / reordering schur: take permutations that keep 2x2 blocks together in real case
        p = sortperm(w, by=real)
        T2, U2 = permuteschur!(copy(T), copy(U), p)
        @test H*U2 ≈ U2*T2
        @test schur2eigvals(T2) ≈ w[p]

        p = sortperm(w, by=abs)
        T2, U2 = permuteschur!(copy(T), copy(U), p)
        @test H*U2 ≈ U2*T2
        @test schur2eigvals(T2) ≈ w[p]
    end
end


text:
export TextConfig, save, load, normalize_text, tokenize
using Unicode

#, language!
# using Languages
# using SnowballStemmer

# const _PUNCTUACTION = """;:,.@#&\\-\"'/:*"""
const _PUNCTUACTION = """;:,.&\\-\"'/:*“”«»"""
const _SYMBOLS = "()[]¿?¡!{}~<>|"
const PUNCTUACTION  = _PUNCTUACTION * _SYMBOLS
# A symbol s in this list will be expanded to BLANK*s if the predecesor of s is neither s nor BLANK
# On changes from s to BLANK or [^s] it will produce also produce an extra BLANK
# Note that enabled del_punc will delete all these symbols without any of the previous expansions

const BLANK_LIST = string(' ', '\t', '\n', '\v', '\r')
const RE_USER = r"""@[^;:,.@#&\\\-\"'/:\*\(\)\[\]\¿\?\¡\!\{\}~\<\>\|\s]+"""
const RE_URL = r"(http|ftp|https)://\S+"
const RE_NUM = r"\d+"
const BLANK = ' '
const PUNCTUACTION_BLANK = string(PUNCTUACTION, BLANK)

# SKIP_WORDS = set(["…", "..", "...", "...."])

mutable struct TextConfig
    del_diac::Bool
    del_dup::Bool
    del_punc::Bool
    group_num::Bool
    group_url::Bool
    group_usr::Bool
    lc::Bool
    qlist::Vector{Int}
    nlist::Vector{Int}
    slist::Vector{Tuple{Int,Int}}

    """
    Initializes a `TextConfig` structure
    """
    function TextConfig(;
        del_diac=true,
        del_dup=false,
        del_punc=false,
        group_num=true,
        group_url=true,
        group_usr=false,
        lc=true,
        qlist=Int[],
        nlist=Int[1],
        slist=Tuple{Int,Int}[]
    )
        new(del_diac, del_dup, del_punc, group_num, group_url, group_usr, lc,
            qlist, nlist, slist)
    end
end

"""
    normalize_text(config::TextConfig, text::String)::Vector{Char}

Normalizes a given text using the specified transformations of `config`

"""
function normalize_text(config::TextConfig, text::String)::Vector{Char}
    if config.lc
        text = lowercase(text)
    end

    if config.group_url
        text = replace(text, RE_URL => "_url")
    end

    if config.group_usr
        text = replace(text, RE_USER => "_usr")
    end

    if config.group_num
        text = replace(text, RE_NUM => "_num")
    end

    L = Char[BLANK]
    prev = BLANK

    @inbounds for u in Unicode.normalize(text, :NFD)
        if config.del_diac
            o = Int(u)
            0x300 <= o && o <= 0x036F && continue
        end

        if u in BLANK_LIST
            u = BLANK
        elseif config.del_dup && prev == u
            continue
        elseif config.del_punc && u in PUNCTUACTION
            prev = u
            continue
        end

        prev = u
        push!(L, u)
    end

    push!(L, BLANK)

    L
end

"""
    push_word!(config::TextConfig, output::Vector{String}, token::String, normalize::Function)

Pushes a word into token list after applying the `normalize` function; it discards empty strings.
"""
function push_word!(config::TextConfig, output::Vector{Symbol}, token::Vector{UInt8}, normalize::Function)
    if length(token) > 0
        t = normalize(String(token))::String

        if length(t) > 0
            push!(output, Symbol(t))
        end
    end
end

"""
    tokenize_words(config::TextConfig, text::Vector{Char}, normalize::Function)

Performs the word tokenization
"""
function tokenize_words(config::TextConfig, text::Vector{Char}, normalize::Function, buff::IOBuffer)
    n = length(text)
    L = Symbol[]
    @inbounds for i in 1:n
        c = text[i]

        if c == BLANK
            push_word!(config, L, take!(buff), normalize)
        elseif i > 1
            if text[i-1] in PUNCTUACTION && !(c in PUNCTUACTION)
                push_word!(config, L, take!(buff), normalize)
                write(buff, c)
                continue
            elseif !(text[i-1] in PUNCTUACTION_BLANK) && c in PUNCTUACTION
                push_word!(config, L, take!(buff), normalize)
                write(buff, c)
                continue
            else
                write(buff, c)
            end
        else
            write(buff, c)
        end
    end

    push_word!(config, L, take!(buff), normalize)

    L
end

"""
    tokenize(config::TextConfig, arr::AbstractVector)::Vector{Symbol}

Tokenizes an array of strings
"""
function tokenize(config::TextConfig, arr::AbstractVector, normalize::Function=identity)::Vector{Symbol}
    L = Symbol[]
    n = length(arr)
    sizehint!(L, (length(config.nlist) + length(config.slist)) * (div(n, 2) + 1) + length(config.qlist) * n)
    for text in arr
        t = normalize_text(config, text)
        tokenize_(config, t, L, normalize)
    end

    L
end

"""
    tokenize(config::TextConfig, text::String)::Vector{Symbol}

Tokenizes a string
"""
function tokenize(config::TextConfig, text::String, normalize::Function=identity)::Vector{Symbol}
    t = normalize_text(config, text)
    n = length(text)
    L = Symbol[]
    sizehint!(L, (length(config.nlist) + length(config.slist)) * (div(n, 2) + 1) + length(config.qlist) * n)
    tokenize_(config, t, L, normalize)
end

"""
    tokenize_(config::TextConfig, text::Vector{Char}, L::Vector{Symbol})::Vector{Symbol}

Tokenizes a vector of characters (internal method)
"""
function tokenize_(config::TextConfig, text::Vector{Char}, L::Vector{Symbol}, normalize::Function)::Vector{Symbol}
    n = length(text)
    buff = IOBuffer(Vector{UInt8}(undef, 64), write=true)
    @inbounds for q in config.qlist
        for i in 1:(n - q + 1)
            last = i + q - 1
            for j in i:last
                # for w in @view text[i:]
                write(buff, text[j])
            end

            push!(L, Symbol(take!(buff)))
        end
    end

    if length(config.nlist) > 0 || length(config.slist) > 0
        ltext = tokenize_words(config, text, normalize, buff)
        n = length(ltext)

        @inbounds for q in config.nlist
            for i in 1:(n - q + 1)
                last = i + q - 1
                for j in i:last-1
                    # for w in @view ltext[i:i+q-1]
                    write(buff, ltext[j])
                    write(buff, BLANK)
                end
                write(buff, ltext[last])
                # w = join(v, BLANK)
                # push!(L, Symbol(w))
                push!(L, Symbol(take!(buff)))
            end
        end

        @inbounds for (qsize, skip) in config.slist
            for start in 1:(n - (qsize + (qsize - 1) * skip) + 1)
                if qsize == 2
                    t = Symbol(ltext[start], BLANK, ltext[start + 1 + skip])
                else
                    t = Symbol(join([ltext[start + i * (1+skip)] for i in 0:(qsize-1)], BLANK))
                end

                push!(L, t)
            end
        end
    end

    L
end

ext: jl


text:
# ---
# title: 140. Word Break II
# id: problem140
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/word-break-ii/description/>
# hidden: true
# ---
#
# Given a **non-empty** string _s_ and a dictionary _wordDict_ containing a list
# of **non-empty** words, add spaces in _s_ to construct a sentence where each
# word is a valid dictionary word. Return all such possible sentences.
#
# **Note:**
#
#   * The same word in the dictionary may be reused multiple times in the segmentation.
#   * You may assume the dictionary does not contain duplicate words.
#
# **Example 1:**
#
#
#
#     Input: s = "catsanddog"
#     wordDict = ["cat", "cats", "and", "sand", "dog"]
#     Output:[
#       "cats and dog",
#       "cat sand dog"
#     ]
#
#
# **Example 2:**
#
#
#
#     Input: s = "pineapplepenapple"
#     wordDict = ["apple", "pen", "applepen", "pine", "pineapple"]
#     Output: [
#       "pine apple pen apple",
#       "pineapple pen apple",
#       "pine applepen apple"
#     ]
#     Explanation: Note that you are allowed to reuse a dictionary word.
#
#
# **Example 3:**
#
#
#
#     Input: s = "catsandog"
#     wordDict = ["cats", "dog", "sand", "and", "cat"]
#     Output: []
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
############################################################################################
include("params.jl")
include("construct.jl")
include("chainrules.jl")

include("types/types.jl")
include("nested/nested.jl")
include("constraints/constraints.jl")

############################################################################################
# Export


text:
# Helper types

struct PVTKArgs
    part::Int
    nparts::Int
    ismain::Bool
    ghost_level::Int
end

struct PVTKFile <: VTKFile
    pvtkargs::PVTKArgs
    xdoc::XMLDocument
    vtk::DatasetFile
    path::String
end

"""
    pvtk_grid(args...; part, nparts, ismain = (part == 1), ghost_level = 0, kwargs...)

Returns a handler representing a parallel vtk file, which can be
eventually written to file with `vtk_save`.

Positional and keyword arguments in `args` and `kwargs`
are passed to `vtk_grid` verbatim (except file names that are augmented with the
corresponding part id).

The extra keyword arguments only apply to parallel VTK file formats.

Mandatory ones are:

- `part` current (1-based) part id,
- `nparts` total number of parts,

Optional ones are:

- `ismain` true if the current part id `part` is the main (the only one that will write the `.pvtk` file),
- `ghost_level` ghost level.
"""
function pvtk_grid(
        filename::AbstractString, args...;
        part, nparts, ismain = (part == 1), ghost_level = 0, kwargs...,
    )
    bname = basename(filename)
    mkpath(filename)
    prefix = joinpath(filename, bname)
    fn = _serial_filename(part, nparts, prefix, "")
    pvtkargs = PVTKArgs(part, nparts, ismain, ghost_level)
    xdoc  = XMLDocument()
    vtk = vtk_grid(fn, args...; kwargs...)
    _, ext = splitext(vtk.path)
    path = filename * ".p" * ext[2:end]
    pvtk = PVTKFile(pvtkargs, xdoc, vtk, path)
    _init_pvtk!(pvtk)
    pvtk
end

# Add point and cell data as usual
function Base.setindex!(
        pvtk::PVTKFile, data, name::AbstractString, args...,
    )
    pvtk.vtk[name, args...] = data
end

# Save as usual
function vtk_save(pvtk::PVTKFile)
    outfiles = String[]
    if isopen(pvtk)
        if pvtk.pvtkargs.ismain
            _update_pvtk!(pvtk)
            save_file(pvtk.xdoc, pvtk.path)
            push!(outfiles, pvtk.path)
        end
        append!(outfiles, vtk_save(pvtk.vtk))
        close(pvtk)
    end
    outfiles
end

# Helper functions
function _serial_filename(part, nparts, prefix, extension)
    p = lpad(part, ndigits(nparts), '0')
    prefix * "_$p" * extension
end

function _init_pvtk!(pvtk::PVTKFile)
    # Recover some data
    vtk = pvtk.vtk
    pvtkargs = pvtk.pvtkargs
    pgrid_type = "P" * vtk.grid_type
    npieces = pvtkargs.nparts
    pref, _ = splitext(pvtk.path)
    _, ext = splitext(vtk.path)
    prefix = joinpath(pref, basename(pref))

    # VTKFile (root) node
    pvtk_root = create_root(pvtk.xdoc, "VTKFile")
    set_attribute(pvtk_root, "type", pgrid_type)
    set_attribute(pvtk_root, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(pvtk_root, "byte_order", "LittleEndian")
    else
        set_attribute(pvtk_root, "byte_order", "BigEndian")
    end

    # Grid node
    pvtk_grid = new_child(pvtk_root, pgrid_type)
    set_attribute(pvtk_grid, "GhostLevel", string(pvtkargs.ghost_level))

    # Pieces (i.e. Pointers to serial files)
    for piece in 1:npieces
        pvtk_piece = new_child(pvtk_grid, "Piece")
        fn = _serial_filename(piece, npieces, prefix, ext)
        set_attribute(pvtk_piece, "Source", fn)
    end

    # Recover point type and number of components
    vtk_root = root(vtk.xdoc)
    vtk_grid = find_element(vtk_root, vtk.grid_type)
    vtk_piece = find_element(vtk_grid, "Piece")
    vtk_points = find_element(vtk_piece, "Points")
    vtk_data_array = find_element(vtk_points, "DataArray")
    point_type = attribute(vtk_data_array, "type")
    Nc = attribute(vtk_data_array, "NumberOfComponents")

    ## PPoints node
    pvtk_ppoints = new_child(pvtk_grid, "PPoints")
    pvtk_pdata_array = new_child(pvtk_ppoints, "PDataArray")
    set_attribute(pvtk_pdata_array, "type", point_type)
    set_attribute(pvtk_pdata_array, "Name", "Points")
    set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)

    pvtk
end

function _update_pvtk!(pvtk::PVTKFile)
    vtk = pvtk.vtk
    vtk_root = root(vtk.xdoc)
    vtk_grid = find_element(vtk_root, vtk.grid_type)
    vtk_piece = find_element(vtk_grid, "Piece")

    pgrid_type = "P" * vtk.grid_type
    pvtk_root = root(pvtk.xdoc)
    pvtk_grid = find_element(pvtk_root, pgrid_type)

    # Generate PPointData
    vtk_point_data = find_element(vtk_piece, "PointData")
    if vtk_point_data !== nothing
        pvtk_ppoint_data = new_child(pvtk_grid, "PPointData")
        for vtk_data_array in child_elements(vtk_point_data)
            t = attribute(vtk_data_array, "type")
            name = attribute(vtk_data_array, "Name")
            Nc = attribute(vtk_data_array, "NumberOfComponents")
            pvtk_pdata_array = new_child(pvtk_ppoint_data, "PDataArray")
            set_attribute(pvtk_pdata_array, "type", t)
            set_attribute(pvtk_pdata_array, "Name", name)
            set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
        end
    end

    # Generate PCellData
    vtk_cell_data = find_element(vtk_piece, "CellData")
    if vtk_cell_data !== nothing
        pvtk_pcell_data = new_child(pvtk_grid, "PCellData")
        for vtk_data_array in child_elements(vtk_cell_data)
            t = attribute(vtk_data_array, "type")
            name = attribute(vtk_data_array, "Name")
            Nc = attribute(vtk_data_array, "NumberOfComponents")
            pvtk_pdata_array = new_child(pvtk_pcell_data, "PDataArray")
            set_attribute(pvtk_pdata_array, "type", t)
            set_attribute(pvtk_pdata_array, "Name", name)
            set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
        end
    end

    pvtk
end


text:
using ChangesOfVariables, InverseFunctions, ArraysOfArrays, Statistics
using Optimisers
using Plots

using EuclidianNormalizingFlows: CenterStretch, JohnsonTrafo, HouseholderTrafo
using EuclidianNormalizingFlows: optimize_whitening

f_fwd_true =
    CenterStretch([4.0], [1.0], [0.0]) ∘
    JohnsonTrafo([10.0], [3.5], [10.0], [1.0])

XW = randn(1, 10^5)
stephist(XW[1,:], nbins = 100; normalize = true)

X = f_fwd_true(XW)
stephist!(X[1,:], nbins = 100; normalize = true)

initial_trafo =
    JohnsonTrafo([0.0], [5.0], [0.0], [5.0]) ∘
    inverse(CenterStretch([0.0], [1.0], [0.0])) ∘
    JohnsonTrafo([0.0], [5.0], [0.0], [5.0]) ∘
    inverse(CenterStretch([0.0], [1.0], [0.0]))

optimizer = ADAGrad()

smpls = nestedview(X)
nbatches = 100
nepochs = 10

r = optimize_whitening(smpls, initial_trafo, optimizer, nbatches = nbatches, nepochs = nepochs)

XW2 = r.result(X)
stephist(XW[1,:], nbins = 100; normalize = true)
stephist!(XW2[1,:], nbins = 100; normalize = true)

using EuclidianNormalizingFlows: mvnormal_negll_trafo, mvnormal_negll_trafograd
batches = collect(Iterators.partition(smpls, round(Int, length(smpls)/nbatches)))
plot(r.negll_history, label = "batch negative log-likelihood")
ref_neg_ll = mean([mvnormal_negll_trafo(inverse(f_fwd_true), flatview(b)) for b in batches])
hline!([ref_neg_ll], label = "initial")
ref_neg_ll = mean([mvnormal_negll_trafo(initial_trafo, flatview(b)) for b in batches])
hline!([ref_neg_ll], label = "target")


text:
function init(args::Dict)
  init()
end

function init()

  parent_folders = [
    "src",
    "test"
  ]

  folders = [
    "functions",
    "constants",
    "macros",
    "types",
    "modules"
  ]

  touch("$(pwd())/defaults.jl")

  touch("$(pwd())/input.jl")
  touch("$(pwd())/test/input.jl")

  remove_default_test()
  remove_default_comments()

  for parent_folder in parent_folders
    add_code_to_main_file(parent_folder)

    for folder in folders
      folder_path = "$(pwd())/$parent_folder/$folder"
      if isdir(folder_path) ; continue ; end

      mkdir(folder_path)
      touch("$folder_path/.keep")
    end
  end

  setup_base_dir()
  setup_config()
  setup_docs()
  setup_lib()

end

function remove_default_comments()
  package_name = replace(rsplit(pwd(), "/"; limit=2)[2], ".jl", "")

  main_file_path = "$(pwd())/src/$package_name.jl"
  main_file_text = readstring(main_file_path)

  line_1 = "# package code goes here"
  line_2 = "end # module"
  default_main = "\n$line_1\n\n$line_2"

  new_file = replace(main_file_text, default_main, "\nend")
  open(main_file_path, "w") do old_file
    write(old_file, new_file)
  end
end

function remove_default_test()
  test_file_path = "$(pwd())/test/runtests.jl"
  test_file_text = readstring(test_file_path)

  line_1 = "# write your own tests here"
  line_2 = "@test 1 == 2"
  default_test = "\n$line_1\n$line_2\n"

  new_file = replace(test_file_text, default_test, "")
  open(test_file_path, "w") do old_file
    write(old_file, new_file)
  end
end

function add_code_to_main_file(parent_folder)
  package_name = replace(rsplit(pwd(), "/"; limit=2)[2], ".jl", "")
  if package_name == "Julz" ; return ; end

  template_dictionary = Dict()

  template_dictionary["src"] = ("package", package_name, "module")
  template_dictionary["test"] = ("runtests", "runtests", "")

  template_name, file_name, buzz_word = template_dictionary[parent_folder]

  file_path = "$(pwd())/$parent_folder/$file_name.jl"

  template = generate_file_template(parent_folder, template_name)
  file = readlines(file_path)

  insert_index = length(file)
  for (index, file_line) in enumerate(file)
    if ismatch(r"(?<!\w)Julz(?!\w)", file_line) ; return ; end
    if contains(file_line, buzz_word) && buzz_word != ""
      if insert_index == length(file)
        insert_index = index
      end
    end
  end

  merge_char = ( template_name == "runtests" ) ? "\n" : ""

  new_text_array = ["$(file[insert_index])$(merge_char)$template"]
  splice!(file, insert_index, new_text_array)

  open(file_path, "w") do old_file
    write(old_file, join(file, merge_char))
  end
end

function setup_base_dir()

  tmp_dir = "$(pwd())/tmp"
  mkdir(tmp_dir)
  touch("$tmp_dir/.keep")

  package_name = replace(rsplit(pwd(), "/"; limit=2)[2], ".jl", "")

  user_name = LibGit2.getconfig("github.user", "")

  cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
  user_email = LibGit2.get(cfg, "user.email", "")

  open("README.md", "a") do readme_file
    write(readme_file, "\n[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://$user_name.github.io/$package_name.jl/stable)\n")
    write(readme_file, "\n[![](https://img.shields.io/badge/docs-latest-blue.svg)](https://$user_name.github.io/$package_name.jl/latest)\n")
  end

  open(".gitignore", "a") do gitignore_file
    write(gitignore_file, "\n.DS_Store")
    write(gitignore_file, "\n/tmp/*")
    write(gitignore_file, "\ndocs/build")
    write(gitignore_file, "\ndocs/site")
    write(gitignore_file, "\n*-checkpoint.ipynb")
    write(gitignore_file, "\noutput.jl")
  end

  open("REQUIRE", "a") do require_file
    write(require_file, "Julz")
  end

  travis_code = "  - julia -e 'cd(Pkg.dir(\"$package_name\")); Pkg.add(\"Documenter\"); using Documenter; include(joinpath(\"docs\", \"make.jl\"))'
matrix:
  allow_failures:
  - julia: nightly
script:
  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
  - julia -e 'Pkg.clone(pwd());'
  - (echo \"y\" && echo \"$package_name Test\" && echo \"$user_email\" && echo \"$package_name Test\" && echo \"n\" && yes && cat) | julia -e 'using PkgDev; PkgDev.config();'
  - julia -e 'Pkg.build(\"$package_name\"); Pkg.test(\"$package_name\"; coverage=true)'"

  open(".travis.yml", "a") do require_file
    write(require_file, travis_code)
  end

end

function setup_config()

  config_dir = "$(pwd())/config"

  mkdir(config_dir)
  mkdir("$config_dir/initializers")
  touch("$config_dir/initializers/.keep")

  file_path = "$config_dir/bootload.jl"
  bootload_file = generate_file_template("config", "bootload")
  open(file_path, "w") do file
    write(file, bootload_file)
  end

end

function setup_docs()

  mkdir("$(pwd())/docs")
  mkdir("$(pwd())/docs/src")

  file_path = "$(pwd())/docs/make.jl"
  make_file = generate_file_template("docs", "make")
  open(file_path, "w") do file
    write(file, make_file)
  end

  doc_files = ["index", "code", "glossary"]

  for doc_file in doc_files
    file_path = "$(pwd())/docs/src/$doc_file.md"
    make_file = generate_file_template("docs", "$doc_file.md")

    open(file_path, "w") do file
      write(file, make_file)
    end
  end

end

function setup_lib()
  lib_folder = "$(pwd())/lib"
  if isdir(lib_folder) ; return ; end

  mkdir("$lib_folder")

  cur_sub_dirs = [ "notebooks" , "input_decks" , "tasks" ]

  for cur_sub_dir in cur_sub_dirs
    mkdir("$lib_folder/$cur_sub_dir")
    touch("$lib_folder/$cur_sub_dir/.keep")
  end
end


text: if true ;;;; else ;;; end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
include("test_tools.jl")
include("test_adt.jl")


text:
using Test
using Random

Random.seed!(0)

@testset "Linear" begin

for R in [Float16, Float32, Float64]
    for shape in [(5,), (3, 5), (3, 4, 5)]
        c = randn(R, shape)
        f = Linear(c)
        predicates_test(f)
        x = randn(R, shape)
        @test gradient(f, x) == (c, f(x))
        call_test(f, x)
        prox_test(f, x, R(0.5)+rand(R))
    end
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Postprocess
using JuliaFEM.Testing

testdata = """\
<?xml version="1.0" encoding="utf-8"?>
<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.1">
  <Domain>
    <Grid CollectionType="Temporal" GridType="Collection" Name="Collection">
      <Grid Name="Grid">
      <Time Value="123"/>
        <Geometry Type="XYZ">
          <DataItem DataType="Float" Dimensions="60" Format="XML" Precision="8">
0.0 0.0 0.0
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
0.5 0.0 0.0
0.5 0.5 0.0
0.0 0.5 0.0
0.0 0.0 0.5
0.5 0.0 0.5
0.0 0.5 0.5
1.0 1.0 1.0
2.0 1.0 1.0
1.0 2.0 1.0
1.0 1.0 2.0
1.5 1.0 1.0
1.5 1.5 1.0
1.0 1.5 1.0
1.0 1.0 1.5
1.5 1.0 1.5
1.0 1.5 1.5
</DataItem>
        </Geometry>
        <Topology TopologyType="Mixed" NumberOfElements="2" >
          <DataItem Format="XML" DataType="Int" Dimensions="22">
38 0 1 2 3 4 5 6 7 8 9
38  10 11 12 13 14 15 16 17 18 19
</DataItem>
        </Topology>
        <Attribute Center="Node" Name="Displacement" Type="Vector">
        <DataItem DataType="Float" Dimensions="60" Format="XML" Precision="8">
0.0 0.0 0.0
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
0.5 0.0 0.0
0.5 0.5 0.0
0.0 0.5 0.0
0.0 0.0 0.5
0.5 0.0 0.5
0.0 0.5 0.5
1.0 1.0 1.0
2.0 1.0 1.0
1.0 2.0 1.0
1.0 1.0 2.0
1.5 1.0 1.0
1.5 1.5 1.0
1.0 1.5 1.0
1.0 1.0 1.5
1.5 1.0 1.5
1.0 1.5 1.5
</DataItem>
        </Attribute>
      </Grid>
    </Grid>
  </Domain>
</Xdmf>
"""

@testset "write simple xmf file" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    u = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [0.0, 0.0],
        3 => [0.5, 1.0],
        4 => [0.0, 0.0])
    n = Dict{Int64, Vector{Float64}}(
        2 => [1.0, 0.0],
        3 => [1.0, 0.0])
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Seg2, [2, 3])
    update!([el1, el2], "geometry", X)
    update!([el1, el2], "displacement", u)
    update!(el2, "normal", n)
    xdmf = XDMF()
    xdmf.dimension = 2
    xdmf_new_result!(xdmf, [el1, el2], 0.0)
    xdmf_save_field!(xdmf, [el1, el2], 0.0, "displacement"; field_type="Vector")
    xdmf_save_field!(xdmf, [el1, el2], 0.0, "normal"; field_type="Vector")
    xdmf_save!(xdmf, "/tmp/test.xmf")
    # TODO: how to test?
end


text:
# Oil well with gas lift, model equations & parameters taken from:
# Krishnamoorthy, D., Fjalestad, K. and Skogestad, S., 2019. Optimal operation of oil and gas production using simple feedback control structures. Control Engineering Practice, 91, p.104107.
# https://doi.org/10.1016/j.conengprac.2019.104107
# taken to Julia by evren, 1 well version

using DifferentialEquations
using Parameters
using Plots
plotlyjs()

## Declaring some constants
@with_kw struct WellPar
    # L, H, D = Length, Height, Diameter [m]

    # well
    L_w = 1500      # length well
    H_w = 1000      # height well
    D_w = 0.121     # diameter well

    # bottom hole
    L_bh = 500
    H_bh = 500
    D_bh = 0.121

    # annulus
    L_a =  L_w
    H_a = H_w
    D_a = 0.189

    rho_o = 800     # density of oil, kg/m3
    mu_oil = 1*0.001

    C_iv = 1e-3     # injection valve characteristic, m2
    C_pc = 2e-3     # choke valve characteristic, m2

    GOR = 0.1       # Gas Oil Ratio
    Press_r = 150   # reservoir pressure
    PI = 2.2        # Productivity index, kg/s/bar
    p_m = 20        # manifold pressure
    T_a = 28 + 273  # K
    T_w = 32  + 273 # K

    Mw = 20e-3 # ? units g/mol

end
# change parameters from the default here
par = WellPar()

# and unpacking all the parameters
@unpack_WellPar par

const A_bh  = π*D_bh^2/4
const LA_bh = L_bh*A_bh
const A_w   = π*D_w^2/4
const LwAw  = L_w*A_w
const g     = 9.81      # m/s2
const R     = 8.314     # J/mol K
const V_a   = L_a*(π*(D_a/2)^2 - π*(D_w/2)^2)

##
function gas_well!(dx,x,p,t)

    #* mass flow rates
    mass_ga = x[1] # mass of gas in annulus
    mass_gt = x[2] # mass of gas in tubing
    mass_ot = x[3] # mass of oil in tubing

    #* Parameters
    w_gl = p
    _GOR = GOR
    #* NL calculations

    #Pressures
    p_ai = 1e-5*(((R*T_a/(V_a*Mw) + 9.81*H_a/V_a)*mass_ga*1e3) + (Mw/(R*T_a)*((R*T_a/(V_a*Mw) + 9.81*H_a/V_a)*mass_ga*1e3))*9.81*H_a) # annulus pressure
    p_wh = 1e-5*(((R*T_w/Mw)*(mass_gt*1e3/(L_w*A_w + L_bh*A_bh - mass_ot*1e3/rho_o))) - ((mass_gt*1e3+mass_ot*1e3 )/(L_w*A_w))*9.81*H_w/2) # wellhead pressure

    # densities
    rho_ai = 1e-2*(Mw/(R*T_a)*p_ai*1e5)  # gas, in annulus
    rho_m = 1e-2*(((mass_gt*1e3 + mass_ot*1e3)*p_wh*1e5*Mw*rho_o)/(mass_ot*1e3*p_wh*1e5*Mw + rho_o*R*T_w*mass_gt*1e3))  # mixture, in tubing

    w_pc = C_pc*sqrt(rho_m*1e2*(p_wh*1e5 - p_m*1e5)) # total flow through choke
    w_pg = (mass_gt*1e3/(mass_gt*1e3+mass_ot*1e3))*w_pc # produced gas rate
    w_po = (mass_ot*1e3/(mass_gt*1e3+mass_ot*1e3))*w_pc # produced oil rate

    # pressure at well injection point
    p_wi = 1e-5*((p_wh*1e5 + 9.81/(A_w*L_w)*(mass_ot*1e3+mass_gt*1e3-rho_o*L_bh*A_bh)*H_w + 128*mu_oil*L_w*w_pc/(3.141*D_w^4*((mass_gt*1e3 + mass_ot*1e3)*p_wh*1e5*Mw*rho_o)/(mass_ot*1e3*p_wh*1e5*Mw + rho_o*R*T_w*mass_gt*1e3))))
    p_bh = 1e-5*(p_wi*1e5 + rho_o*9.81*H_bh + 128*mu_oil*L_bh*w_po/(3.14*D_bh^4*rho_o))
    w_iv = C_iv*sqrt(rho_ai*1e2*(p_ai*1e5 - p_wi*1e5))
    w_ro = (PI)*1e-6*(Press_r*1e5 - p_bh*1e5)
    w_rg = 1e1*_GOR*w_ro

    #* Differential Equations
    dx[1] = (w_gl - w_iv)*1e-3
    dx[2] = (w_iv + w_rg*1e-1 - w_pg)*1e-3
    dx[3] = (w_ro - w_po)*1e-3
end

#* Initial conditions

x0 = zeros(3)
#
x0[1:3]= [1.32, 0.8, 6.0]

#* Parameters
w_gl = 4.0 # kg/s

p = w_gl # GOR can also be an uncertain parameter

tspan = (0.0, 2000.0)
prob = ODEProblem(gas_well!, x0, tspan, p)
sol = solve(prob)

plot(sol, label = ["gas annulus" "gas tubing" "oil tubing"])

max_stars_repo_path: ODE in Julia_corrected_evren.jl

max_issues_repo_path: ODE in Julia_corrected_evren.jl

max_forks_repo_path: ODE in Julia_corrected_evren.jl


text:
mutable struct Feature{T}
    fixation_time::Float64
    visible::Bool
	value::T
end

Feature(;fixation_time=0.0, visible=false, value) = Feature(fixation_time, visible, value)

function populate_features(features, values)
	vals = [Feature(value=v) for v in values]
	return NamedTuple{features}(vals)
end

mutable struct VisualObject{F,T} <:AbstractVisualObject
	features::F
    attended::Bool
    visible::Bool
	width::Float64
	angular_size::Float64
	activation::T
	bottomup_activation::Float64
	topdown_activation::Float64
	location::Vector{Float64}
	attend_time::Float64
	target::Bool
end

function VisualObject(;features, attended=false, visible=false, width=0.0, angular_size=0.0, location=[0.0,0.0],
		attend_time=0.0, target=false, activation=0.0)
	return VisualObject(features, attended, visible, width, 0.0, activation, 0.0, 0.0, location, attend_time,
		target)
end

"""
* `viewing_distance`: distance between model and screen in inches (30 inches)
* `current_time`: current model processing time during a trial
* `topdown_weight`: a weight for the influence of top-down activation (1.1)
* `bottomup_weight`: a weight for the influence of bottom-up activation (0.4)
* `noise`: noise added to visual activation
* `persistence`: time during which visible visual objects can stay in iconic memory
* `acuity`: NamedTuple of acuity parameters
* `n_finst`: number of visual objects in finst
* `finst_span`: the duration of finst
* `τₐ`: activation threshold for terminating search
*  `Δτ`: activation threshold increment following a fixation to a distractor
"""
mutable struct Parm{A,T,T1} <: AbstractParms
    viewing_distance::Float64
	topdown_weight::T1
	bottomup_weight::Float64
	noise::Bool
	rnd_time::Bool
	σ::Float64
	persistence::Float64
	acuity::A
	n_finst::Int64
	finst_span::Float64
	β₀exe::Float64
	Δexe::Float64
	τₐ::Float64
	Δτ::Float64
	misc::T
end

function Parm(;viewing_distance=30.0, topdown_weight=.66, bottomup_weight=1.1, noise=false, rnd_time=false,
	σ=.2*π/sqrt(3), persistence=4.0, a_color=.104, b_color=.85, a_shape=.142, b_shape=.96, n_finst=4,
	finst_span=3.0, β₀exe=.02, Δexe=.002, τₐ=0.0, Δτ=0.39, args...)
	acuity = (color = (a=a_color,b=b_color), shape = (a=a_shape,b=b_shape))
	return Parm(viewing_distance, topdown_weight, bottomup_weight,
		noise, rnd_time, σ, persistence, acuity, n_finst, finst_span, β₀exe, Δexe, τₐ, Δτ, NamedTuple(args))
end

"""
**ACTRV**

ACTR model object
- `declarative`: declarative memory module
- `imaginal`: imaginal memory module
- `visual`: visual module
- `goal`: goal module
- `visual_location`: visual location module
- `parms`: model parameters
- `time`: model time

Constructor
````julia
ACTRV(;declarative=Declarative(), imaginal=Imaginal(), goal = Goal(),
    scheduler=nothing, visual=nothing, visual_location=nothing, time=0.0, parms...)
````
"""
@concrete mutable struct ACTRV <: AbstractACTR
    declarative
    imaginal
    visual
    visual_location
    goal
	visicon
    parms
    time
end

function ACTRV(;declarative=Declarative(), imaginal=Imaginal(), goal = Goal(),
    scheduler=nothing, visual=nothing, visual_location=nothing, visicon=nothing, time=0.0, parms...)
    parms′ = Parm(;parms...)
    ACTRV(declarative, imaginal, visual, visual_location, goal, visicon, parms′, time)
end

mutable struct Fixation
	target_color::Symbol
	target_shape::Symbol
	attend_time::Float64
	idx::Int
	stop::Bool
end

Fixation(;color, shape, attend_time, idx, stop) = Fixation(color, shape, attend_time, idx, stop)

mutable struct Data
	target_present::Symbol
	target_color::Symbol
	target_shape::Symbol
	trial_type::Symbol
	response::Symbol
	rt::Float64
end

Data() = Data(fill(:_,5)..., 0.0)

"""
** Experiment **

* `array_width`: with of visual array in pixels
* `n_cells`: number of cells in visual array grid
* `n_trials`: number of trials in simulation
* `cell_width`: width of cell that contains ≤ 1 visual objects
* `object_width`: width of visual object in pixels
* `n_color_distractors`: number of distractors for color
* `n_shape_distractors`: number of distractors for shape
* `set_size`: number of elements in visual array
* `colors`: tuple of colors
* `shapes`: tuple of shapes
* `base_rate`: probability that target is present
* `data`: Array of Data for all trials
* `trial_data`: Data for current trial
* `fixations`: Array of fixations for all trials
* `trial_fixations`: fixations for current trial
* `trace`: displays trace if true
* `window`: GUI window
* `canvas`: GUI canvas
* `visible`: displays GUI if true
* `speed`: how quickly to simulate the model if trace is on
* `populate_visicon`: a function that populates the visicon
"""
mutable struct Experiment{T1,T2,F<:Function}
	array_width::Float64
	n_cells::Int64
	n_trials::Int64
	cell_width::Float64
	object_width::Float64
	n_color_distractors::Int64
	n_shape_distractors::Int64
	set_size::Int64
	colors::Vector{Symbol}
	shapes::Vector{Symbol}
	base_rate::Float64
	data::Vector{Data}
	trial_data::Data
	fixations::Vector{Vector{Fixation}}
	trial_fixations::Vector{Fixation}
	trace::Bool
	window::T1
	canvas::T2
	visible::Bool
	speed::Float64
	populate_visicon::F
end

function Experiment(;array_width=428., object_width=32.0, n_cells=8, n_trials=20,
	n_color_distractors=20, set_size=40, n_shape_distractors=20, shapes=[:p,:q], colors=[:red,:blue],
	base_rate=.50, data=Data[], trial_data=Data(), fixations=Vector{Vector{Fixation}}(), trial_fixations=Fixation[],
	trace=false, window=nothing, canvas=nothing, visible=false, speed=1.0, populate_visicon=conjunctive_ratio)
	cell_width = array_width/n_cells
	@argcheck  cell_width > object_width
	@argcheck n_color_distractors + n_shape_distractors + 1 <= n_cells^2
	visible ? ((canvas,window) = setup_window(array_width)) : nothing
    visible ? Gtk.showall(window) : nothing
	return Experiment(array_width, n_cells, n_trials, cell_width, object_width,
		n_color_distractors, n_shape_distractors, set_size, colors, shapes, base_rate,
		data, trial_data, fixations, trial_fixations, trace, window, canvas, visible, speed, populate_visicon)
end

function import_gui()
    path = @__DIR__
    include(path*"/draw.jl")
end


text:
# ------------------------------------------------------------------------------------------
# ## Learning with a single neuron
#
# In this notebook, we'll build a neuron that classifies an image as an apple or as a banana
# using multiple features from the image. We'll **train** our neuron using data from many
# images that have already been correctly categorised; our neuron will thereby **learn**
# what parameters to use, by minimizing a loss function using gradient descent.
#
# We'll do this with almost the simplest possible neuron, namely one that takes just two
# inputs:
# ------------------------------------------------------------------------------------------

include("draw_neural_net.jl")
number_inputs, number_neurons = 2, 1
draw_network([number_inputs, number_neurons])

# ------------------------------------------------------------------------------------------
# To do this, we need to work with and *clean* some real data. Let's get started!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Loading in some data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's load in some real data! We'll use data that we have prepared from photos of apples
# and bananas; it turns out to be stored on disk in data files as "tab-separated values". We
# can read this data in with the `CSV.jl` package, as follows.
# ------------------------------------------------------------------------------------------

;head data/Apple_Golden_1.dat

;head data/bananas.dat

using CSV
using TextParse

applecols, applecolnames = TextParse.csvread("data/Apple_Golden_1.dat", '\t')
bananacols, bananacolnames = TextParse.csvread("data/bananas.dat", '\t');

# ------------------------------------------------------------------------------------------
# Next, we want to use `DataFrames` to store the information from our CSV files.
# ------------------------------------------------------------------------------------------

using DataFrames

# ------------------------------------------------------------------------------------------
# One way we can create a `DataFrame` is to pass a dictionary that contain arrays as values
# and descriptive names for those arrays as keys to the `DataFrame` function:
# ------------------------------------------------------------------------------------------

apples =  DataFrame(Dict(strip(name) => col for (name, col) in zip(applecolnames, applecols)))
bananas = DataFrame(Dict(strip(name) => col for (name, col) in zip(bananacolnames, bananacols)))

# ------------------------------------------------------------------------------------------
# Above, we used a "dictionary comprehension" to create each `DataFrame`. Here is some code
# to create a dictionary, `appledict`, via a dictionary comprehension:
#
# ```julia
# appledict = Dict(strip(name)=>col for (name, col) in zip(applecolnames, applecols))
# ```
# For now, don't worry about the exact syntax used here. Just note that the dictionary
# created associates names with arrays, and that we've used a dictionary like this to create
# a `DataFrame` with named columns!
#
# So for now, each of the two data sets is stored in a `DataFrame` (from the `DataFrames.jl`
# package).
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Roadmap
#
# To use a neuron with two inputs, we will use just two of the data (numbers) for each
# image, say columns 3 and 4, the mean amount of red and the mean amount of green; each data
# point will then be a 2-dimensional vector, and the data points lie on a two-dimensional
# plane. We will have many input data, labelled with an index $i$. We will denote the $i$th
# data point as  $\mathbf{x}^{(i)}$.
#
# The goal is that our neuron will take a single point on the two-dimensional plane as its
# input, and should return a single output that **classifies** it as an apple ($0$) or a
# banana ($1$).
#
# To do so, it must "**learn**" the correct values of its parameters $\mathbf{w}$ and $b$.
# For this learning to take place, we'll need **labels** for each of the input data, which
# identify them as either an apple (0) or as a banana (1).
#
# These labels will, in turn, allow us to create a loss function, which will allow our
# algorithm to learn to determine if a given choice of parameters does a good or a poor job
# of classifying our fruit images.
#
# **So what do we have left to do?**
#
# The above might sound complex, but luckily we can break it down into a series of
# actionable steps:
#
# 1. Clean our input data (amounts of red and green) to get it into a usable format;
# 2. Create a sequence of labels that we can use to identify correct and incorrect
# classifications;
# 3. Define a loss function that contains parameters;
# 4. Implement an algorithm to pick parameters for our neuron by minimizing the loss
# function with respect to the parameters;
# 5. Use all of the above to train our neuron how to classify images as apples or bananas!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Note:
#
# Note that *in general we cannot expect that a single neuron will be adequate for
# classification.*
#
# If a single neuron struggles to classify our images, we may need to use a more complicated
# neural network structure (which corresponds to using a more complicated function).
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Cleaning the data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Usually it will be necessary to "clean" the data in some way, i.e. pre-process it, before
# it can be used for whichever task you are interested in.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Our next *meta*-exercise will be to collect all the data from columns 3 and 4 into a
# *single* Julia vector `x` (of which each entry is itself a vector), and the labels into a
# single vector `y`. Let's do this in a series of steps!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# First, let's practice grabbing a single column of a `DataFrame`. To grab a column, you can
# index into the `DataFrame` with the name of the column you want, passed as a symbol. In
# Julia, symbols are names preceded by a `:`. For example, we could grab the "height" column
# from `apples` by indexing into `apples` with the symbol, `:height`:
#
# ```julia
# apples[:height]
# ```
#
# Index into `apples` to get the "red" column. What is the type of the object returned? How
# many entries does it have?
#
# A) Array, 5

# B) DataArray, 5

# C) Array, 64

# D) DataArray, 64

# E) Array, 492

# F) DataArray, 492
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# We can grab an individual entry of a `DataFrame` by specifying the row index of the entry
# and the column symbol. For example, to access the height of the 4th image of an apple, we
# would execute
#
# ```julia
# apples[4, :height]
# ```
#
# How much red is there in the 63rd image of a banana?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# We want to reorganize data from the 3rd and 4th columns of `apples` and `bananas` to put
# that data in a single array. Let's start by organizing the data from the 3rd and 4th
# columns of `apples` into a single array, `x_apples`. Create `x_apples` such that there is
# a single element in `x_apples` for each image in `apples`. The $i^\mathrm{th}$ element in
# `x_apples` should be a `Vector`, i.e. a 1D `Array`, with two elements - the amount of red
# and the amount of blue in the $i^\mathrm{th}$ image from `apples`.
#
# Similarly create the `Array` `x_bananas` using data from `bananas`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# Next we want to combine the elements of `x_apples` and `x_bananas` into a single array,
# `xs`. `xs` should contain, first, all the elements of `x_apples`, and then all the
# elements of `x_bananas`. Use the `vcat` function to create `xs`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# If you've gotten this far, our data is in the format we want for learning. Now we need
# labels! We want to store a label (either a `0` or a `1` for every apple or banana image in
# our data set in an array called `ys`. Recall that "0" refers to an apple and "1" refers to
# a banana.
#
# Create an array `ys` where the $i^\mathrm{th}$ element of `ys` is a `0` if the
# $i^\mathrm{th}$ element of `xs` is an apple, and where the $i^\mathrm{th}$ element of `ys`
# is a `1` if the $i^\mathrm{th}$ element of `xs` is a banana.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# Add data points for all apple and banana images in our data sets to a plot using
# `scatter`. Plot data points for apples in one color and use a different color for banana
# data points.
#
# Hint: You may want to use the `first` and `last` functions.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## "Learning" by hand
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Intuitively, looking at the plot of the data, we see that it should be "easy" to find a
# function that separates the data into bananas on one side and apples on the other: we just
# need to draw a straight line that divides the two clouds of data. We can do this "by
# hand", as follows.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In the code below, the neuron will learn a function of the form $\sigma(\mathbf{w} \cdot
# \mathbf{x} + b)$. Since $\sigma$ looks like a smooth version of a step function, we can
# think of $\sigma$ classifying based on whether the value of its output argument is less
# than `0.5` or greater than `0.5`.
#
# **Game**: Use the interactive visualization to find suitable values of $\mathbf{w}$ and
# $b$ such that the hyperplane $\sigma(w_1 x_1 + w_2 x_2 + b) = 0.5$ divides the data. This
# is the same as the hyperplane for which $w_1 x_1 + w_2 x_2 + b = 0$ ! (Note that there are
# many such values!)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can solve for $x_2$ via
#
# $$x_2 = -(w_1 x_1 + b) / w_2,$$
#
# and use this to draw the corresponding hyperplane.
# ------------------------------------------------------------------------------------------

using Interact

@manipulate for w1 in -2:0.01:3, w2 in -2:0.01:3, b in -2:0.01:3

    scatter(first.(x_apples), last.(x_apples), m=:cross, label="apples")
    scatter!(first.(x_bananas), last.(x_bananas), label="bananas")

    ylims!(0.3, 0.66)
    xlims!(0.45, 0.75)

    plot!(x -> -(w1*x + b) / w2)
end

# ------------------------------------------------------------------------------------------
# ## How can the neuron *learn* to classify the data?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We are now ready for our first experience of **machine learning**: we will let the neuron
# learn automatically by processing data and tuning model parameters accordingly (the
# process we call "learning")!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# For given values of the parameters $w_1$, $w_2$ and $b$, the function $f_{\mathbf{w}, b}$
# maps a vector of length $2$ to a number between $0$ and $1$ (due to the definition of
# $\sigma$). Now we want to have a neuron *learn* suitable values of these parameters.
#
# We want to discover (learn!) the parameters such that $f$ models the relationship between
# the data we explored above about the fruit images, and outputs a classification of the
# fruit: $0$ if it corresponds to an apple, and $1$ if it is a banana.
#
# So the neuron's input will be a vector of 2 pieces of information about an image; let's
# call the data about the $i$th image $\mathbf{x}^{(i)}$.
# We also are given the label that says which type of fruit it is: $0$ for an apple, and $1$
# for a banana; let's call this *desired* output number $y^{(i)}$. When we feed in the $i$th
# data, $\mathbf{x}^{(i)}$, we want the neuron to give an output that is *as close as
# possible* to the desired output $y^{(i)}$; i.e. it should **minimize** the mean-squared
# distance
#
# $$L_i = [f_{\mathbf{w}, b}(\mathbf{x}^{(i)}) - y^{(i)} ]^2.$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# However, now we see a key difference from what we did previously: the neuron should vary
# its parameters in such a way that it manages to minimize this distance for *all* of the
# input data, simultaneously!
#
# How can we express this mathematically? We once again define a loss function,
# $L(\mathbf{w}, b)$, which tells us "how wrong" we are when the parameters take on the
# given values, and then **minimize** this loss function with respect to all of its
# parameters.
#
# One way to take account of all the data at once is to use the "mean-squared error" loss
# function, which is the mean (squared) over all the differences between the output of the
# network, $f_{\mathbf{w}, b}(\mathbf{x}^{(i)})$ on the $i$th data, and the desired output
# $y^{(i)}$:
#
# $$L_\mathrm{total}(\mathbf{w}, b) = \frac{1}{N} \sum_i L_i = \frac{1}{N} \sum_i
# [f_{\mathbf{w}, b}(\mathbf{x}^{(i)}) - y^{(i)} ]^2,$$
#
# where $N$ is the total number of data in the training set.
#
# Why do we choose this particular loss function? Because the minimum possible value of this
# loss function is $0$ (since it is a sum of squares), and this is reached only when the
# neural network perfectly predicts the output. If we can find a way to minimize this loss
# function, we will get as close as possible to this perfect prediction. (In general,
# though, we won't be able to get an exact prediction.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Minimizing the loss function: *stochastic* gradient descent
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We already know how to minimise loss functions on a computer: we just calculate the
# gradient, and do gradient descent! But here we hit a problem: the function
# $L_\mathrm{total}$ usually has a *lot* of terms, and so calculating the gradient of that
# function will be very time-consuming.
#
# Instead, we will use a variant of gradient descent, called *stochastic* gradient descent.
# Here, the idea is that we will not use the complete loss function; instead, at each step
# we will choose a random data point, number $i$, and do a step of gradient descent for the
# partial loss function $L_i$ *corresponding to only that data point*.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# **Exercise 7:**
#
# Write functions for the partial loss function `L(w, b, x, y)`.
#
# To do this, recall
#
# $$
# \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix};
# \qquad
# \mathbf{w} = \begin{pmatrix} w_1 \\ w_2 \end{pmatrix};
# \qquad
# f_{\mathbf{w}, b}(\mathbf{x}) = \sigma(\mathbf{w} \cdot \mathbf{x} + b),$$
#
# and declare `f(x, w, b)` as in notebook 8.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 8
#
# Write a function for the gradient of `L`, i.e. `∇L(w, b, x, y)`, with respect to the
# parameters $(w_1, w_2, b)$, using finite differences. $∇L$ will be a vector with one
# component per parameter:
#
# $$∇L = \left( \frac{\partial L}{\partial w_1}, \frac{\partial L}{\partial w_2},
# \frac{\partial L}{\partial b} \right).$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 9
#
# Implement stochastic gradient descent in the function `stochastic_gradient_descent(L, w,
# b, xs, ys, N = 1000)`. Use this function to minimize the function $L_\mathrm{total}$.
#
# The algorithm: For each of `N` steps, randomly select an index $i$ into the vector `xs`
# storing your image data. Calculate the gradient of the cost function, $L_i$, for this
# image and update each of the parameters, $p_j$, of $L_i$ according to
#
# $$p_j = p_j - 0.01 * ∇L_j$$
#
# (Here, $j$ signifies the $j^{th}$ parameter of $L$ and similarly the $j^{th}$ component of
# $∇L$.)
#
# `stochastic_gradient_descent` should return the updated values for vector $\mathbf{w}$ and
# scalar $b$.
#
# Optional: Keep track of the value of $L_\mathrm{total}$ over time if you want to visualize
# the learning process.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 10
#
# Use the values of `w` and `b` from the last exercise to see how `f` is classifying a
# couple of the images in out dataset.
#
# In particular, calculate `f` using the 1st and the 900th image in `xs`. For which image is
# the output of `f` closer to the value of its label?
#
# A) The output of `f` for the 1st image in `xs` is closer to its label
#
# B) The output of `f` for the 900th image in `xs` is closer to its label.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 11
#
# Use the `maximum` function to determine the maximum squared distance of the prediction
# from the true value. (For each image, this formula is $y_i - f_{w, b}(x_i)$.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 12
#
# Use `w` and `b` from stochastic gradient descent to draw the function that the network has
# learned, as before, as the hyperplane $w_1 x + w_2 y + b = 0$. Overlay this with the data.
#
#
# Does this hyperplane correctly separate the data? (I.e., is the data for all apples on one
# side of the line, and is the data for all bananas on the other side of the line?)
#
# A) Yes
# B) No
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/12. Tools - Learning with a single neuron.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# file for functions related to interpolating from solution points to
# faces

"""
  Holds all the temporary arrays needed to interpolate the metrics
"""
mutable struct Interpolation{Tmsh, Tdim}
  dxdxi_el::Array{Tmsh, 3}
  dxdxi_el_bar::Array{Tmsh, 3}

  dxdxi_elface::Array{Tmsh, 3}
  dxdxi_elface_bar::Array{Tmsh, 3}

  dxidx_node::Array{Tmsh, 2}
  dxidx_node_bar::Array{Tmsh, 2}

  dxdxi_node::Array{Tmsh, 2}
  dxdxi_node_bar::Array{Tmsh, 2}

  bndry_arr::Array{Boundary, 1}
end

function Interpolation(mesh::PumiMesh2D{Tmsh}) where Tmsh
  dxdxi_el = zeros(Tmsh, 4, mesh.numNodesPerElement, 1)
  dxdxi_el_bar = zeros(dxdxi_el)

  dxdxi_elface = zeros(Tmsh, 4, mesh.sbpface.numnodes, 1)
  dxdxi_elface_bar = zeros(dxdxi_elface)

  dxidx_node = zeros(Tmsh, 2, 2)
  dxidx_node_bar = zeros(dxidx_node)

  dxdxi_node = zeros(Tmsh, 2, 2)
  dxdxi_node_bar = zeros(dxdxi_node)

  bndry_arr = Array{Boundary}(1)
  return Interpolation{Tmsh, 2}(dxdxi_el, dxdxi_el_bar,
                                dxdxi_elface, dxdxi_elface_bar,
                                dxidx_node, dxidx_node_bar,
                                dxdxi_node, dxdxi_node_bar,
                                bndry_arr)
end

function Interpolation(mesh::PumiMesh2D, ::Type{Tmsh}) where Tmsh
  dxdxi_el = zeros(Tmsh, 4, mesh.numNodesPerElement, 1)
  dxdxi_el_bar = zeros(dxdxi_el)

  dxdxi_elface = zeros(Tmsh, 4, mesh.sbpface.numnodes, 1)
  dxdxi_elface_bar = zeros(dxdxi_elface)

  dxidx_node = zeros(Tmsh, 2, 2)
  dxidx_node_bar = zeros(dxidx_node)

  dxdxi_node = zeros(Tmsh, 2, 2)
  dxdxi_node_bar = zeros(dxdxi_node)

  bndry_arr = Array{Boundary}(1)
  return Interpolation{Tmsh, 2}(dxdxi_el, dxdxi_el_bar,
                                dxdxi_elface, dxdxi_elface_bar,
                                dxidx_node, dxidx_node_bar,
                                dxdxi_node, dxdxi_node_bar,
                                bndry_arr)
end

function Interpolation(mesh::PumiMesh3D{Tmsh}) where Tmsh
  dxdxi_el = zeros(Tmsh, 9, mesh.numNodesPerElement, 1)
  dxdxi_el_bar = zeros(dxdxi_el)

  dxdxi_elface = zeros(Tmsh, 9, mesh.sbpface.numnodes, 1)
  dxdxi_elface_bar = zeros(dxdxi_elface)

  dxidx_node = zeros(Tmsh, 3, 3)
  dxidx_node_bar = zeros(dxidx_node)

  dxdxi_node = zeros(Tmsh, 3, 3)
  dxdxi_node_bar = zeros(dxdxi_node)

  bndry_arr = Array{Boundary}(1)
  return Interpolation{Tmsh, 3}(dxdxi_el, dxdxi_el_bar,
                                dxdxi_elface, dxdxi_elface_bar,
                                dxidx_node, dxidx_node_bar,
                                dxdxi_node, dxdxi_node_bar,
                                bndry_arr)
end

function Interpolation(mesh::PumiMesh3D, ::Type{Tmsh}) where Tmsh
  dxdxi_el = zeros(Tmsh, 9, mesh.numNodesPerElement, 1)
  dxdxi_el_bar = zeros(dxdxi_el)

  dxdxi_elface = zeros(Tmsh, 9, mesh.sbpface.numnodes, 1)
  dxdxi_elface_bar = zeros(dxdxi_elface)

  dxidx_node = zeros(Tmsh, 3, 3)
  dxidx_node_bar = zeros(dxidx_node)

  dxdxi_node = zeros(Tmsh, 3, 3)
  dxdxi_node_bar = zeros(dxdxi_node)

  bndry_arr = Array{Boundary}(1)
  return Interpolation{Tmsh, 3}(dxdxi_el, dxdxi_el_bar,
                                dxdxi_elface, dxdxi_elface_bar,
                                dxidx_node, dxidx_node_bar,
                                dxdxi_node, dxdxi_node_bar,
                                bndry_arr)
end

"""
  Interpolates dxidx and jac to the face nodes, faces shared in parallel and
  boundary faces, as well as gets the coordinates of the nodes for
  the shared faces and boundaryfaces.  All the functions use smart allocators,
  so the fields of the mesh object will be allocated if needed, but not if
  they were previously allocated.

  This is the main entry point for everything in this file.
"""
function interpolateCoordinatesAndMetrics(mesh::PumiMeshDG)

  # interpolate metric terms
  interpolateMapping(mesh)

  #  allocate coordinate arrays if needed
  allocateFaceCoordinates(mesh)

  # get the coordinates
  getBndryCoordinates(mesh, mesh.bndryfaces, mesh.coords_bndry)
  getInterfaceCoordinates(mesh, mesh.interfaces, mesh.coords_interface)
  for i=1:mesh.npeers
    getBndryCoordinates(mesh, mesh.bndries_local[i],
                        mesh.coords_sharedface[i])
  end

  return nothing
end

"""
  Allocates the dxidx and jac arrays for the face nodes/ shared face
  data.  Non-curvilinear meshes only.  Specifically:
    mesh.dxidx_face, mesh.jac_face
    mesh.dxidx_bndry, mesh.jac_bndry
    mesh.dxidx_sharedface, mesh.jac_shareface (and their internal arrays)

    Also the corresponding bar variables
      mesh.dxidx_face_bar, mesh.jac_face_bar
      mesh.dxidx_bndry_bar, mesh.jac_bndry_bar
      mesh.dxidx_sharedface_bar, mesh.jac_sharedface_bar
"""
function allocateInterpolatedMetrics(mesh::PumiMeshDG{Tmsh}) where Tmsh

  dim = mesh.dim
  sbpface = mesh.sbpface
   if !isFieldDefined(mesh, :dxidx_face, :jac_face, :dxidx_sharedface,
                           :jac_sharedface, :dxidx_bndry, :jac_bndry)
    # interior faces
    mesh.dxidx_face = zeros(Tmsh, dim, dim, sbpface.numnodes, mesh.numInterfaces)
    mesh.jac_face = zeros(Tmsh, sbpface.numnodes, mesh.numInterfaces)

    mesh.dxidx_face_bar = zeros(mesh.dxidx_face)
    mesh.jac_face_bar = zeros(mesh.jac_face)

    # boundary faces
    mesh.dxidx_bndry = zeros(Tmsh, dim, dim, sbpface.numnodes,
                                   mesh.numBoundaryFaces)
    mesh.jac_bndry = zeros(Tmsh, sbpface.numnodes, mesh.numBoundaryFaces)
    mesh.dxidx_bndry_bar = zeros(mesh.dxidx_bndry)
    mesh.jac_bndry_bar = zeros(mesh.jac_bndry)

    # parallel shared faces
    mesh.dxidx_sharedface = Array{Array{Tmsh, 4}}(mesh.npeers)
    mesh.jac_sharedface = Array{Array{Tmsh, 2}}(mesh.npeers)
    mesh.dxidx_sharedface_bar = Array{Array{Tmsh, 4}}(mesh.npeers)
    mesh.jac_sharedface_bar = Array{Array{Tmsh, 2}}(mesh.npeers)
    for i=1:mesh.npeers
      mesh.dxidx_sharedface[i] = zeros(Tmsh, dim, dim, sbpface.numnodes,
                                       mesh.peer_face_counts[i])

      mesh.jac_sharedface[i] = Array{Tmsh}(sbpface.numnodes, mesh.peer_face_counts[i])
      mesh.dxidx_sharedface_bar[i] = zeros(mesh.dxidx_sharedface[i])
      mesh.jac_sharedface_bar[i] = zeros(mesh.jac_sharedface[i])

    end
  else
    # zero out metrics?
    fill!(mesh.dxidx_face, 0.0)
    fill!(mesh.jac_face, 0.0)
    fill!(mesh.dxidx_bndry, 0.0)
    fill!(mesh.jac_bndry, 0.0)
    for i=1:mesh.npeers
      fill!(mesh.dxidx_sharedface[i], 0.0)
      fill!(mesh.jac_sharedface[i], 0.0)
    end
  end

  return nothing
end

"""
  Allocates coordinate arays for the boundary and shared face nodes.
  Specifically:
    mesh.coords_bndry
    mesh.coords_sharedface (and its internal arrays)
"""
function allocateFaceCoordinates(mesh::PumiMeshDG{Tmsh}) where Tmsh

  sbpface = mesh.sbpface

  if !isFieldDefined(mesh, :coords_bndry, :coords_sharedface)
    mesh.coords_bndry = zeros(Tmsh, mesh.dim, sbpface.numnodes,
                                  mesh.numBoundaryFaces)
    mesh.coords_bndry_bar = zeros(Tmsh, mesh.dim, sbpface.numnodes,
                                  mesh.numBoundaryFaces)
    mesh.coords_interface = zeros(Tmsh, mesh.dim, sbpface.numnodes, mesh.numInterfaces)
    mesh.coords_sharedface = Array{Array{Tmsh, 3}}(mesh.npeers)
    for i=1:mesh.npeers
      mesh.coords_sharedface[i] = zeros(Tmsh, mesh.dim, sbpface.numnodes,
                                             mesh.peer_face_counts[i])
    end
  else
    fill!(mesh.coords_bndry, 0.0)
    fill!(mesh.coords_interface, 0.0)
    for i=1:mesh.npeers
      fill!(mesh.coords_sharedface[i], 0.0)
    end

  end

  return nothing
end

"""
  Interpolates dxidx, jac to the face nodes of local interfaces, shared
  faces, and boundary faces.  This function uses a smart allocator to
  allocate the relevent fields of the mesh object if they have been been
  allocated previously and to not allocate them if they have been previously.

  By convention, the metrics of elementL of an interface are interpolated
  (recall elementL is always the local elements for shared faces).

  See allocateInterpolatedMetrics for the fields that are populated by this
  function.
"""
function interpolateMapping(mesh::PumiMeshDG{Tmsh}) where Tmsh

  allocateInterpolatedMetrics(mesh)

  sbpface = mesh.sbpface
  dim = mesh.dim

  # pull out the needed arrays
  dxidx_face = mesh.dxidx_face
  dxidx_sharedface = mesh.dxidx_sharedface
  jac_face = mesh.jac_face

  dxidx_bndry = mesh.dxidx_bndry
  jac_bndry = mesh.jac_bndry
  jac_sharedface = mesh.jac_sharedface

  # temporary storage
  interp_data = Interpolation(mesh)

  for i=1:mesh.numInterfaces
    dxidx_i = sview(dxidx_face, :, :, :, i)
    jac_i = sview(jac_face, :, i)

    interface_i = mesh.interfaces[i]
    el = interface_i.elementL
    face = interface_i.faceL
    interp_data.bndry_arr[1] =  Boundary(1, face)

    dxidx_in = sview(mesh.dxidx, :, :, :, el)
    jac_in = sview(mesh.jac, :, el)

    interpolateFace(interp_data, mesh.sbpface, dxidx_in, jac_in, dxidx_i, jac_i)

  end
  # now do shared edges
  for peer = 1:mesh.npeers
    dxidx_p = dxidx_sharedface[peer]
    jac_p = jac_sharedface[peer]
    interfaces_p = mesh.shared_interfaces[peer]
    for i=1:mesh.peer_face_counts[peer]
      dxidx_i = sview(dxidx_p, :, :, :, i)
      jac_i = sview(jac_p, :, i)

      interface_i = interfaces_p[i]
      el = interface_i.elementL
      face = interface_i.faceL

      interp_data.bndry_arr[1] = Boundary(1, face)

      dxidx_in = sview(mesh.dxidx, :, :, :, el)
      jac_in = sview(mesh.jac, :, el)

      interpolateFace(interp_data, mesh.sbpface, dxidx_in, jac_in, dxidx_i, jac_i)
    end
  end

  # now do boundary
  for i=1:mesh.numBoundaryFaces
    bndry = mesh.bndryfaces[i]

    dxidx_i = sview(dxidx_bndry, :, :, :, i)
    jac_i = sview(jac_bndry, :, i)

    el = bndry.element
    dxidx_in = sview(mesh.dxidx, :, :, :, el)
    jac_in = sview(mesh.jac, :, el)
    interp_data.bndry_arr[1] = bndry
    interpolateFace(interp_data, mesh.sbpface, dxidx_in, jac_in, dxidx_i, jac_i)
  end

  return nothing

end  # end function

"""
  Reverse mode of interpolateMapping.  This function propigates the adjoint
  parts of the metrics at the interfaces, boundary faces, and shared faces
  back to the volume nodes.

  This function back propigates both dxidx and jac, although for 2D it turns
  out jac has no effect on the computation.

"""
function interpolateMapping_rev(mesh::PumiMeshDG{Tmsh}) where Tmsh

  sbpface = mesh.sbpface
  dim = mesh.dim

  # pull out the needed arrays
  dxidx_face = mesh.dxidx_face
  dxidx_face_bar = mesh.dxidx_face_bar

  dxidx_sharedface = mesh.dxidx_sharedface
  dxidx_sharedface_bar = mesh.dxidx_sharedface_bar

  dxidx_bndry = mesh.dxidx_bndry
  dxidx_bndry_bar = mesh.dxidx_bndry_bar

  jac_face = mesh.jac_face
  jac_face_bar = mesh.jac_face_bar

  jac_bndry = mesh.jac_bndry
  jac_bndry_bar = mesh.jac_bndry_bar

  jac_sharedface = mesh.jac_sharedface
  jac_sharedface_bar = mesh.jac_sharedface_bar

  # temporary storage
  interp_data = Interpolation(mesh)

  for i=1:mesh.numInterfaces
    dxidx_i = sview(dxidx_face, :, :, :, i)
    dxidx_bar_i = sview(dxidx_face_bar, :, :, :, i)
    jac_i = sview(jac_face, :, i)
    jac_bar_i = sview(jac_face_bar, :, i)

    interface_i = mesh.interfaces[i]
    el = interface_i.elementL
    face = interface_i.faceL
    interp_data.bndry_arr[1] =  Boundary(1, face)

    dxidx_in = sview(mesh.dxidx, :, :, :, el)
    dxidx_bar_in = sview(mesh.dxidx_bar, :, :, :, el)
    jac_in = sview(mesh.jac, :, el)
    jac_bar_in = sview(mesh.jac_bar, :, el)

    interpolateFace_rev(interp_data, mesh.sbpface,
                         dxidx_in, dxidx_bar_in, jac_in, jac_bar_in,
                         dxidx_i, dxidx_bar_i, jac_i, jac_bar_i)

  end
  # now do shared edges
  for peer = 1:mesh.npeers
    dxidx_p = dxidx_sharedface[peer]
    dxidx_bar_p = dxidx_sharedface_bar[peer]
    jac_p = jac_sharedface[peer]
    jac_bar_p = jac_sharedface_bar[perr]

    interfaces_p = mesh.shared_interfaces[peer]
    for i=1:mesh.peer_face_counts[peer]
      dxidx_i = sview(dxidx_p, :, :, :, i)
      dxidx_bar_i = sview(dxidx_bar_p, :, :, :, i)
      jac_i = sview(jac_p, :, i)
      jac_bar_i = sview(jac_bar_p, :, i)

      interface_i = interfaces_p[i]
      el = interface_i.elementL
      face = interface_i.faceL

      interp_data.bndry_arr[1] = Boundary(1, face)

      dxidx_in = sview(mesh.dxidx, :, :, :, el)
      dxidx_bar_in = sview(mesh.dxidx_bar, :, :, :, el)
      jac_in = sview(mesh.jac, :, el)
      jac_bar_in = sview(mesh.jac_bar, :, el)

      interpolateFace_rev(interp_data, mesh.sbpface,
                      dxidx_in, dxidx_bar_in, jac_in, jac_bar_in,
                      dxidx_i, dxidx_bar_i, jac_i, jac_bar_i)
    end
  end

  # now do boundary
  for i=1:mesh.numBoundaryFaces
    bndry = mesh.bndryfaces[i]

    dxidx_i = sview(dxidx_bndry, :, :, :, i)
    dxidx_bar_i = sview(dxidx_bndry_bar, :, :, :, i)
    jac_i = sview(jac_bndry, :, i)
    jac_bar_i = sview(jac_bndry_bar, :, i)

    el = bndry.element

    dxidx_in = sview(mesh.dxidx, :, :, :, el)
    dxidx_bar_in = sview(mesh.dxidx_bar, :, :, :, el)
    jac_in = sview(mesh.jac, :, el)
    jac_bar_in = sview(mesh.jac_bar, :, el)

    interp_data.bndry_arr[1] = bndry
    interpolateFace_rev(interp_data, mesh.sbpface,
                    dxidx_in, dxidx_bar_in, jac_in, jac_bar_in,
                    dxidx_i, dxidx_bar_i, jac_i, jac_bar_i)
  end

  return nothing
end

# 2D version
function interpolateFace(interp_data::Interpolation{Tmsh, 2}, sbpface,
                         dxidx_hat_in::AbstractArray{Tmsh, 3},
                         jac_in::AbstractVector{Tmsh},
                         dxidx_hat_out::AbstractArray{Tmsh, 3},
                         jac_out::AbstractVector{Tmsh}) where Tmsh

  # unpack argumetns
  dxdxi_el = interp_data.dxdxi_el
  dxdxi_elface = interp_data.dxdxi_elface
  dxdxi_node = interp_data.dxdxi_node
  dxidx_node = interp_data.dxidx_node

  dim = 2
  numNodesPerElement = size(dxidx_hat_in, 3)
  # get the data
  for j=1:numNodesPerElement
    dxidx_hat = sview(dxidx_hat_in, :, :, j)
    detJ = jac_in[j]

    adjugate2(dxidx_hat, dxidx_node)
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_el[pos,j,1] = dxidx_node[k, p]
      end
    end
  end

  # interpolate to the face
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)
  boundaryinterpolate!(sbpface, interp_data.bndry_arr, dxdxi_el, dxdxi_elface)

  # now store dxidx, |J| at the boundary nodesa
  for j=1:sbpface.numnodes
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det2(dxdxi_node)
    adjugate2(dxdxi_node, dxidx_node)
    scale!(dxidx_node, 1./det_dxdxi)

    detJ = det2(dxidx_node)
    fac = 1./detJ
    for k=1:dim
      for p=1:dim
        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
    jac_out[j] = detJ
  end

  return nothing

end

"""
  Reverse mode of interpolateFace.

  Propigates the adjoint parts of the inputs to the adjoint parts of the
  outputs.  This function also verifies the primal parts have not been modified
  (ie. it computes values dxidx_hat_out and verifies they are the same as
  the ones already present in the array).

  Inputs (see interpolateFace):
    interp_data
    sbpface
    dxidx_hat_in
    jac_in
    dxidx_hat_out
    dxidx_hat_out_bar: the adjoint part of dxidx_hat_out
    jac_out
    jac_out_bar: the adjoint part of jac_out

  Outputs:
    dxidx_hat_in_bar: the adjoint part of dxidx_hat_in
    jac_in_bar: the adjoint part of jac_in

  Aliasing restrictions: please don't
"""
function interpolateFace_rev(interp_data::Interpolation{Tmsh, 2},
                         sbpface,
                         dxidx_hat_in::AbstractArray{Tmsh, 3},
                         dxidx_hat_in_bar::AbstractArray{Tmsh, 3},
                         jac_in::AbstractVector{Tmsh},
                         jac_in_bar::AbstractVector{Tmsh},
                         dxidx_hat_out::AbstractArray{Tmsh, 3},
                         dxidx_hat_out_bar::AbstractArray{Tmsh, 3},
                         jac_out::AbstractVector{Tmsh},
                         jac_out_bar::AbstractVector{Tmsh}) where Tmsh

  #----------------------------------------------------------------------------
  # forward sweep
  # unpack arguments
  dxdxi_el = interp_data.dxdxi_el
  dxdxi_elface = interp_data.dxdxi_elface
  dxdxi_node = interp_data.dxdxi_node
  dxidx_node = interp_data.dxidx_node

  fill!(dxdxi_el, 0.0)
  fill!(dxdxi_elface, 0.0)

  dim = 2
  numNodesPerElement = size(dxidx_hat_in, 3)

  # get the data
  for j=1:numNodesPerElement
    dxidx_hat = sview(dxidx_hat_in, :, :, j)
#    detJ = jac_in[j]

    adjugate2(dxidx_hat, dxidx_node)
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_el[pos,j,1] = dxidx_node[k, p]
      end
    end
  end

  # interpolate to the face
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)
  boundaryinterpolate!(sbpface, interp_data.bndry_arr, dxdxi_el, dxdxi_elface)

# everything in here is local to the loop, no need to compute it here
  for j=1:sbpface.numnodes
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det2(dxdxi_node)
    adjugate2(dxdxi_node, dxidx_node)
    for i=1:length(dxidx_node)
       dxidx_node[i] = dxidx_node[i] / det_dxdxi
    end

    detJ = det2(dxidx_node)
    fac = 1./detJ

    for k=1:dim
      for p=1:dim
        # verify this forward sweep is consistent with the original function
        @assert (dxidx_hat_out[p, k, j] - dxidx_node[p, k]*fac) < 1e-13
#        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
    jac_out[j] = detJ
  end

  #----------------------------------------------------------------------------
  # reverse sweep

  # unpack arguments
  dxdxi_el_bar = interp_data.dxdxi_el_bar
  dxdxi_elface_bar = interp_data.dxdxi_elface_bar
  dxdxi_node_bar = interp_data.dxdxi_node_bar
  dxidx_node_bar = interp_data.dxidx_node_bar

  # zero them out
  fill!(dxdxi_el_bar, 0.0)
  fill!(dxdxi_elface_bar, 0.0)

  # each interation of this loop is independent - combine forward and reverse
  for j=sbpface.numnodes:-1:1

    #--------------------------------------------------------------------------
    # forward sweep
    fill!(dxdxi_node_bar, 0.0)
    fill!(dxidx_node_bar, 0.0)

    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det2(dxdxi_node)
    adjugate2(dxdxi_node, dxidx_node)
    for i=1:length(dxidx_node)
       dxidx_node[i] = dxidx_node[i] / det_dxdxi
    end

    detJ = det2(dxidx_node)
    fac = 1./detJ
#=
    for k=1:dim
      for p=1:dim
        @assert (dxidx_hat_out[p, k, j] - dxidx_node[p, k]*fac) < 1e-13
        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
    jac_out[j] = detJ
=#
    #--------------------------------------------------------------------------
    # reverse sweep

    det_dxdxi_bar = zero(det_dxdxi)
    fac_bar = zero(eltype(dxidx_hat_out))

    detJ_bar = jac_out_bar[j]*1

    for k=dim:-1:1
      for p=dim:-1:1
        dxidx_node_bar[p, k] += dxidx_hat_out_bar[p, k, j]*fac

        fac_bar += dxidx_hat_out_bar[p, k, j]*dxidx_node[p, k]
      end
    end

    detJ_bar -= fac_bar*1/(detJ*detJ)
    det2_rev(dxidx_node, dxidx_node_bar, detJ_bar)

    # uncomment this to introdue the bug
    for i=1:length(dxidx_node)
       dxidx_node[i] *= det_dxdxi  # need to reverse the primal variables too
       det_dxdxi_bar += dxidx_node_bar[i]*dxidx_node[i]*(-1/(det_dxdxi*det_dxdxi))
       dxidx_node_bar[i] = dxidx_node_bar[i]/det_dxdxi
    end

    adjugate2_rev(dxdxi_node, dxdxi_node_bar, dxidx_node, dxidx_node_bar)
    det2_rev(dxdxi_node, dxdxi_node_bar, det_dxdxi_bar)

    for k=dim:-1:1
      for p=dim:-1:1
        pos = p + dim*(k-1)
        dxdxi_elface_bar[pos, j, 1] += dxdxi_node_bar[k, p]
      end
    end

  end  # end loop j

  # interpolate to the face
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)
  boundaryinterpolate_rev!(sbpface, interp_data.bndry_arr, dxdxi_el_bar, dxdxi_elface_bar)

  # propagate back to the inputs
  for j=numNodesPerElement:-1:1
    fill!(dxidx_node_bar, 0.0)
    fill!(dxdxi_node_bar, 0.0)

    dxidx_hat = sview(dxidx_hat_in, :, :, j)
    dxidx_hat_bar = sview(dxidx_hat_in_bar, :, :, j)

    for k=dim:-1:1
      for p=dim:-1:1
        pos = p + dim*(k-1)
        dxidx_node_bar[k, p] += dxdxi_el_bar[pos, j, 1]
      end
    end

    adjugate2_rev(dxidx_hat, dxidx_hat_bar, dxidx_node, dxidx_node_bar)

    # it turns out jac_in isn't used for anything
#    jac_in_bar[j] = detJ_bar
  end  # end loop j

  return nothing

end

# 3d version
"""
  Interpolates the metric terms from the volume nodes of an element to
  the nodes of a particular face.

  Inputs:
    interp_data:  Interpolation object containing the temporary arrays.
                  Which face is interpolated to depends
                  interp_data.bndry_arr.face.  The element field does
                  not matter.  This array is overwritten.

    sbpface: an SBPFace
    dxidx_hat_in: 3D array holding dxidx/|J| at the volume nodes,
                   2 x 2 x numNodesPerElement in 2D or
                   3 x 3 x numNodesPerElemente in 3D.
    jac_hat_in: determanent of the dxidx at the volume nodes, vector of
                length numNodesPerElement

  Outputs:
    dxidx_i: array to be populated with dxidx/|J| at the face nodes,
             2 x 2 x numNodesPerFace in 2 or
             3 x 3 x numNodesPerFace in 3D

    jac_i: array to be populated with |J|, vector of length numNodesPerFAce

  Aliasing restrictions: do not alias
"""
function interpolateFace(interp_data::Interpolation{Tmsh, 3}, sbpface,
                         dxidx_hat_in::AbstractArray{Tmsh, 3},
                         jac_in::AbstractVector{Tmsh},
                         dxidx_hat_out::AbstractArray{Tmsh, 3},
                         jac_out::AbstractVector{Tmsh}) where Tmsh
# herein, dxidx_hat is dxidx/|J|, and dxidx is the true dxidx

  dxdxi_el = interp_data.dxdxi_el
  dxdxi_elface = interp_data.dxdxi_elface
  dxdxi_node = interp_data.dxdxi_node
  dxidx_node = interp_data.dxidx_node

  dim = 3
  numNodesPerElement = size(dxidx_hat_in, 3)
  # get the data
  for j=1:numNodesPerElement
    dxidx_hat = sview(dxidx_hat_in, :, :, j)
    detJ = jac_in[j]

    adjugate3(dxidx_hat, dxidx_node)
    for k=1:dim
      for p=1:dim
        dxidx_node[k, p] *= detJ
        pos = p + dim*(k-1)
        dxdxi_el[pos,j,1] = dxidx_node[k, p]
      end
    end
  end

  # interpolate to the face
#  face = bndry.face
#  bndry_arr = [Boundary(1, face)]
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)

  boundaryinterpolate!(sbpface, interp_data.bndry_arr, dxdxi_el, dxdxi_elface)

  # now store dxidx, |J| at the boundary nodesa
  for j=1:sbpface.numnodes
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det3(dxdxi_node)
    adjugate3(dxdxi_node, dxidx_node)
    scale!(dxidx_node, 1./det_dxdxi)

    detJ = det3(dxidx_node)
    fac = 1./detJ
    for k=1:dim
      for p=1:dim
        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
    jac_out[j] = detJ
  end

  return nothing

end

function interpolateFace_rev(interp_data::Interpolation{Tmsh, 3},
                         sbpface,
                         dxidx_hat_in::AbstractArray{Tmsh, 3},
                         dxidx_hat_in_bar::AbstractArray{Tmsh, 3},
                         jac_in::AbstractVector{Tmsh},
                         jac_in_bar::AbstractVector{Tmsh},
                         dxidx_hat_out::AbstractArray{Tmsh, 3},
                         dxidx_hat_out_bar::AbstractArray{Tmsh, 3},
                         jac_out::AbstractVector{Tmsh},
                         jac_out_bar::AbstractVector{Tmsh}) where Tmsh

  #----------------------------------------------------------------------------
  # forward sweep
  # unpack arguments
  dxdxi_el = interp_data.dxdxi_el
  dxdxi_elface = interp_data.dxdxi_elface
  dxdxi_node = interp_data.dxdxi_node
  dxidx_node = interp_data.dxidx_node

  fill!(dxdxi_el, 0.0)
  fill!(dxdxi_elface, 0.0)

  dim = 3
  numNodesPerElement = size(dxidx_hat_in, 3)

  # get the data
  for j=1:numNodesPerElement
    dxidx_hat = sview(dxidx_hat_in, :, :, j)
    detJ = jac_in[j]

    adjugate3(dxidx_hat, dxidx_node)
    for k=1:dim
      for p=1:dim
        dxidx_node[k, p] *= detJ
        pos = p + dim*(k-1)
        dxdxi_el[pos,j,1] = dxidx_node[k, p]
      end
    end
  end

  # interpolate to the face
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)
  boundaryinterpolate!(sbpface, interp_data.bndry_arr, dxdxi_el, dxdxi_elface)

  # everything in here is local to the loop, no need to compute it here
  for j=1:sbpface.numnodes
    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det3(dxdxi_node)
    adjugate3(dxdxi_node, dxidx_node)
    for i=1:length(dxidx_node)
       dxidx_node[i] = dxidx_node[i] / det_dxdxi
    end

    detJ = det3(dxidx_node)
    fac = 1./detJ

    for k=1:dim
      for p=1:dim
        # verify this forward sweep is consistent with the original function
        @assert (dxidx_hat_out[p, k, j] - dxidx_node[p, k]*fac) < 1e-13
#        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
#    jac_out[j] = detJ
  end

  #----------------------------------------------------------------------------
  # reverse sweep

  # unpack arguments
  dxdxi_el_bar = interp_data.dxdxi_el_bar
  dxdxi_elface_bar = interp_data.dxdxi_elface_bar
  dxdxi_node_bar = interp_data.dxdxi_node_bar
  dxidx_node_bar = interp_data.dxidx_node_bar

  # zero them out
  fill!(dxdxi_el_bar, 0.0)
  fill!(dxdxi_elface_bar, 0.0)

  # each interation of this loop is independent - combine forward and reverse
  for j=sbpface.numnodes:-1:1

    #--------------------------------------------------------------------------
    # forward sweep
    fill!(dxdxi_node_bar, 0.0)
    fill!(dxidx_node_bar, 0.0)

    for k=1:dim
      for p=1:dim
        pos = p + dim*(k-1)
        dxdxi_node[k,p] = dxdxi_elface[pos, j, 1]
      end
    end

    # inv(A) = adj(A)/|A|
    det_dxdxi = det3(dxdxi_node)
    adjugate3(dxdxi_node, dxidx_node)
    for i=1:length(dxidx_node)
       dxidx_node[i] = dxidx_node[i] / det_dxdxi
    end

    detJ = det3(dxidx_node)
    fac = 1./detJ
#=
    for k=1:dim
      for p=1:dim
        @assert (dxidx_hat_out[p, k, j] - dxidx_node[p, k]*fac) < 1e-13
        dxidx_hat_out[p, k, j] = dxidx_node[p, k]*fac
      end
    end
    jac_out[j] = detJ
=#
    #--------------------------------------------------------------------------
    # reverse sweep
    det_dxdxi_bar = zero(det_dxdxi)
    fac_bar = zero(eltype(dxidx_hat_out))

    detJ_bar = jac_out_bar[j]*1

    for k=dim:-1:1
      for p=dim:-1:1
        dxidx_node_bar[p, k] += dxidx_hat_out_bar[p, k, j]*fac

        fac_bar += dxidx_hat_out_bar[p, k, j]*dxidx_node[p, k]
      end
    end

    detJ_bar -= fac_bar*1/(detJ*detJ)
    det3_rev(dxidx_node, dxidx_node_bar, detJ_bar)

    # uncomment this to introdue the bug
    for i=1:length(dxidx_node)
       dxidx_node[i] *= det_dxdxi  # need to reverse the primal variables too
       det_dxdxi_bar += dxidx_node_bar[i]*dxidx_node[i]*(-1/(det_dxdxi*det_dxdxi))
       dxidx_node_bar[i] = dxidx_node_bar[i]/det_dxdxi
    end

    adjugate3_rev(dxdxi_node, dxdxi_node_bar, dxidx_node, dxidx_node_bar)
    det3_rev(dxdxi_node, dxdxi_node_bar, det_dxdxi_bar)

    for k=dim:-1:1
      for p=dim:-1:1
        pos = p + dim*(k-1)
        dxdxi_elface_bar[pos, j, 1] += dxdxi_node_bar[k, p]
      end
    end

  end  # end loop j

  # interpolate to the face
  face = interp_data.bndry_arr[1].face
  interp_data.bndry_arr[1] = Boundary(1, face)
  boundaryinterpolate_rev!(sbpface, interp_data.bndry_arr, dxdxi_el_bar, dxdxi_elface_bar)

  # propagate back to the inputs
  for j=numNodesPerElement:-1:1
    fill!(dxidx_node_bar, 0.0)
    fill!(dxdxi_node_bar, 0.0)
    detJ = jac_in[j]
    detJ_bar = zero(detJ)

    dxidx_hat = sview(dxidx_hat_in, :, :, j)
    dxidx_hat_bar = sview(dxidx_hat_in_bar, :, :, j)

    for k=dim:-1:1
      for p=dim:-1:1
        pos = p + dim*(k-1)
        dxidx_node_bar[k, p] += dxdxi_el_bar[pos, j, 1]

        dxidx_node[k, p] /= detJ
        detJ_bar += dxidx_node_bar[k, p]*dxidx_node[k, p]
        dxidx_node_bar[k, p] = dxidx_node_bar[k, p]*detJ
      end
    end

    adjugate3_rev(dxidx_hat, dxidx_hat_bar, dxidx_node, dxidx_node_bar)

    jac_in_bar[j] = detJ_bar
  end  # end loop j

  return nothing

end

#TODO: eorder these operations for column major arrays
function det2(A::AbstractMatrix)
# determinet of 2x2 matrix

  return A[1,1]*A[2,2] - A[1,2]*A[2,1]
end

"""
  Computes the reverse mode of det2, returning only the derivative and not
  the value of the determinant.

  Inputs:
    A: The matrix
    y_bar: the adjoint part of the determinant (the primal part is not required
           for the reverse calculation)

  Inputs/Outputs:
    A_bar: the adjoint part of A, which is incremented, not overwritten

"""
function det2_rev(A::AbstractMatrix, A_bar::AbstractMatrix, y_bar::Number)

  A_bar[1,1] += y_bar*A[2,2]
  A_bar[1,2] += -y_bar*A[2,1]
  A_bar[2,1] += -y_bar*A[1,2]
  A_bar[2,2] += y_bar*A[1,1]

  return nothing
end

function det3(A::AbstractMatrix)
# determinent of 3x3 matrix
# cofactor expansion about first column
  val1 = A[1,1]*(A[2,2]*A[3,3] - A[2,3]*A[3,2])
  val2 = -A[2,1]*(A[1,2]*A[3,3] - A[1,3]*A[3,2])
  val3 = A[3,1]*(A[1,2]*A[2,3] - A[1,3]*A[2,2])
  return val1 + val2 + val3
end

function det3_rev(A::AbstractMatrix, A_bar::AbstractMatrix, y_bar::Number)

  # initialize intermediate adjoints to zero
  val1_bar = y_bar
  val2_bar = y_bar
  val3_bar = y_bar

  # val3 expression
  # for true alias compatability, need to apply reverse mode to each line
  # right to left, rather than left to right
  A_bar[1,1] += val3_bar*(A[2,2]*A[3,3] - A[2,3]*A[3,2])
  A_bar[2,2] += val3_bar*(A[1,1]*A[3,3])
  A_bar[3,3] += val3_bar*(A[1,1]*A[2,2])

  A_bar[2,3] += -val3_bar*A[1,1]*A[3,2]
  A_bar[3,2] += -val3_bar*A[1,1]*A[2,3]

  # val2 expression
  A_bar[2,1] += -val2_bar*(A[1,2]*A[3,3] - A[1,3]*A[3,2])
  A_bar[1,2] += -val2_bar*A[2,1]*A[3,3]
  A_bar[3,3] += -val2_bar*A[2,1]*A[1,2]

  A_bar[1,3] += val2_bar*A[2,1]*A[3,2]
  A_bar[3,2] += val2_bar*A[2,1]*A[1,3]

  # val1 expression
  A_bar[3,1] += val1_bar*(A[1,2]*A[2,3] - A[1,3]*A[2,2])
  A_bar[1,2] += val1_bar*A[3,1]*A[2,3]
  A_bar[2,3] += val1_bar*A[3,1]*A[1,2]

  A_bar[1,3] += -val1_bar*A[3,1]*A[2,2]
  A_bar[2,2] += -val1_bar*A[3,1]*A[1,3]

  return nothing
end

function adjugate2(A::AbstractMatrix, B::AbstractMatrix)
  B[1,1] = A[2,2]
  B[1,2] = -A[1,2]
  B[2,1] = -A[2,1]
  B[2,2] = A[1,1]
end

"""
  Computes the reverse mode of adjugate2, without modifying A or B

  Inputs:
    A: The original matrix
    B: The adjuage matrix (presumably computed by adjugate2)
    B_bar: the adjoint part of B

  Inputs/Outputs:
    A_bar: The adjoint part of A, which is incremented, not overwritten

  Aliasing restrictions: I'm not sure what happens if things alias, so don't
                         do it.
"""
function adjugate2_rev(A::AbstractMatrix, A_bar::AbstractMatrix,
                       B::AbstractMatrix, B_bar::AbstractMatrix)

  A_bar[2,2] += B_bar[1,1]
  A_bar[1,2] += -B_bar[1,2]
  A_bar[2,1] += -B_bar[2,1]
  A_bar[1,1] += B_bar[2,2]

end

function adjugate3(A::AbstractMatrix, B::AbstractMatrix)

  B[1,1] = A[2,2]*A[3,3] - A[2,3]*A[3,2]
  B[1,2] = -(A[1,2]*A[3,3] - A[1,3]*A[3,2])
  B[1,3] = A[1,2]*A[2,3] - A[1,3]*A[2,2]

  B[2,1] = -(A[2,1]*A[3,3] - A[2,3]*A[3,1])
  B[2,2] = A[1,1]*A[3,3] - A[1,3]*A[3,1]
  B[2,3] = -(A[1,1]*A[2,3] - A[1,3]*A[2,1])

  B[3,1] = A[2,1]*A[3,2] - A[2,2]*A[3,1]
  B[3,2] = -(A[1,1]*A[3,2] - A[1,2]*A[3,1])
  B[3,3] = A[1,1]*A[2,2] - A[1,2]*A[2,1]
end

"""
  Computes the reverse mode of adjugate2, without modifying A or B

  Inputs:
    A: the original matrix
    B: the adjugate matrix (presumably computed by adjugate3)
    B_bar: the adjoint part of B

  Inputs/Outputs
    A_bar: the adjoint part of A, which is incremented, not overwritten

  Aliasing restrictions: Please don't
"""
function adjugate3_rev(A::AbstractMatrix, A_bar::AbstractMatrix,
                       B::AbstractMatrix, B_bar::AbstractMatrix)

  # B11
  B11_bar = B_bar[1,1]
  A_bar[2,2] += B11_bar*A[3,3]
  A_bar[3,3] += B11_bar*A[2,2]
  A_bar[2,3] += -B11_bar*A[3,2]
  A_bar[3,2] += -B11_bar*A[2,3]

  bar = B_bar[1,2]
  A_bar[1,2] += -bar*A[3,3]
  A_bar[3,3] += -bar*A[1,2]
  A_bar[1,3] += bar*A[3,2]
  A_bar[3,2] += bar*A[1,3]

  bar = B_bar[1,3]
  A_bar[1,2] += bar*A[2,3]
  A_bar[2,3] += bar*A[1,2]
  A_bar[1,3] += -bar*A[2,2]
  A_bar[2,2] += -bar*A[1,3]

  bar = B_bar[2,1]
  A_bar[2,1] += -bar*A[3,3]
  A_bar[3,3] += -bar*A[2,1]
  A_bar[2,3] += bar*A[3,1]
  A_bar[3,1] += bar*A[2,3]

  bar = B_bar[2,2]
  A_bar[1,1] += bar*A[3,3]
  A_bar[3,3] += bar*A[1,1]
  A_bar[1,3] += -bar*A[3,1]
  A_bar[3,1] += -bar*A[1,3]

  bar = B_bar[2,3]
  A_bar[1,1] += -bar*A[2,3]
  A_bar[2,3] += -bar*A[1,1]
  A_bar[1,3] += bar*A[2,1]
  A_bar[2,1] += bar*A[1,3]

  bar = B_bar[3,1]
  A_bar[2,1] += bar*A[3,2]
  A_bar[3,2] += bar*A[2,1]
  A_bar[2,2] += -bar*A[3,1]
  A_bar[3,1] += -bar*A[2,2]

  bar = B_bar[3,2]
  A_bar[1,1] += -bar*A[3,2]
  A_bar[3,2] += -bar*A[1,1]
  A_bar[1,2] += bar*A[3,1]
  A_bar[3,1] += bar*A[1,2]

  bar = B_bar[3,3]
  A_bar[1,1] += bar*A[2,2]
  A_bar[2,2] += bar*A[1,1]
  A_bar[1,2] += -bar*A[2,1]
  A_bar[2,1] += -bar*A[1,2]

  return nothing
end


text:
import LinearAlgebra: eigvals

function orig_to_line( c1::Number, c2::Number )
    if isapprox(c1,c2)
        return abs(c1)
    else
        x1 = real(c1);
        y1 = imag(c1);
        x2 = real(c2);
        y2 = imag(c2);
        return abs((x2-x1)*y1-x1*(y2-y1))/sqrt((x2-x1)^2+(y2-y1)^2);
    end
end

function dist_to_neighbour( ar )
    d = length(ar);
    #return abs(ar - ar[mod((0:d-1)-1,d)+1]);
    neighbour_index = broadcast(mod, (1:d) .-2, d) .+ 1
    return abs.(ar - ar[neighbour_index]);
end

"""
worstfidelity(U,V)

For two unitaries `U` and `V`, `worstfidelity(U,V)` is the worst case
Jozsa fidelity between `U*a` and `V*a`, where `a` is a complex vector
with norm 1. That is, `worstfidelity(U,V)` is equal to the minimum
of `|a'*U'*V*a|^2` over all complex vectors `a` with unit norm.
"""
function worstfidelity(u::Matrix, v::Matrix)
    if size(u) != size(v)
        error("Input matrices do not have the same size")
    end

    #local variables
    e  = eigvals(u'*v) # eigenvalues of u'*v
    shifted_e = map(x -> angle(x) + pi, e)
    es = sortslices(hcat(e,shifted_e),dims=1,by=x->real(x[2]))
    d  = es[:,1]
    f  = 0

    #es(:,2) = es(:,2)-pi;
    n = length(d)

    # find the minimum distance from the convex hull
    # of the distinct eigenvalues to the origin. if
    # the origin is contained in the convex hull,
    # that distance is defined as 0.
    if n==1
        f = 1
    else
        if n==2
            f = orig_to_line(d[1],d[2])
        else
            dn = dist_to_neighbour(angle.(d))
            dn[1] = 2*pi-sum(dn[2:n])  # the sum of the angular separations is 2*pi
            # and the boundary cases are funny
            # so it is best to calc it this way
            dn = findall(dn .> pi)
            if length(dn)==1
                f = orig_to_line(d[dn[1]],d[mod(dn[1]-2,n)+1])
            end
        end
    end
    return abs(f)^2
end

"""
ddist(U,V)

Diamond norm distance between two linear CPTP superoperators.
  Equivalent to dnrom(E-F), but under the assumption that `E` and `F`
  are CPTP supeoperators, but the `ddist` call is
  much more accurate and faster due to properties of CPTP
  superoperators.

See `ddistu` for a similar implementation of `dnorm` customized
to the difference between unitary operations.

"""
ddist(E::Matrix,F::Matrix) = dnormcptp2(E,F)

"""
ddist(U,V)

Diamond norm distance between two unitary operations.
  Equivalent to `dnorm(liou(U)-liou(V))`, under the assumption `U` and
  `V` are unitary matrices. However the `ddist` call is much more
  accurate and much faster due to properties of unitary matrices and
  the corresponding superoperators.

  **Note:** for this particular case, the matrices in question are not
    the superoperators corresponding to the unitary operation, but
    rather the unitary operations themselves.

See `ddist` for a similar implementation of `dnorm` customized
to the difference between CPTP maps.

"""
function ddistu(U::Matrix,V::Matrix)
    w = worstfidelity(U,V)
    return 2*sqrt(1-w)
end


text:
using testJuliaPackage
using Test

#without export
#testJuliaPackage.my_f(2,1)
# by using export, this is available, otherwith an UndefVarError is thrown
#my_f(2,1)

@testset "testJuliaPackage.jl" begin
    #2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end


text:
"""

# `backdoor_paths`

$(SIGNATURES)

Internal.
"""
function backdoor_paths(d::DAG, paths::Vector{Vector{Symbol}}, f::Symbol)
  backdoors = Vector{Symbol}[]
  for path in paths
    if d.a[f, path[2]] == 1
      append!(backdoors, [path])
    end
  end
  backdoors
end

export
  backdoor_paths


text:
## use is diffrent than plots
## using DThree, Plots # no using Plots; plotly()
## plot(sin, 0, 2pi)
## figure()
## plot(cos, 0, 2pi)
## plot!(sin)

using Plots
plotly()
PlotlyBackend = Plots.PlotlyBackend

function html_head(plt::Plots.Plot{PlotlyBackend})
    "<script src=\"$(Pkg.dir("Plots","deps","plotly-latest.min.js"))\"></script>\n" *
    "<script src=\"https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.16/d3.min.js\"></script>"
end

# ----------------------------------------------------------------

Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", p::Plots.Plot{PlotlyBackend}) = display(_display, p)

function Base.writemime(io::IO, ::MIME"text/html", plt::Plots.Plot{PlotlyBackend})
  write(io, html_head(plt) * Plots.html_body(plt))
end