text:
## Defines all cipher functionality
## As usual, check out
#http://www.lysator.liu.se/~nisse/nettle/nettle.html#Cipher-functions

import Base: show
export CipherType, get_cipher_types
export gen_key32_iv16, add_padding_PKCS5, trim_padding_PKCS5
export Encryptor, Decryptor, decrypt, decrypt!, encrypt, encrypt!

# This is a mirror of the nettle-meta.h:nettle_cipher struct
struct NettleCipher
    name::Ptr{UInt8}
    context_size::Cuint
    block_size::Cuint
    key_size::Cuint
    set_encrypt_key::Ptr{Cvoid}
    set_decrypt_key::Ptr{Cvoid}
    encrypt::Ptr{Cvoid}
    decrypt::Ptr{Cvoid}
end

# For much the same reasons as in hash_common.jl, we define a separate, more "Julia friendly" type
struct CipherType
    name::AbstractString
    context_size::Cuint
    block_size::Cuint
    key_size::Cuint
    set_encrypt_key::Ptr{Cvoid}
    set_decrypt_key::Ptr{Cvoid}
    encrypt::Ptr{Cvoid}
    decrypt::Ptr{Cvoid}
end

# These are the user-facing types that are used to actually {en,de}cipher stuff
struct Encryptor
    cipher_type::CipherType
    state::Vector{UInt8}
end
struct Decryptor
    cipher_type::CipherType
    state::Vector{UInt8}
end

# The function that maps from a NettleCipher to a CipherType
function CipherType(nc::NettleCipher)
    CipherType( uppercase(unsafe_string(nc.name)),
                nc.context_size, nc.block_size, nc.key_size,
                nc.set_encrypt_key, nc.set_decrypt_key, nc.encrypt, nc.decrypt)
end

# The global dictionary of hash types we know how to construct
const _cipher_types = Dict{AbstractString,CipherType}()

# We're going to load in each NettleCipher struct individually, deriving
# HashAlgorithm types off of the names we find, and storing the output
# and context size from the data members in the C structures
function get_cipher_types()
    # If we have already gotten the hash types from libnettle, don't query again
    if isempty(_cipher_types)
        cipher_idx = 1
        # nettle_ciphers is an array of pointers ended by a NULL pointer, continue reading hash types until we hit it
        while( true )
            ncptr = unsafe_load(cglobal(("nettle_ciphers",libnettle),Ptr{Ptr{Cvoid}}),cipher_idx)
            if ncptr == C_NULL
                break
            end
            cipher_idx += 1
            nc = unsafe_load(convert(Ptr{NettleCipher}, ncptr))
            cipher_type = CipherType(nc)
            _cipher_types[cipher_type.name] = cipher_type
        end
    end
    return _cipher_types
end

const _cipher_suites = [:CBC] # [:CBC, :GCM, :CCM]

function gen_key32_iv16(pw::Vector{UInt8}, salt::Vector{UInt8})
    s1 = digest("MD5", [pw; salt])
    s2 = digest("MD5", [s1; pw; salt])
    s3 = digest("MD5", [s2; pw; salt])
    return ([s1; s2], s3)
end

function add_padding_PKCS5(data::Vector{UInt8}, block_size::Int)
  padlen = block_size - (sizeof(data) % block_size)
  return [data; map(i -> UInt8(padlen), 1:padlen)]
end

function trim_padding_PKCS5(data::Vector{UInt8})
  padlen = data[sizeof(data)]
  return data[1:sizeof(data)-padlen]
end

function Encryptor(name::AbstractString, key)
    cipher_types = get_cipher_types()
    name = uppercase(name)
    if !haskey(cipher_types, name)
        throw(ArgumentError("Invalid cipher type $name: call Nettle.get_cipher_types() to see available list"))
    end
    cipher_type = cipher_types[name]

    if sizeof(key) != cipher_type.key_size
        throw(ArgumentError("Key must be $(cipher_type.key_size) bytes long"))
    end

    state = Vector{UInt8}(undef, cipher_type.context_size)
    ccall( cipher_type.set_encrypt_key, Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), state, pointer(key))

    return Encryptor(cipher_type, state)
end

function Decryptor(name::AbstractString, key)
    cipher_types = get_cipher_types()
    name = uppercase(name)
    if !haskey(cipher_types, name)
        throw(ArgumentError("Invalid cipher type $name: call Nettle.get_cipher_types() to see available list"))
    end
    cipher_type = cipher_types[name]

    if sizeof(key) != cipher_type.key_size
        throw(ArgumentError("Key must be $(cipher_type.key_size) bytes long"))
    end

    state = Vector{UInt8}(undef, cipher_type.context_size)
    ccall( cipher_type.set_decrypt_key, Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), state, pointer(key))

    return Decryptor(cipher_type, state)
end

