text:
export ProdFinal

import Base: eltype

"""
    ProdFinal{T}

The `ProdFinal` is a wrapper around a distribution. By passing it as a message along an edge of the graph the corresponding marginal is calculated as the distribution of the `ProdFinal`.
In a sense, the `ProdFinal` ignores any further prod with any other distribution for calculating the marginal and only check for variate types of two distributions. Trying to prod two instances of `ProdFinal` will result in an error.
Note: `ProdFinal` is not a prod strategy, as opposed to `ProdAnalytical` and `ProdGeneric`.

See also: [`BIFM`]
"""
struct ProdFinal{T}
    dist :: T
end

ProdFinal(prod::ProdFinal) = prod

Base.show(io::IO, prod::ProdFinal) = print(io, "ProdFinal(", getdist(prod), ")")
Base.eltype(prod::ProdFinal)       = Base.eltype(getdist(prod))

getdist(dist::ProdFinal) = dist.dist

MacroHelpers.@proxy_methods ProdFinal getdist [
    Distributions.mean,
    Distributions.median,
    Distributions.mode,
    Distributions.shape,
    Distributions.scale,
    Distributions.rate,
    Distributions.var,
    Distributions.std,
    Distributions.cov,
    Distributions.invcov,
    Distributions.logdetcov,
    Distributions.entropy,
    Distributions.params,
    Base.precision,
    Base.length,
    Base.ndims,
    Base.size,
    mean_cov,
    mean_var,
    mean_invcov,
    mean_precision,
    weightedmean_cov,
    weightedmean_var,
    weightedmean_invcov,
    weightedmean_precision,
    probvec,
    weightedmean,
    inversemean,
    logmean,
    meanlogmean,
    mirroredlogmean,
    loggammamean
]

Distributions.pdf(prod::ProdFinal, x)    = Distributions.pdf(getdist(prod), x)
Distributions.logpdf(prod::ProdFinal, x) = Distributions.logpdf(getdist(prod), x)

custom_isapprox(dist1::ProdFinal, dist2::ProdFinal; kwargs...)  = custom_isapprox(getdist(dist1), getdist(dist2); kwargs...)
custom_isapprox(dist1::Any, dist2::ProdFinal; kwargs...)        = custom_isapprox(dist1, getdist(dist2); kwargs...)
custom_isapprox(dist1::ProdFinal, dist2::Any; kwargs...)        = custom_isapprox(getdist(dist1), dist2; kwargs...)

# This function is unsafe and uses internal fields of Julia types, but should be used only in tests
convert_eltype(::Type{ <: ProdFinal }, ::Type{T}, prod::ProdFinal{D}) where { T, D } = ProdFinal(convert_eltype(D.name.wrapper, T, getdist(prod)))

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

# product of distribution message with `ProdFinal` always returns the same `ProdFinal` object directly
prod_final_check_variate_types(::Type{ T }, ::Type{ T },   result) where { T <: Distributions.VariateForm } = result
prod_final_check_variate_types(::Type{ T1 }, ::Type{ T2 }, result) where { T1, T2 } = error("Different variate types in a prod with `ProdFinal`: $(T1) × $(T2)")

prod(::ProdAnalytical, left::ProdFinal, right) = prod_final_check_variate_types(variate_form(getdist(left)), variate_form(right), left)
prod(::ProdAnalytical, left, right::ProdFinal) = prod_final_check_variate_types(variate_form(left), variate_form(getdist(right)), right)

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

prod(::ProdAnalytical, left::ProdFinal, right::ProdFinal) = error("Invalid product: `ProdFinal` × `ProdFinal`")


text:
using Printf                    #

mutable struct LinkedInts
    next :: Vector{Int}
    prev :: Vector{Int}

    free_ptr :: Int
    free_cap :: Int
    root     :: Int

    block :: Vector{Int}
    size  :: Vector{Int}
end

LinkedInts(capacity=128) =
    LinkedInts(Int[],Int[],
               0,0,0,
               Int[],
               Int[])

allocatedlist(s::LinkedInts) = findall(s.block .> 0)
allocated(s::LinkedInts, id :: Int) = id > 0 && id <= length(s.block) && s.block[id] > 0
blocksize(s::LinkedInts, id :: Int) = s.size[id]
Base.length(s::LinkedInts) = length(s.next)
numblocks(s::LinkedInts) = length(s.block)

function Base.show(f::IO, s::LinkedInts)
    print(f,"LinkedInts(\n")
    @printf(f,"  Number of blocks: %d\n", length(s.block))
    @printf(f,"  Number of elements: %d\n", length(s.next))
    print(f,"  Blocks:\n")
    for i in 1:length(s.block)
        if s.block[i] > 0
            idxs = getindexes(s,i)
            print(f,"    #$i: $idxs\n")
        end
    end
    p = s.free_ptr
    freelst = Int[]
    while p > 0
        push!(freelst,p)
        p = s.prev[p]
    end
    print(f,"  Free: $freelst\n")

    println(f,"  free_ptr = $(s.free_ptr)")
    println(f,"  root     = $(s.root)")
    println(f,"  next     = $(s.next)")
    println(f,"  prev     = $(s.prev)")

    print(f,")")
end

"""
    ensurefree(s::LinkedInts, N :: Int)

Ensure that there are at least `N` elements free, and allocate as necessary.
"""
function ensurefree(s::LinkedInts, N :: Int)
    if s.free_cap < N
        num = N - s.free_cap

        cap = length(s.next)
        first = cap+1
        last  = cap+num

        append!(s.next, Int[i+1 for i in first:last])
        append!(s.prev, Int[i-1 for i in first:last])

        s.next[last] = 0
        s.prev[first] = s.free_ptr
        if s.prev[first] > 0
            s.next[s.prev[first]] = first
        end
        s.free_ptr = last
        s.free_cap += num

        return num
    else
        return 0
    end
end

"""
    newblock(s::LinkedInts, N :: Int)

Add a new block in list `idx`
"""
function newblock(s::LinkedInts, N :: Int) :: Int
    @assert(N > 0)
    ensurefree(s, N)
    # remove from free list
    ptre = s.free_ptr
    ptrb = ptre
    for i = 1:N-1
        ptrb = s.prev[ptrb]
    end

    prev = s.prev[ptrb]

    if prev > 0
        s.next[prev] = 0
    end

    s.free_ptr = s.prev[ptrb]
    s.free_cap -= N

    # insert into list `idx`
    s.prev[ptrb] = s.root
    if s.root > 0
        s.next[s.root] = ptrb
    end
    s.root = ptre
    push!(s.block, ptrb)
    push!(s.size, N)

    id = length(s.block)

    #if ! checkconsistency(s)
    #    println("List = ",s)
    #    assert(false)
    #end

    return id
end

"""
Move a block to the free list.
"""
function deleteblock(s::LinkedInts, id :: Int)
    if s.size[id] > 0
        ptrb = s.block[id]
        N = s.size[id]
        ptre = ptrb
        for i in 2:N
            ptre = s.next[ptre]
        end
        prev = s.prev[ptrb]
        next = s.next[ptre]

        # remove from list and clear the block id
        if s.root == ptre s.root = prev end
        if prev > 0 s.next[prev] = next end
        if next > 0 s.prev[next] = prev end

        s.size[id]  = 0
        s.block[id] = 0

        # add to free list
        if s.free_ptr > 0
            s.next[s.free_ptr] = ptrb
        end
        s.prev[ptrb] = s.free_ptr
        s.free_ptr = ptre
        s.next[ptre] = 0

        s.free_cap += N
    end
end

"""
    getindex(s::LinkedInts, id::Int)

Shortcut for `getindexes(s, id)[1]` when `s.size[id]` is 1.
"""
function getindex(s::LinkedInts, id::Int)
    @assert s.size[id] == 1
    return s.block[id]
end

"""
    getindexes(s::LinkedInts, id :: Int)

Return the vector of indices for the block `id`.
"""
function getindexes(s::LinkedInts, id :: Int)
    N = s.size[id]
    r = Vector{Int}(undef, N)
    p = s.block[id]
    for i in 1:N
        r[i] = p
        p = s.next[p]
    end
    return r
end

function getindexes(s::LinkedInts, id::Int, target::Vector{Int}, offset::Int)
    N = s.size[id]
    p = s.block[id]
    for i in 1:N
        target[i+offset-1] = p
        p = s.next[p]
    end
    return N
end

function getindexes(s::LinkedInts, ids::Vector{Int})
    N = sum(map(id -> s.size[id], ids))
    r = Vector{Int}(undef,N)
    offset = 1
    for id in ids
        offset += getindexes(s, id, r, offset)
    end
    return r
end

function getoneindex(s::LinkedInts, id :: Int)
    N = s.size[id]
    if N < 1
        error("No values at id")
    end

    s.block[i]
end

"""
Get a list if the currently free elements.
"""
function getfreeindexes(s::LinkedInts)
    N = s.free_cap
    r = Array{Int}(undef,N)
    ptr = s.free_ptr
    for i in 1:N
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Get a list if the currently used elements.
"""
function getusedindexes(s::LinkedInts)
    N = length(s.next) - s.free_cap
    r = Array{Int}(undef,N)
    ptr = s.root
    for i in 1:N
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Check consistency of the internal structures.
"""
function checkconsistency(s::LinkedInts) :: Bool
    if length(s.prev) != length(s.next)
        return false
    end

    N = length(s.prev)

    if ! (all(i -> s.prev[i] == 0 || s.next[s.prev[i]] == i, 1:N) &&
          all(i -> s.next[i] == 0 || s.prev[s.next[i]] == i, 1:N))
        @assert(false)
    end

    mark = fill(false,length(s.prev))

    p = s.free_ptr
    while p != 0
        mark[p] = true
        p = s.prev[p]
    end

    p = s.root
    while p != 0
        @assert(!mark[p])
        mark[p] = true
        p = s.prev[p]
    end

    if !all(mark)
        println(s)
        println(mark)
        @assert(all(mark))
    end

    return true
end

ext: jl


text:
module KyotoCabinet

include("c.jl")

# Basic methods
import Base: open, close

# Iteration
import Base: start, next, done

# Generic collections
import Base: isempty, empty!, length

# Indexed collections
import Base: getindex, setindex!

# Dict
import Base: Associative, haskey, getkey, get, get!, delete!, pop!

using .c

export
  # Types
  Bytes, Db, KyotoCabinetException,

  # Db methods
  get, set!, path, cas, bulkset!, bulkdelete!,
  pack, unpack

typealias Bytes Array{Uint8,1}

type Db{K,V} <: Associative{K,V}
  ptr :: Ptr{Void}

  function Db(null=false)
    if (null)
      new(C_NULL)
    else
      ptr = kcdbnew()
      self = new(ptr)
      finalizer(self, destroy)
      self
    end
  end
end

type Cursor{K,V}
  ptr :: Ptr{Void}
  db :: Db{K,V} # to prevent DB GCed before cursor

  function Cursor()
    new(C_NULL, Db{K,V}(true))
  end

  function Cursor(db::Db{K,V})
    ptr = kcdbcursor(db.ptr)
    self = Cursor{K,V}()
    self.ptr = ptr
    self.db = db
    finalizer(self, destroy)
    self
  end
end

immutable RecordIterator{K,V}
  cursor :: Cursor{K,V}
end

immutable KyotoCabinetException <: Exception
  code :: Int32
  message :: String
end

# Pack value into byte array (Array{Uint8,1})
pack(v::Bytes) = v

# Unpack byte array (Array{Uint8,1}) into value of type T.
# buf is not GCed and will be freed right after unpack
# use copy() to own
unpack(T::Type{Bytes}, buf::Bytes) = copy(buf)

# Generic collections

isempty(db::Db) = (length(db) == 0)

function length(db::Db)
  count = kcdbcount(db.ptr)
  if (count == -1) throw(kcexception(db)) end
  count
end

function empty!(db::Db)
  ok = kcdbclear(db.ptr)
  if (ok == 0) throw(kcexception(db)) end
  db
end

# Iterable interface for Db

function start{K,V}(db::Db{K,V})
  cur = Cursor{K,V}(db)
  _start!(cur) ? RecordIterator{K,V}(cur) : RecordIterator{K,V}(Cursor{K,V}())
end

function next{K,V}(db::Db{K,V}, it::RecordIterator{K,V})
  if done(db, it)
    throw(KyotoCabinetException(KCENOREC, "Can not move forward"))
  end
  kv = _record(it.cursor)
  if (_next!(it.cursor))
    (kv, it)
  else
    destroy(it.cursor)
    (kv, RecordIterator{K,V}(Cursor{K,V}()))
  end
end

function done(db::Db, it::RecordIterator)
  it.cursor.ptr == C_NULL
end

# Db methods

function open{K,V}(db::Db{K,V}, file::String, mode::String)
  open(db, file, _mode(mode))
end

function open{K,V}(f::Function, db::Db{K,V}, file::String, mode::String)
  open(f, db, file, _mode(mode))
end

function open{K,V}(db::Db{K,V}, file::String, mode::Uint)
  ok = kcdbopen(db.ptr, bytestring(file), mode)
  if (ok == 0) throw(kcexception(db)) end
  db
end

function open{K,V}(f::Function, db::Db{K,V}, file::String, mode::Uint)
  db = open(db, file, mode)
  try
    f(db)
  finally
    close(db)
    destroy(db)
  end
end

function close(db::Db)
  ok = kcdbclose(db.ptr)
  if (ok == 0) throw(kcexception(db)) end
end

function destroy(db::Db)
  if db.ptr == C_NULL
    return
  end
  kcdbclose(db.ptr)
  kcdbdel(db.ptr)
  db.ptr = C_NULL
end

function cas{K,V}(db::Db{K,V}, key::K, old::V, new::V)
  kbuf = pack(key)
  ovbuf = pack(old)
  nvbuf = pack(new)
  ok, code = throw_if(db, 0, KCELOGIC) do
    kcdbcas(db.ptr, pointer(kbuf), length(kbuf), pointer(ovbuf), length(ovbuf), pointer(nvbuf), length(nvbuf))
  end
  code == KCESUCCESS
end

# To resolve conflict with V=nothing
cas{K}(db::Db{K,Nothing}, key::K, old::Nothing, new::Nothing) = KCESUCCESS

function cas{K,V}(db::Db{K,V}, key::K, old::Nothing, new::V)
  kbuf = pack(key)
  nvbuf = pack(new)
  ok, code = throw_if(db, 0, KCELOGIC) do
    kcdbcas(db.ptr, pointer(kbuf), length(kbuf), C_NULL, 0, pointer(nvbuf), length(nvbuf))
  end
  code == KCESUCCESS
end

function cas{K,V}(db::Db{K,V}, key::K, old::V, new::Nothing)
  kbuf = pack(key)
  ovbuf = pack(old)
  ok, code = throw_if(db, 0, KCELOGIC) do
    kcdbcas(db.ptr, pointer(kbuf), length(kbuf), pointer(ovbuf), length(ovbuf), C_NULL, 0)
  end
  code == KCESUCCESS
end

# Dict methods
function haskey{K,V}(db::Db{K,V}, k::K)
  kbuf = pack(k)
  v, code = throw_if(db, -1, KCENOREC) do
    kcdbcheck(db.ptr, pointer(kbuf), length(kbuf))
  end
  code != KCENOREC
end

function getkey{K,V}(db::Db{K,V}, key::K, default::K)
  haskey(db, key) ? key : default
end

function set!{K,V}(db::Db{K,V}, k::K, v::V)
  kbuf = pack(k)
  vbuf = pack(v)
  ok = kcdbset(db.ptr, pointer(kbuf), length(kbuf), pointer(vbuf), length(vbuf))
  if (ok == 0) throw(kcexception(db)) end
  v
end

function bulkset!{K,V}(db::Db{K,V}, kvs::Dict{K,V}, atomic::Bool)
  # make a copy to prevent GC
  recbuf = [(pack(k), pack(v)) for (k, v) in kvs]
  recs = [KCREC(KCSTR(k, length(k)), KCSTR(v, length(v))) for (k, v) in recbuf]
  c = kcdbsetbulk(db.ptr, recs, length(recs), atomic ? 1 : 0)
  if (c == -1) throw(kcexception(db)) end
  c
end

function bulkdelete!{K,V}(db::Db{K,V}, keys::Array{K,1}, atomic::Bool)
  keybuf = [pack(k) for k in keys]
  ks = [KCSTR(k, length(k)) for k in keybuf]
  c = kcdbremovebulk(db.ptr, ks, length(ks), atomic ? 1 : 0)
  if (c == -1) throw(kcexception(db)) end
  c
end

function get{K,V}(db::Db{K,V}, k::K)
  kbuf = pack(k)
  vsize = Csize_t[0]
  pv = kcdbget(db.ptr, pointer(kbuf), length(kbuf), pointer(vsize))
  if (pv == C_NULL) throw(kcexception(db)) end
  _unpack(V, pv, int(vsize[1]))
end

get(db::Db, k, default) = get(()->default, db, k)

function get{K,V}(default::Function, db::Db{K,V}, k::K)
  kbuf = pack(k)
  vsize = Csize_t[0]
  pv, code = throw_if(db, C_NULL, KCENOREC) do
    kcdbget(db.ptr, pointer(kbuf), length(kbuf), pointer(vsize))
  end
  code == KCENOREC ? default() : _unpack(V, pv, int(vsize[1]))
end

get!(db::Db, k, default) = get!(()->default, db, k)

function get!{K,V}(default::Function, db::Db{K,V}, k::K)
  kbuf = pack(k)
  vsize = Csize_t[0]
  pv, code = throw_if(db, C_NULL, KCENOREC) do
    kcdbget(db.ptr, pointer(kbuf), length(kbuf), pointer(vsize))
  end
  code == KCENOREC ? set!(db, k, default()) : _unpack(V, pv, int(vsize[1]))
end

function delete!{K,V}(db::Db{K,V}, k::K)
  kbuf = pack(k)
  ok = kcdbremove(db.ptr, pointer(kbuf), length(kbuf))
  if (ok == 0) throw(kcexception(db)) end
  db
end

function pop!{K,V}(db::Db{K,V}, k::K)
  kbuf = pack(k)
  vsize = Csize_t[0]
  pv = kcdbseize(db.ptr, pointer(kbuf), length(kbuf), pointer(vsize))
  if (pv == C_NULL) throw(kcexception(db)) end
  _unpack(V, pv, int(vsize[1]))
end

function pop!{K,V}(db::Db{K,V}, k::K, default::V)
  kbuf = pack(k)
  vsize = Csize_t[0]
  pv, code = throw_if(db, C_NULL, KCENOREC) do
    kcdbseize(db.ptr, pointer(kbuf), length(kbuf), pointer(vsize))
  end
  code == KCENOREC ? default : _unpack(V, pv, int(vsize[1]))
end

# Indexable collection
getindex(db::Db, k) = get(db, k)
setindex!(db::Db, v, k) = set!(db, k, v)

# Cursor

function _start!(cursor::Cursor)
  f(cursor::Cursor) = kccurjump(cursor.ptr)
  _move!(cursor, f)
end

# Move to the next record. Return false if there no next record.
function _next!(cursor::Cursor)
  f(cursor::Cursor) = kccurstep(cursor.ptr)
  _move!(cursor, f)
end

function _move!(cursor::Cursor, f)
  ok, code = throw_if(cursor, 0, KCENOREC) do
    f(cursor)
  end
  code != KCENOREC
end

function _record{K,V}(cursor::Cursor{K,V})
  pkSize = Csize_t[0]
  pvSize = Csize_t[0]
  pv = CString[0]
  pk = kccurget(cursor.ptr, pointer(pkSize), pointer(pv), pointer(pvSize), 0)
  if (pk == C_NULL) throw(kcexception(cursor)) end

  res = (_unpack(K, pk, int(pkSize[1]), false), _unpack(V, pv[1], int(pvSize[1]), false))
  ok = kcfree(pk)
  if (ok == 0) throw(kcexception(cursor)) end

  res
end

function path(db::Db)
  p = kcdbpath(db.ptr)
  v = bytestring(p)
  ok = kcfree(p)
  if (ok == 0) throw(KyotoCabinetException(KCESYSTEM, "Can not free memory")) end
  v
end

# KyotoCabinet exceptions
function throw_if(f::Function, db::Db, result_invalid, ecode_valid)
  result = f()
  if (result == result_invalid)
    code = kcdbecode(db.ptr)
    if (code == ecode_valid)
      return (result, code)
    else
      message = bytestring(kcdbemsg(db.ptr))
      throw(KyotoCabinetException(code, message))
    end
  end
  (result, KCESUCCESS)
end

function throw_if(f::Function, cursor::Cursor, result_invalid, ecode_valid)
  result = f()
  if (result == result_invalid)
    code = kccurecode(cursor.ptr)
    if (code == ecode_valid)
      return (result, code)
    else
      message = bytestring(kccuremsg(cursor.ptr))
      throw(KyotoCabinetException(code, message))
    end
  end
  (result, KCESUCCESS)
end

function kcexception(db::Db)
  @assert db.ptr != C_NULL

  code = kcdbecode(db.ptr)
  message = bytestring(kcdbemsg(db.ptr))

  KyotoCabinetException(code, message)
end

function kcexception(cur::Cursor)
  @assert cur.ptr != C_NULL

  code = kccurecode(cur.ptr)
  message = bytestring(kccuremsg(cur.ptr))

  KyotoCabinetException(code, message)
end

function destroy(cursor::Cursor)
  if cursor.ptr == C_NULL
    return
  end
  kccurdel(cursor.ptr)
  cursor.ptr = C_NULL
end

function _unpack(T, p::Ptr{Uint8}, length, free=true)
  v = unpack(T, pointer_to_array(p, length))
  if free
    ok = kcfree(p)
    if (ok == 0) throw(KyotoCabinetException(KCESYSTEM, "Can not free memory")) end
  end
  v
end

_modes = {
  "r" => KCOREADER,
  "w" => KCOWRITER,
  "w+" => KCOWRITER | KCOCREATE
}
_mode(mode::String) = get(_modes, mode, KCOREADER)

end # module KyotoCabinet


text:
using Test
using Distances, LinearAlgebra
using KernelFunctions

A = rand(10,5)
B = rand(20,5)
@testset "Distance" begin
    @testset "Dot Product" begin
        d = KernelFunctions.DotProduct()
        @test diag(pairwise(d,A,dims=2)) == [dot(A[:,i],A[:,i]) for i in 1:size(A,2)]
        @test_throws DimensionMismatch d(rand(3),rand(4))
        @test d(3.0,2.0) == 6.0
    end
    @testset "Delta" begin
        d = KernelFunctions.Delta()
        @test pairwise(d,A,dims=1) == Matrix(I,size(A,1),size(A,1))
        @test pairwise(d,A,B,dims=1) == zeros(size(A,1),size(B,1))
        @test d(1,2) == 0
        @test d(1,1) == 1
        @test_throws DimensionMismatch d(rand(3),rand(4))
    end
end


text:
module TwoDNavierStokes

export
  Problem,
  set_ζ!,
  updatevars!,

  energy,
  energy_dissipation_hyperviscosity,
  energy_dissipation_hypoviscosity,
  energy_work,
  enstrophy,
  enstrophy_dissipation_hyperviscosity,
  enstrophy_dissipation_hypoviscosity,
  enstrophy_work

using
  CUDA,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

nothingfunction(args...) = nothing

"""
    Problem(dev::Device=CPU();
                nx = 256,
                ny = nx,
                Lx = 2π,
                Ly = Lx,
                 ν = 0,
                nν = 1,
                 μ = 0,
                nμ = 0,
                dt = 0.01,
           stepper = "RK4",
             calcF = nothingfunction,
        stochastic = false,
  aliased_fraction = 1/3,
                 T = Float64)

Construct a two-dimensional Navier-Stokes `problem` on device `dev`.

Keyword arguments
=================
  - `dev`: (required) `CPU()` or `GPU()`; computer architecture used to time-step `problem`.
  - `nx`: Number of grid points in ``x``-domain.
  - `ny`: Number of grid points in ``y``-domain.
  - `Lx`: Extent of the ``x``-domain.
  - `Ly`: Extent of the ``y``-domain.
  - `ν`: Small-scale (hyper)-viscosity coefficient.
  - `nν`: (Hyper)-viscosity order, `nν```≥ 1``.
  - `μ`: Large-scale (hypo)-viscosity coefficient.
  - `nμ`: (Hypo)-viscosity order, `nμ```≤ 0``.
  - `dt`: Time-step.
  - `stepper`: Time-stepping method.
  - `calcF`: Function that calculates the Fourier transform of the forcing, ``F̂``.
  - `stochastic`: `true` or `false`; boolean denoting whether `calcF` is temporally stochastic.
  - `aliased_fraction`: the fraction of high-wavenumbers that are zero-ed out by `dealias!()`.
  - `T`: `Float32` or `Float64`; floating point type used for `problem` data.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
                nx = 256,
                ny = nx,
                Lx = 2π,
                Ly = Lx,
  # Drag and/or hyper-/hypo-viscosity
                 ν = 0,
                nν = 1,
                 μ = 0,
                nμ = 0,
  # Timestepper and equation options
                dt = 0.01,
           stepper = "RK4",
             calcF = nothingfunction,
        stochastic = false,
  # Float type and dealiasing
  aliased_fraction = 1/3,
                 T = Float64)

  grid = TwoDGrid(dev, nx, Lx, ny, Ly; aliased_fraction=aliased_fraction, T=T)

  params = Params{T}(ν, nν, μ, nμ, calcF)

  vars = calcF == nothingfunction ? DecayingVars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid))

  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

"""
    Params{T}(ν, nν, μ, nμ, calcF!)