function decrypt!(state::Decryptor, e::Symbol, iv::Vector{UInt8}, result, data)
    if sizeof(result) < sizeof(data)
        throw(ArgumentError("Output array of length $(sizeof(result)) insufficient for input data length ($(sizeof(data)))"))
    end
    if sizeof(result) % state.cipher_type.block_size > 0
        throw(ArgumentError("Output array of length $(sizeof(result)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(data) % state.cipher_type.block_size > 0
        throw(ArgumentError("Input array of length $(sizeof(data)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(iv) != state.cipher_type.block_size
        throw(ArgumentError("Iv must be $(state.cipher_type.block_size) bytes long"))
    end
    if ! (Symbol(uppercase(string(e))) in _cipher_suites)
        throw(ArgumentError("now supports $(_cipher_suites) only but ':$(e)'"))
    end
hdl = Libdl.dlopen_e(libnettle)
s = Symbol("nettle_", lowercase(string(e)), "_decrypt")
c = Libdl.dlsym(hdl, s)
if c == C_NULL
    throw(ArgumentError("not found function '$(s)' for ':$(e)'"))
end
# c points (:nettle_***_decrypt, nettle) may be loaded as another instance
iiv = copy(iv)
ccall(c, Cvoid, (
    Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, Ptr{UInt8},
    Csize_t, Ptr{UInt8}, Ptr{UInt8}),
    state.state, state.cipher_type.decrypt, sizeof(iiv), iiv,
    sizeof(data), pointer(result), pointer(data))
Libdl.dlclose(hdl)
    return result
end

function decrypt!(state::Decryptor, result, data)
    if sizeof(result) < sizeof(data)
        throw(ArgumentError("Output array of length $(sizeof(result)) insufficient for input data length ($(sizeof(data)))"))
    end
    if sizeof(result) % state.cipher_type.block_size > 0
        throw(ArgumentError("Output array of length $(sizeof(result)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(data) % state.cipher_type.block_size > 0
        throw(ArgumentError("Input array of length $(sizeof(data)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    ccall(state.cipher_type.decrypt, Cvoid, (Ptr{Cvoid},Csize_t,Ptr{UInt8},Ptr{UInt8}),
        state.state, sizeof(data), pointer(result), pointer(data))
    return result
end

function decrypt(state::Decryptor, e::Symbol, iv::Vector{UInt8}, data)
    result = Vector{UInt8}(undef, sizeof(data))
    decrypt!(state, e, iv, result, data)
    return result
end

function decrypt(state::Decryptor, data)
    result = Vector{UInt8}(undef, sizeof(data))
    decrypt!(state, result, data)
    return result
end

function encrypt!(state::Encryptor, e::Symbol, iv::Vector{UInt8}, result, data)
    if sizeof(result) < sizeof(data)
        throw(ArgumentError("Output array of length $(sizeof(result)) insufficient for input data length ($(sizeof(data)))"))
    end
    if sizeof(result) % state.cipher_type.block_size > 0
        throw(ArgumentError("Output array of length $(sizeof(result)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(data) % state.cipher_type.block_size > 0
        throw(ArgumentError("Input array of length $(sizeof(data)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(iv) != state.cipher_type.block_size
        throw(ArgumentError("Iv must be $(state.cipher_type.block_size) bytes long"))
    end
    if ! (Symbol(uppercase(string(e))) in _cipher_suites)
        throw(ArgumentError("now supports $(_cipher_suites) only but ':$(e)'"))
    end
hdl = Libdl.dlopen_e(libnettle)
s = Symbol("nettle_", lowercase(string(e)), "_encrypt")
c = Libdl.dlsym(hdl, s)
if c == C_NULL
    throw(ArgumentError("not found function '$(s)' for ':$(e)'"))
end
# c points (:nettle_***_encrypt, nettle) may be loaded as another instance
iiv = copy(iv)
ccall(c, Cvoid, (
    Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, Ptr{UInt8},
    Csize_t, Ptr{UInt8}, Ptr{UInt8}),
    state.state, state.cipher_type.encrypt, sizeof(iiv), iiv,
    sizeof(data), pointer(result), pointer(data))
Libdl.dlclose(hdl)
    return result
end

function encrypt!(state::Encryptor, result, data)
    if sizeof(result) < sizeof(data)
        throw(ArgumentError("Output array of length $(sizeof(result)) insufficient for input data length ($(sizeof(data)))"))
    end
    if sizeof(result) % state.cipher_type.block_size > 0
        throw(ArgumentError("Output array of length $(sizeof(result)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    if sizeof(data) % state.cipher_type.block_size > 0
        throw(ArgumentError("Input array of length $(sizeof(data)) must be N times $(state.cipher_type.block_size) bytes long"))
    end
    ccall(state.cipher_type.encrypt, Cvoid, (Ptr{Cvoid},Csize_t,Ptr{UInt8},Ptr{UInt8}),
        state.state, sizeof(data), pointer(result), pointer(data))
    return result
end

function encrypt(state::Encryptor, e::Symbol, iv::Vector{UInt8}, data)
    result = Vector{UInt8}(undef, sizeof(data))
    encrypt!(state, e, iv, result, data)
    return result
end

function encrypt(state::Encryptor, data)
    result = Vector{UInt8}(undef, sizeof(data))
    encrypt!(state, result, data)
    return result
end

# The one-shot functions that make this whole thing so easy
decrypt(name::AbstractString, key, data) = decrypt(Decryptor(name, key), data)
encrypt(name::AbstractString, key, data) = encrypt(Encryptor(name, key), data)

decrypt(name::AbstractString, e::Symbol, iv::Vector{UInt8}, key, data) = decrypt(Decryptor(name, key), e, iv, data)
encrypt(name::AbstractString, e::Symbol, iv::Vector{UInt8}, key, data) = encrypt(Encryptor(name, key), e, iv, data)

# Custom show overrides make this package have a little more pizzaz!
function show(io::IO, x::CipherType)
    write(io, "$(x.name) Cipher\n")
    write(io, "  Context size: $(x.context_size) bytes\n")
    write(io, "  Block size: $(x.block_size) bytes\n")
    write(io, "  Key size: $(x.key_size) bytes")
end
show(io::IO, x::Encryptor) = write(io, "$(x.cipher_type.name) Encryption state")
show(io::IO, x::Decryptor) = write(io, "$(x.cipher_type.name) Decryption state")


text:
remove_vk_prefix(name) = replace(name, r"^(?:vk|Vk|VK_)" => "")
remove_vk_prefix(name::Symbol) = Symbol(remove_vk_prefix(string(name)))

function remove_vk_prefix(ex::Expr)
    postwalk(ex) do ex
        ex isa Symbol ? Symbol(remove_vk_prefix(string(ex))) : ex
    end
end

const convention_exceptions = Dict(
    :textureCompressionASTC_LDR => :texture_compression_astc_ldr,
    :textureCompressionASTC_HDR => :texture_compression_astc_hdr,
    :formatA4R4G4B4 => :format_a4r4g4b4,
    :formatA4B4G4R4 => :format_a4b4g4r4,
)

ext: jl


text:
"""
$(TYPEDEF)

Defines a time dependent Hamiltonian object with custom function.

# Fields

$(FIELDS)
"""
struct CustomDenseHamiltonian{T<:Number,in_place} <: AbstractDenseHamiltonian{T}
    """Function for the Hamiltonian `H(s)`"""
    f::Any
    """Internal cache"""
    u_cache::Any
    """Size"""
    size::Any
end

function hamiltonian_from_function(func; in_place = false)
    hmat = func(0.0)
    CustomDenseHamiltonian{eltype(hmat),in_place}(func, nothing, size(hmat))
end

get_cache(H::CustomDenseHamiltonian) = zeros(eltype(H), size(H))
(H::CustomDenseHamiltonian)(s) = H.f(s)


text:
for N in [Float64, Rational{Int}, Float32]
    # compact set
    @test Ball1 <: CompactSet && !(Ball1 <: NonCompactSet)
    b = Ball1(N[0, 1], N(1))
    @test b isa CompactSet && !(b isa NonCompactSet)

    # non-compact set
    @test Line <: NonCompactSet && !(Line <: CompactSet)
    l = Line(N[0, 1], N(1))
    @test l isa NonCompactSet && !(l isa CompactSet)

    # lazy operations are neither compact nor non-compact by default
    @test !(MinkowskiSum <: CompactSet) && !(MinkowskiSum <: NonCompactSet)
    ms = b + l
    @test !(ms isa CompactSet) && !(ms isa NonCompactSet)
end


text:
@testset "sdfilewriter" begin

@testset "molblock" begin
    demomol = joinpath(dirname(@__FILE__), "..", "assets", "test", "demo.mol")
    mol = sdftomol(demomol)
    mol2 = sdftomol(split(printv2mol(mol), "\n"))
    @test nodecount(mol2) == 37

    smol = smilestomol("CCC1CC(C=O)CCC1N")
    smol2 = sdftomol(split(printv2mol(smol), "\n"))
    @test nodecount(smol2) == 11
end

end #sdfilewriter


text:
using ImageView, Images, FileIO;
using StatsBase;
using Luxor;
using Colors;
using Plots;

img = load("./img/landscape.png");
tgsize = 35;

training_size = 200;

for split in 1:training_size
	fx, fy = size(img).-(2*tgsize, 2*tgsize)
	tx = sample((tgsize):fx)
	ty = sample((tgsize):fy)
  trial = img[(tx-tgsize):(tx),(ty-tgsize):(ty)]
  save("splits/train_"*string(split)*".jpg", trial)
end

max_stars_repo_name: tpoisot/WaterFinder

max_issues_repo_name: tpoisot/WaterFinder

max_forks_repo_name: tpoisot/WaterFinder


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    monoisotopicmass, nominalmass, exactmass, standardweight,
    isotopiccomposition, massspecpeaks, simulatemassspec

const DEFAULT_WEIGHT_DIGITS = 2
const DEFAULT_MASS_DIGITS = 6

function molecularmass(mol::GraphMol, massfunc, implicithfunc)
    mass = 0.0
    unc = 0.0
    hm, hu = implicithfunc(:H)
    for (atom, hcnt) in zip(nodeattrs(mol), implicithcount(mol))
        m, u = massfunc(atom)
        mass += m + hm * hcnt
        unc += u + hu * hcnt
    end
    return (mass, unc)
end

"""
    monoisotopicmass(atomsymbol::Symbol) -> Tuple{Float64,Float64}
    monoisotopicmass(atom::Atom) -> Tuple{Float64,Float64}
    monoisotopicmass(mol::GraphMol) -> Tuple{Float64,Float64}

Return a tuple of monoisotopic mass of the atom/molecule and its uncertainty.

Monoisotopic mass is the relative atomic mass of the most abundant isotope. Even if there is specific `Atom.mass` value, it will be ignored.
"""
function monoisotopicmass(atomsymbol::Symbol)
    num = atomnumber(atomsymbol)
    mass = ATOMTABLE[num]["Monoisotopic"]
    unc = ATOMTABLE[num]["MonoisotopicUncertainty"]
    return (mass, unc)
end

monoisotopicmass(atom::Atom) = monoisotopicmass(atom.symbol)
monoisotopicmass(mol::GraphMol
    ) = molecularmass(mol, monoisotopicmass, monoisotopicmass)

"""
    monoisotopicmass(::Type{Float64}, atomsymbol::Symbol) -> Float64
    monoisotopicmass(::Type{Float64}, atom::Atom) -> Float64
    monoisotopicmass(::Type{Float64}, mol::GraphMol) -> Float64

Return monoisotopic mass of the atom/molecule rounded to `digits=6`.
"""
function monoisotopicmass(::Type{Float64}, atomsymbol::Symbol)
    mass, unc = monoisotopicmass(atomsymbol)
    mass === NaN && return NaN
    return round(mass, digits=DEFAULT_MASS_DIGITS)
end

monoisotopicmass(::Type{Float64}, atom::Atom
    ) = monoisotopicmass(Float64, atom.symbol)

function monoisotopicmass(::Type{Float64}, mol::GraphMol)
    mass, unc = molecularmass(mol, monoisotopicmass, monoisotopicmass)
    mass === NaN && return NaN
    return round(mass, digits=DEFAULT_MASS_DIGITS)
end

"""
    monoisotopicmass(::Type{T}, atomsymbol::Symbol) where {T<:Integer} -> T
    monoisotopicmass(::Type{T}, atom::Atom) where {T<:Integer} -> T
    monoisotopicmass(::Type{T}, mol::GraphMol) where {T<:Integer} -> T
    nominalmass(atomsymbol::Symbol) -> Int
    nominalmass(atom::Atom) -> Int
    nominalmass(mol::GraphMol) -> Int

Return nominal mass of the atom/molecule.
"""
monoisotopicmass(::Type{T}, atomsymbol::Symbol
    ) where {T<:Integer} = round(T, monoisotopicmass(Float64, atomsymbol))

monoisotopicmass(::Type{T}, atom::Atom
    ) where {T<:Integer} = round(T, monoisotopicmass(Float64, atom.symbol))

monoisotopicmass(::Type{T}, mol::GraphMol
    ) where {T<:Integer} = round(T, monoisotopicmass(Float64, mol))

nominalmass(atomsymbol::Symbol) = monoisotopicmass(Int, atomsymbol)
nominalmass(atom::Atom) = monoisotopicmass(Int, atom)
nominalmass(mol::GraphMol) = monoisotopicmass(Int, mol)

"""
    exactmass(atomsymbol::Symbol) -> Tuple{Float64,Float64}
    exactmass(atomsymbol::Symbol, number::Int) -> Tuple{Float64,Float64}
    exactmass(atom::Atom) -> Tuple{Float64,Float64}
    exactmass(mol::GraphMol) -> Tuple{Float64,Float64}

Return a tuple of calculated exact mass and its uncertainty.

If `number` is not given or `Atom.mass` is not specified, monoisotopic mass will be used instead.
"""
exactmass(atomsymbol::Symbol) = monoisotopicmass(atomsymbol)

function exactmass(atomsymbol::Symbol, number::Int)
    pred = rcd -> rcd["Number"] == number
    iso = ATOMTABLE[atomnumber(atomsymbol)]["Isotopes"]
    k = findfirst(pred.(iso))
    mass = iso[k]["Mass"]
    unc = iso[k]["MassUncertainty"]
    return (mass, unc)
end

function exactmass(atom::Atom)
    atom.mass === nothing || return exactmass(atom.symbol, atom.mass)
    return monoisotopicmass(atom.symbol)
end

exactmass(mol::GraphMol) = molecularmass(mol, exactmass, monoisotopicmass)

"""
    exactmass(::Type{Float64}, atomsymbol::Symbol) -> Float64
    exactmass(::Type{Float64}, atomsymbol::Symbol, number::Int) -> Float64
    exactmass(::Type{Float64}, atom::Atom) -> Float64
    exactmass(::Type{Float64}, mol::GraphMol) -> Float64

Return calculated exact mass rounded to `digit=6`.
"""
exactmass(::Type{Float64}, atomsymbol::Symbol
    ) = monoisotopicmass(Float64, atomsymbol)

function exactmass(::Type{Float64}, atomsymbol::Symbol, number::Int)
    mass, unc = exactmass(atomsymbol, number)
    mass === NaN && return NaN
    return round(mass, digits=DEFAULT_MASS_DIGITS)
end

function exactmass(::Type{Float64}, atom::Atom)
    atom.mass === nothing || return exactmass(Float64, atom.symbol, atom.mass)
    return monoisotopicmass(Float64, atom.symbol)
end

function exactmass(::Type{Float64}, mol::GraphMol)
    mass, unc = molecularmass(mol, exactmass, monoisotopicmass)
    mass === NaN && return NaN
    return round(mass, digits=DEFAULT_MASS_DIGITS)
end

"""
    standardweight(atomsymbol::Symbol) -> Tuple{Float64,Float64}
    standardweight(atom::Atom) -> Tuple{Float64,Float64}
    standardweight(mol::GraphMol) -> Tuple{Float64,Float64}

Return a tuple of standard atomic weight (or molecular weight) and its uncertainty.

If `Atom.mass` is specified, calculated exact mass of the atom will be used instead.
"""
function standardweight(atomsymbol::Symbol)
    num = atomnumber(atomsymbol)
    wt = ATOMTABLE[num]["Weight"]
    unctype = ATOMTABLE[num]["WeightType"]
    if unctype == "Interval"
        lower = ATOMTABLE[num]["WeightLower"]
        higher = ATOMTABLE[num]["WeightHigher"]
        unc = max(higher - wt, wt - lower)
    elseif unctype == "Uncertainty"
        unc = ATOMTABLE[num]["WeightUncertainty"]
    else
        unc = NaN  # unctype in ("MostStable", "NotAvailable")
    end
    return (wt, unc)
end

function standardweight(atom::Atom)
    atom.mass === nothing || return exactmass(atom.symbol, atom.mass)
    return standardweight(atom.symbol)
end

standardweight(mol::GraphMol
    ) = molecularmass(mol, standardweight, standardweight)

"""
    standardweight(::Type{Float64}, atomsymbol::Symbol) -> Float64
    standardweight(::Type{Float64}, atom::Atom) -> Float64
    standardweight(::Type{Float64}, mol::GraphMol) -> Float64

Return standard atomic weight (or molecular weight) rounded to `digit=2`.
"""
function standardweight(::Type{Float64}, atomsymbol::Symbol)
    wt, unc = standardweight(atomsymbol)
    wt === NaN && return NaN
    return round(wt, digits=DEFAULT_WEIGHT_DIGITS)
end

function standardweight(::Type{Float64}, atom::Atom)
    atom.mass === nothing || return exactmass(Float64, atom.symbol, atom.mass)
    return standardweight(Float64, atom.symbol)
end

function standardweight(::Type{Float64}, mol::GraphMol)
    wt, unc = molecularmass(mol, standardweight, standardweight)
    wt === NaN && return NaN
    return round(wt, digits=DEFAULT_WEIGHT_DIGITS)
end

"""
    isotopiccomposition(atomsymbol::Symbol, number::Int; threshold=0.001
        ) -> Vector{Tuple{Float64,Float64}}
    isotopiccomposition(mol::GraphMol; threshold=0.001
        ) -> Vector{Tuple{Float64,Float64}}

Return isotopic composition of the atoms/molecule as a vector of tuples of mass and composition.

Records that have lower abundance than the given threshold will be filtered out (default 0.001 = 0.1%)
"""
function isotopiccomposition(atomsymbol::Symbol, number::Int; threshold=0.001)
    z = atomnumber(atomsymbol)
    isotopes = []
    for iso in ATOMTABLE[z]["Isotopes"]
        cmp = iso["Composition"]
        cmp === NaN && continue
        push!(isotopes, iso)
    end
    isocnt = length(isotopes)
    data = Tuple{Float64,Float64}[]
    for seps in combinations(number + isocnt - 1, isocnt - 1)
        nums = Int[]
        prev = 0
        for sep in seps
            push!(nums, sep - prev - 1)
            prev = sep
        end
        push!(nums, number + isocnt - prev - 1)
        mass = 0.0
        dups = [logfactorial(n) for n in nums]
        comp = exp(reduce(-, dups; init=logfactorial(number)))
        for (i, iso) in enumerate(isotopes)
            mass += isotopes[i]["Mass"] * nums[i]
            comp *= isotopes[i]["Composition"] ^ nums[i]
        end
        comp > threshold && push!(data, (mass, comp))
    end
    return data
end

function isotopiccomposition(mol::GraphMol; threshold=0.001)
    data = Tuple{Float64,Float64}[]
    for tup in Iterators.product(
            (isotopiccomposition(sym, cnt; threshold=threshold)
            for (sym, cnt) in countatoms(mol))...)
        mass = 0.0
        comp = 1.0
        for t in tup
            mass += t[1]
            comp *= t[2]
        end
        comp > threshold && push!(data, (mass, comp))
    end
    return sort(data, by=x->x[1])
end

"""
    massspecpeaks(mol::GraphMol; threshold=0.001) -> Matrix{Float64}

Return a vector of tuples of each isotopic masses and their relative intensity in the simulated mass spectrum (base peak intensity = 100).

Records that have lower abundance (not peak intensity) than the given threshold will be filtered out (default 0.001 = 0.1%)
"""
function massspecpeaks(mol::GraphMol; threshold=0.001)
    isocomp = isotopiccomposition(mol; threshold=threshold)
    bp = maximum(x->x[2], isocomp)
    data = Tuple{Float64,Float64}[]
    for (mass, comp) in isocomp
        mass = round(mass, digits=DEFAULT_MASS_DIGITS)
        rint = round(comp / bp * 100, digits=ceil(Int, -log(10, threshold))-2)
        push!(data, (mass, rint))
    end
    return data
end

function gaussianpeak(mass::Float64, intensity::Float64, resolution::Int)
    fwhm = mass / resolution
    s = fwhm / (2 * sqrt(2 * log(2)))
    return x -> exp(- (x - mass)^2 / (2 * s^2)) * intensity
end

"""
    simulatemassspec(peaks::Vector{Tuple{Float64,Float64}};
        resolution=10000, rate=0.01) -> Matrix{Float64}
    simulatemassspec(mol::GraphMol;
        threshold=0.001, resolution=10000, rate=0.01) -> Matrix{Float64}

Return a matrix of simulate mass spectrum (dim 1: datapoints, dim 2: mass and intensity).

# Usage (with Plot.jl)
```julia
using MolecularGraph
using Plots
gr()
Plots.GRBackend()

mol = smilestomol("CCO")
data = simulatemassspec(mol)
plot(
    data[:, 1], data[:, 2],
    leg=false, xlabel = "Mass", ylabel = "Intensity"
)
```
"""
function simulatemassspec(
        peaks::Vector{Tuple{Float64,Float64}}; resolution=10000, rate=0.001)
    fs = Function[]
    for (mass, intensity) in peaks
        push!(fs, gaussianpeak(mass, intensity, resolution))
    end
    superposed(x) = reduce(+, [f(x) for f in fs]; init=0.0)
    xleft = round(Int, minimum(x->x[1], peaks) - 1)
    xright = round(Int, maximum(x->x[1], peaks) + 1)
    arr = collect(xleft:rate:xright)
    return hcat(arr, superposed.(arr))
end

simulatemassspec(mol::GraphMol; kwargs...
    ) =  simulatemassspec(massspecpeaks(mol; kwargs...); kwargs...)


text:
using ViscosityDrop
using Test

@testset "Grids equal points" begin
    grid = make_grid(32,32,2,2)
    @test grid.Nx == 1
    @test grid.Ny == 32
    @test grid.Nz == 32
    @test grid.Ly == 2.0
    @test grid.Lz == 2.0
end

@testset "Grids unequal points" begin
    grid = make_grid(16,32,2,3)
    @test grid.Nx == 1
    @test grid.Ny == 16
    @test grid.Nz == 32
    @test grid.Ly == 2.0
    @test grid.Lz == 3.0
end

@testset "Boundary conditions" begin
    grid32 = make_grid(16,32,2,3)
    v_bc, w_bc = no_slip_bc(grid32)

    @test v_bc.x.left.condition == nothing
    @test v_bc.x.right.condition == nothing
    @test v_bc.y.left.condition == 0.0
    @test v_bc.y.right.condition == 0.0
    @test v_bc.z.left.condition == 0.0
    @test v_bc.z.right.condition == 0.0

    @test w_bc.x.left.condition == nothing
    @test w_bc.x.right.condition == nothing
    @test w_bc.y.left.condition == 0.0
    @test w_bc.y.right.condition == 0.0
    @test w_bc.z.left.condition == 0.0
    @test w_bc.z.right.condition == 0.0
end

@testset "Making models (default)" begin
    grid = make_grid(16,32,2,3)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, 1e-2)

    @test model.closure.ν == 1e-2
    @test Core.Typeof(model.timestepper).name.name == :RungeKutta3TimeStepper
    @test Core.Typeof(model.advection).name.name == :WENO5
    @test model.clock.iteration == 0
end

@testset "Making models (other options)" begin
    using Oceananigans.Advection
    grid = make_grid(16,32,2,3)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, 2e-2; timestepper = :QuasiAdamsBashforth2,
                                     advection = CenteredSecondOrder())

    @test model.closure.ν == 2e-2
    @test Core.Typeof(model.timestepper).name.name == :QuasiAdamsBashforth2TimeStepper
    @test Core.Typeof(model.advection).name.name == :CenteredSecondOrder
    @test model.clock.iteration == 0
end

@testset "Adding buoyancy" begin
    grid = make_grid(32,32,2,2)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, 1e-2)

    @test all(model.tracers.b[1,1:32,1:32] .== 0.0)

    set_buoyancy!(model)
    @test !all(model.tracers.b[1,1:32,1:32] .== 0.0)
    @test model.tracers.b[1,15,25] == -1e-1
    @test model.tracers.b[1,5,25] == 0.0

    set_buoyancy!(model, z_position=0.6, width=0.3, amplitude=-2e-1)
    @test !all(model.tracers.b[1,1:32,1:32] .== 0.0)
    @test model.tracers.b[1,15,25] == -2e-1
    @test model.tracers.b[1,5,25] == 0.0
end

@testset "Simulation setup" begin
    filename = "tmp"
    dir = "."
    grid = make_grid(32,32,2,2)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, 1e-2)
    set_buoyancy!(model)

    sim = make_simulation(model, dir, filename, output_fields=:all)
    @test sim.iteration_interval == 1000
    @test sim.run_time == 0.0
    @test sim.stop_time == 40.0
    @test sim.Δt == 0.002
    @test sim.output_writers.vals[1].schedule.interval == 0.5
    @test typeof(sim.output_writers.vals[1]).parameters[1].parameters[1] == (:b, :v, :w, :ω, :s)
    @test isfile("./tmp.jld2")

    if isfile("./tmp.jld2")
        rm("./tmp.jld2")
    end

    sim2 = make_simulation(model, dir, filename, Δt=0.001,
                                                 stop_time=20.0,
                                                 iteration_interval=2000,
                                                 schedule_interval=0.4,
                                                 output_fields=:b)
    @test sim2.iteration_interval == 2000
    @test sim2.run_time == 0.0
    @test sim2.stop_time == 20.0
    @test sim2.Δt == 0.001
    @test sim2.output_writers.vals[1].schedule.interval == 0.4
    @test isfile("./tmp.jld2")

    if isfile("./tmp.jld2")
        rm("./tmp.jld2")
    end
end

@testset "Running simulation" begin
    filename = "tmp"
    dir = "."
    grid = make_grid(16,16,2,2)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, 1e-2)
    set_buoyancy!(model)

    sim = make_simulation(model, dir, filename, output_fields=:all, stop_time=2.0)
    run!(sim)

    @test sim.model.clock.time == 2.0
    @test isfile("./tmp.jld2")

    if isfile("./tmp.jld2")
        rm("./tmp.jld2")
    end
end

@testset "Wrapper" begin
    filename = "tmp"
    dir = "."

    do_all(filename, dir, stop_time=2)
    @test isfile("./tmp.jld2")

    using JLD2
    file = jldopen("./tmp.jld2")
    @test typeof(file["timeseries/t"]) == JLD2.Group{JLD2.JLDFile{JLD2.MmapIO}}

    last_key = last(keys(file["timeseries/t"]))
    last_time = "timeseries/t/" * last_key
    @test file[last_time] >= 2

    if isfile("./tmp.jld2")
        rm("./tmp.jld2")
    end
end

@testset "Post-processing" begin
    filename = "tmp"
    dir = "."
    do_all(filename, dir, stop_time=2)
    out_results = get_results(dir, dir, "tmp.jld2"; Ncourse=11, verbose=false)

    @test isfile("./tmp.txt")
    @test size(out_results, 2) == 5

    using JLD2
    file = jldopen("./tmp.jld2")
    Nt = length(keys(file["timeseries/t"]))

    @test size(out_results, 1) == Nt

    if isfile("./tmp.jld2")
        rm("./tmp.jld2")
    end
    if isfile("./tmp.txt")
        rm("./tmp.txt")
    end
end


text:
#=
Implement an LRU (Least Recently Used) cache. It should be able to be initialized with a cache size n, and contain the following methods:

    set(key, value): sets key to value. If there are already n items in the cache and we are adding a new item, then it should also remove the least recently used item.
    get(key): gets the value at key. If no such key exists, return null.

Each operation should run in O(1) time.
=#

mutable struct LRU_cache
    cache::Dict{String, Int}
    n::Int
    key::Array{String, 1}
end

function LRU_cache(n::Int)
    LRU_cache(Dict(), n, [])
end

function set(key::String,value::Int)
    LRU.cache[key] = value
    push!(LRU.key, key)
    if length(LRU.cache) > LRU.n
        key_to_remove = LRU.key[1]
        delete!(LRU.cache, key_to_remove)
        popfirst!(LRU.key)
    end
end

function get(key::String)
    try
        return LRU.cache[key]
    catch
        return nothing
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module TransitiveDependence

import ..DebugMsg
DebugMsg.init()

using CompilerTools
using CompilerTools.Helper
using CompilerTools.AstWalker
using CompilerTools.LambdaHandling
using CompilerTools.LivenessAnalysis
using CompilerTools.Helper

import Base.show

function show(io::IO, d :: Dict{LHSVar, Set{LHSVar}})
    for ventry in d
        vdef    = ventry[1]
        vdepset = ventry[2]
        println(vdef, " => ", vdepset)
    end
end

type cddata
    ret :: Dict{LHSVar, Set{LHSVar}}
    defs
    uses
    cb
    data
end

type CallbackResult
    ret :: Dict{LHSVar, Set{LHSVar}}
    defs
    uses
    read_exprs
    write_exprs
end

function mergeDepSet(merge_into :: Dict{LHSVar, Set{LHSVar}},
                     merge_from :: Dict{LHSVar, Set{LHSVar}})
    for cbitem in merge_from
        vdef    = cbitem[1]
        vdefset = cbitem[2]

        if haskey(merge_into, vdef)
            cur_set = merge_into[vdef]
        else
            cur_set = Set{LHSVar}()
        end

        # Compute the new immediate dependency set for this variable.
        cur_set = union(cur_set, vdefset)
        # Remember the new immediate dependency set for this variable.
        merge_into[vdef] = cur_set
    end

    return merge_into
end

function finishStatement(data :: cddata)
    for cur_def in data.defs
        # Start with an empty use set if we haven't processed a statement that defines this variable before.
        # Else, start with the immediate dependency set as previously computed for this variable.
        if haskey(data.ret, cur_def)
            cur_set = data.ret[cur_def]
        else
            cur_set = Set{LHSVar}()
        end

        # Compute the new immediate dependency set for this variable.
        cur_set = union(cur_set, data.uses)
        # Remember the new immediate dependency set for this variable.
        data.ret[cur_def] = cur_set
    end
    data.defs = Set{LHSVar}()
    data.uses = Set{LHSVar}()
end

function cdwalk(x :: ANY, data :: cddata, top_level_number, is_top_level, read)
    if is_top_level
        finishStatement(data)
    end

    if data.cb != nothing
        cbres = data.cb(x, data.data)
        if cbres != nothing
            assert(isa(cbres, CallbackResult))
            mergeDepSet(data.ret, cbres.ret)
            data.uses = union(data.uses, cbres.uses)
            data.defs = union(data.defs, cbres.defs)
            for r in cbres.read_exprs
                CompilerTools.AstWalker.from_expr(r, 1, cdwalk, data, 1, false, true)
            end
            for w in cbres.write_exprs
                CompilerTools.AstWalker.from_expr(w, 1, cdwalk, data, 1, false, false)
            end
            return x
        end
    end

    if isa(x, RHSVar)
        lhsvar = toLHSVar(x)
        if read
            push!(data.uses, lhsvar)
        else
            push!(data.defs, lhsvar)
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function computeDependenciesAST(body :: ANY, callback=not_handled, cbdata :: ANY=nothing)
    data = cddata(Dict{LHSVar,Set{LHSVar}}(), Set{LHSVar}(), Set{LHSVar}(), callback, cbdata)
    AstWalk(body, cdwalk, data)
    finishStatement(data)
    @dprintln(3, "Before transitive in computeDependenciesAST")
    @dprintln(3, data.ret)
    return transitive(data.ret)
end

function transitive(ret :: Dict{LHSVar, Set{LHSVar}})
    change = true

    while change
        change = false

        for ventry in ret
            vdef    = ventry[1]
            vdepset = ventry[2]

            new_depset = vdepset

            for dep in vdepset
                if haskey(ret, dep)
                    new_depset = union(new_depset, ret[dep])
                    @dprintln(3, vdef, " after merging ", dep, " => ", new_depset)
                end
            end

            if new_depset != vdepset
                change = true
                ret[vdef] = new_depset
            end
        end
    end

    return ret
end

function computeDependencies(bl :: CompilerTools.LivenessAnalysis.BlockLiveness)
    @dprintln(3,"computeDependencies bl = ", bl)

    ret = Dict{LHSVar, Set{LHSVar}}()

    # Make one pass over all the basic blocks.
    for bbentry in bl.basic_blocks
        bb = bbentry[2]

        # For each statement in every basic block.
        for stmt in bb.statements
            # For each variable defined by this statement.
            for cur_def in stmt.def
                # Start with an empty use set if we haven't processed a statement that defines this variable before.
                # Else, start with the immediate dependency set as previously computed for this variable.
                if haskey(ret, cur_def)
                    cur_set = ret[cur_def]
                else
                    cur_set = Set{LHSVar}()
                end

                # Compute the new immediate dependency set for this variable.
                cur_set = union(cur_set, stmt.use)
                # Remember the new immediate dependency set for this variable.
                ret[cur_def] = cur_set
            end
        end
    end

    @dprintln(3, "Initial vdepset based on statements.")
    @dprintln(3, ret)

    return transitive(ret)
end

end


text:
#function show(io::IO, M::ResidueRingPolyMap)
#  println(io, "Map $(domain(M)) -> $(codomain(M)): gen -> $(M.gen_image)")
#  if isdefined(M, :coeff_map)
#    println(io, "using coeff map: ", M.coeff_map)
#  end
#end
#
#function ResidueRingPolyMap_image(M::ResidueRingPolyMap{D, C}, a::Generic.Res{T}) where {D, C, T}
#  #a should be in the domain of M...
#  if isdefined(M, :coeff_map)
#    I = codomain(M)(0)
#    for i=degree(a.data):-1:0
#      I = I*M.gen_image + M.coeff_map(coeff(a.data, i))
#    end
#  else
#    I = codomain(M)(0)
#    for i=degree(a.data):-1:0
#      I = I*M.gen_image + coeff(a.data, i)
#    end
#  end
#  return I
#end
#
#function ResidueRingPolyMap_preimage(M::ResidueRingPolyMap{D, C}, a::Generic.Res{T}) where {D, C, T}
#  #a should be in the codomain of M...
#  #
#  # a good 1st attempt. If many preimages are needed, the matrix Mt
#  # and the rref or whatever solve is used should be cached and reused.
#  #
#  R = codomain(M)
#  parent(a) == R || throw("mixed rings in preimage")
#  g = gens(domain(M))
#  im_gen = map(x->M(x), g) ## possibly should be cached and stored
#  ## need to write the elements in a matrix, solve the eqn for a
#  Mt = zero_matrix(base_ring(base_ring(R)), degree(R.modulus), degree(R.modulus))
#  for i=1:degree(R.modulus)
#    elem_to_mat_row!(Mt, i, im_gen[i])
#  end
#  b = zero_matrix(base_ring(base_ring(R)), 1, degree(R.modulus))
#  elem_to_mat_row!(b, 1, a)
#  s = solve_rational(Mt', b') # why, oh why is solve operating on columns????
#  if isa(s, Tuple) ## again, why, oh why is solve doing things differently
#                   ## over rings than fields?
#    s = s[1] * inv(s[2]) # all rings here (type) are actually fields (math)
#  end
#  ## need a version returning a bool to test is solution exists
#  ## without an error!
#  if isdefined(M, :coeff_map)
#    s = sum([preimage(M.coeff_map, s[i,1])*g[i] for i=1:length(g)])
#  else
#    s = sum([s[i,1]*g[i] for i=1:length(g)])
#  end
#  @assert parent(s) == domain(M)
#  return s
#end
#


text:
#=---------------------------------------------------------------------

ransacdemo - Functions demonstrating the use of RANSAC to fit lines,
planes, fundamental matrices and homographies.

Copyright (c) 2016 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK March     2016
   September 2018 Updates for v0.7/v1.0

---------------------------------------------------------------------=#

export fitlinedemo, fitplanedemo
export fitfunddemo, fithomogdemo

using ImageProjectiveGeometry, PyPlot, Printf, FileIO

#-----------------------------------------------------------------------
"""
fitlinedemo - Demonstrates RANSAC line fitting.

Function generates a noisy set of data points with outliers and uses
RANSAC to fit a line.

```
Usage: fitlinedemo(outliers, sigma, t, feedback=false)

Arguments:
              outliers - Fraction specifying how many points are to be
                         outliers.
              sigma    - Standard deviation of inlying points from the
                         true line.
              t        - Distance threshold to be used by the RANSAC
                         algorithm for deciding whether a point is an
                         inlier.
              feedback - Optional Boolean flag turn on RANSAC feedback
                         information.

Try using:  fitlinedemo(0.3, 0.05, 0.05)
```
See also: ransacfitplane(), fitplane()
"""
function fitlinedemo(outliers, sigma, t::Real, feedback::Bool = false)
    # August 2006  testfitline created from testfitplane
    #              author: Felix Duvallet

    # Hard wire some constants - vary these as you wish
    npts = 100  # Number of 3D data points

    # Define a line:
    #    Y = m*X
    #    Z = n*X + Y + b
    # This definition needs fixing, but it works for now
    m = 6
    n = -3
    b = -4

    outsigma = 30*sigma  # outlying points have a distribution that is
                         # 30 times as spread as the inlying points

    vpts = round(Int, (1-outliers)*npts)  # No of valid points
    opts = npts - vpts                    # No of outlying points

    # Generate npts points in the line
    X = rand(Float64, 1,npts)

    Y = m*X
    Z = n*X .+ Y .+ b
    Z = zeros(size(Y))

    XYZ =  [X
    	    Y
    	    Z]

    # Add uniform noise of +/-sigma
    XYZ = XYZ .+ (2*rand(Float64, size(XYZ)) .- 1)*sigma

    # Generate opts random outliers
    n = size(XYZ,2)
    ind = randperm(n)[1:opts]  # get a random set of point indices of length opts

    # Add uniform noise of outsigma to the points chosen to be outliers.
    XYZ[:,ind] = XYZ[:,ind] .+ sign.(rand(Float64, 3,opts).-0.5).*(rand(Float64, 3,opts).+1)*outsigma

    # Perform RANSAC fitting of the line
    (V, P, inliers) = ransacfitline(XYZ, t, feedback)

    if feedback
        @printf("Number of Inliers: %d\n", length(inliers))
    end

    # Plot the inlier points blue, with the outlier points in red.
    # Display the cloud of outlier points
    figure(1); clf()
    outliers = setdiff(1:npts, inliers)
    plot3D(vec(XYZ[1,outliers]), vec(XYZ[2, outliers]), vec(XYZ[3, outliers]), "r*")
#    hold(true)

    # Plot the inliers as blue points
    plot3D(vec(XYZ[1,inliers]), vec(XYZ[2, inliers]), vec(XYZ[3, inliers]), "b*")

    # Display the line formed by the 2 points that gave the
    # line of maximum consensus as a green line
    plot3D(vec(P[1,:]), vec(P[2,:]), vec(P[3,:]), "g-", linewidth=4)

    #Display the line formed by the covariance fitting in magenta
    plot3D(vec(V[1,:]), vec(V[2, :]), vec(V[3,:]), "m-", linewidth=4)
    box("on"), grid("on")
#    hold(false)

end
#-----------------------------------------------------------------------
"""
fitplanedemo - Demonstrates RANSAC plane fitting

Function generates a noisy set of data points with outliers and uses
RANSAC to fit a plane.

```
Usage: fitplanedemo(outliers, sigma, t, feedback)

Arguments:
              outliers - Fraction specifying how many points are to be
                         outliers.
              sigma    - Standard deviation of inlying points from the
                         true plane.
              t        - Distance threshold to be used by the RANSAC
                         algorithm for deciding whether a point is an
                         inlier.
              feedback - Optional flag 0 or 1 to turn on RANSAC feedback
                         information.

Try using:  fitplanedemo(0.3, 0.02, 0.05)
```
See also: ransacfitplane(), fitplane()
"""
function fitplanedemo(outliers, sigma, t, feedback::Bool = false)

    # Hard wire some constants - vary these as you wish
    npts = 100  # Number of 3D data points

    # Define a plane  ax + by + cz + d = 0
    a = 10; b = -3; c = 5; d = 1

    B = [a, b, c, d]

    outsigma = 30*sigma  # outlying points have a distribution that is
                         # 30 times as spread as the inlying points

    vpts = round(Int, (1-outliers)*npts)  # No of valid points
    opts = npts - vpts                    # No of outlying points

    # Generate npts points in the plane
    X = rand(Float64, npts)'
    Y = rand(Float64, npts)'
    Z = (-a*X .- b*Y .-d)/c

    XYZ =  [X
	    Y
	    Z]

    # Add uniform noise of +/-sigma
    XYZ = XYZ .+ (2*rand(Float64, size(XYZ)) .- 1)*sigma

    # Generate opts random outliers

    n = size(XYZ,2)
    ind = randperm(n)[1:opts]  # get a random set of point indices of length opts

    # Add uniform noise of outsigma to the points chosen to be outliers.
    #  XYZ(:,ind) = XYZ(:,ind)  + (2*rand(3,opts)-1)*outsigma
    XYZ[:,ind] = XYZ[:,ind] .+ sign.(rand(Float64, 3,opts).-0.5).*(rand(Float64, 3,opts).+1)*outsigma

    # Display the cloud of points (vec() needed for 0.4)
    figure(1); clf(); plot3D(vec(XYZ[1,:]),vec(XYZ[2,:]),vec(XYZ[3,:]), "r*")
    box("on")

    # Perform RANSAC fitting of the plane
    (Bfitted, P, inliers) = ransacfitplane(XYZ, t, feedback)

    Bfitted = Bfitted/Bfitted[4]
    @printf("Original plane coefficients:  ")
    @printf("%8.3f %8.3f %8.3f %8.3f \n",B[1], B[2], B[3], B[4])
    @printf("Fitted plane coefficients:    ")
    @printf("%8.3f %8.3f %8.3f %8.3f \n",Bfitted[1], Bfitted[2], Bfitted[3], Bfitted[4])

    # Display the triangular patch formed by the 3 points that gave the
    # plane of maximum consensus
#    hold(true)
    pts = [P P[:,1]]'
    plot3D(pts[:,1], pts[:,2], pts[:,3], "k-")
#    hold(false)

    @printf("\nRotate image so that the triangular patch is seen edge on\n")
    @printf("These are the points that form the plane of max consensus.\n\n")

end

#-----------------------------------------------------------------------
"""
fitfunddemo - Example of fundamental matrix computation

Demonstration of feature matching via simple correlation, and then using
RANSAC to estimate the fundamental matrix and at the same time identify
(mostly) inlying matches
```
Usage:  fitfunddemo           - Demonstrates fundamental matrix calculation
                                on two default images.
        funddemo(img1,img2)   - Computes fundamental matrix on two supplied
                                images.
```
Edit code as necessary to tweak parameters

See also: ransacfitfundmatrix(), fundmatrix()
"""
function fitfunddemo(img1=[], img2=[])

    if isempty(img1)
	#img1 = PyPlot.imread("img02.jpg")
	#img2 = PyPlot.imread("img03.jpg")
        img1 = float(load("img02.jpg"))
	img2 = float(load("img03.jpg"))
    end

    (rows1,cols1) = size(img1)
    (rows2,cols2) = size(img2)

    nonmaxrad = 3  # Non-maximal suppression radius
    dmax = 100     # Maximum search distance for matching
    w = 11         # Window size for correlation matching

    # Find Shi Tomasi corners in image1 and image2
    (cim1, r1, c1) = shi_tomasi(img1, 1, radius=nonmaxrad, N=200, img=img1, fig=1)
    (cim2, r2, c2) = shi_tomasi(img2, 1, radius=nonmaxrad, N=200, img=img2, fig=2)
    figure(1); axis("off");
    figure(2); axis("off");
    keypause()

    @printf("Matching features...\n")
    (m1,m2) = matchbycorrelation(copy(img1), [r1';c1'], copy(img2), [r2';c2'], w, dmax)

    # Display putative matches
    figure(3); clf(); imshow(img1); axis("off")  # hold(true)
    nMatch = size(m1,2)
    for n = 1:nMatch
	plot([m1[2,n], m2[2,n]], [m1[1,n], m2[1,n]],"r-")
    end
    title("Putative matches")
#    hold(false)

    # Assemble homogeneous feature coordinates for fitting of the
    # fundamental matrix, note that [x,y] corresponds to [col, row]
    x1 = [m1[2:2,:]; m1[1:1,:]; ones(1,nMatch)]
    x2 = [m2[2:2,:]; m2[1:1,:]; ones(1,nMatch)]

    t = .002  # Distance threshold for deciding outliers

    # Change the commenting on the lines below to switch between the 8
    # point fundamental matrix solution or the affine fundamental
    # matrix solution.
    (F, inliers) = ransacfitfundmatrix(x1, x2, t, true)
#   (F, inliers) = ransacfitaffinefundmatrix(x1, x2, t, true)

    @printf("Number of inliers was %d (%d%%) \n",
	    length(inliers),round(Int, 100*length(inliers)/nMatch))
    @printf("Number of putative matches was %d \n", nMatch)

    # Display both images overlayed with inlying matched feature points
    figure(4); clf();  imshow(img1);  axis("off") # hold(true)
    plot(m1[2,inliers], m1[1,inliers],"r+")
    plot(m2[2,inliers], m2[1,inliers],"g+")
    for n = inliers
	plot([m1[2,n], m2[2,n]], [m1[1,n], m2[1,n]], "b-")
    end
    title("Inlying matches")
#    hold(false)

    @printf("Step through each epipolar line [y/n]?\n")
    response = readline()
    if response[1] == 'n'
	return
    end

    # Step through each matched pair of points and display the
    # corresponding epipolar lines on the two images.
    l2 = F*x1    # Epipolar lines in image2
    l1 = F'*x2   # Epipolar lines in image1

    # Solve for epipoles
    (U,D,V) = svd(F)
    e1 = hnormalise(V[:,3])
    e2 = hnormalise(U[:,3])

    for n = inliers
	figure(1); clf(); imshow(img1)
        axis([1,cols1,rows1,1])
#        hold(true)
        plot(x1[1,n],x1[2,n],"r+",markersize=8)
	hline(l1[:,n]);
        plot(e1[1], e1[2], "g*");
#        hold(false)

	figure(2); clf(); imshow(img2);
        axis([1,cols2,rows2,1])
 #       hold(true);
        plot(x2[1,n],x2[2,n],"r+",markersize=8)
	hline(l2[:,n]);
        plot(e2[1], e2[2], "g*")
#        hold(false)
	keypause()
    end

    @printf("                                         \n")

end

#-----------------------------------------------------------------------
"""
fithomogdemo - Example of finding a homography

Demonstration of feature matching via simple correlation, and then using
RANSAC to estimate the homography between two images and at the same time
identify (mostly) inlying matches
```
Usage:  fithomogdemo            - Demonstrates homography calculation on two
                                  default images
        fithomogdemo(img1,img2) - Computes homography on two supplied images

Edit code as necessary to tweak parameters
```
See also: ransacfithomography(), homography2d()
"""
function fithomogdemo(img1=[], img2=[])

    if isempty(img1)
	#img1 = PyPlot.imread("boats.jpg")
	#img2 = PyPlot.imread("boatsrot.jpg")
        img1 = float(load("boats.jpg"))
	img2 = float(load("boatsrot.jpg"))
    end

    nonmaxrad = 3  # Non-maximal suppression radius
    dmax = 100     # Maximum search distance for matching
    w = 11         # Window size for correlation matching

    # Find Shi Tomasi corners in image1 and image2
    (cim1, r1, c1) = shi_tomasi(img1, 1, radius=nonmaxrad, N=100, img=img1, fig=1)
    (cim2, r2, c2) = shi_tomasi(img2, 1, radius=nonmaxrad, N=100, img=img2, fig=2)
    keypause()
    @printf("Matching features...\n")
    (m1,m2) = matchbycorrelation(img1, [r1';c1'], img2, [r2';c2'], w, dmax)

    # Display putative matches
    figure(3); clf();  imshow(img1);  # hold(true)
    nMatch = size(m1,2)
    for n = 1:nMatch
	plot([m1[2,n], m2[2,n]], [m1[1,n], m2[1,n]],"r-")
    end
    title("Putative matches")
#    hold(false)

    # Assemble homogeneous feature coordinates for fitting of the
    # fundamental matrix, note that [x,y] corresponds to [col, row]
    x1 = [m1[2:2,:]; m1[1:1,:]; ones(1,nMatch)]
    x2 = [m2[2:2,:]; m2[1:1,:]; ones(1,nMatch)]

    t = .001  # Distance threshold for deciding outliers
    (H, inliers) = ransacfithomography(x1, x2, t)

    @printf("Number of inliers was %d (%d%%) \n",
	    length(inliers),round(Int, 100*length(inliers)/nMatch))
    @printf("Number of putative matches was %d \n", nMatch)

    # Display both images overlayed with inlying matched feature points
    figure(4); clf();  imshow(img1);  # hold(true)
    plot(m1[2,inliers], m1[1,inliers],"r+")
    plot(m2[2,inliers], m2[1,inliers],"g+")
    for n = inliers
	plot([m1[2,n], m2[2,n]], [m1[1,n], m2[1,n]], "b-")
    end

    title("Inlying matches")
#    hold(false)
end


text:
#### Coordinates

export Coordinates, coords
export h_extreme, boundary_pts, over_points
export DataLocation, Center1D, Vertex1D
export CellBoundary, NodeBoundary
export All, Ghost, Boundary, Interior

"""
    DataLocation

Subtypes are used for dispatching
on the data location
"""
abstract type DataLocation end

""" 1D cell center """
struct Center1D <: DataLocation end

""" 1D cell vertex """
struct Vertex1D <: DataLocation end

abstract type CollocatedCoordinates{DL,FT} end

struct Coordinates1D{DL,FT,VFT,IT} <: CollocatedCoordinates{DL,FT}
  s::IT
  h::VFT
  Δh::VFT
end

function Base.show(io::IO, c::CollocatedCoordinates)
  print(io, "\n")
  println(io, "  s     = ",c.s)
  println(io, "  h     = ",c.h)
  println(io, "  Δh    = ",c.Δh)
end

struct SparseMat{D,U}
  diag::D
  upper::U
end

"""
    Coordinates

Cell node (`n`) and center (`c`)
1D coordinates.
"""
struct Coordinates{N, C, I, D}
  n::N
  c::C
  interp::I
  deriv::D
end

coords(c::Coordinates, ::Center1D) = c.c
coords(c::Coordinates, ::Vertex1D) = c.n

Base.length(c::Coordinates, ::Center1D) = c.c.s
Base.length(c::Coordinates, ::Vertex1D) = c.n.s

include("extremum_funcs.jl")

abstract type BoundaryPoints end

""" Cell center boundary type, containing ghost and first interior indices """
struct CellBoundary{IT<:Integer} <: BoundaryPoints
  "ghost index"
  g::IT
  "interior index"
  i::IT
end

""" Node boundary type, containing ghost, boundary, and first interior indices """
struct NodeBoundary{IT<:Integer} <: BoundaryPoints
  "ghost index"
  g::IT
  "boundary index"
  b::IT
  "interior index"
  i::IT
end

boundary_pts(::Coordinates1D{Center1D}, ::Min) = CellBoundary(1, 2)
boundary_pts(::Coordinates1D{Vertex1D}, ::Min) = NodeBoundary(1, 2, 3)

boundary_pts(c::Coordinates1D{Center1D}, ::Max) = CellBoundary(c.s, c.s-1)
boundary_pts(n::Coordinates1D{Vertex1D}, ::Max) = NodeBoundary(n.s, n.s-1, n.s-2)

h_extreme(n::Coordinates1D{Vertex1D}, e::Extremum) = n.h[boundary_pts(n,e).b]

abstract type DomainDecomp end
""" Type for dispatching on all points """
struct All <: DomainDecomp end
""" Type for dispatching on boundary points """
struct Boundary <: DomainDecomp end
""" Type for dispatching on interior points """
struct Interior <: DomainDecomp end
""" Type for dispatching on ghost points """
struct Ghost <: DomainDecomp end

function over_points(c::Coordinates, ::All     , dl::DataLocation)
  c = coords(c, dl)
  boundary_pts(c,Min()).g:boundary_pts(c,Max()).g
end
function over_points(c::Coordinates, ::Boundary, dl::DataLocation)
  c = coords(c, dl)
  boundary_pts(c,Min()).b:boundary_pts(c,Max()).b
end
function over_points(c::Coordinates, ::Interior, dl::DataLocation)
  c = coords(c, dl)
  boundary_pts(c,Min()).i:boundary_pts(c,Max()).i
end
function over_points(c::Coordinates, ::Ghost   , dl::DataLocation)
  c = coords(c, dl)
  b_pts_zmin = boundary_pts(c,Min())
  b_pts_zmax = boundary_pts(c,Max())
  return (b_pts_zmin.g:b_pts_zmin.g,
          b_pts_zmax.g:b_pts_zmax.g)
end

function pad!(h::Vector{FT}) where FT
  pushfirst!(h, h[1] - (h[2]-h[1]))
  push!(h, h[end] + (h[end]-h[end-1]))
end

function Coordinates(a::FT, b::FT, n::IT; warpfun::F=nothing, args=nothing) where {FT, IT<:Integer, F}
  sn = n+3

  if !(warpfun==nothing)
    hn = warpfun(a, b, n, args...)
  else
    hn = uniform(a,b,n)
  end
  pad!(hn)
  @assert length(hn) == n+3
  Δhn = FT[hn[i+1]-hn[i] for i in 1:sn-1]
  cn = Coordinates1D{Vertex1D,FT,typeof(hn),typeof(sn)}(sn, hn, Δhn)

  sc = n+2
  hc = FT[hn[i]+Δhn[i]/2 for i in 1:sc]
  @assert length(hc) == n+2
  Δhc = FT[hc[i+1]-hc[i] for i in 1:sc-1]
  cc = Coordinates1D{Center1D,FT,typeof(hc),typeof(sc)}(sc, hc, Δhc)

  if !(all(isfinite.(hn)) && all(isfinite.(hc)))
    throw(ArgumentError("Coordinates are not finite."))
  end

  interp = SparseMat(interp_diag(hn, hc, sc), interp_upper_diag(hn, hc, sc))
  deriv = (
    c2v=SparseMat(deriv_diag_C2V(Δhc, sc), deriv_upper_diag_C2V(Δhc, sc)),
    v2c=SparseMat(deriv_diag_V2C(Δhn, sn), deriv_upper_diag_V2C(Δhn, sn)),
    )
  return Coordinates(cn, cc, interp, deriv)
end

interp_diag(hn, hc, sc) =
    [(hn[i+1] - hc[i])/(hc[i+1] - hc[i]) for i in 1:sc-1]

interp_upper_diag(hn, hc, sc) = 1 .- interp_diag(hn, hc, sc)

deriv_diag_C2V(Δh, sc) = [-1/Δh[i] for i in 1:sc-1]
deriv_upper_diag_C2V(Δh, sc) = [1/Δh[i] for i in 1:sc-1]

deriv_diag_V2C(Δh, sc) = [-1/Δh[i] for i in 1:sc-1]
deriv_upper_diag_V2C(Δh, sc) = [1/Δh[i] for i in 1:sc-1]


text:
using DoseCalculators

include "medication_formulas.jl"

dose_calculator_app(amoxicillin_po, "Amoxicillin PO Dosing", "Amoxicillin dosage(mg)")