A struct containing the parameters for the two-dimensional Navier-Stokes. Included are:

$(TYPEDFIELDS)
"""
struct Params{T} <: AbstractParams
    "small-scale (hyper)-viscosity coefficient"
       ν :: T
    "(hyper)-viscosity order, `nν```≥ 1``"
      nν :: Int
    "large-scale (hypo)-viscosity coefficient"
       μ :: T
    "(hypo)-viscosity order, `nμ```≤ 0``"
      nμ :: Int
    "function that calculates the Fourier transform of the forcing, ``F̂``"
  calcF! :: Function
end

Params(ν, nν) = Params(ν, nν, typeof(ν)(0), 0, nothingfunction)

# ---------
# Equations
# ---------

"""
    Equation(params, grid)

Return the `equation` for two-dimensional Navier-Stokes with `params` and `grid`. The linear
operator ``L`` includes (hyper)-viscosity of order ``n_ν`` with coefficient ``ν`` and
hypo-viscocity of order ``n_μ`` with coefficient ``μ``,

```math
L = - ν |𝐤|^{2 n_ν} - μ |𝐤|^{2 n_μ} .
```

Plain old viscocity corresponds to ``n_ν=1`` while ``n_μ=0`` corresponds to linear drag.

The nonlinear term is computed via the function `calcN!`.
"""
function Equation(params::Params, grid::AbstractGrid)
  L = @. - params.ν * grid.Krsq^params.nν - params.μ * grid.Krsq^params.nμ
  CUDA.@allowscalar L[1, 1] = 0

  return FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

abstract type TwoDNavierStokesVars <: AbstractVars end

"""
    Vars{Aphys, Atrans, F, P}(ζ, u, v, ζh, uh, vh, Fh, prevsol)

The variables for two-dimensional Navier-Stokes:

$(FIELDS)
"""
struct Vars{Aphys, Atrans, F, P} <: TwoDNavierStokesVars
    "relative vorticity"
        ζ :: Aphys
    "x-component of velocity"
        u :: Aphys
    "y-component of velocity"
        v :: Aphys
    "Fourier transform of relative vorticity"
       ζh :: Atrans
    "Fourier transform of x-component of velocity"
       uh :: Atrans
    "Fourier transform of y-component of velocity"
       vh :: Atrans
    "Fourier transform of forcing"
       Fh :: F
    "`sol` at previous time-step"
  prevsol :: P
end

const DecayingVars = Vars{<:AbstractArray, <:AbstractArray, Nothing, Nothing}
const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    DecayingVars(dev, grid)

Return the `vars` for unforced two-dimensional Navier-Stokes problem on device `dev` and
with `grid`.
"""
function DecayingVars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) ζ u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) ζh uh vh

  return Vars(ζ, u, v, ζh, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Return the `vars` for forced two-dimensional Navier-Stokes on device `dev` and with `grid`.
"""
function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) ζ u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) ζh uh vh Fh

  return Vars(ζ, u, v, ζh, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Return the `vars` for stochastically forced two-dimensional Navier-Stokes on device `dev` and
with `grid`.
"""
function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny) ζ u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) ζh uh vh Fh prevsol

  return Vars(ζ, u, v, ζh, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, ζ)`` in conservative form,
i.e., ``- ∂_x[(∂_y ψ)ζ] - ∂_y[(∂_x ψ)ζ]`` and store it in `N`:

```math
N = - \\widehat{𝖩(ψ, ζ)} = - i k_x \\widehat{u ζ} - i k_y \\widehat{v ζ} .
```
"""
function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.uh =   im * grid.l  * grid.invKrsq * sol
  @. vars.vh = - im * grid.kr * grid.invKrsq * sol
  @. vars.ζh = sol

  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)
  ldiv!(vars.ζ, grid.rfftplan, vars.ζh)

  uζ = vars.u                  # use vars.u as scratch variable
  @. uζ *= vars.ζ              # u*ζ
  vζ = vars.v                  # use vars.v as scratch variable
  @. vζ *= vars.ζ              # v*ζ

  uζh = vars.uh                # use vars.uh as scratch variable
  mul!(uζh, grid.rfftplan, uζ) # \hat{u*ζ}
  vζh = vars.vh                # use vars.vh as scratch variable
  mul!(vζh, grid.rfftplan, vζ) # \hat{v*ζ}

  @. N = - im * grid.kr * uζh - im * grid.l * vζh

  return nothing
end

"""
    calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term, that is the advection term and the forcing,

```math
N = - \\widehat{𝖩(ψ, ζ)} + F̂ .
```
"""
function calcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  calcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

"""
    addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term ``F̂`` and add it to the
nonlinear term ``N``.
"""
addforcing!(N, sol, t, clock, vars::DecayingVars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)

  @. N += vars.Fh

  return nothing
end
function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end
  @. N += vars.Fh

  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update variables in `vars` with solution in `sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  dealias!(sol, grid)

  @. vars.ζh = sol
  @. vars.uh =   im * grid.l  * grid.invKrsq * sol
  @. vars.vh = - im * grid.kr * grid.invKrsq * sol

  ldiv!(vars.ζ, grid.rfftplan, deepcopy(vars.ζh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh)) # deepcopy() since inverse real-fft destroys its input

  return nothing
end

"""
    set_ζ!(prob, ζ)

Set the solution `sol` as the transform of `ζ` and then update variables in `vars`.
"""
function set_ζ!(prob, ζ)
  mul!(prob.sol, prob.grid.rfftplan, ζ)

  CUDA.@allowscalar prob.sol[1, 1] = 0 # enforce zero domain average

  updatevars!(prob)

  return nothing
end

"""
    energy(prob)

Return the domain-averaged kinetic energy. Since ``u² + v² = |{\\bf ∇} ψ|²``, the domain-averaged
kinetic energy is

```math
\\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² .
```
"""
@inline function energy(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  energyh = vars.uh # use vars.uh as scratch variable

  @. energyh = 1 / 2 * grid.invKrsq * abs2(sol)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energyh, grid)
end

"""
    enstrophy(prob)

Returns the domain-averaged enstrophy,

```math
\\int \\frac1{2} ζ² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |ζ̂|² .
```
"""
@inline function enstrophy(prob)
  sol, grid = prob.sol, prob.grid
  return 1 / (2 * grid.Lx * grid.Ly) * parsevalsum(abs2.(sol), grid)
end

"""
    energy_dissipation(prob, ξ, νξ)

Return the domain-averaged energy dissipation rate done by the viscous term,

```math
- ξ (-1)^{n_ξ+1} \\int ψ ∇^{2n_ξ} ζ \\frac{𝖽x 𝖽y}{L_x L_y} = - ξ \\sum_{𝐤} |𝐤|^{2(n_ξ-1)} |ζ̂|² ,
```
where ``ξ`` and ``nξ`` could be either the (hyper)-viscosity coefficient ``ν`` and its order
``n_ν``, or the hypo-viscocity coefficient ``μ`` and its order ``n_μ``.
"""
@inline function energy_dissipation(prob, ξ, nξ)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  energy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. energy_dissipationh = - ξ * grid.Krsq^(nξ - 1) * abs2(sol)
  CUDA.@allowscalar energy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_dissipationh, grid)
end

"""
    energy_dissipation_hyperviscosity(prob, ξ, νξ)

Return the domain-averaged energy dissipation rate done by the ``ν`` (hyper)-viscosity.
"""
energy_dissipation_hyperviscosity(prob) = energy_dissipation(prob, prob.params.ν, prob.params.nν)

"""
    energy_dissipation_hypoviscosity(prob, ξ, νξ)

Return the domain-averaged energy dissipation rate done by the ``μ`` (hypo)-viscosity.
"""
energy_dissipation_hypoviscosity(prob) = energy_dissipation(prob, prob.params.μ, prob.params.nμ)

"""
    enstrophy_dissipation(prob, ξ, νξ)

Return the domain-averaged enstrophy dissipation rate done by the viscous term,

```math
ξ (-1)^{n_ξ+1} \\int ζ ∇^{2n_ξ} ζ \\frac{𝖽x 𝖽y}{L_x L_y} = - ξ \\sum_{𝐤} |𝐤|^{2n_ξ} |ζ̂|² ,
```

where ``ξ`` and ``nξ`` could be either the (hyper)-viscosity coefficient ``ν`` and its order
``n_ν``, or the hypo-viscocity coefficient ``μ`` and its order ``n_μ``.
"""
@inline function enstrophy_dissipation(prob, ξ, nξ)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  enstrophy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_dissipationh = - ξ * grid.Krsq^nξ * abs2(sol)
  CUDA.@allowscalar enstrophy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_dissipationh, grid)
end

"""
    enstrophy_dissipation_hyperviscosity(prob, ξ, νξ)

Return the domain-averaged enstrophy dissipation rate done by the ``ν`` (hyper)-viscosity.
"""
enstrophy_dissipation_hyperviscosity(prob) = enstrophy_dissipation(prob, prob.params.ν, prob.params.nν)

"""
    enstrophy_dissipation_hypoviscosity(prob, ξ, νξ)

Return the domain-averaged enstrophy dissipation rate done by the ``μ`` (hypo)-viscosity.
"""
enstrophy_dissipation_hypoviscosity(prob) = enstrophy_dissipation(prob, prob.params.μ, prob.params.nμ)

"""
    energy_work(prob)
    energy_work(sol, vars, grid)

Return the domain-averaged rate of work of energy by the forcing ``F``,

```math
- \\int ψ F \\frac{𝖽x 𝖽y}{L_x L_y} = - \\sum_{𝐤} ψ̂ F̂^* .
```
"""
@inline function energy_work(sol, vars::ForcedVars, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid)
end

@inline function energy_work(sol, vars::StochasticForcedVars, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid)
end

@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.grid)

"""
    enstrophy_work(prob)
    enstrophy_work(sol, vars, grid)

Return the domain-averaged rate of work of enstrophy by the forcing ``F``,

```math
\\int ζ F \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} ζ̂ F̂^* .
```
"""
@inline function enstrophy_work(sol, vars::ForcedVars, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline function enstrophy_work(sol, vars::StochasticForcedVars, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.grid)

end # module


text:
@testset "Alphabet" begin

    @testset "Creation, Iteration and getindex" begin

        # Creation & Iteration
        @test [ i for i in UngappedAlphabet() ] == Int[ i for i in 1:20 ]
        @test [ i for i in GappedAlphabet() ]   == Int[ i for i in 1:21 ]
        @test [ i for i in ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP") ]  ==
           Int[ i for i in 1:8 ]
    end

    @testset "getindex and names" begin

        reduced = ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP")
        strings = ["AILMV", "NQST", "RHK", "DE", "FWY", "C", "G", "P"]

        @test length(reduced) == length(strings)
        @test names(reduced) == strings

        for i in 1:length(reduced)
            @test reduced[strings[i]] == i
        end

        strings_gapped = ["A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S",
            "T","W","Y","V","-"]

        @test names(GappedAlphabet()) == strings_gapped
        @test names(UngappedAlphabet()) == strings_gapped[1:end-1]

        for i in 1:20
            @test UngappedAlphabet()[strings_gapped[i]] == i
            @test GappedAlphabet()[strings_gapped[i]] == i
        end

        @test GappedAlphabet()["-"] == 21
        @test UngappedAlphabet()["-"] == 22
    end

    @testset "in" begin

        # Creation & Iteration
        @test in(Residue('A'), UngappedAlphabet())
        @test in(Residue('A'), GappedAlphabet())
        @test in(Residue('A'), ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

        @test !in(GAP, UngappedAlphabet())
        @test  in(GAP, GappedAlphabet())
        @test !in(GAP, ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

        @test !in(XAA, UngappedAlphabet())
        @test !in(XAA, GappedAlphabet())
        @test !in(XAA, ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
    end

    @testset "Show" begin

        tmp = IOBuffer()
        show(tmp, UngappedAlphabet())
        @test String(take!(tmp)) ==
        "UngappedAlphabet of length 20. Residues : res\"ARNDCQEGHILKMFPSTWYV\""

        show(tmp, GappedAlphabet())
        @test String(take!(tmp)) ==
        "GappedAlphabet of length 21. Residues : res\"ARNDCQEGHILKMFPSTWYV-\""

        show(tmp, ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
        @test String(take!(tmp)) ==
        "ReducedAlphabet of length 8 : \"(AILMV)(NQST)(RHK)(DE)(FWY)CGP\""
    end
end


text:
using SparsityOperators

# test FourierOperators
function testFT(N=16)
  # random image
  x = zeros(ComplexF64,N,N)
  for i=1:N,j=1:N
    x[i,j] = rand()
  end

  # FourierMatrix
  idx = CartesianIndices((N,N))[collect(1:N^2)]
  F = [ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1))/N) for j=1:N^2, k=1:N^2 ]
  F_adj = F'

  # Operators
  tr = CartesianTrajectory(Float64,N,N)
  F_nfft = NFFTOp((N,N),tr,symmetrize=false)
  F_exp = ExplicitOp((N,N),tr,zeros(ComplexF64,N,N),symmetrize=false)

  # test agains FourierOperators
  y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N)) )
  y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N)) )

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  y_exp = F_exp*vec(x)
  y_adj_exp = adjoint(F_exp) * vec(x)

  @test y     ≈ y_nfft      rtol = 1e-2
  @test y     ≈ y_exp       rtol = 1e-2
  @test y_adj ≈ y_adj_nfft  rtol = 1e-2
  @test y_adj ≈ y_adj_exp   rtol = 1e-2

  # test AHA w/o Toeplitz
  F_nfft.toeplitz = false
  AHA = SparsityOperators.normalOperator(F_nfft)
  y_AHA_nfft = AHA * vec(x)
  y_AHA = F' * F * vec(x)
  @test y_AHA ≈ y_AHA_nfft   rtol = 1e-2

  # test AHA with Toeplitz
  F_nfft.toeplitz = true
  AHA = SparsityOperators.normalOperator(F_nfft)
  y_AHA_nfft = AHA * vec(x)
  y_AHA_nfft = adjoint(F_nfft) * F_nfft * vec(x)
  y_AHA = F' * F * vec(x)
  @test y_AHA ≈ y_AHA_nfft   rtol = 1e-2

  # test type stability;
  # TODO: Ensure type stability for Trajectory objects and test here
  nodes = Float32.(tr.nodes)
  F_nfft = NFFTOp((N,N),nodes,symmetrize=false)

  y_nfft = F_nfft * vec(ComplexF32.(x))
  y_adj_nfft = adjoint(F_nfft) * vec(ComplexF32.(x))

  @test Complex{eltype(nodes)} === eltype(y_nfft)
  @test Complex{eltype(nodes)} === eltype(y_adj_nfft)
end

function testFT3d(N=12)
  # random image
  x = zeros(ComplexF64,N,N,N)
  for i=1:N,j=1:N,k=1:N
    x[i,j,k] = rand()
  end

  # FourierMatrix
  idx = CartesianIndices((N,N,N))[collect(1:N^3)]
  F = [ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1)+(idx[j][3]-1)*(idx[k][3]-1))/N) for j=1:N^3, k=1:N^3 ]
  F_adj = F'

  # Operators
  tr = CartesianTrajectory3D(Float64,N,N,numSlices=N)
  F_nfft = NFFTOp((N,N,N),tr,symmetrize=false)
  F_exp = ExplicitOp((N,N,N),tr,zeros(ComplexF64,N,N,N),symmetrize=false)

  # test agains FourierOperators
  y = vec( ifftshift(reshape(F*vec(fftshift(x)),N,N,N)) )
  y_adj = vec( ifftshift(reshape(F_adj*vec(fftshift(x)),N,N,N)) )

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  y_exp = F_exp*vec(x)
  y_adj_exp = adjoint(F_exp) * vec(x)

  @test  y     ≈ y_exp      rtol = 1e-2
  @test  y     ≈ y_nfft     rtol = 1e-2
  @test  y_adj ≈ y_adj_exp  rtol = 1e-2
  @test  y_adj ≈ y_adj_nfft rtol = 1e-2

  # test AHA w/o Toeplitz
  F_nfft.toeplitz = false
  AHA = SparsityOperators.normalOperator(F_nfft)
  y_AHA_nfft = AHA * vec(x)
  y_AHA = F' * F * vec(x)
  @test y_AHA ≈ y_AHA_nfft   rtol = 1e-2

  # test AHA with Toeplitz
  F_nfft.toeplitz = true
  AHA = SparsityOperators.normalOperator(F_nfft)
  y_AHA_nfft = AHA * vec(x)
  y_AHA_nfft = adjoint(F_nfft) * F_nfft * vec(x)
  y_AHA = F' * F * vec(x)
  @test y_AHA ≈ y_AHA_nfft   rtol = 1e-2
end

## test FieldmapNFFTOp
function testFieldmapFT(N=16)
  # random image
  x = zeros(ComplexF64,N,N)
  for i=1:N,j=1:N
    x[i,j] = rand()
  end

  tr = CartesianTrajectory(Float64,N,N;TE=0.0,AQ=0.01)
  times = readoutTimes(tr)
  nodes = kspaceNodes(tr)
  cmap = im*quadraticFieldmap(N,N)[:,:,1]

  # FourierMatrix
  idx = CartesianIndices((N,N))[collect(1:N^2)]
  F = [exp(-2*1im*pi*(nodes[1,k]*(idx[l][1]-size(x,1)/2-1)+nodes[2,k]*(idx[l][2]-size(x,2)/2-1))-cmap[idx[l][1],idx[l][2]]*times[k]) for k=1:size(nodes,2), l=1:length(x)]
  F_adj = F'

  # Operators
  F_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false)

  # test agains FourierOperators
  y = F*vec(x)
  y_adj = F_adj*vec(x)

  y_nfft = F_nfft*vec(x)
  y_adj_nfft = adjoint(F_nfft) * vec(x)

  @test (norm(y-y_nfft)/norm(y)) < 1e-2
  @test (norm(y_adj-y_adj_nfft)/norm(y_adj)) < 1e-2
end

function testOperators()
  @testset "Linear Operator" begin
    testFT()
    testFT3d()
    testFieldmapFT()
  end
end

testOperators()


text:
#!/usr/bin/env julia

using Eigenbroetler

squarified_sinusoid(x::Number, y::Number) =
    exp(im * (x ^ 4 + y ^ 4) ^ (1 / 4))

eb = Eigenbrot(squarified_sinusoid, 512, 512)
save("squarified_sinusoid.png", eb, ImageSetting(RealPart, LinearScale))
save("squarified_sinusoid_dft.png", fft(eb), ImageSetting(Magn, LogScale))


text:
struct HCCache{F1<:AbstractMatrix, F2<:AbstractMatrix, V<:AbstractVector} <: AbstractCache
    q::F1
    X::F1
    x::F2
    v::V
    w::V
    η::V
    u::V
end

function HCCache(X::AbstractMatrix{T1}) where T1<:Real
    n, p = size(X)
    HCCache(similar(X), X, Array{T1, 2}(undef, p, p),
             Array{T1, 1}(undef, n), Array{T1, 1}(undef, n),
             Array{T1, 1}(undef, n), Array{T1, 1}(undef, n))
end

function emptyHCCache()
    V = Array{WFLOAT, 1}(undef, 0)
    M = Array{WFLOAT, 2}(undef, 0, 0)
    HCCache(M, M, M, V, V, V, V)
end


text:
using KernelAbstractions
using KernelAbstractions.NDIteration
using Test

function nditeration_testsuite()
@testset "iteration" begin
    let ndrange = NDRange{2, DynamicSize, DynamicSize}(CartesianIndices((256, 256)), CartesianIndices((32, 32)));
        @test length(ndrange) == 256*256
        @test all(p->p[1]==p[2], zip(ndrange, CartesianIndices((256, 256))))
    end
    let ndrange = NDRange{2, StaticSize{(256,256)}, DynamicSize}(nothing, CartesianIndices((32, 32)));
        @test length(ndrange) == 256*256
        @test all(p->p[1]==p[2], zip(ndrange, CartesianIndices((256, 256))))
    end
end

# GPU scenario where we get a linear index into workitems/blocks
function linear_iteration(ndrange)
    idx = Array{CartesianIndex{2}}(undef, length(blocks(ndrange)) * length(workitems(ndrange)))
    for i in 1:length(blocks(ndrange))
        for j in 1:length(workitems(ndrange))
            I = j + (i-1) * length(workitems(ndrange))
            idx[I] = expand(ndrange, i, j)
        end
    end
    return idx
end

function check(idx, offset, offset_x, offset_y, Dim_x, Dim_y)
    N = Dim_x * Dim_y
    all(p->p[1]==p[2], zip(idx[offset*N  .+ 1:N], CartesianIndices((offset_x*Dim_x .+ 1:Dim_x, offset_y*Dim_y .+ 1:Dim_y))))
end

@testset "linear_iteration" begin
    Dim_x = 32
    Dim_y = 32
    let ndrange = NDRange{2, StaticSize{(4,4)}, StaticSize{(Dim_x,Dim_y)}}();
        idx = linear_iteration(ndrange)
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end
    let ndrange = NDRange{2, DynamicSize, DynamicSize}(CartesianIndices((4,4)), CartesianIndices((Dim_x, Dim_y)));
        idx = linear_iteration(ndrange)
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end

    Dim_x = 32
    Dim_y = 1
    let ndrange = NDRange{2, StaticSize{(4,4*32)}, StaticSize{(Dim_x,Dim_y)}}();
        idx = linear_iteration(ndrange)
        N = length(workitems(ndrange))
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end
    let ndrange = NDRange{2, DynamicSize, DynamicSize}(CartesianIndices((4,4*32)), CartesianIndices((Dim_x, Dim_y)));
        idx = linear_iteration(ndrange)
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end

    Dim_x = 1
    Dim_y = 32
    let ndrange = NDRange{2, StaticSize{(4*32,4)}, StaticSize{(Dim_x,Dim_y)}}();
        idx = linear_iteration(ndrange)
        N = length(workitems(ndrange))
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end
    let ndrange = NDRange{2, DynamicSize, DynamicSize}(CartesianIndices((4*32,4)), CartesianIndices((Dim_x, Dim_y)));
        idx = linear_iteration(ndrange)
        for (i, I) in zip(1:length(blocks(ndrange)), blocks(ndrange))
            I = Tuple(I)
            @test check(idx, i-1, ntuple(i->I[i]-1, length(I))..., Dim_x, Dim_y)
        end
    end
end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1NodeStatus <: SwaggerModel
    addresses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NodeAddress} } # spec name: addresses
    allocatable::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: allocatable
    capacity::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: capacity
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NodeCondition} } # spec name: conditions
    config::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigStatus } # spec name: config
    daemonEndpoints::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeDaemonEndpoints } # spec name: daemonEndpoints
    images::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerImage} } # spec name: images
    nodeInfo::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSystemInfo } # spec name: nodeInfo
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase
    volumesAttached::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1AttachedVolume} } # spec name: volumesAttached
    volumesInUse::Any # spec type: Union{ Nothing, Vector{String} } # spec name: volumesInUse

    function IoK8sApiCoreV1NodeStatus(;addresses=nothing, allocatable=nothing, capacity=nothing, conditions=nothing, config=nothing, daemonEndpoints=nothing, images=nothing, nodeInfo=nothing, phase=nothing, volumesAttached=nothing, volumesInUse=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("addresses"), addresses)
        setfield!(o, Symbol("addresses"), addresses)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("allocatable"), allocatable)
        setfield!(o, Symbol("allocatable"), allocatable)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("config"), config)
        setfield!(o, Symbol("config"), config)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("daemonEndpoints"), daemonEndpoints)
        setfield!(o, Symbol("daemonEndpoints"), daemonEndpoints)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("images"), images)
        setfield!(o, Symbol("images"), images)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("nodeInfo"), nodeInfo)
        setfield!(o, Symbol("nodeInfo"), nodeInfo)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("volumesAttached"), volumesAttached)
        setfield!(o, Symbol("volumesAttached"), volumesAttached)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("volumesInUse"), volumesInUse)
        setfield!(o, Symbol("volumesInUse"), volumesInUse)
        o
    end
end # type IoK8sApiCoreV1NodeStatus

const _property_map_IoK8sApiCoreV1NodeStatus = Dict{Symbol,Symbol}(Symbol("addresses")=>Symbol("addresses"), Symbol("allocatable")=>Symbol("allocatable"), Symbol("capacity")=>Symbol("capacity"), Symbol("conditions")=>Symbol("conditions"), Symbol("config")=>Symbol("config"), Symbol("daemonEndpoints")=>Symbol("daemonEndpoints"), Symbol("images")=>Symbol("images"), Symbol("nodeInfo")=>Symbol("nodeInfo"), Symbol("phase")=>Symbol("phase"), Symbol("volumesAttached")=>Symbol("volumesAttached"), Symbol("volumesInUse")=>Symbol("volumesInUse"))
const _property_types_IoK8sApiCoreV1NodeStatus = Dict{Symbol,String}(Symbol("addresses")=>"Vector{IoK8sApiCoreV1NodeAddress}", Symbol("allocatable")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("capacity")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("conditions")=>"Vector{IoK8sApiCoreV1NodeCondition}", Symbol("config")=>"IoK8sApiCoreV1NodeConfigStatus", Symbol("daemonEndpoints")=>"IoK8sApiCoreV1NodeDaemonEndpoints", Symbol("images")=>"Vector{IoK8sApiCoreV1ContainerImage}", Symbol("nodeInfo")=>"IoK8sApiCoreV1NodeSystemInfo", Symbol("phase")=>"String", Symbol("volumesAttached")=>"Vector{IoK8sApiCoreV1AttachedVolume}", Symbol("volumesInUse")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeStatus }) = collect(keys(_property_map_IoK8sApiCoreV1NodeStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeStatus[property_name]

function check_required(o::IoK8sApiCoreV1NodeStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NodeStatus }, name::Symbol, val)
end


text:
abstract type BasisFunction end

#1s Gaussian-type function
abstract type Gaussian1s <: BasisFunction end
    #Gaussian orbital exponent
    alpha :: Float64
    #nuclear coordinates
    center::Float64
end

#Slater Type Orbital fit with N primative gausians (STO-NG) type basis
type STONG <: BasisFunction
    n::Integer
    #contraction coeffiecents
    d::Array{Float64}
    #primative gaussians
    g::Array{Gaussian1s}
end

#STO-3G basis for hydrogen
function sto3g_hydrogen(center)
    sto3g(center, 1.24)
end

#STO-3G basis for helium
function sto3g_helium(center)
    sto3g(center, 2.0925)
end

#Builds a STO-3G basis that best approximates a single slater type
#orbital with Slater orbital exponent zeta
function sto3g(center, zeta)
    scaling = zeta^2
    STONG(3,[0.444635, 0.535328, 0.154329],
            [Gaussian1s(scaling*.109818, center),
             Gaussian1s(scaling*.405771, center),
             Gaussian1s(scaling*2.22766, center)])
end

#The overlap integrals describe how the basis functions overlap
#as the atom centered gaussian basis functions are non-orthognal
#they have a non-zero overlap. The integral has the following form:
#S_{ij} = \int \phi_i(r-R_a) \phi_j(r-R_b) \mathrm{d}r
function overlap_integral(b1::STONG, b2::STONG)
    two_center_contraction(b1, b2, overlap_integral)
end

#This function calculates the overlap integral for 1s Gaussian
#orbitals using a closed-form expression.
function overlap_integral(g1::Gaussian1s, g2::Gaussian1s)
    alpha = g1.alpha
    beta = g2.alpha
    Ra = g1.center
    Rb = g2.center

    #normalization constant
    n = (2*alpha/pi)^(3/4) * (2*beta/pi)^(3/4)

    I  = n * (pi/(alpha+beta))^(3/2)
    I *= exp(-alpha*beta/(alpha+beta) * abs(Ra-Rb)^2)

    return I
end

#
function nuclear_attraction_integral(Zc::Int, Rc::Float64, b1::STONG,
                                     b2::STONG)
    integral(g1, g2) = nuclear_attraction_integral(Zc, Rc, g1, g2)
    two_center_contraction(b1, b2, integral)
end

function nuclear_attraction_integral(Zc::Int, Rc::Float64, g1::Gaussian1s, g2::Gaussian1s)
    alpha = g1.alpha
    beta  = g2.alpha
    Ra = g1.center
    Rb = g2.center
    Rp = (alpha*Ra + beta*Rb)/(alpha + beta)

    n = (2*alpha/pi)^(3/4) * (2*beta/pi)^(3/4)
    matrix_element  = n*-2*pi/(alpha+beta)*Zc
    matrix_element *= exp(-alpha*beta/(alpha+beta)*abs(Ra-Rb)^2)

    t = (alpha+beta)*abs(Rp-Rc)^2
    if abs(t) < 1e-8
        return matrix_element
    end

    matrix_element *= 0.5 * sqrt(pi/t) * erf(sqrt(t))
    return matrix_element
end

function kinetic_energy_integral(b1::STONG, b2::STONG)
    two_center_contraction(b1, b2, kinetic_energy_integral)
end

function kinetic_energy_integral(g1::Gaussian1s, g2::Gaussian1s)
    alpha = g1.alpha
    beta = g2.alpha
    Ra = g1.center
    Rb = g2.center

    n = (2*alpha/pi)^(3/4) * (2*beta/pi)^(3/4)

    matrix_element  = n * alpha*beta/(alpha+beta)
    matrix_element *= (3-2*alpha*beta/((alpha+beta)/abs(Ra-Rb)^2 ))
    matrix_element *= (pi/(alpha+beta))^(3/2)
    matrix_element *= exp(-alpha*beta/(alpha+beta) * abs(Ra-Rb)^2)
end

function two_electron_integral(g1::STONG, g2::STONG, g3::STONG, g4::STONG)
    four_center_contraction(g1, g2, g3, g4, two_electron_integral)
end

function two_electron_integral(g1::Gaussian1s, g2::Gaussian1s, g3::Gaussian1s,
                               g4::Gaussian1s)
    alpha = g1.alpha
    beta  = g2.alpha
    gamma = g3.alpha
    delta = g4.alpha
    Ra = g1.center
    Rb = g2.center
    Rc = g3.center
    Rd = g4.center
    Rp = (alpha*Ra + beta*Rb)/(alpha + beta)
    Rq = (gamma*Rc + delta*Rd)/(gamma + delta)

    n  = (2*alpha/pi)^(3/4) * (2*beta/pi)^(3/4)
    n *= (2*gamma/pi)^(3/4) * (2*delta/pi)^(3/4)

    matrix_element  = n*2*pi^(5/2)
    matrix_element /= ((alpha+beta)*(gamma+delta)*sqrt(alpha+beta+gamma+delta))
    matrix_element *= exp(-alpha*beta/(alpha+beta)*abs(Ra-Rb)^2 - gamma*delta/(gamma+delta)*abs(Rc-Rd)^2)
    t = (alpha+beta)*(gamma+delta)/(alpha+beta+gamma+delta)*abs(Rp-Rq)^2
    if abs(t) < 1e-8
        return matrix_element
    end

    matrix_element *= 0.5 * sqrt(pi/t) * erf(sqrt(t))
    return matrix_element
end

function two_center_contraction(b1::STONG, b2::STONG, integral::Function)
    total = 0.0
    for p = 1:b1.n
        for q = 1:b2.n
            d1 = b1.d[p]
            d2 = b2.d[q]
            total += d1*d2*integral(b1.g[p], b2.g[q])
        end
    end
    return total
end

function four_center_contraction(b1::STONG, b2::STONG, b3::STONG, b4::STONG,
                                 integral::Function)
    total = 0.0
    for p = 1:b1.n
        for q = 1:b2.n
            for r = 1:b3.n
                for s = 1:b4.n
                    dp = b1.d[p]
                    dq = b2.d[q]
                    dr = b3.d[r]
                    ds = b4.d[s]
                    total += dp*dq*dr*ds*integral(b1.g[p], b2.g[q], b3.g[r], b4.g[s])
                end
            end
        end
    end
    total
 end

function hartree_fock(R, Z)
    #println("constructing basis set")
    phi = Array(BasisFunction, length(Z))
    for A = 1:length(Z)
        if Z[A] == 1
            phi[A] = sto3g_hydrogen(R[A])
        elseif Z[A] == 2
            phi[A] = sto3g_helium(R[A])
        end
    end
    #calculate the overlap matrix S
    #the matrix should be symmetric with diagonal entries equal to one
    #println("building overlap matrix")
    S = eye(length(phi))
    for i = 1:length(phi)
        for j = (i+1):length(phi)
            S[i,j] = S[j,i] = overlap_integral(phi[i], phi[j])
        end
    end

    #println("S: ", S)

    #calculate the kinetic energy matrix T
    #println("building kinetic energy matrix")
    T = zeros(length(phi), length(phi))
    for i = 1:length(phi)
        for j = i:length(phi)
            T[i,j] = T[j,i] = kinetic_energy_integral(phi[i], phi[j])
        end
    end

    #println("T: ", T)

    #calculate nuclear attraction matrices V_i
    #println("building nuclear attraction matrices")
    V = zeros(length(Z), length(Z))
    for A = 1:length(Z)
        for i = 1:length(phi)
            for j = i:length(phi)
                v = nuclear_attraction_integral(Z[A], R[A], phi[i], phi[j])
                V[i,j] += v
                if i != j
                    V[j,i] += v
                end
            end
        end
    end

    #println("V: ", V)

    #build core-Hamiltonian matrix
    #println("building core-Hamiltonian matrix")
    Hcore = T + V

    #println("Hcore: ", Hcore)

    #diagonalize overlap matrix to get transformation matrix X
    #println("diagonalizing overlap matrix")
    s, U = eig(S)
    #println("building transformation matrix")
    X = U*diagm(s.^(-1/2))*U'
    #println("X: ", X)

    #calculate all of the two-electron integrals
    K = length(phi)
    two_electron = zeros(K,K,K,K)
    for mu = 1:K
        for v = 1:K
            for lambda = 1:K
                for sigma = 1:K
                    coulomb  = two_electron_integral(phi[mu], phi[v],
                                                     phi[sigma], phi[lambda])
                    two_electron[mu,v,sigma,lambda] = coulomb
                    #println("coulomb  ($mu $v | $sigma $lambda): $coulomb")
                    exchange = two_electron_integral(phi[mu], phi[lambda],
                                                     phi[sigma], phi[v])
                    #println("exchange ($mu $lambda | $sigma $v): $exchange")
                    two_electron[mu,lambda,sigma,v] = exchange
                end
            end
        end
    end

    P = zeros(K,K)

    total_energy = 0.0
    old_energy = 0.0
    electronic_energy = 0.0

    printf("%4s %13s de\n", "iter", "total energy")
    for scf_iter = 1:100
        #calculate the two electron part of the Fock matrix
        G = zeros(size(Hcore))
        K = length(phi)
        for mu = 1:K
            for v = 1:K
                for lambda = 1:K
                    for sigma = 1:K
                        coulomb  = two_electron[mu,v,sigma,lambda]
                        exchange = two_electron[mu,lambda,sigma,v]
                        G[mu,v] += P[lambda,sigma]*(coulomb - 0.5*exchange)
                    end
                end
            end
        end

        #println("G: ", G)
        F = Hcore + G

        nuclear_energy = 0.0
        for A = 1:length(Z)
            for B = (A+1):length(Z)
                nuclear_energy += Z[A]*Z[B]/abs(R[A]-R[B])
            end
        end
        #println("E_nclr: $nuclear_energy")

        electronic_energy = 0.0
        for mu = 1:length(phi)
            for v = 1:length(phi)
                electronic_energy += P[v,mu]*(Hcore[mu,v]+F[mu,v])
            end
        end
        electronic_energy *= 0.5
        #println("E_elec: $electronic_energy")
        total_energy = electronic_energy + nuclear_energy
        printf("%3i %12.8f %12.4e\n", scf_iter, total_energy,
               total_energy - old_energy)

        if scf_iter > 2 && abs(old_energy - total_energy) < 1e-6
            break
        end

        #println("F: ", F)
        Fprime = X' * F * X
        #println("F': $Fprime")
        epsilon, Cprime = eig(Fprime)
        #println("epsilon: ", epsilon)
        #println("C': ", Cprime)
        C = real(X*Cprime)
        #println("C: ", C)

        P = zeros(size(Hcore))
        for mu = 1:length(phi)
            for v = 1:length(phi)
                P[mu,v] = 2*C[mu,1]*C[v,1]
            end
        end
        #println("P: ", P)

        old_energy = total_energy
    end

    return total_energy, electronic_energy
end

function test_h2()
    println("TESTING H2")
    total_energy, electronic_energy = hartree_fock([0., 1.4], [1, 1])
    szabo_energy = -1.8310
    if abs(electronic_energy - szabo_energy) > 1e-6
        println("TEST FAILED")
    else
        println("TEST PASSED")
    end

end

function test_heh()
    println("TESTING HEH+")
    total_energy, electronic_energy = hartree_fock([0., 1.4632], [2, 1])
    szabo_energy = -4.227529
    if abs(electronic_energy - szabo_energy) > 1e-6
        println("TEST FAILED")
    else
        println("TEST PASSED")
    end

end

function heh_pes()
    file = open("heh_pes.dat", "w")
    energy_he, energy_he = hartree_fock([0.0], [2])
    for r in linspace(0.7, 3.5, 25)
        total_energy, electronic_energy = hartree_fock([0., r], [2, 1])
        write(file, "$r $(total_energy - energy_he)\n")
    end
    close(file)
end

function tests()
    test_h2()
    test_heh()
end

tests()
heh_pes()


text:
# Executable to search PubMed for articles to a related term
# Date: May 6, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using ArgParse
using PubMedMiner
using SQLite
using BioMedQuery.UMLS: Credentials
using BioMedQuery.Entrez

function main(args)

    # initialize the settings (the description is for the help screen)
   s = ArgParseSettings()
   s.description = "This is pubMedSearch.jl"
   s.version = "1.0"

   @add_arg_table s begin
    "--clean_db"
        help = "Flag indicating wether to empty database"
        action = :store_true
    "--email"
         help = "your email address - required by NCBI"
         arg_type = ASCIIString
         required = true
    "--max_articles"
         help = "Maximum number of articles"
         arg_type = Int
         default = typemax(Int64)
    "--search_term"
         help = "Term to search"
         default = "obesity"
    "--verbose"
         help = "Verbose. Store temp files with server response"
         action = :store_true
    "mysql"
        help = "Use MySQL backend"
        action = :command
    "sqlite"
         help = "Use SQLite backend"
         action = :command
   end

   @add_arg_table s["mysql"] begin
    "--host"
        help = "Host where you database lives"
        arg_type = ASCIIString
        default = "localhost"
    "--dbname"
        help = "Database name"
        arg_type = ASCIIString
        required = true
    "--username"
        help = "MySQL username"
        arg_type = ASCIIString
        default = "root"
    "--password"
        help = "MySQL password"
        arg_type = ASCIIString
        default = ""
   end

    @add_arg_table s["sqlite"] begin
    "--db_path"
         help = "Path to SQLite database file to store results"
         arg_type = ASCIIString
         required = true
    end

   parsed_args = parse_args(s) # the result is a Dict{String,Any}
   println("-------------------------------------------------------------")
   println(s.description)
   println("Parsed args:")
   for (key,val) in parsed_args
       println("  $key  =>  $(repr(val))")
   end
   println("-------------------------------------------------------------")

   email = parsed_args["email"]
   search_term = parsed_args["search_term"]
   max_articles = parsed_args["max_articles"]
   verbose= parsed_args["verbose"]
   db_config = Dict()
   save_func = nothing

   if haskey(parsed_args, "mysql")
       db_config = Dict(:host=>parsed_args["mysql"]["host"],
                        :dbname=>parsed_args["mysql"]["dbname"],
                        :username=>parsed_args["mysql"]["username"],
                        :pswd=>parsed_args["mysql"]["password"],
                        :overwrite=>parsed_args["clean_db"])
       save_func = save_efetch_mysql
   elseif haskey(parsed_args, "sqlite")
      db_config = Dict(:db_path=>parsed_args["sqlite"]["db_path"],
                       :overwrite=>parsed_args["clean_db"])
      save_func = save_efetch_sqlite
    else
      error("Unsupported database backend")
  end

   @time begin
       db = PubMedMiner.pubmed_search(email, search_term, max_articles,
       save_func, db_config, parsed_args["verbose"])
   end

end

main(ARGS)


text:
import TrajectoryOptimization

urdf_folder = joinpath(@__DIR__, "urdfs")
urdf_doublecartpole = joinpath(urdf_folder, "doublecartpole.urdf")

doublecartpole_model = TrajectoryOptimization.Model(urdf_doublecartpole, [1., 0., 0.])
doublecartpole_model_full = TrajectoryOptimization.Model(urdf_doublecartpole, [1., 1., 1.])


text:
using Butterfly
using CSV
using Random

Random.seed!(777)

simdata = generate_simdata(100)
CSV.write(joinpath(@__DIR__, "benchmark_data.csv"), simdata)

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
struct Problem
    C::Vector{Int64} # Vecteur des pofits de chaque variables
    A::Union{Vector{Vector{Int64}}, Array{Int64, 2}} # Matrice creuse/Liste d'adgacence des contraintes

    #Donées annexes

    listOccur::Array{Int64, 2}
    occurVar :: Vector{Int64} # Nombre de variable que chaque variable penalise
    occurConst::Vector{Int64} #Nombre de fois qu'une variable apparait  dans une contraintes
end

struct Solution
    x::Vector{Int64}
    z::Int64

	# Données annexes

	var1::Vector{Int64} # Indice des variables misent à 1
	var0::Vector{Int64} # Indice des variables misent à 0
end

function Solution(x::Vector{Int64}, z::Int64)
	num1 = sum(x)
	var1 = Vector{Int64}(undef, num1)
	var0 = Vector{Int64}(undef, length(x)-num1)

	num0 = 0
	num1 = 0
	for ind = 1:length(x)
		if x[ind] == 0
			num0 += 1
			var0[num0] = ind
		else
			num1 += 1
			var1[num1] = ind
		end
	end
	return Solution(x, z, var1, var0)
end

function loadSPP(fname)
    f=open(fname)
    # lecture du nbre de contraintes (m) et de variables (n)
    m, n = parse.(Int, split(readline(f)) )
    # lecture des n coefficients de la fonction economique et cree le vecteur d'entiers c
    C = parse.(Int, split(readline(f)) )
    # lecture des m contraintes et reconstruction de la matrice binaire A
    #A=Vector{Vector{Int64}}(undef, m )
    # A = Array{Int64, 2}(undef, m, n)
    A = zeros(Int64, m, n)
    listOccur = zeros(Int64, n, n)
    occurVar = zeros(Int64, n)
    occurConst = zeros(Int64, n)

    for i=1:m
        # lecture du nombre d'elements non nuls sur la contrainte i (non utilise)
        nbVarConst = parse.(Int, split(readline(f)))
        # A[i] = Vector{Int64}(undef, nbVarConst[1])
        # lecture des indices des elements non nuls sur la contrainte i
        line = parse.(Int, split(readline(f)))
        # for ind in 1:nbVarConst
        #     A[i][ind] = line[ind]
        # end
        for var in line
            A[i, var] = 1
            occurConst[var] += 1
            for var2 in line
                listOccur[var, var2] = 1
                occurVar[var] += 1
            end
        end
    end
    close(f)
    return Problem(C, A, listOccur, occurVar, occurConst)
end

function calculUtil(var::Int64, problem::Problem, x::Vector{Int64})
	util = undef
	if x[var] != -1
		util = -Inf
	else
		prof::Int64 = problem.C[var]
		nbOccur::Int64 = problem.occurConst[var]
		util = prof/(nbOccur)^0.5
	end
	return util
end

function findVariable(x::Vector{Int64}, problem::Problem)::Int64
    var0 = -1
    utilMax = -Inf
    for var = 1:length(x)
        util = calculUtil(var, problem, x)
        if util > utilMax
            var0 = var
            utilMax = util
        end
    end
    return var0
end

function gloutonBasique(fname)
    problem = loadSPP(fname)
    nbVar::Int64 = length(problem.C)

    x::Vector{Int64} = zeros(Int64, nbVar)
    for ind = 1:nbVar
    	x[ind] = -1
    end
    z::Int64 = 0

	nbVariablesLibre = nbVar

    while nbVariablesLibre > 0
        x_i = findVariable(x, problem)
        x[x_i] = 1
        nbVariablesLibre -= 1
        z += problem.C[x_i]
        filterList = []
        for ind = 1:nbVar
            if problem.listOccur[x_i, ind] == 1 && x[ind] == -1
                x[ind] = 0
                nbVariablesLibre -= 1
            end
        end
    end
    return Solution(x, z)
end

function getfname(target)
    # target : string := chemin + nom du repertoire ou se trouve les instances

    # positionne le currentdirectory dans le repertoire cible
    cd(joinpath(homedir(),target))

    # recupere tous les fichiers se trouvant dans le repertoire data
    allfiles = readdir()

    # vecteur booleen qui marque les noms de fichiers valides
    flag = trues(size(allfiles))

    k=1
    for f in allfiles
        # traite chaque fichier du repertoire
        if f[1] == "."
            # fichier cache => supprimer
            flag[k] = false
        end
        k = k+1
    end

    # extrait les noms valides et retourne le vecteur correspondant
    finstances = allfiles[flag]
    return finstances
end

function main()
    gloutonBasique("Data/test.txt")
    fnames = getfname("/home/lucas/Bureau/Dm/Data")
	for fname in fnames
		print("File : "*fname*" -> ")
		sol = gloutonBasique(fname)
		println(sol.z)
	end
end

function exchangeK1(x::Vector{Int64}, z::Int64, var0::Vector{Int64}, var1::Vector{Int64}, prob::Problem, k::Int64)
	#A chaque appel de cette fonction p = 1
	varPossible = Vector{Int64}(undef, length(x))
	indPossible = 1

	permut::Vector{Tuple{Int64, Int64, Int64}} = Vector(undef, length(var0) * length(var1)*(length(var1)-1))
	indPermut = 1

	for ind0 in var0
		indPossible = 1
		for var = 1:length(x)
			if prob.listOccur[ind0,var] == 1
				varPossible[indPossible] = var
				indPossible += 1
			end
		end
		for ind = 1:indPossible
			for indBis = 1:indPossible
				if prob.listOccur[ind, indBis] == 0
					permut[indPermut] = [ind, indBis, ind0]
				end
			end
		end
	end
	return permut
end

function amelioration(sol::Solution, prob::Problem)
	chgt::Bool = true
	chgt21::Bool = true
	chgt11::Bool = true
	chgt10::Bool = true
	x, z, var0, var1 = sol.x, sol.z, sol.var0, sol.var1
	nbchgt21 = 0
	nbchgt11 = 0
	nbchgt10 = 0

	while chgt
		chgt = false
		# nbchgt21 = 0
		while chgt21
			chgt21 = true
			newX, newZ = exchangeKp(x, z, var0, var1, prob, 2,1)
			if z!=newZ
				x = newX
				z = newZ
			else
				chgt21 = false
			end
			# nbchgt21+=1
		end
		nbchgt11 = 0
		while chgt11
			chgt11 = true
			newX, newZ = exchangeKp(x, z, var0, var1, prob, 1, 1)
			if z!=newZ
				x = newX
				z = newZ
			else
				chgt11 = false
			end
			nbchgt11 +=1
		end
		nbchgt10 = 0
		while chgt10
			chgt10 = true
			newX, newZ = exchangeKp(x, z, var0, var1, prob, 1, 0)
			if z!=newZ
				x = newX
				z = newZ
			else
				chgt10 = false
			end
			nbchgt10 += 1
		end

		if nbchgt11 + nbchgt10 == 2
			chgt = false
		end
	end
	return Solution(x, z)
end


text:
using Simulacrum

# all to all
n = rand(collect(1:1000))
@test length(indices(AllToAll(n))) == n * (n-1) # n-1 hyperedges (excluding loops)
@test_throws ErrorException AllToAll(1)
@test_throws ErrorException AllToAll(-1)

# cycle
n = rand((1:1000))
@test length(indices(Cycle(n))) == n
@test_throws ErrorException Cycle(0)
@test isa(Cycle(4), Simulacrum.AbstractPattern)
@test indices(Cycle(1)) == [(1,1)] # a cycle of length 1 is a loop

# path
n = rand((1:1000))
@test length(indices(Path(n))) == n-1
@test_throws ErrorException Path(0)


text:
module Pilot

include("data.jl")
include("model/embedding.jl")

end # module Pilot


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real}
  # Make sure to balance the p4est before creating any containers
  # in case someone has tampered with the p4est after creating the mesh
  balance!(mesh)

  elements   = init_elements(mesh, equations, dg.basis, uEltype)
  interfaces = init_interfaces(mesh, equations, dg.basis, elements)
  boundaries = init_boundaries(mesh, equations, dg.basis, elements)
  mortars    = init_mortars(mesh, equations, dg.basis, elements)

  cache = (; elements, interfaces, boundaries, mortars)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

# Extract outward-pointing normal direction
# (contravariant vector ±Ja^i, i = index)
# Note that this vector is not normalized
@inline function get_normal_direction(direction, contravariant_vectors, indices...)

  orientation = (direction + 1) >> 1
  normal = get_contravariant_vector(orientation, contravariant_vectors, indices...)

  # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards
  if isodd(direction)
    return -normal
  else
    return normal
  end
end

include("containers.jl")
include("dg_2d.jl")
include("dg_3d.jl")

end # @muladd


text: include("dfsane.jl")


text:
function rk4step(p::VDPTagMDP, pos::Vec2)
    h = p.dt
    k1 = vdp_dynamics(p.mu, pos)
    k2 = vdp_dynamics(p.mu, pos+h/2*k1)
    k3 = vdp_dynamics(p.mu, pos+h/2*k2)
    k4 = vdp_dynamics(p.mu, pos+h*k3)
    return pos + h/6*(k1+2*k2+2*k3+k4)
end

function vdp_dynamics(mu::Float64, pos::Vec2)
    x = pos[1]
    y = pos[2]
    return SVector(mu*(x-x^3/3-y), x/mu)
end


text:
using LinQuadOptInterface

const LQOI = LinQuadOptInterface
const MOI  = LQOI.MathOptInterface

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear
    LQOI.SinVar
    LQOI.Quad
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecVar, MOI.SOS1{Float64}),
    (LQOI.VecVar, MOI.SOS2{Float64}),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros),
    (LQOI.Quad, LQOI.EQ),
    (LQOI.Quad, LQOI.LE),
    (LQOI.Quad, LQOI.GE)
]

mutable struct Optimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    env::Union{Nothing, Env}
    params::Dict{String, Any}

    """
        Optimizer(env = nothing; kwargs...)

    Create a new Optimizer object.

    You can share CPLEX `Env`s between models by passing an instance of `Env` as
    the first argument. By default, a new environment is created for every
    model.
    """
    function Optimizer(env::Union{Nothing, Env} = nothing; kwargs...)
        model = new()
        model.env = env
        # For consistency with MPB, output logs to stdout by default.
        model.params = Dict{String, Any}("CPX_PARAM_SCRIND" => 1)
        for (name, value) in kwargs
            model.params[string(name)] = value
        end
        MOI.empty!(model)
        return model
    end
end

# The existing env is `Nothing`, so create a new one.
LQOI.LinearQuadraticModel(::Type{Optimizer}, ::Nothing) = Model(Env())
# The existing env is `Env`, so pass it through.
LQOI.LinearQuadraticModel(::Type{Optimizer}, env::Env) = Model(env)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model, model.env)
    for (name, value) in model.params
        set_param!(model.inner.env, name, value)
    end
    return
end

MOI.get(::Optimizer, ::MOI.SolverName) = "CPLEX"

LQOI.supported_constraints(::Optimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(::Optimizer)  = SUPPORTED_OBJECTIVES

LQOI.backend_type(model::Optimizer, ::MOI.EqualTo{Float64})     = Cchar('E')
LQOI.backend_type(model::Optimizer, ::MOI.LessThan{Float64})    = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.GreaterThan{Float64}) = Cchar('G')

LQOI.backend_type(model::Optimizer, ::MOI.Zeros)                = Cchar('E')
LQOI.backend_type(model::Optimizer, ::MOI.Nonpositives)         = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.Nonnegatives)         = Cchar('G')

function LQOI.change_variable_bounds!(model::Optimizer,
        columns::Vector{Int}, values::Vector{Float64}, senses::Vector{Cchar})
    c_api_chgbds(model.inner, ivec(columns), senses, values)
    return
end

function LQOI.get_variable_lowerbound(model::Optimizer, column::Int)
    return c_api_getlb(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_variable_upperbound(model::Optimizer, column::Int)
    return c_api_getub(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_number_linear_constraints(model::Optimizer)
    return c_api_getnumrows(model.inner)
end

function LQOI.add_linear_constraints!(model::Optimizer,
        A::LQOI.CSRMatrix{Float64}, sense::Vector{Cchar}, rhs::Vector{Float64})
    c_api_addrows(model.inner, ivec(A.row_pointers), ivec(A.columns),
        A.coefficients, sense, rhs)
    return
end

function LQOI.get_rhs(model::Optimizer, row::Int)
    rhs = Vector{Cdouble}(undef, 1)
    c_api_getrhs(model.inner, rhs, Cint(row), Cint(row))
    return rhs[1]
end

function LQOI.get_linear_constraint(model::Optimizer, row::Int)
    (nzcnt, rmatbeg, rmatind, rmatval) =
            c_api_getrows(model.inner, Cint(row), Cint(row))
    return rmatind[1:nzcnt], rmatval[1:nzcnt]
end

function LQOI.change_matrix_coefficient!(
        model::Optimizer, row::Int, col::Int, coef::Float64)
    c_api_chgcoef(model.inner, Cint(row), Cint(col), coef)
    return
end

function LQOI.change_objective_coefficient!(
        model::Optimizer, col::Int, coef::Float64)
    c_api_chgobj(model.inner, [Cint(col)], [coef])
    return
end

function LQOI.change_rhs_coefficient!(
        model::Optimizer, row::Int, coef::Float64)
    c_api_chgrhs(model.inner, [Cint(row)], [coef])
    return
end

function LQOI.delete_linear_constraints!(
        model::Optimizer, first_row::Int, last_row::Int)
    c_api_delrows(model.inner, Cint(first_row), Cint(last_row))
    return
end

function LQOI.change_variable_types!(
        model::Optimizer, columns::Vector{Int}, vtypes::Vector{Cchar})
    c_api_chgctype(model.inner, ivec(columns), vtypes)
    return
end

function LQOI.change_linear_constraint_sense!(
        model::Optimizer, rows::Vector{Int}, senses::Vector{Cchar})
    c_api_chgsense(model.inner, ivec(rows), senses)
    return
end

function LQOI.set_linear_objective!(
        model::Optimizer, columns::Vector{Int}, coefficients::Vector{Float64})
    n = num_var(model.inner)
    all_coefs = zeros(Float64, n)
    for (col, coef) in zip(columns, coefficients)
        all_coefs[col] += coef
    end
    c_api_chgobj(model.inner, Cint[1:n;], all_coefs)
    return
end

function LQOI.change_objective_sense!(model::Optimizer, symbol)
    if symbol == :min
        c_api_chgobjsen(model.inner, Cint(1))
    else
        @assert symbol == :max
        c_api_chgobjsen(model.inner, Cint(-1))
    end
    return
end

function LQOI.get_linear_objective!(model::Optimizer, dest)
    c_api_getobj(model.inner, dest, Cint(1), c_api_getnumcols(model.inner))
    return
end

function LQOI.get_objectivesense(model::Optimizer)
    s = c_api_getobjsen(model.inner)
    if s == 1
        return MOI.MIN_SENSE
    else
        @assert s == -1
        return MOI.MAX_SENSE
    end
end

function LQOI.get_number_variables(model::Optimizer)
    return c_api_getnumcols(model.inner)
end

function LQOI.add_variables!(model::Optimizer, N::Int)
    add_vars!(model.inner, zeros(Float64, N), fill(-Inf, N), fill(Inf, N))
    return
end

function LQOI.delete_variables!(model::Optimizer, first_col::Int, last_col::Int)
    c_api_delcols(model.inner, Cint(first_col), Cint(last_col))
    return
end

function LQOI.solve_mip_problem!(model::Optimizer)
    model.inner.has_int = true
    LQOI.solve_linear_problem!(model)
    return
end

function LQOI.solve_linear_problem!(model::Optimizer)
    optimize!(model.inner)
    return
end

function LQOI.get_termination_status(model::Optimizer)
    stat = c_api_getstat(model.inner)

    if stat == 1 # CPX_STAT_OPTIMAL
        return MOI.OPTIMAL
    elseif stat == 3 # CPX_STAT_INFEASIBLE
        return MOI.INFEASIBLE
    elseif stat == 4 # CPX_STAT_INForUNBD
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif stat == 2 # CPX_STAT_UNBOUNDED
        return MOI.DUAL_INFEASIBLE
    elseif stat in (12, 21, 22, 36) # CPX_STAT_*ABORT*_OBJ_LIM
        return MOI.OBJECTIVE_LIMIT
    elseif stat in (10, 34) # CPX_STAT_*ABORT_IT_LIM
        return MOI.ITERATION_LIMIT
    elseif stat == 53 # CPX_STAT_CONFLICT_ABORT_NODE_LIM
        return MOI.NODE_LIMIT
    elseif stat in (11, 25, 33, 39) # CPX_STAT_*ABORT*TIME_LIM
        return MOI.TIME_LIMIT
    elseif stat == 5 # CPX_STAT_OPTIMAL_INFEAS
        return MOI.NUMERICAL_ERROR
    # MIP STATUS
    elseif stat in (101, 102) # CPXMIP_OPTIMAL, CPXMIP_OPTIMAL_TOL
        return MOI.OPTIMAL
    elseif stat == 103 # CPXMIP_INFEASIBLE
        return MOI.INFEASIBLE
    elseif stat == 119 # CPXMIP_INForUNBD
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif stat == 118 # CPXMIP_UNBOUNDED
        return MOI.DUAL_INFEASIBLE
    elseif stat in (105, 106) # CPXMIP_NODE_LIM*
        return MOI.NODE_LIMIT
    elseif stat in (107, 108, 131, 132) # CPXMIP_*TIME_LIM*
        return MOI.TIME_LIMIT
    else
        return MOI.OTHER_ERROR
    end
end

function LQOI.get_primal_status(model::Optimizer)
    soln_method, soln_type, primal_stat, dual_stat = c_api_solninfo(model.inner)
    if primal_stat == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function LQOI.get_dual_status(model::Optimizer)
    if model.inner.has_int
        return MOI.NO_SOLUTION
    end
    soln_method, soln_type, primal_stat, dual_stat = c_api_solninfo(model.inner)
    if dual_stat == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function LQOI.get_variable_primal_solution!(model::Optimizer, dest)
    c_api_getx(model.inner, dest)
    return
end

function LQOI.get_linear_primal_solution!(model::Optimizer, dest)
    c_api_getax(model.inner, dest)
    return
end

function LQOI.get_variable_dual_solution!(model::Optimizer, dest)
    c_api_getdj(model.inner, dest)
    return
end

function LQOI.get_linear_dual_solution!(model::Optimizer, dest)
    c_api_getpi(model.inner, dest)
    return
end

function LQOI.get_objective_value(model::Optimizer)
    return c_api_getobjval(model.inner)
end

function LQOI.get_objective_bound(model::Optimizer)
    return get_best_bound(model.inner)
end

function LQOI.get_farkas_dual!(model::Optimizer, dest)
    copy!(dest, get_infeasibility_ray(model.inner))
    return
end

function LQOI.get_unbounded_ray!(model::Optimizer, dest)
    copy!(dest, get_unbounded_ray(model.inner))
    return
end

function LQOI.add_sos_constraint!(model::Optimizer, columns::Vector{Int},
        weights::Vector{Float64}, sos_type::Symbol)
    add_sos!(model.inner, sos_type, columns, weights)
    return
end

function LQOI.get_sos_constraint(model::Optimizer, index::Int)
    (cols, weights, sos_type) = c_api_getsos(model.inner, index - 1)
    if sos_type == Cchar('1')
        return (cols .+ 1, weights, :SOS1)
    else
        @assert sos_type == Cchar('2')
        return (cols .+ 1, weights, :SOS2)
    end
    return
end

function LQOI.delete_sos!(model::Optimizer, begin_index::Int, end_index::Int)
    c_api_delsos(model.inner, begin_index - 1, end_index - 1)
    return
end

function scalediagonal!(V, I, J, scale)
    #  LQOI assumes 0.5 x' Q x, but CPLEX requires the list of terms, e.g.,
    #  2x^2 + xy + y^2, so we multiply the diagonal of V by 0.5. We don't
    #  multiply the off-diagonal terms since we assume they are symmetric and we
    #  only need to give one.
    #
    #  We also need to make sure that after adding the constraint we un-scale
    #  the vector because we can't modify user-data.
    for i in 1:length(I)
        if I[i] == J[i]
            V[i] *= scale
        end
    end
    return
end

function LQOI.set_quadratic_objective!(model::Optimizer, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})
    @assert length(I) == length(J) == length(V)
    CPLEX.add_qpterms!(model.inner, I, J, V)
    return
end

function LQOI.solve_quadratic_problem!(model::Optimizer)
    model.inner.has_qc = true
    LQOI.solve_linear_problem!(model)
end

function LQOI.get_quadratic_primal_solution!(model::Optimizer, dest)
    c_api_getxqxax(model.inner, dest)
    return
end

function LQOI.get_quadratic_dual_solution!(model::Optimizer, dest)
    c_api_getqconstrslack(model.inner, dest)
    return
end

function LQOI.add_quadratic_constraint!(model::Optimizer,
        affine_columns::Vector{Int}, affine_coefficients::Vector{Float64},
        rhs::Float64, sense::Cchar,
        I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})
    @assert length(I) == length(J) == length(V)
    scalediagonal!(V, I, J, 0.5)
    add_qconstr!(model.inner, Cint.(affine_columns), affine_coefficients,
                 Cint.(I), Cint.(J), V, sense, rhs)
    scalediagonal!(V, I, J, 2.0)
    return
end

function LQOI.get_quadratic_constraint(model::Optimizer, row::Int)
    affine_cols, affine_coefficients, I, J, V, _, _ = c_api_getqconstr(model.inner, row)
    for i in 1:length(I)
        if I[i] != J[i]
            V[i] *= 0.5  # Account for the 0.5 term in 0.5 x' Q x.
        end
    end
    # Convert 0-based indices into 1-based indices for the variables.
    return Int.(affine_cols .+ 1), affine_coefficients, sparse(I .+ 1, J .+ 1, V)
end

function LQOI.get_quadratic_rhs(model::Optimizer, row::Int)
    _, _, _, _, _, _, rhs = c_api_getqconstr(model.inner, row)
    return rhs
end

function LQOI.get_number_quadratic_constraints(model::Optimizer)
    return CPLEX.num_qconstr(model.inner)
end

 function LQOI.get_quadratic_terms_objective(model::Optimizer)
    qmatbeg, qmatind, qmatval = c_api_getquad(model.inner)
    qmatind .+= 1
    qmatcol = fill(length(qmatbeg), length(qmatind))
    # qmatbeg[i] stores the initial element (0-indexed) in qmatind and qmatval
    # for the i'th variable. In the next loop, we exclude the last variable
    # because it is implicitly set via the call to `fill`.
    for i in 1:length(qmatbeg) - 1
        start_index = qmatbeg[i] + 1  # +1 converts to 1-based.
        stop_index = qmatbeg[i + 1]
        for j in start_index:stop_index
        qmatcol[j] = i
        end
    end
    for i in 1:length(qmatval)
        if qmatind[i] != qmatcol[i]
            qmatval[i] *= 0.5  # Account for the 0.5 term in 0.5 x' Q x.
        end
    end
    # qmatind is ::Vector{Cint}, so we convert back to ::Vector{Int}.
    return sparse(Int.(qmatind), qmatcol, qmatval)
 end


text:
using MultivariateStats

# Simply calling Julia's PCA routine for testing purposes
function fit_pca_default(X::Array; max_num_components::Integer=5)
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  M = fit(PCA,X,maxoutdim=max_num_components,pratio=0.999999)
  num_components = size(projection(M),2)
  #Y = transform(M,X)
  #Xr = reconstruct(M,Y)
  #save("pca.jld", "num_lambda", num_lambda, "num_spectra", num_spectra, "num_components", num_components, "X", X, "M", M, "Y", Y, "Xr", Xr)
  return M
end

# Compute the PCA component with the largest eigenvalue
# X is data, r is vector of random numbers, s is preallocated memory; r & s  are of same length as each data point
function compute_pca_component!{T}(X::AbstractArray{T,2}, r::AbstractArray{T,1}, s::AbstractArray{T,1}; tol::Float64=1e-8, max_it::Int64 = 10 )
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  @assert(length(r)==num_lambda)
  #rand!(r)                                               # assume r is already randomized
  last_mag_s = 0.0
  #s = zeros(T,num_lambda)
  for c in 1:max_it
	s[:] = zero(T)
	for i in 1:num_spectra
	   BLAS.axpy!(dot(view(X,:,i),r),view(X,:,i),s)         # s += dot(X[:,i],r)*X[:,i]
	end
	mag_s = sqrt(sum(abs2,s))
	r[:]  = s/mag_s
	if abs(mag_s-last_mag_s) < tol*mag_s break end
	last_mag_s = mag_s
  end
  return r
end

# Compute the PCA component with the largest eigenvalue
function compute_pca_component{T}(X::AbstractArray{T,2}; tol::Float64=1e-12, max_it::Int64 = 10)
  num_lambda = size(X,1)
  r = rand(T,num_lambda)
  s = zeros(T,num_lambda)
  compute_pca_component!(X,r,s,tol=tol,max_it=max_it)
end

# Compute first num_components basis vectors for PCA, only return the basis vectors
function fit_pca_eford_proj_only{T}(X::Array{T,2}; num_components::Integer=4, tol::Float64=1e-12, max_it::Int64 = 10)
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  M = rand(T, (num_lambda,num_components) )  # random initialization is part of algorithm (i.e., not zeros)
  s = zeros(T,num_lambda)                    # pre-allocated memory for compute_pca_component
  Xtmp = X.-mean(X,2)                        # perform PCA after subtracting off mean
  compute_pca_component!(Xtmp,view(M,:,1),s,tol=tol)  # compute first component separately, since generalized algorithm will need to
  for j in 2:num_components
     for i in 1:num_spectra
       score = dot(view(Xtmp,:,i),view(M,:,j-1))
       Xtmp[:,i] .-= score*view(M,:,j-1)
	   end
	  compute_pca_component!(Xtmp,view(M,:,j),s,tol=tol,max_it=max_it)
  end
  return M
end

# Compute first num_components basis vectors for PCA, also return scores and print even more info
function fit_pca_eford{T}(X::Array{T,2}; num_components::Integer=4, tol::Float64=1e-12, max_it::Int64 = 10)
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  M = rand(T, (num_lambda,num_components) )  # random initialization is part of algorithm (i.e., not zeros)
  s = zeros(T,num_lambda)                    # pre-allocated memory for compute_pca_component
  scores = zeros(num_spectra,num_components)
  mu = vec(mean(X,2))
  Xtmp = X.-mu                               # perform PCA after subtracting off mean
  totalvar = sum(abs2,Xtmp)
  #println("# j = ", 0, " sum(abs2,Xtmp) = ", totalvar)
  compute_pca_component!(Xtmp,view(M,:,1),s,tol=tol)  # compute first component separately, since generalized algorithm will need to
  fracvar = zeros(num_components)
  compute_pca_component!(Xtmp,view(M,:,1),s,tol=tol)
  for j in 2:num_components
	  for i in 1:num_spectra
      scores[i,j-1] = dot(view(Xtmp,:,i),view(M,:,j-1)) #/sum(abs2,view(M,:,j-1))
      Xtmp[:,i] .-= scores[i,j-1]*view(M,:,j-1)
	  end
    fracvar[j-1] = sum(abs2,Xtmp)/totalvar
    println("# j = ", j-1, " sum(abs2,Xtmp) = ", sum(abs2,Xtmp), " frac_var_remain= ", fracvar[j-1] )
	  compute_pca_component!(Xtmp,view(M,:,j),s,tol=tol,max_it=max_it)
  end
	for i in 1:num_spectra
       scores[i,num_components] = dot(view(Xtmp,:,i),view(M,:,num_components))
       Xtmp[:,i] .-= scores[i,num_components]*view(M,:,num_components)
	end
  fracvar[num_components] = sum(abs2,Xtmp)/totalvar
  println("# j = ", num_components, " sum(abs2,Xtmp) = ", sum(abs2,Xtmp), " fracvar= ", fracvar[num_components] )
	return (mu,M,scores)
end

# test that eford's PCA routine is giving accurate results for dataset X and min(num_components,number of components kept by Julia's PCA routine)
function test_pca_eford{T}(obs::Array{T,2}; num_components::Integer=4 )
  pca_out = fit_pca_default(obs)
  (pca_eford_mu, pca_eford_out, pca_eford_scores) = fit_pca_eford(obs,num_components=num_components)
  for i in 1:min(size(projection(pca_out),2),size(pca_eford_out,2))
     if ! (projection(pca_out)[1,i] * pca_eford_out[1,i] > 0.0)
	    pca_eford_out[:,i] = -pca_eford_out[:,i]
        pca_eford_scores[:,i] = -pca_eford_scores[:,i]
	 end
  end
  maxabsdif = maximum(abs(projection(pca_out).-pca_eford_out[:,1:size(projection(pca_out),2)]))
  return maxabsdif
end

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
__precompile__()

module MP3

deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(deps) ? include(deps) : error("MP3 is not properly installed. Please run: Pkg.build(\"MP3\")")

# package code goes here
using SampledSignals
using FileIO
using Libdl

# methods to override
import SampledSignals: nchannels, nframes, samplerate, unsafe_read!, unsafe_write
import FileIO: load, save

# re-export
export load, save
export Hz, kHz, s

struct MP3INFO
    nframes::Int64
    nchannels::Int32
    samplerate::Int32
    datatype::DataType
end

"""create an MP3INFO object from given audio buffer"""
function MP3INFO(buf::SampleBuf{T}) where {T}
    MP3INFO(nframes(buf), nchannels(buf), samplerate(buf), T)
end

include("readers.jl")
include("writers.jl")

function __init__()
    initialize_readers()
    initialize_writers()

    # MP3 files with ID3v1 (or no tags) and ID3v2 tags have different headers
    magic = (UInt8[0xff, 0xfb], UInt8[0x49, 0x44, 0x33])
    add_format(format"MP3", magic, [".mp3"], [:MP3])
end

end # module


text:
module NaiveVariationalBayesTest

using Test
using ForneyLab
using ForneyLab: SoftFactor, generateId, addNode!, associate!, inferUpdateRule!, outboundType, isApplicable, setTargets!, messagePassingSchedule
using ForneyLab: VBGaussianMeanVarianceOut, VBGaussianMeanPrecisionM, SPEqualityGaussian

# Integration helper
mutable struct MockNode <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Array{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@naiveVariationalRule(:node_type     => MockNode,
                      :outbound_type => Message{PointMass},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBMockOut)

@testset "@naiveVariationalRule" begin
    @test VBMockOut <: NaiveVariationalRule{MockNode}
    @test isApplicable(VBMockOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBMockOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
end

@testset "inferUpdateRule!" begin
    FactorGraph()
    nd = MockNode([Variable(), constant(0.0), constant(0.0)])

    entry = ScheduleEntry(nd.i[1], NaiveVariationalRule{MockNode})
    inferUpdateRule!(entry, entry.message_update_rule, Dict{Interface, Type}())

    @test entry.message_update_rule == VBMockOut
end

@testset "messagePassingSchedule" begin
    g = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    w = Variable()
    nd_w = Gamma(w, constant(1.0), constant(1.0))
    y = Variable[]
    nd_y = FactorNode[]
    for i = 1:3
        y_i = Variable()
        nd_y_i = GaussianMeanPrecision(y_i, m, w)
        placeholder(y_i, :y, index=i)
        push!(y, y_i)
        push!(nd_y, nd_y_i)
    end

    pfz = PosteriorFactorization()
    q_m = PosteriorFactor(m)

    setTargets!(q_m, pfz, external_targets=true)
    schedule = messagePassingSchedule(q_m)

    @test length(schedule) == 6
    @test ScheduleEntry(nd_m.i[:out], VBGaussianMeanVarianceOut) in schedule
    @test ScheduleEntry(nd_y[2].i[:m], VBGaussianMeanPrecisionM) in schedule
    @test ScheduleEntry(nd_y[3].i[:m], VBGaussianMeanPrecisionM) in schedule
    @test ScheduleEntry(nd_m.i[:out].partner.node.i[3].partner, SPEqualityGaussian) in schedule
    @test ScheduleEntry(nd_y[1].i[:m], VBGaussianMeanPrecisionM) in schedule
    @test ScheduleEntry(nd_m.i[:out].partner, SPEqualityGaussian) in schedule
end

@testset "messagePassingAlgorithm" begin
    fg = FactorGraph()
    m = Variable()
    nd_m = GaussianMeanVariance(m, constant(0.0), constant(1.0))
    w = Variable()
    nd_w = Gamma(w, constant(1.0), constant(1.0))
    y = Variable[]
    nd_y = FactorNode[]
    for i = 1:3
        y_i = Variable()
        nd_y_i = GaussianMeanPrecision(y_i, m, w)
        placeholder(y_i, :y, index=i)
        push!(y, y_i)
        push!(nd_y, nd_y_i)
    end

    pfz = PosteriorFactorization(m, w)
    algo = messagePassingAlgorithm(m)

    @test isa(algo, InferenceAlgorithm)
end

@testset "messagePassingAlgorithm with custom target" begin
    fg = FactorGraph()

    # Define model
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    @RV y = 2.0*x
    @RV z = y + 1.0
    @RV w ~ GaussianMeanPrecision(z, 1.0)
    placeholder(w, :w)

    # Derive algorithm
    pfz = PosteriorFactorization(y)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)

    @test occursin("marginals[:y]", code)
end

end # module