text:
# This is the high level interface for building dependencies using the declarative BinDeps Interface

# A dependency provider, if succcessfully exectued will satisfy the dependency
abstract DependencyProvider

# A library helper may be used by `DependencyProvider`s but will by iteself not provide the library
abstract DependencyHelper

type PackageContext
    do_install::Bool
    dir::String
    package::String
    deps::Vector{Any}
end

type LibraryDependency
    name::String
    context::PackageContext
    providers::Vector{(DependencyProvider,Dict{Symbol,Any})}
    helpers::Vector{(DependencyHelper,Dict{Symbol,Any})}
    properties::Dict{Symbol,Any}
    libvalidate::Function
end

type LibraryGroup
    name::String
    deps::Vector{LibraryDependency}
end

import Base: show

function show(io::IO, dep::LibraryDependency)
    print(io,"LibraryDependency(",dep.name,")")
end

# Default directory organization
pkgdir(dep) = dep.context.dir
depsdir(dep) = joinpath(pkgdir(dep),"deps")
usrdir(dep) = joinpath(depsdir(dep),"usr")
libdir(dep) = joinpath(usrdir(dep),"lib")
bindir(dep) = joinpath(usrdir(dep),"bin")
includedir(dep) = joinpath(usrdir(dep),"include")
builddir(dep) = joinpath(depsdir(dep),"builds")
downloadsdir(dep) = joinpath(depsdir(dep),"downloads")
srcdir(dep) = joinpath(depsdir(dep),"src")
libdir(provider, dep) = [libdir(dep), libdir(dep)*"32", libdir(dep)*"64"]
bindir(provider, dep) = bindir(dep)

successful_validate(l,p) = true

function _library_dependency(context::PackageContext, name; properties...)
    validate = successful_validate
    group = nothing
    for i in 1:length(properties)
        k,v = properties[i]
        if k == :validate
            validate = v
            splice!(properties,i)
        end
        if k == :group
            group = v
        end
    end
    r = LibraryDependency(name,context,Array((DependencyProvider,Dict{Symbol,Any}),0),Array((DependencyHelper,Dict{Symbol,Any}),0),(Symbol=>Any)[name => value for (name,value) in properties],validate)
    if group !== nothing
        push!(group.deps,r)
    else
        push!(context.deps,r)
    end
    r
end

function _library_group(context,name)
    r = LibraryGroup(name,LibraryDependency[])
    push!(context.deps,r)
    r
end

# This macro expects to be the first thing run. It attempts to deduce the package name and initializes the context
macro setup()
    dir = normpath(joinpath(pwd(),".."))
    package = basename(dir)
    esc(quote
        if length(ARGS) > 0 && isa(ARGS[1],BinDeps.PackageContext)
            bindeps_context = ARGS[1]
        else
            bindeps_context = BinDeps.PackageContext(true,$dir,$package,{})
        end
        library_group(args...) = BinDeps._library_group(bindeps_context,args...)
        library_dependency(args...; properties...) = BinDeps._library_dependency(bindeps_context,args...;properties...)
    end)
end

export library_dependency, bindir, srcdir, usrdir, libdir

library_dependency(args...; properties...) = error("No context provided. Did you forget `@BinDeps.setup`?")

abstract PackageManager <: DependencyProvider

DEBIAN_VERSION_REGEX = r"^
    ([0-9]+\:)?                                           # epoch
    (?:(?:([0-9][a-z0-9.\-+:~]*)-([0-9][a-z0-9.+~]*)) |   # upstream version + debian revision
          ([0-9][a-z0-9.+:~]*))                           # upstream version
"ix

const has_apt = try success(`apt-get -v`) catch e false end
type AptGet <: PackageManager
    package::String
end
can_use(::Type{AptGet}) = has_apt && OS_NAME == :Linux
package_available(p::AptGet) = can_use(AptGet) && !isempty(available_versions(p))
function available_versions(p::AptGet)
    vers = ASCIIString[]
    lookfor_version = false
    for l in eachline(`apt-cache showpkg $(p.package)`)
        if beginswith(l,"Version:")
            try
                vs = l[(1+length("Version: ")):end]
                push!(vers, vs)
            end
        elseif lookfor_version && (m = match(DEBIAN_VERSION_REGEX, l)) != nothing
            m.captures[2] != nothing ? push!(vers, m.captures[2]) :
                                       push!(vers, m.captures[4])
        elseif beginswith(l, "Versions:")
            lookfor_version = true
        elseif beginswith(l, "Reverse Depends:")
            lookfor_version = false
        end
    end
    return vers
end
function available_version(p::AptGet)
    vers = available_versions(p)
    isempty(vers) && error("apt-cache did not return version information. This shouldn't happen. Please file a bug!")
    length(vers) > 1 && warn("Multiple versions of $(p.package) are available.  Use BinDeps.available_versions to get all versions.")
    return vers[end]
end
pkg_name(a::AptGet) = a.package

libdir(p::AptGet,dep) = ["/usr/lib", "/usr/lib64", "/usr/lib32", "/usr/lib/x86_64-linux-gnu", "/usr/lib/i386-linux-gnu"]

const has_yum = try success(`yum --version`) catch e false end
type Yum <: PackageManager
    package::String
end
can_use(::Type{Yum}) = has_yum && OS_NAME == :Linux
package_available(y::Yum) = can_use(Yum) && success(`yum list $(y.package)`)
function available_version(y::Yum)
    uname = readchomp(`uname -m`)
    found_uname = false
    found_version = false
    for l in eachline(`yum info $(y.package)`)
        l = chomp(l)
        if !found_uname
            # On 64-bit systems, we may have multiple arches installed
            # this makes sure we get the right one
            found_uname = endswith(l, uname)
            continue
        end
        if beginswith(l, "Version")
            return convert(VersionNumber, split(l)[end])
        end
    end
    error("yum did not return version information.  This shouldn't happen. Please file a bug!")
end
pkg_name(y::Yum) = y.package

# Pacman/Yaourt are package managers for Arch Linux.

# Note that `pacman --version` has an unreliable return value.
const has_pacman = try success(`pacman -Qq`) catch e false end
type Pacman <: PackageManager
    package::String
end
can_use(::Type{Pacman}) = has_pacman && OS_NAME == :Linux
package_available(p::Pacman) = can_use(Pacman) && success(`pacman -Si $(p.package)`)
# Only one version is usually available via pacman, hence no `available_versions`.
function available_version(p::Pacman)
    for l in eachline(`/usr/bin/pacman -Si $(p.package)`) # To circumvent alias problems
        if beginswith(l, "Version")
            # The following isn't perfect, but it's hopefully less brittle than
            # writing a regex for pacman's nonexistent version-string standard.

            # This also strips away the sometimes leading epoch as in ffmpeg's
            # Version        : 1:2.3.3-1
            versionstr = strip(split(l, ":")[end])
            try
                return convert(VersionNumber, versionstr)
            catch e
                # For too long versions like imagemagick's 6.8.9.6-1, give it
                # a second try just discarding superfluous stuff.
                return convert(VersionNumber, join(split(versionstr, '.')[1:3], '.'))
            end
        end
    end
    error("pacman did not return version information. This shouldn't happen. Please file a bug!")
end
pkg_name(p::Pacman) = p.package

libdir(p::Pacman,dep) = ["/usr/lib", "/usr/lib32"]

# The nix pacakge manager is available on most UNIX platforms

const has_nix = try success(`nix-env --version`) catch e false end
type Nix <: PackageManager
    package::String
end
can_use(::Type{Nix}) = has_nix && OS_NAME == :Linux
package_available(p::Nix) = can_use(Nix) && !isempty(available_versions(p))
function available_versions(p::Nix)
    vers = ASCIIString[]
    for l in eachline(`nix-env -qa` |> `grep $(p.package)`)
		parts = split(l, "-")
        if length(parts)>1
			push!(vers, parts[2])
		end
    end
    return vers
end
function available_version(p::Nix)
    vers = available_versions(p)
    isempty(vers) && error("nix-env did not return version information. This shouldn't happen. Please file a bug!")
    length(vers) > 1 && warn("Multiple versions of $(p.package) are available.  Use BinDeps.available_versions to get all versions.")
    return vers[1]  # to be consistent with nix-env behaviour we install the first encountered version
end
pkg_name(a::Nix) = a.package

libdir(p::Nix,dep) = ["$(homedir())/.nix-profile/lib", "$(homedir())/.nix-profil/lib64"]

# Can use everything else without restriction by default
can_use(::Type) = true

abstract Sources <: DependencyHelper
abstract Binaries <: DependencyProvider

#
# A dummy provider checked for every library that
# indicates the library was found somewhere on the
# system using dlopen.
#
immutable SystemPaths <: Binaries; end

show(io::IO, ::SystemPaths) = print(io,"System Paths")

using URIParser
export URI

type NetworkSource <: Sources
    uri::URI
end

srcdir(s::Sources, dep::LibraryDependency) = srcdir(dep,s,(Symbol=>Any)[])
function srcdir( dep::LibraryDependency, s::NetworkSource,opts)
    joinpath(srcdir(dep),get(opts,:unpacked_dir,splittarpath(basename(s.uri.path))[1]))
end

type RemoteBinaries <: Binaries
    uri::URI
end

type CustomPathBinaries <: Binaries
    path::String
end

libdir(p::CustomPathBinaries,dep) = p.path

abstract BuildProcess <: DependencyProvider

type SimpleBuild <: BuildProcess
    steps
end

type Autotools <: BuildProcess
    source
    opts
end

type GetSources <: BuildStep
    dep::LibraryDependency
end

lower(x::GetSources,collection) = push!(collection,generate_steps(x.dep,gethelper(x.dep,Sources)...))

Autotools(;opts...) = Autotools(nothing,{k => v for (k,v) in opts})

export AptGet, Yum, Pacman, Nix, Sources, Binaries, provides, BuildProcess, Autotools, GetSources, SimpleBuild, available_version

provider{T<:PackageManager}(::Type{T},package::String; opts...) = T(package)
provider(::Type{Sources},uri::URI; opts...) = NetworkSource(uri)
provider(::Type{Binaries},uri::URI; opts...) = RemoteBinaries(uri)
provider(::Type{Binaries},path::String; opts...) = CustomPathBinaries(path)
provider(::Type{SimpleBuild},steps; opts...) = SimpleBuild(steps)
provider{T<:BuildProcess}(::Type{BuildProcess},p::T; opts...) = provider(T,p; opts...)
provider(::Type{BuildProcess},steps::Union(BuildStep,SynchronousStepCollection); opts...) = provider(SimpleBuild,steps; opts...)
provider(::Type{Autotools},a::Autotools; opts...) = a

provides(provider::DependencyProvider,dep::LibraryDependency; opts...) = push!(dep.providers,(provider,(Symbol=>Any)[k=>v for (k,v) in opts]))
provides(helper::DependencyHelper,dep::LibraryDependency; opts...) = push!(dep.helpers,(helper,(Symbol=>Any)[k=>v for (k,v) in opts]))
provides{T}(::Type{T},p,dep::LibraryDependency; opts...) = provides(provider(T,p; opts...),dep; opts...)
function provides{T}(::Type{T},ps,deps::Vector{LibraryDependency}; opts...)
    p = provider(T,ps; opts...)
    for dep in deps
        provides(p,dep; opts...)
    end
end

function provides{T}(::Type{T},providers::Dict; opts...)
    for (k,v) in providers
        provides(T,k,v;opts...)
    end
end

generate_steps(h::DependencyProvider,dep::LibraryDependency) = error("Must also pass provider options")
generate_steps(h::BuildProcess,dep::LibraryDependency,opts) = h.steps
function generate_steps(dep::LibraryDependency,h::AptGet,opts)
    if get(opts,:force_rebuild,false)
        error("Will not force apt-get to rebuild dependency \"$(dep.name)\".\n"*
              "Please make any necessary adjustments manually (This might just be a version upgrade)")
    end
    @build_steps begin
        println("Installing dependency $(h.package) via `sudo apt-get install $(h.package)`:")
        `sudo apt-get install $(h.package)`
        ()->(ccall(:jl_read_sonames,Void,()))
    end
end
function generate_steps(dep::LibraryDependency,h::Yum,opts)
    if get(opts,:force_rebuild,false)
        error("Will not force yum to rebuild dependency \"$(dep.name)\".\n"*
              "Please make any necessary adjustments manually (This might just be a version upgrade)")
    end

    @build_steps begin
        println("Installing dependency $(h.package) via `sudo yum install $(h.package)`:")
        `sudo yum install $(h.package)`
        ()->(ccall(:jl_read_sonames,Void,()))
    end
end
function generate_steps(dep::LibraryDependency,h::Pacman,opts)
    if get(opts,:force_rebuild,false)
        error("Will not force pacman to rebuild dependency \"$(dep.name)\".\n"*
              "Please make any necessary adjustments manually (This might just be a version upgrade)")
    end
    @build_steps begin
        println("Installing dependency $(h.package) via `sudo pacman -S --needed $(h.package)`:")
        `sudo pacman -S --needed $(h.package)`
        ()->(ccall(:jl_read_sonames,Void,()))
    end
end
function generate_steps(dep::LibraryDependency,h::Nix,opts)
    if get(opts,:force_rebuild,false)
        error("Will not force nix to rebuild dependency \"$(dep.name)\".\n"*
              "Please make any necessary adjustments manually (This might just be a version upgrade)")
    end
    @build_steps begin
        println("Installing dependency $(h.package) via `nix-env -i $(h.package)`:")
        `nix-env -i $(h.package)`
        ()->(ccall(:jl_read_sonames,Void,()))
    end
end
function generate_steps(dep::LibraryDependency,h::NetworkSource,opts)
    localfile = joinpath(downloadsdir(dep),basename(h.uri.path))
    @build_steps begin
        FileDownloader(string(h.uri),localfile)
        ChecksumValidator(get(opts,:SHA,""),localfile)
        CreateDirectory(srcdir(dep))
        FileUnpacker(localfile,srcdir(dep),srcdir(dep,h,opts))
    end
end
function generate_steps(dep::LibraryDependency,h::RemoteBinaries,opts)
    get(opts,:force_rebuild,false) && error("Force rebuild not allowed for binaries. Use a different download location instead.")
    localfile = joinpath(downloadsdir(dep),basename(h.uri.path))
    steps = @build_steps begin
        FileDownloader(string(h.uri),localfile)
        ChecksumValidator(get(opts,:SHA,""),localfile)
        FileUnpacker(localfile,depsdir(dep),get(opts,:unpacked_dir,"usr"))
    end
end
generate_steps(dep::LibraryDependency,h::SimpleBuild,opts) = h.steps

function getoneprovider(dep::LibraryDependency,method)
    for (p,opts) = dep.providers
        if typeof(p) <: method && can_use(typeof(p))
            return (p,opts)
        end
    end
    return (nothing,nothing)
end

function getallproviders(dep::LibraryDependency,method)
    ret = {}
    for (p,opts) = dep.providers
        if typeof(p) <: method && can_use(typeof(p))
            push!(ret,(p,opts))
        end
    end
    ret
end

function gethelper(dep::LibraryDependency,method)
    for (p,opts) = dep.helpers
        if typeof(p) <: method
            return (p,opts)
        end
    end
    return (nothing,nothing)
end

function generate_steps(dep::LibraryDependency,method)
    (p,opts) = getoneprovider(dep,method)
    !is(p,nothing) && return generate_steps(p,dep,opts)
    (p,hopts) = gethelper(dep,method)
    !is(p,nothing) && return generate_steps(p,dep,hopts)
    error("No provider or helper for method $method found for dependency $(dep.name)")
end

function generate_steps(dep::LibraryDependency, h::Autotools,  provider_opts)
    if is(h.source, nothing)
        h.source = gethelper(dep,Sources)
    end
    if isa(h.source,Sources)
        h.source = (h.source,(Symbol=>Any)[])
    end
    is(h.source[1], nothing) && error("Could not obtain sources for dependency $(dep.name)")
    steps = lower(generate_steps(dep,h.source...))
    opts = {:srcdir=>srcdir(dep,h.source...), :prefix=>usrdir(dep), :builddir=>joinpath(builddir(dep),dep.name)}
    merge!(opts,h.opts)
    if haskey(opts,:installed_libname)
        !haskey(opts,:installed_libpath) || error("Can't specify both installed_libpath and installed_libname")
        opts[:installed_libpath] = ByteString[joinpath(libdir(dep),opts[:installed_libname])]
        delete!(opts, :installed_libname)
    elseif !haskey(opts,:installed_libpath)
        opts[:installed_libpath] = ByteString[joinpath(libdir(dep),x)*"."*dlext for x in get(dep.properties,:aliases,ByteString[])]
    end
    if !haskey(opts,:libtarget) && haskey(dep.properties,:aliases)
        opts[:libtarget] = ByteString[x*"."*dlext for x in dep.properties[:aliases]]
    end
    if !haskey(opts,:include_dirs)
        opts[:include_dirs] = String[]
    end
    if !haskey(opts,:lib_dirs)
        opts[:lib_dirs] = String[]
    end
    if !haskey(opts,:pkg_config_dirs)
        opts[:pkg_config_dirs] = String[]
    end
    if !haskey(opts,:rpath_dirs)
        opts[:rpath_dirs] = String[]
    end
    if haskey(opts,:configure_subdir)
        opts[:srcdir] = joinpath(opts[:srcdir],opts[:configure_subdir])
        delete!(opts, :configure_subdir)
    end
    unshift!(opts[:include_dirs],includedir(dep))
    unshift!(opts[:lib_dirs],libdir(dep))
    unshift!(opts[:rpath_dirs],libdir(dep))
    unshift!(opts[:pkg_config_dirs],joinpath(libdir(dep),"pkgconfig"))
    env = Dict{ByteString,ByteString}()
    env["PKG_CONFIG_PATH"] = join(opts[:pkg_config_dirs],":")
    delete!(opts,:pkg_config_dirs)
    @unix_only env["PATH"] = bindir(dep)*":"*ENV["PATH"]
    @windows_only env["PATH"] = bindir(dep)*";"*ENV["PATH"]
    haskey(opts,:env) && merge!(env,opts[:env])
    opts[:env] = env
    if get(provider_opts,:force_rebuild,false)
        opts[:force_rebuild] = true
    end
    steps |= AutotoolsDependency(;opts...)
    steps
end

@osx_only const EXTENSIONS = ["",".dylib"]
@windows_only const EXTENSIONS = ["", ".dll"]
@linux_only const EXTENSIONS = ["",".so"]

#
# Finds all copies of the library on the system, listed in preference order.
# Return value is an array of tuples of the provider and the path where it is found
#
function _find_library(dep::LibraryDependency; provider = Any)
    ret = Any[]
    # Same as find_library, but with extra check defined by dep
    libnames = [dep.name;get(dep.properties,:aliases,ASCIIString[])]
    # Make sure we keep the defaults first, but also look in the other directories
    providers = unique([reduce(vcat,[getallproviders(dep,p) for p in defaults]),dep.providers])
    for (p,opts) in providers
        (p != nothing && can_use(typeof(p)) && can_provide(p,opts,dep)) || continue
        paths = String[]

        # Allow user to override installation path
        if haskey(opts,:installed_libpath) && isdir(opts[:installed_libpath])
            unshift!(paths,opts[:installed_libpath])
        end

        ppaths = libdir(p,dep)
        append!(paths,isa(ppaths,Array) ? ppaths : [ppaths])

        if haskey(opts,:unpacked_dir) && isdir(joinpath(depsdir(dep),opts[:unpacked_dir]))
            push!(paths,joinpath(depsdir(dep),opts[:unpacked_dir]))
        end

        # Windows, do you know what `lib` stands for???
        @windows_only push!(paths,bindir(p,dep))
        (isempty(paths) || all(map(isempty,paths))) && continue
        for lib in libnames, path in paths
            l = joinpath(path, lib)
            h = dlopen_e(l, RTLD_LAZY)
            if h != C_NULL
                works = dep.libvalidate(l,h)
                if VERSION >= v"0.3-"
                    l = Sys.dlpath(h)
                end
                dlclose(h)
                if works
                    push!(ret, ((p, opts), l))
                else
                    # We tried to load this providers' library, but it didn't satisfy
                    # the requirements, so tell it to force a rebuild since the requirements
                    # have most likely changed
                    opts[:force_rebuild] = true
                end
            end
        end
    end
    # 0.2 compatibility
    if VERSION < v"0.3-"
        for lib in libnames
            p = dlopen_e(lib, RTLD_LAZY)
            if p != C_NULL
                works = dep.libvalidate(lib,p)
                dlclose(p)
                if works
                    push!(ret,((SystemPaths(),(Any=>Any)[]),lib))
                end
            end
        end
    else
        # Now check system libraries
        for lib in libnames
            # We don't want to use regular dlopen, because we want to get at
            # system libraries even if one of our providers is higher in the
            # DL_LOAD_PATH
            for path in Base.DL_LOAD_PATH
                for ext in EXTENSIONS
                    opath = string(joinpath(path,lib),ext)
                    check_path!(ret,dep,opath)
                end
            end
            for ext in EXTENSIONS
                opath = string(lib,ext)
                check_path!(ret,dep,opath)
            end
            @linux_only begin
                soname = ccall(:jl_lookup_soname, Ptr{Uint8}, (Ptr{Uint8}, Csize_t), lib, sizeof(lib))
                soname != C_NULL && check_path!(ret,dep,bytestring(soname))
            end
        end
    end
    return ret
end

if VERSION >= v"0.3-"
    function check_path!(ret,dep,opath)
        flags = RTLD_LAZY
        handle = c_malloc(2*sizeof(Ptr{Void}))
        err = ccall(:jl_uv_dlopen,Cint,(Ptr{Uint8},Ptr{Void},Cuint),opath,handle,flags)
        if err == 0
            check_system_handle!(ret,dep,handle)
            dlclose(handle)
            c_free(handle)
        end
    end

    function check_system_handle!(ret,dep,handle)
        if handle != C_NULL
            libpath = Sys.dlpath(handle)
            # Check that this is not a duplicate
            for p in ret
                try
                    if realpath(p[2]) == realpath(libpath)
                        return
                    end
                catch
                    warn("""
                        Found a library that does not exist.
                        This may happen if the library has an active open handle.
                        Please quit julia and try again.
                        """)
                    return
                end
            end
            works = dep.libvalidate(libpath,handle)
            if works
                push!(ret, ((SystemPaths(),(Any=>Any)[]), libpath))
            end
        end
    end
end

# Default installation method
if OS_NAME == :Darwin
    defaults = [Binaries,PackageManager,SystemPaths,BuildProcess]
elseif OS_NAME == :Linux
    defaults = [PackageManager,SystemPaths,BuildProcess]
elseif OS_NAME == :Windows
    defaults = [Binaries,PackageManager,SystemPaths]
else
    defaults = [SystemPaths,BuildProcess]
end

function applicable(dep::LibraryDependency)
    if haskey(dep.properties,:os)
        if (dep.properties[:os] != OS_NAME && dep.properties[:os] != :Unix) || (dep.properties[:os] == :Unix && !Base.is_unix(OS_NAME))
            return false
        end
    elseif haskey(dep.properties,:runtime) && dep.properties[:runtime] == false
        return false
    end
    return true
end

applicable(deps::LibraryGroup) = any([applicable(dep) for dep in deps.deps])

function can_provide(p,opts,dep)
    if p === nothing || (haskey(opts,:os) && opts[:os] != OS_NAME && (opts[:os] != :Unix || !Base.is_unix(OS_NAME)))
        return false
    end
    if !haskey(opts,:validate)
        return true
    elseif isa(opts[:validate],Bool)
        return opts[:validate]
    else
        return opts[:validate](p,dep)
    end
end

function can_provide(p::PackageManager,opts,dep)
    if p === nothing || (haskey(opts,:os) && opts[:os] != OS_NAME && (opts[:os] != :Unix || !Base.is_unix(OS_NAME)))
        return false
    end
    if !package_available(p)
        return false
    end
    if !haskey(opts,:validate)
        return true
    elseif isa(opts[:validate],Bool)
        return opts[:validate]
    else
        return opts[:validate](p,dep)
    end
end

issatisfied(dep::LibraryDependency) = !isempty(_find_library(dep))

allf(deps) = [dep => _find_library(dep) for dep in deps.deps]
function satisfied_providers(deps::LibraryGroup, allfl = allf(deps))
    viable_providers = nothing
    for dep in deps.deps
        if !applicable(dep)
            continue
        end
        providers = map(x->typeof(x[1][1]),allfl[dep])
        if viable_providers == nothing
            viable_providers = providers
        else
            viable_providers = intersect(viable_providers,providers)
        end
    end
    viable_providers
end

function viable_providers(deps::LibraryGroup)
    vp = nothing
    for dep in deps.deps
        if !applicable(dep)
            continue
        end
        providers = map(x->typeof(x[1]),dep.providers)
        if vp == nothing
            vp = providers
        else
            vp = intersect(vp,providers)
        end
    end
    vp
end

#
# We need to make sure all libraries are satisfied with the
# additional constraint that all of them are satisfied by the
# same provider.
#
issatisfied(deps::LibraryGroup) = !isempty(satisfied_providers(deps))

function _find_library(deps::LibraryGroup, allfl = allf(deps); provider = Any)
    providers = satisfied_providers(deps,allfl)
    p = nothing
    if isempty(providers)
        return (Any=>Any)[]
    else
        for p2 in providers
            if p2 <: provider
                p = p2
            end
        end
    end
    p === nothing && error("Given provider does not satisfy the library group")
    [dep => begin
        thisfl = allfl[dep]
        ret = nothing
        for fl in thisfl
            if isa(fl[1][1],p)
                ret = fl
                break
            end
        end
        @assert ret != nothing
        ret
    end for dep in filter(applicable,deps.deps)]
end

function satisfy!(deps::LibraryGroup, methods = defaults)
    sp = satisfied_providers(deps)
    if !isempty(sp)
        for m in methods
            for s in sp
                if s <: m
                    return s
                end
            end
        end
    end
    if !applicable(deps)
        return Any
    end
    vp = viable_providers(deps)
    didsatisfy = false
    for method in methods
        for p in vp
            if !(p <: method) || !can_use(p)
                continue
            end
            skip = false
            for dep in deps.deps
                !applicable(dep) && continue
                hasany = false
                for (p2,opts) in getallproviders(dep,p)
                    can_provide(p2, opts, dep) && (hasany = true)
                end
                if !hasany
                    skip = true
                    break
                end
            end
            if skip
                continue
            end
            for dep in deps.deps
                satisfy!(dep,[p])
            end
            return p
        end
    end
    error("""
        None of the selected providers could satisfy library group $(deps.name)
        Use BinDeps.debug(package_name) to see available providers
        """)
end

function satisfy!(dep::LibraryDependency, methods = defaults)
    sp = map(x->typeof(x[1][1]),_find_library(dep))
    if !isempty(sp)
        for m in methods
            for s in sp
                if s <: m
                    return s
                end
            end
        end
    end
    if !applicable(dep)
        return
    end
    for method in methods
        for (p,opts) in getallproviders(dep,method)
            can_provide(p,opts,dep) || continue
            if haskey(opts,:force_depends)
                for (dmethod,ddep) in opts[:force_depends]
                    (dp,dopts) = getallproviders(ddep,dmethod)[1]
                    run(lower(generate_steps(ddep,dp,dopts)))
                end
            end
            run(lower(generate_steps(dep,p,opts)))
            !issatisfied(dep) && error("Provider $method failed to satisfy dependency $(dep.name)")
            return p
        end
    end
    error("""
        None of the selected providers can install dependency $(dep.name).
        Use BinDeps.debug(package_name) to see available providers
        """)
end

execute(dep::LibraryDependency,method) = run(lower(generate_steps(dep,method)))

macro install (_libmaps...)
    if length(_libmaps) == 0
        return esc(quote
            if bindeps_context.do_install
                for d in bindeps_context.deps
                    BinDeps.satisfy!(d)
                end
            end
        end)
    else
        libmaps = eval(_libmaps[1])
        load_cache = gensym()
        ret = Expr(:block)
        push!(ret.args,
            esc(quote
                    load_cache = Dict()
                    pre_hooks = Set{String}()
                    load_hooks = Set{String}()
                    if bindeps_context.do_install
                        for d in bindeps_context.deps
                            p = BinDeps.satisfy!(d)
                            libs = BinDeps._find_library(d; provider = p)
                            if isa(d, BinDeps.LibraryGroup)
                                if !isempty(libs)
                                    for dep in d.deps
                                        !BinDeps.applicable(dep) && continue
                                        load_cache[dep.name] = libs[dep][2]
                                        opts = libs[dep][1][2]
                                        haskey(opts, :preload) && push!(pre_hooks,opts[:preload])
                                        haskey(opts, :onload) && push!(load_hooks,opts[:onload])
                                    end
                                end
                            else
                                for (k,v) in libs
                                    load_cache[d.name] = v
                                    opts = k[2]
                                    haskey(opts, :preload) && push!(pre_hooks,opts[:preload])
                                    haskey(opts, :onload) && push!(load_hooks,opts[:onload])
                                end
                            end
                        end
                    end

                    # Generate "deps.jl" file for runtime loading
                    depsfile = open(joinpath(splitdir(Base.source_path())[1],"deps.jl"), "w")
                    println(depsfile, join(pre_hooks, "\n"))
                    println(depsfile,
                        """
                        macro checked_lib(libname, path)
                            (dlopen_e(path) == C_NULL) && error("Unable to load \\n\\n\$libname (\$path)\\n\\nPlease re-run Pkg.build(package), and restart Julia.")
                            quote const \$(esc(libname)) = \$path end
                        end""")
                    for libkey in keys($libmaps)
                        ((cached = get(load_cache,string(libkey),nothing)) === nothing) && continue
                        println(depsfile, "@checked_lib ", $libmaps[libkey], " \"", escape_string(cached), "\"")
                    end
                    println(depsfile, join(load_hooks,"\n"))
                    close(depsfile)
                end))
        if !(typeof(libmaps) <: Associative)
            warn("Incorrect mapping in BinDeps.@install call. No dependencies will be cached.")
        end
        ret
    end
end

# Usage: @load_dependencies [file] [filter]
#
# Load dependencies as declared in `file` (default to "../deps/build.jl")
#
# This will also assign global variables in the current module that are the result of the symbol lookup
# The filter argument determines which libaries are being loaded and can be used to modify
# the name of the global variable to be set to the result of the lookup.
# The second argument may be as follows:
#
#  1. Vector{Symbol} or Vector{T <: String}
#       Only load that are declared whose name is listed in the Array
#       E.g. @load_dependencies "file.jl" [:cairo, :tk]
#
#  2. Associative{S<:Union(Symbol,String),S<:Union(Symbol,String)}
#       Only loads libraries whose name matches a key in the Associative collection, but assigns it
#       to the name matiching the corresponsing value
#       E.g. @load_dependencies "file.jl" [:cairo=>:libcairo, :tk=>:libtk]
#       will assign the result of the lookup for :cairo and :tk to the variables `libcairo` and `libtk`
#       respectively.
#
#  3. Function
#       A filter function
#       E.g. @load_dependencies "file.jl" x->x=="tk"
#
#
macro load_dependencies(args...)
    dir = dirname(normpath(joinpath(dirname(Base.source_path()),"..")))
    arg1 = nothing
    file = "../deps/build.jl"
    if length(args) == 1
        if isa(args[1],Expr)
            arg1 = eval(args[1])
        elseif typeof(args[1]) <: String
            file = args[1]
            dir = dirname(normpath(joinpath(dirname(file),"..")))
        elseif typeof(args[1]) <: Associative || isa(args[1],Vector)
            arg1 = args[1]
        else
            error("Type $(typeof(args[1])) not recognized for argument 1. See usage instructions!")
        end
    elseif length(args) == 2
        file = args[1]
        arg1 = typeof(args[2]) <: Associative || isa(args[2],Vector) ? args[2] : eval(args[2])
    elseif length(args) != 0
        error("No version of @load_dependencies takes $(length(args)) arguments. See usage instructions!")
    end
    pkg = ""
    r = search(dir,Pkg.Dir.path())
    if r != 0:-1
        s = search(dir,"/",last(r)+2)
        if s != 0:-1
            pkg = dir[(last(r)+2):(first(s)-1)]
        else
            pkg = dir[(last(r)+2):end]
        end
    end
    context = BinDeps.PackageContext(false,dir,pkg,{})
    m = Module(:__anon__)
    body = Expr(:toplevel,:(ARGS=[$context]),:(include($file)))
    eval(m,body)
    ret = Expr(:block)
    for dep in context.deps
        if !applicable(dep)
            continue
        end
        name = sym = dep.name
        if arg1 !== nothing
            if (typeof(arg1) <: Associative) && all(map(x->(x == Symbol || x <: String),eltype(arg1)))
                found = false
                for need in keys(arg1)
                    found = (dep.name == string(need))
                    if found
                        sym = arg1[need]
                        delete!(arg1,need)
                        break
                    end
                end
                if !found
                    continue
                end
            elseif isa(arg1,Vector) && ((eltype(arg1) == Symbol) || (eltype(arg1) <: String))
                found = false
                for i = 1:length(args)
                    found = (dep.name == string(arg1[i]))
                    if found
                        sym = arg1[i]
                        splice!(arg1,i)
                        break
                    end
                end
                if !found
                    continue
                end
            elseif isa(arg1,Function)
                if !f(name)
                    continue
                end
            else
                error("Can't deal with argument type $(typeof(arg1)). See usage instructions!")
            end
        end
        s = symbol(sym)
        errorcase = Expr(:block)
        push!(errorcase.args,:(error("Could not load library "*$(dep.name)*". Try running Pkg.build() to install missing dependencies!")))
        push!(ret.args,quote
            const $(esc(s)) = BinDeps._find_library($dep)
            if isempty($(esc(s)))
                $errorcase
            end
        end)
    end
    if arg1 != nothing && !isa(arg1,Function)
        if !isempty(arg1)
            errrormsg = "The following required libraries were not declared in build.jl:\n"
            for k in (isa(arg1,Vector) ? arg1 : keys(arg1))
                errrormsg *= " - $k\n"
            end
            error(errrormsg)
        end
    end
    ret
end

function build(pkg::String, method; dep::String="", force=false)
    dir = Pkg.dir(pkg)
    file = joinpath(dir,"deps/build.jl")
    context = BinDeps.PackageContext(false,dir,pkg,{})
    m = Module(:__anon__)
    body = Expr(:toplevel,:(ARGS=[$context]),:(include($file)))
    eval(m,body)
    for d in context.deps
        BinDeps.satisfy!(d,[method])
    end
end

# Calculate the SHA-512 hash of a file
using SHA

function sha_check(path, sha)
    calc_sha = sha256(open(path))
    if calc_sha != sha
        error("Checksum mismatch!  Expected:\n$sha\nCalculated:\n$calc_sha\nDelete $path and try again")
    end
end


text:
using Base: LogicalIndex
@testset "event" begin

    for platform in testplatforms

        #Execute Division without problems
        @test_nowarn begin
            m = @agent(
                0,

                [active,id1]::Identity,
                loc::Local,

                EventDivision =
                begin
                    if active == 1
                        active_1 = 2
                        active_2 = 1
                    end
                end
            )

            m = compile(m,platform=platform)
        end

        #Linear growth
        m = @agent(
            0,

            [active,id1]::Identity,
            loc::Local,

            EventDivision =
            begin
                if active == 1
                    active_1 = 2
                    active_2 = 1
                end
            end
        )

        m = compile(m,platform=platform,debug=false)
        # println(m.program)

        @test_nowarn begin

            com = Community(m,N=1)
            n = com.N

            com.id1 .= 4
            com.active .= 1
            com.loc .= 67.

            comt = m.evolve(com,dt=1,tMax=10,nMax=15)

            if comt[end].N != n*comt[end].t+1 error() end
            if !prod(comt[end].loc .== 67.) error() end
            if !prod(comt[end].id1 .== 4) error() end
            if maximum(x->isnan(x) ? -Inf : x, comt.id) != 2*length(comt)-1 error() end

        end

        #Exponential growth
        m = @agent(
            0,

            [active,id1,d]::Identity,
            loc::Local,

            EventDivision =
            begin
                if active == 2
                    id1_1 = id1 + 1
                    id1_2 = id1 + 1
                end
            end
        )

        m = compile(m,platform=platform)
        #println(m.program)

        @test_nowarn begin

            com = Community(m,N=1)
            n = com.N

            com.id1 .= 0
            com.active .= 2
            com.loc .= 67.

            comt = m.evolve(com,dt=1,tMax=6,nMax=100)

            if comt[end].N != n*2^comt[end].t error() end
            if !prod(comt[end].id1 .== comt[end].t) error() end

        end

        #Execute Death without problems
        @test_nowarn begin
            m = @agent(
                0,

                [active,id1]::Identity,
                loc::Local,

                EventDeath = loc < t
            )

            m = compile(m,platform=platform)
        end

        m = @agent(
            1,

            [active,id1]::Identity,
            [loc]::Local,

            UpdateLocal = begin
                id1 = 3
                x = x
            end,

            EventDeath = loc <= t
        )

        m = compile(m,platform=platform)
        #println(m.program)

        @test_nowarn begin

            com = Community(m,N=100)
            n = com.N

            com.loc .= 0:99
            com.x .= 1:100

            comt = m.evolve(com,dt=1,tMax=105,nMax=100)

            if !Bool(prod(comt.N[1:101] .== Array(100:-1:0))) error() end
            if !Bool(prod(comt.N[101:end] .== 0)) error() end
            for i in 1:100
                if !Bool(prod(comt[i].id .== comt[i].x)) error() end
            end

        end

    end
end

ext: jl

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
## ---------------  H Production

abstract type AbstractHcProdSwitches end
abstract type AbstractHcProdSet <: ModelObject end
abstract type AbstractHcProdFct end

Base.show(io :: IO, h :: AbstractHcProdSet) =
    print(io, typeof(h));

## -----------  Access routines

StructLH.describe(h :: AbstractHcProdSwitches) = settings_table(h);

settings_table(h :: AbstractHcProdSet) = settings_table(h.switches);

n_colleges(h :: AbstractHcProdSet) = h.nc;

cal_delta_h!(h :: AbstractHcProdSwitches) = (h.calDeltaH = true);
cal_delta_h(h :: AbstractHcProdSwitches) = h.calDeltaH;
cal_delta_h(h :: AbstractHcProdSet) = cal_delta_h(h.switches);
delta_h(h :: AbstractHcProdSwitches) = h.deltaH;
delta_h(h :: AbstractHcProdSet) = h.deltaH;
delta_h(h :: AbstractHcProdFct) = h.deltaH;

time_exp(h :: AbstractHcProdSet) = h.timeExp;
time_exp(h :: AbstractHcProdFct) = h.timeExp;
h_exp(h :: AbstractHcProdFct) = h.hExp;
set_hexp_lb!(switches :: AbstractHcProdSwitches, hExpLb :: Double) =
    switches.hExpLb = hExpLb;

# Not defined for bounded case
has_tfp(h :: AbstractHcProdFct) = true;
tfp(h :: AbstractHcProdFct) = h.tfp;
set_tfp(h :: AbstractHcProdFct, tfp :: Double) = (h.tfp = tfp);
tfp(h :: AbstractHcProdSet, ic) = values(h.tfpV)[ic];

function h_exp(h :: AbstractHcProdSet)
    if h.switches.sameExponents
        return h.timeExp;
    else
        return h.hExp;
    end
end

function same_exponents!(s :: AbstractHcProdSwitches)
    s.sameExponents = true;
    s.calHExp = false;
end

function separate_exponents!(s :: AbstractHcProdSwitches)
    s.sameExponents = false;
    s.calTimeExp = true;
    s.calHExp = true;
end

same_exponents(h :: AbstractHcProdSwitches) = h.sameExponents;

min_time_per_course(h :: AbstractHcProdFct) = h.minTimePerCourse;

## ----------  Generic methods

# Study time per course, from total study time.
function study_time_per_course(hS :: AbstractHcProdFct, sTimeTotalV, nTriedV)
    return max.(min_time_per_course(hS),  study_time_per_course_no_min(hS, sTimeTotalV, nTriedV))
end

study_time_per_course_no_min(hS :: AbstractHcProdFct, sTimeTotalV, nTriedV) =
    sTimeTotalV ./ nTriedV .- hS.timePerCourse;

# # Study time total, from study time per course.
# function study_time_total(hS :: AbstractHcProdFct, sTimePerCourseV, nTriedV)
#     sTimeV = (max.(min_time_per_course(hS), sTimePerCourseV) .+ hS.timePerCourse) .* nTriedV;
#     return sTimeV
# end

# Study time per course from choice of study time per course (not taking out the fixed time per course).
# function study_time_per_course_from_choice_per_course(hS :: AbstractHcProdFct,
#     sTimeV, nTriedV)

#     sTimePerCourseV = max.(hS.minTimePerCourse, sTimeV .- hS.timePerCourse);
#     return sTimePerCourseV;
# end

"""
	$(SIGNATURES)

h next period. h is restricted to never decline. Study time is total.
"""
hprime(hS :: AbstractHcProdFct, abil, h, h0, sTime, nTried) =
    max.(h, h .* (1.0 .- delta_h(hS)) .+ dh(hS, abil, h, h0, sTime, nTried));

# Path h[t] for given endowments and decisions
function h_path(hS :: AbstractHcProdFct, abil :: Double, h0 :: Double,
    sTimeV, nTriedV)

    T = length(sTimeV);
    hV = zeros(Double, T+1);
    hV[1] = h0;
    for t = 1 : T
        hV[t+1] = hprime(hS, abil, hV[t], h0, sTimeV[t], nTriedV[t]);
    end
    return hV
end

## ----------  Construction

function init_h_depreciation(switches :: AbstractHcProdSwitches)
    deltaH = delta_h(switches);
    pDeltaH = Param(:deltaH, ldescription(:ddh), lsymbol(:ddh),
        deltaH, deltaH, 0.0, 0.1, cal_delta_h(switches));
    return pDeltaH
end

# Fixed time per course that gets subtracted from study time
# Study time per course = study time / nCourses - fixed time per course.
function init_time_per_course()
    # Bounds are arbitrary
    mCourses = data_to_model_courses(1);
    # Time endowment: about 120 hours per week.
    # Normal load: 10 courses per year (5 courses at a time).
    # 1 hour per course => 5 hours per week => 5/120 of time endowment.
    timePerCourseLb = hours_per_week_to_mtu(0.1 / mCourses);
    timePerCourseUb = hours_per_week_to_mtu(2.0 / mCourses);
    timePerCourse = hours_per_week_to_mtu(0.2 / mCourses);

    return Param(:timePerCourse,  ldescription(:sTimeFixed), lsymbol(:sTimeFixed),
        timePerCourse, timePerCourse,
        timePerCourseLb, timePerCourseUb, true)
end

# ----------


text:
include(joinpath(pwd(),"src/models/cartpole.jl"))
include(joinpath(pwd(),"src/constraints/friction.jl"))

optimize = true

# Model
μ0 = 0.1 # coefficient of friction

model_nominal = CartpoleFriction(n, m, d, 1.0, 0.2, 0.5, 9.81, 0.0)
model_friction = CartpoleFriction(n, m, d, 1.0, 0.2, 0.5, 9.81, μ0)

# Horizon
T = 51
tf = 5.0
h = tf / (T-1)

# Bounds
ul_friction = zeros(7)
ul_friction[1] = -10.0
uu_friction = Inf * ones(7)
uu_friction[1] = 10.0

ul, uu = control_bounds(model, T, ul_friction, uu_friction)

# Initial and final states
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [0.0; π; 0.0; 0.0]

xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
    [(t < T ? Diagonal(ones(model_nominal.n))
        : Diagonal(zeros(model_nominal.n))) for t = 1:T],
    [Diagonal([0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model_nominal.m) for t = 1:T])

penalty_obj = PenaltyObjective(1000.0, 7)

multi_obj = MultiObjective([obj, penalty_obj])

# Constraints
n_stage = 5

con_friction = FrictionConstraints(n_stage * (T - 1),
    vcat([(t - 1) * n_stage .+ (3:5) for t = 1:T-1]...),
    n_stage)

# Problem
prob_nominal = trajectory_optimization_problem(model_nominal,
                    multi_obj,
                    T,
                    xl = xl,
                    xu = xu,
                    ul = ul,
                    uu = uu,
                    h = h,
                    con = con_friction)

prob_friction = trajectory_optimization_problem(model_friction,
                    multi_obj,
                    T,
                    xl = xl,
                    xu = xu,
                    ul = ul,
                    uu = uu,
                    h = h,
                    con = con_friction)

# Trajectory initialization
X0 = linear_interp(x1, xT, T) # linear interpolation on state
U0 = [ones(model_nominal.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob_nominal)

# Solve nominal problem
if optimize
    @time Z̄_nominal = solve(prob_nominal, copy(Z0), tol = 1.0e-5, c_tol = 1.0e-5)
    @time Z̄_friction = solve(prob_friction, copy(Z0), tol = 1.0e-5, c_tol = 1.0e-5)

    @save joinpath(@__DIR__, "sol.jld2") Z̄_nominal Z̄_friction
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol.jld2") Z̄_nominal Z̄_friction
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# ---
# title: ICVI Simple Example
# id: icvi_example
# cover: ../assets/logo.png
# date: 2021-12-6
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo is a simple example of how to use a CVI in batch mode.
# ---

# ## Overview

# This demo is a simple example of how to use CVIs incrementally.
# Here, we load a simple dataset and run a basic clustering algorithm to prescribe a set of clusters to the features.
# We will take advantage of the fact that we can compute a criterion value at every step by running the ICVI alongside an online clustering algorithm.
# This simple example demonstrates the usage of a single ICVI, but it may be substituted for any other ICVI in the `ClusterValidityIndices.jl` package.

# ## Online Clustering

# ### Data Setup

# First, we must load all of our dependencies.
# We will load the `ClusterValidityIndices.jl` along with some data utilities and the Julia `Clustering.jl` package to cluster that data.
using ClusterValidityIndices    # CVI/ICVI
using AdaptiveResonance         # DDVFA
using MLDatasets                # Iris dataset
using MLDataUtils               # Shuffling and splitting
using Printf                    # Formatted number printing
using Plots

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# ### ART Online Clustering

# Adaptive Resonance Theory (ART) is a neurocognitive theory that is the basis of a class of online clustering algorithms.
# Because these clustering algorithms run online, we can both cluster and compute a new criterion value at every step.
# For more on these ART algorithms, see [`AdaptiveResonance.jl`](https://github.com/AP6YC/AdaptiveResonance.jl).

## Create a Distributed Dual-Vigilance Fuzzy ART (DDVFA) module with default options
art = DDVFA()
typeof(art)

# Because we are streaming clustering, we must setup the internal data setup of the DDVFA module.
# This is akin to doing some data preprocessing and communicating the dimension of the data, bounds, etc. to the module beforehand.
## Setup the data configuration for the module
data_setup!(art.config, features)
## Verify that the data is setup
art.config.setup

# We can now cluster and get the criterion values online.
# We will do this by creating an ICVI object, setting up containers for the iterations, and then iterating.

## Create an ICVI object
icvi = CH()

## Setup the online/streaming clustering
n_samples = length(labels)          # Number of samples
c_labels = zeros(Int, n_samples)    # Clustering labels
criterion_values = zeros(n_samples) # ICVI outputs

## Iterate over all samples
for ix = 1:n_samples
    ## Extract one sample
    sample = features[:, ix]
    ## Cluster the sample online
    c_labels[ix] = train!(art, sample)
    ## Get the new criterion value (ICVI output)
    criterion_values[ix] = get_icvi!(icvi, sample, c_labels[ix])
end

## See the list of criterion values
criterion_values


text:
################################################################################
#
#       GrpAb/Map.jl : Functions for maps between
#                      finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# 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.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

export haspreimage, has_image, hom, kernel, cokernel, image, is_injective, is_surjective,
       is_bijective

################################################################################
#
#  Check for existence of (pre)image
#
################################################################################

@doc Markdown.doc"""
    haspreimage(M::GrpAbFinGenMap, a::GrpAbFinGenElem) -> Bool, GrpAbFinGenElem

Returns whether $a$ is in the image of $M$. If so, the second return value is
an element $b$ with $M(b) = a$.
"""
function haspreimage(M::GrpAbFinGenMap, a::GrpAbFinGenElem)
  if isdefined(M, :imap)
    return true, preimage(M, a)
  end

  m = vcat(M.map, rels(codomain(M)))
  fl, p = can_solve_with_solution(m, a.coeff, side = :left)

  if fl
    return true, GrpAbFinGenElem(domain(M), view(p, 1:1, 1:ngens(domain(M))))
  else
    return false, id(domain(M))
  end
end

function haspreimage(M::GrpAbFinGenMap, a::Vector{GrpAbFinGenElem})
  if isdefined(M, :imap)
    return true, map(x->preimage(M, x), a)
  end

  m = vcat(M.map, rels(codomain(M)))
  G = domain(M)
  if isdefined(G, :exponent) && fits(Int, G.exponent) && is_prime(G.exponent)
    e = G.exponent
    RR = GF(Int(e))
    fl, p = can_solve_with_solution(map_entries(RR, m), map_entries(RR, vcat([x.coeff for x = a])), side = :left)
    p = map_entries(lift, p)
  else
    fl, p = can_solve_with_solution(m, vcat([x.coeff for x = a]), side = :left)
  end

  if fl
    s = [GrpAbFinGenElem(domain(M), p[i, 1:ngens(domain(M))]) for i=1:length(a)]
    # @assert all(i->M(s[i]) == a[i], 1:length(a))
    return fl, s
  else
    return false, GrpAbFinGenElem[id(domain(M))]
  end
end

# Note that a map can be a partial function. The following function
# checks if an element is in the domain of definition.
#
# Here is an example:
# S, mS = sub(...), so ms: S -> G
# h = pseudo_inv(mS)
# Now h is a partial function on G with domain of definition the image of mS.
# Then has_image(h, x) would check if x is in the image of mS.
function has_image(M::GrpAbFinGenMap, a::GrpAbFinGenElem)
  if isdefined(M, :map)
    return true, image(M, a)
  end

  m = vcat(M.imap, rels(domain(M)))
  fl, p = can_solve_with_solution(m, a.coeff, side = :left)
  if fl
    return true, codomain(M)(view(p, 1:1, 1:ngens(codomain(M))))
  else
    return false, codomain(M)[1]
  end
end

################################################################################
#
#  Homomorphisms specified on arbitrary sets of generators
#
################################################################################

id_hom(G::GrpAbFinGen) = hom(G, G, identity_matrix(FlintZZ, ngens(G)), identity_matrix(FlintZZ, ngens(G)), check = false)

@doc Markdown.doc"""
    hom(A::Vector{GrpAbFinGenElem}, B::Vector{GrpAbFinGenElem}) -> Map

Creates the homomorphism $A[i] \mapsto B[i]$.
"""
function hom(A::Vector{GrpAbFinGenElem}, B::Vector{GrpAbFinGenElem}; check::Bool = true)
  GA = parent(A[1])
  GB = parent(B[1])
  @assert length(B) == length(A)
  @assert length(A) > 0
  #=
  if (check)
    m = vcat(fmpz_mat[x.coeff for x in A])
    m = vcat(m, rels(parent(A[1])))
    i, T = nullspace(transpose(m))
    T = transpose(T)
    T = sub(T, 1:nrows(T), 1:length(A))
    n = vcat(fmpz_mat[x.coeff for x in B])
    n = T*n
    if !can_solve_with_solution(rels(parent(B[1])), n, side = :left)[1]
      error("Data does not define a homomorphism")
    end
  end
  =#
  if ngens(GB) == 0
    return hom(GA, GB, matrix(FlintZZ, ngens(GA), 0, fmpz[]), check = check)
  end

  M = vcat([hcat(A[i].coeff, B[i].coeff) for i = 1:length(A)])
  RA = rels(GA)
  M = vcat(M, hcat(RA, zero_matrix(FlintZZ, nrows(RA), ncols(B[1].coeff))))
  if isdefined(GB, :exponent) && nrows(M) >= ncols(M)
    H = hnf_modular_eldiv(M, exponent(GB))
  else
    H = hnf(M)
  end
  H = sub(H, 1:ngens(GA), ngens(GA)+1:ngens(GA)+ngens(GB))
  h = hom(GA, GB, H, check = check)
  return h
end

@doc Markdown.doc"""
    hom(G::GrpAbFinGen, B::Vector{GrpAbFinGenElem}) -> Map

Creates the homomorphism which maps $G[i]$ to $B[i]$.
"""
function hom(G::GrpAbFinGen, B::Vector{GrpAbFinGenElem}; check::Bool = true)
  GB = parent(B[1])
  @assert length(B) == ngens(G)
  M = vcat([B[i].coeff for i = 1:length(B)])
  h = hom(G, GB, M, check = check)
  return h
end

function hom(G::GrpAbFinGen, H::GrpAbFinGen, B::Vector{GrpAbFinGenElem}; check::Bool = true)
  @assert length(B) == ngens(G)
  @assert all(i -> parent(i) == H, B)
  M = zero_matrix(FlintZZ, ngens(G), ngens(H))
  for i = 1:ngens(G)
    for j = 1:ngens(H)
      M[i, j] = B[i][j]
    end
  end
  h = hom(G, H, M, check = check)
  return h
end

function check_mat(A::GrpAbFinGen, B::GrpAbFinGen, M::fmpz_mat)
  #we have A = X/Y  and B = U/V (generators and relations)
  #        M defines a map (hom) from X -> U
  #        Y -> X is the canonical embedding
  #        V -> U dito
  #                          M
  # need to check if Y -> X --> U lands in V
  # if Y -> X -> B is zero.
  R = rels(A) * M
  return all(x -> iszero(GrpAbFinGenElem(B, R[x, :])), 1:nrows(R))
end

function hom(A::GrpAbFinGen, B::GrpAbFinGen, M::fmpz_mat; check::Bool = true)
  if check
    check_mat(A, B, M) || error("Matrix does not define a morphism of abelian groups")
  end

  return GrpAbFinGenMap(A, B, M)::GrpAbFinGenMap
end

function hom(A::GrpAbFinGen, B::GrpAbFinGen, M::fmpz_mat, Minv; check::Bool = true)
  if check
    check_mat(A, B, M) || error("Matrix does not define a morphism of abelian groups")
    check_mat(B, A, Minv) || error("Matrix does not define a morphism of abelian groups")
    h = GrpAbFinGenMap(A, B, M, Minv)
    ph = pseudo_inv(h)
    all(x -> x == ph(h(x)), gens(A)) || error("Matrix does not define a morphism of abelian groups")
    return h::GrpAbFinGenMap
  end

  return GrpAbFinGenMap(A, B, M, Minv)::GrpAbFinGenMap
end

==(f::GrpAbFinGenMap, g::GrpAbFinGenMap) = domain(f) === domain(g) && codomain(f) === codomain(g) && all(x -> f(x) == g(x), gens(domain(f)))
################################################################################
#
#  Inverse of a map
#
################################################################################

function inv(f::GrpAbFinGenMap)
  if isdefined(f, :imap)
    return hom(codomain(f), domain(f), f.imap, f.map, check = false)
  end
  if !is_injective(f)
    error("The map is not invertible")
  end
  gB = gens(codomain(f))
  imgs = Vector{GrpAbFinGenElem}(undef, length(gB))
  for i = 1:length(imgs)
    fl, el = haspreimage(f, gB[i])
    if !fl
      error("The map is not invertible")
    end
    imgs[i] = el
  end
  return hom(codomain(f),domain(f), imgs, check = false)
end

################################################################################
#
#  Kernel, image, cokernel
#
################################################################################

#TODO: store and reuse on map. Maybe need to change map

@doc Markdown.doc"""
    kernel(h::GrpAbFinGenMap) -> GrpAbFinGen, Map

Let $G$ be the domain of $h$. This function returns an abelian group $A$ and an
injective morphism $f \colon A \to G$, such that the image of $f$ is the kernel
of $h$.
"""
function kernel(h::GrpAbFinGenMap, add_to_lattice::Bool = true)
  G = domain(h)
  H = codomain(h)
  m=zero_matrix(FlintZZ, nrows(h.map)+nrows(rels(H)), ncols(h.map))
  for i=1:nrows(h.map)
    for j=1:ncols(h.map)
      m[i,j]=h.map[i,j]
    end
  end
  if !is_snf(H)
    for i=1:nrels(H)
      for j=1:ngens(H)
        m[nrows(h.map)+i,j]=H.rels[i,j]
      end
    end
  else
    for i=1:length(H.snf)
      m[nrows(h.map)+i,i]=H.snf[i]
    end
  end
  hn, t = hnf_with_transform(m)
  for i = 1:nrows(hn)
    if is_zero_row(hn, i)
      return sub(G, sub(t, i:nrows(t), 1:ngens(G)), add_to_lattice)
    end
  end
  if nrows(hn) == 0
    return sub(G, elem_type(G)[], add_to_lattice)
  end
  error("Something went terribly wrong in kernel computation")
end

@doc Markdown.doc"""
    image(h::GrpAbFinGenMap) -> GrpAbFinGen, Map

Let $G$ be the codomain of $h$. This functions returns an abelian group $A$ and
an injective morphism $f \colon A \to G$, such that the image of $f$ is the
image of $h$.
"""
function image(h::GrpAbFinGenMap, add_to_lattice::Bool = true)
  G = domain(h)
  H = codomain(h)
  hn = hnf(vcat(h.map, rels(H)))
  im = GrpAbFinGenElem[]
  for i = 1:nrows(hn)
    if !is_zero_row(hn, i)
      push!(im, GrpAbFinGenElem(H, sub(hn, i:i, 1:ngens(H))))
    else
      break
    end
  end
  return sub(H, im, add_to_lattice)  # too much, this is sub in hnf....
end

@doc Markdown.doc"""
    cokernel(h::GrpAbFinGenMap) -> GrpAbFinGen, Map

Let $G$ be the codomain of $h$. This function returns an abelian group $A$ and
a morphism $f \colon G \to A$, such that $A$ is the quotient of $G$ with
respect to the image of $h$.
"""
function cokernel(h::GrpAbFinGenMap, add_to_lattice::Bool = true)
  S, mS = image(h, false)
  return quo(codomain(h), GrpAbFinGenElem[mS(g) for g in gens(S)], add_to_lattice)
end

################################################################################
#
#  Surjectivity
#
################################################################################

@doc Markdown.doc"""
    is_surjective(h::GrpAbFinGenMap) -> Bool

Returns whether $h$ is surjective.
"""
function is_surjective(A::GrpAbFinGenMap)
  if isfinite(codomain(A))
    H, mH = image(A)
    return (order(codomain(A)) == order(H))::Bool
  else
    return (order(cokernel(A)[1]) == 1)::Bool
  end
end

################################################################################
#
#  Is zero, Is identity
#
################################################################################

function iszero(h::T) where T <: Map{<:GrpAbFinGen, <:GrpAbFinGen}
  return all(x -> iszero(h(x)), gens(domain(h)))
end

function is_identity(h::T) where T <: Map{<:GrpAbFinGen, <:GrpAbFinGen}
  @assert domain(h) === codomain(h)
  return all(x -> iszero(h(x)-x), gens(domain(h)))
end
################################################################################
#
#  Injectivity
#
################################################################################

#TODO: Improve in the finite case
@doc Markdown.doc"""
    is_injective(h::GrpAbFinGenMap) -> Bool

Returns whether $h$ is injective.
"""
function is_injective(A::GrpAbFinGenMap)
  K = kernel(A, false)[1]
  return isfinite(K) && isone(order(K))
end

################################################################################
#
#  Bijectivity
#
################################################################################

#TODO: Improve in the finite case
@doc Markdown.doc"""
    is_bijective(h::GrpAbFinGenMap) -> Bool

Returns whether $h$ is bijective.
"""
function is_bijective(A::GrpAbFinGenMap)
  return is_injective(A) && is_surjective(A)
end

###############################################################################
#
#  Compose and Squash for abelian group maps
#
###############################################################################

function compose(f::GrpAbFinGenMap, g::GrpAbFinGenMap)
  @assert domain(g) == codomain(f)
  C = codomain(g)
  if isdefined(C, :exponent)
    if fits(Int, C.exponent)
      RR = ResidueRing(FlintZZ, Int(C.exponent), cached = false)
      fRR = map_entries(RR, f.map)
      gRR = map_entries(RR, g.map)
      MRR = fRR*gRR
      M = lift(MRR)
    else
      R = ResidueRing(FlintZZ, C.exponent, cached = false)
      fR = map_entries(R, f.map)
      gR = map_entries(R, g.map)
      MR = fR*gR
      M = map_entries(lift, MR)
    end
  else
    M = f.map*g.map
  end
  if is_snf(C)
    reduce_mod_snf!(M, C.snf)
  else
    assure_has_hnf(C)
    reduce_mod_hnf_ur!(M, C.hnf)
  end
  return hom(domain(f), codomain(g), M, check = false)

end

###############################################################################
struct MapParent
  dom
  codom
  typ::String
end

elem_type(::Type{MapParent}) = Map

function show(io::IO, MP::MapParent)
  print(io, "Set of all $(MP.typ) from $(MP.dom) to $(MP.codom)")
end

parent(f::GrpAbFinGenMap) = MapParent(domain(f), codomain(f), "homomorphisms")

function cyclic_hom(a::fmpz, b::fmpz)
  #hom from Z/a -> Z/b
  if iszero(a)
    return (b, fmpz(1))
  end
  if iszero(b)
    return (fmpz(1), fmpz(1))
  end
  g = gcd(a, b)
  return (g, divexact(b, g))
end

@doc Markdown.doc"""
    hom(G::GrpAbFinGen, H::GrpAbFinGen; task::Symbol = :map) -> GrpAbFinGen, Map

Computes the group of all homomorpisms from $G$ to $H$ as an abstract group.
If `task` is ":map", then a map $\phi$ is computed that can be used
to obtain actual homomorphisms. This map also allows preimages.
Set `task` to ":none" to not compute the map.
"""
function hom(G::GrpAbFinGen, H::GrpAbFinGen; task::Symbol = :map)
  @assert task in [:map, :none]
  sG, mG = snf(G)  # mG : sG -> G
  sH, mH = snf(H)  # mH : sH -> G
  n = ngens(sG)
  m = ngens(sH)
  r = [cyclic_hom(x, y) for x = sG.snf for y = sH.snf]
  R = GrpAbFinGen([x[1] for x = r])
  set_attribute!(R, :hom=>(G, H), :show => show_hom)
  if task == :none
    return R
  end

  c = [x[2] for x = r]

  function phi(r::GrpAbFinGenElem)
    return GrpAbFinGenMap(inv(mG) * hom(sG, sH, matrix(FlintZZ, n, m, [r[i]*c[i] for i=1:ngens(R)]), check = true) * mH)
  end

  function ihp(r::GrpAbFinGenMap)
    #transpose is due to linear indexing being wrong order
    if ngens(sG) == 0
      return R[0]
    end
    local m = transpose(vcat([preimage(mH, r(mG(sG[i]))).coeff for i = 1:ngens(sG)]))
    return R([divexact(m[i], c[i]) for i = 1:ngens(R)])
  end
  return R, MapFromFunc(phi, ihp, R, MapParent(G, H, "homomorphisms"))
end


text:
#!/usr/bin/env julia

using SpectralElements, OrdinaryDiffEq, Flux
using LinearAlgebra,Plots
#----------------------------------#
Ex = 2; nr1 = 8;
Ey = 2; ns1 = 8;

m1 = Mesh(nr1,ns1,Ex,Ey,[true,true])
bc = ['D','D','D','D']
M  = generateMask(bc,m1)
opM(u) = m1.Bi .* u
opB(u) = gatherScatter(mass(u,m1),m1)
#----------------------------------#
function setIC!(u,x,y,t)
    u = 0.0
    return
end

function setBC!(ub,x,y,t)
    ub = 0+0*x
    return
end

function setForcing!(f,x,y,t)
    f = 1+0*x
    return
end

setVisc!(ν,x,y,t) = 1+0*x

#----------------------------------#
it = 0
function cond(u,t,integrator)
    global it
    cond = (it % 100) == 0
    it += 1
    return cond
end

function affect!(integrator)
    global it
    u = integrator.u
    p = integrator.p
    t = integrator.t
    dt = integrator.dt

    ut = utrue(m1.x,m1.y,t)
    er = norm(ut-u,Inf)
    println("Step=$it, Time=$t, dt=$dt, er=$er")
    plt = meshplt(u,m1)
    plt = plot!(zlims=(-1,1))
    display(plt)
    return
end

cb = DiscreteCallback(cond,affect!,save_positions=(false,false))

#----------------------------------#
function dudt(u,p,t)
    f = setForcing.(m1.x,m1.y,t)

    rhs = -lapl(u,m1) + mass(f,m1)
    rhs = gatherScatter(rhs,m1)

    dudt = pcg(rhs,opB;opM=opM,mult=m1.mult,ifv=false)

    return dudt
end

u0   = setIC(m1.x,m1.y)
dt   = 0.01
tspn = (0.0,1.0)

prob = ODEProblem(dudt,u0,tspn)
sol  = solve(prob,RK4();saveat=0.1,callback=cb)

err = sol.u - [utrue.(m1.x,m1.y,sol.t[i]) for i=axes(sol.t,1)]
ee = maximum.(abs.(err[i]) for i=axes(sol.t,1))
#----------------------------------#
return ee


text:
mutable struct TServerBase
    srvr_t::TServerTransport
    processor::TProcessor
    in_t::Function
    in_p::Function
    out_t::Function
    out_p::Function
end

function serve_accepted(client::TTransport, s::TServerBase)
    itrans = s.in_t(client)
    otrans = s.out_t(client)
    iprot = s.in_p(itrans)
    oprot = s.out_p(otrans)

    try
        while true
            process(s.processor, iprot, oprot)
        end
    catch ex
        if !isa(ex, EOFError)
            println(ex)
            Base.show_backtrace(stderr, catch_backtrace())
        end
    end
    close(itrans)
    close(otrans)
end

close(srvr::TServer) = close(srvr.base.srvr_t)

##
# Blocking server. Requests are processed in the main task.
mutable struct TSimpleServer <: TServer
    base::TServerBase
    TSimpleServer(srvr_t::TServerTransport, processor::TProcessor, in_t::Function, in_p::Function, out_t::Function, out_p::Function) = new(TServerBase(srvr_t, processor, in_t, in_p, out_t, out_p))
end

function serve(ss::TSimpleServer)
    s = ss.base
    listen(s.srvr_t)

    while true
        client = accept(s.srvr_t)
        serve_accepted(client, s)
    end
end

##
# Task server. Tasks are spawned for each connection.
mutable struct TTaskServer <: TServer
    base::TServerBase
    TTaskServer(srvr_t::TServerTransport, processor::TProcessor, in_t::Function, in_p::Function, out_t::Function, out_p::Function) = new(TServerBase(srvr_t, processor, in_t, in_p, out_t, out_p))
end

##
# Process Pool Server
mutable struct TProcessPoolServer <: TServer
    base::TServerBase
    function TProcessPoolServer(srvr_t::TServerTransport, processor::TProcessor, in_t::Function, in_p::Function, out_t::Function, out_p::Function)
        distribute(processor)
        new(TServerBase(srvr_t, processor, in_t, in_p, out_t, out_p))
    end
end

const TAsyncServer = Union{TTaskServer, TProcessPoolServer}

function serve(ss::TAsyncServer)
    s = ss.base
    listen(s.srvr_t)

    while true
        client = accept(s.srvr_t)
        @async serve_accepted(client, s)
    end
end


text:
module Invest

include("option.jl")

export AmericanCall, AmericanPut
export EuropeanCall, EuropeanPut
export Δ
export Γ
export Θ
export isat
export price
export dprice

include("interest.jl")

export aer
export cci
export find_effective_rate
export irr
export npv
export pv
export present_value
export yearly_compounded_interest

end


text:
const alias_dict = Dict(
    Symbol(alias["name"]) => Symbol(alias["alias"]) for alias ∈ findall("//*[@alias and not(@feature)]", xroot)
)

"""
Whether this type is an alias for another name.
"""
isalias(name) = name ∈ keys(alias_dict)

"""
Whether an alias was built from this name.
"""
hasalias(name) = name ∈ values(alias_dict)

alias_verts = unique(vcat(collect(keys(alias_dict)), collect(values(alias_dict))))

function build_alias_graph(alias_verts, alias_dict)
    g = SimpleDiGraph(length(alias_verts))

    for (j, (src, dst)) ∈ enumerate(alias_dict)
        i = findfirst(==(dst), alias_verts)
        add_edge!(g, i, j)
    end

    g
end

const alias_graph = build_alias_graph(alias_verts, alias_dict)

aliases(alias_graph::SimpleDiGraph, index) = getindex.(Ref(alias_verts), outneighbors(alias_graph, index))

function aliases(name)
    index = findfirst(==(name), alias_verts)
    isnothing(index) ? Symbol[] : aliases(alias_graph, index)
end

function follow_alias(alias_graph::SimpleDiGraph, index)
    indices = inneighbors(alias_graph, index)
    if isempty(indices)
        alias_verts[index]
    elseif length(indices) > 1
        error(
            "More than one indices returned for $(alias_verts[index]) when following alias $(getindex.(Ref(alias_verts), indices))",
        )
    else
        i = first(indices)
        if i == index
            alias_verts[i]
        else
            follow_alias(alias_graph, i)
        end
    end
end

function follow_alias(name)
    index = findfirst(==(name), alias_verts)
    isnothing(index) ? name : follow_alias(alias_graph, index)
end


text:
@recipe function f(prefit::PreFit, xs; labels=("pedestal", "spe"))
    xguide := "charge [A.U.]"
    yguide := "counts"
    yscale := :log10
    ylims := (1, 1e5)
    @series begin
        linewidth := 2
        label := labels[1]
        xs, gauss(xs, [prefit.μₚ, prefit.σₚ, prefit.Aₚ])
    end
    @series begin
        linewidth := 2
        label := labels[2]
        xs, gauss(xs, [prefit.μₛ, prefit.σₛ, prefit.Aₛ])
    end
end

@recipe function f(pmtrespfit::PMTRespFit, xs; show_single_gaussians=false, n_gaussians=5)
    xguide := "charge [A.U.]"
    yguide := "counts"
    yscale := :log10
    ylims := (1, 1e5)
    @series begin
        linewidth := 2
        xs, pmtresp(xs, [pmtrespfit.μₚ, pmtrespfit.σₚ,
                         pmtrespfit.μₛ, pmtrespfit.σₛ, pmtrespfit.nₚₑ,
                         pmtrespfit.A])
    end
    if show_single_gaussians
        poisson = Poisson(abs(pmtrespfit.nₚₑ))
        for i in 1:n_gaussians
            @series begin
                linewidth := 2
                linestyle := :dash
                seriescolor := :black
                xs, pmtrespfit.A .* gauss(xs,
                                        [i * pmtrespfit.μₛ,
                                         sqrt(i) * pmtrespfit.σₛ,
                                         pdf(poisson, i)])
            end
        end
    end
end

@recipe function f(pmtrespfit::PMTRespUapFit, xs; show_single_gaussians=false)
    xguide := "charge [A.U.]"
    yguide := "counts"
    yscale := :log10
    ylims := (1, 1e5)
    @series begin
        linewidth := 2
        xs, pmtresp_uap(xs, [pmtrespfit.μₚ, pmtrespfit.σₚ,
                         pmtrespfit.μₛ, pmtrespfit.σₛ, pmtrespfit.nₚₑ,
                         pmtrespfit.A, pmtrespfit.μᵤₐₚ, pmtrespfit.σᵤₐₚ,
                         pmtrespfit.Aᵤₐₚ])
    end
    if show_single_gaussians
        println("fdsagf")
        poisson = Poisson(abs(pmtrespfit.nₚₑ))
        for i in 1:n_gaussians
            @series begin
                linewidth := 2
                linestyle := :dash
                seriescolor := :black
                xs, pmtrespfit.A .* gauss(xs,
                                        [i * pmtrespfit.μₛ,
                                         sqrt(i) * pmtrespfit.σₛ,
                                         pdf(poisson, i)])
            end
        end
    end
end


text:
module Align

using Rematch, Dates
using LinearAlgebra
using ProgressMeter

using Base.Threads: @spawn, @threads

using ..Utility: read_paf, enforce_cutoff!
using ..Blocks
using ..Paths: replace!
using ..Nodes
using ..Graphs
using ..Mash

import ..Minimap

export align

# ------------------------------------------------------------------------
# helper functions

function log(msg...)
    println(stderr, now(), " ", join(msg, " ")...)
    flush(stderr)
end

# ------------------------------------------------------------------------
# guide tree for order of pairwise comparison for multiple alignments

# TODO: distance?
"""
	mutable struct Clade
		name   :: String
		parent :: Union{Clade,Nothing}
		left   :: Union{Clade,Nothing}
		right  :: Union{Clade,Nothing}
		graph  :: Channel{Graph}
	end

Clade is a node (internal or leaf) of a binary guide tree used to order pairwise alignments
associated to a multiple genome alignment in progress.
`name` is only non-empty for leaf nodes.
`parent` is `nothing` for the root node.
`graph` is a 0-sized channel that is used as a message passing primitive in alignment.
"""
mutable struct Clade
    name   :: String
    parent :: Union{Clade,Nothing}
    left   :: Union{Clade,Nothing}
    right  :: Union{Clade,Nothing}
    graph  :: Channel{Graph}
end

# ---------------------------
# constructors

"""
	Clade()

Generate an empty, disconnected clade.
"""
Clade()     = Clade("",nothing,nothing,nothing,Channel{Graph}(2))
"""
	Clade(name)

Generate an empty, disconnected clade with name `name`.
"""
Clade(name) = Clade(name,nothing,nothing,nothing,Channel{Graph}(2))
"""
	Clade(left::Clade, right::Clade)

Generate an nameless clade with `left` and `right` children.
"""
function Clade(left::Clade, right::Clade)
    parent = Clade("",nothing,left,right,Channel{Graph}(2))

    left.parent = parent
    right.parent = parent

    parent
end

"""
	Clade(distance, names; algo=:nj)

Generate a tree from a matrix of pairwise distances `distance`.
The names of leafs are given by an array of strings `names`.
`algo` dictates the algorithm used to transform the distance matrix into a tree.
Currently on neighbor joining (:nj) is supported.
"""
function Clade(distance, names; algo=:nj)
    @match algo begin
        :nj => return nj(distance, names)
        _ => error("unrecognized algorithm $(algo)")
    end
end

"""
	nj(distance, names)

Lower-level function.
Generate a tree from a matrix of pairwise distances `distance`.
The names of leafs are given by an array of strings `names`.
Uses neighbor joining.
"""
function nj(distance, names)
    nodes = [Clade(names[i]) for i in 1:size(distance,1)]

    # internally scoped functions
    Q = function(D)
        n = size(D,1)
        q = (n-2)*D .- sum(D,dims=1) .- sum(D,dims=2)
        q[diagind(q)] .= Inf
        return q
    end

    pair = function(Q)
        ι  = argmin(Q)
        q₀ = Q[ι]
        if ι.I[1] > ι.I[2]
            return ι.I[2], ι.I[1], q₀
        end
        return ι.I[1], ι.I[2], q₀
    end

    dist = function(D, i, j)
        n  = size(D,1)
        d₁ = .5*D[i,j] + 1/(2*(n-2))*(sum(D[i,:]) - sum(D[j,:]))
        d₂ = D[i,j] - d₁

        if d₁ < 0
            d₂ -= d₁
            d₁ = 0
        end

        if d₂ < 0
            d₁ -= d₂
            d₂ = 0
        end

        dₙ = .5*(D[i,:] .+ D[j,:] .+ D[i,j])

        return d₁, d₂, dₙ
    end

    join! = function(D)
        q = Q(D)
        i, j, q₀ = pair(q)

        node = Clade(nodes[i], nodes[j])
        nodes[i].parent = node
        nodes[j].parent = node

        d₁, d₂, dₙ = dist(D, i, j)
        D[i,:] .= dₙ
        D[:,i] .= dₙ
        D[i,i]  = 0

        nodes[i] = node

        D = D[1:end .!= j, 1:end .!= j]
        deleteat!(nodes, j)

        return D
    end

    # body of function
    @assert length(names) == length(Set(names))

    while length(nodes) > 2
        distance = join!(distance)
    end

    return Clade(nodes[1], nodes[2])
end

# ---------------------------
# operators

"""
	isleaf(c::Clade)

Return if Clade `c` is a terminal node, i.e. a leaf.
"""
isleaf(c::Clade) = isnothing(c.left) && isnothing(c.right)

# serialization to newick format
function Base.show(io::IO, c::Clade)
    if isleaf(c)
        print(io, c.name)
    else
        print(io, "(")
        show(io, c.left)
        print(io, ",")
        show(io, c.right)
        print(io, ")")
    end
end

"""
	leaves(root::Clade)

Return all terminal nodes that have `root` as their common ancestor.
"""
function leaves(root::Clade)
    itr = Channel{Clade}(0)
    function traverse(node::Clade)
        if isleaf(node)
            put!(itr, node)
        else
            traverse(node.left)
            traverse(node.right)
        end
    end

    @async begin
        traverse(root)
        close(itr)
    end

    return itr
end

"""
	postorder(root::Clade)

Return an postorder iterator over descendents of `root`.
"""
function postorder(root::Clade)
    itr = Channel{Clade}(0)
    function traverse(node::Clade)
        if isleaf(node)
            put!(itr, node)
        else
            traverse(node.left)
            traverse(node.right)
            put!(itr,node)
        end
    end

    @async begin
        traverse(root)
        close(itr)
    end

    return itr
end

"""
	preorder(root::Clade)

Return an preorder iterator over descendents of `root`.
"""
function preorder(root::Clade)
    itr = Channel{Clade}(0)
    function traverse(node::Clade)
        if isleaf(node)
            put!(itr, node)
        else
            put!(itr,node)
            traverse(node.left)
            traverse(node.right)
        end
    end

    @async begin
        traverse(root)
        close(itr)
    end

    return itr
end

function bisect(tips)
    length(tips) > 1 || return tips[1]

    m = length(tips) ÷ 2
    l = bisect(tips[1:m])
    r = bisect(tips[(m+1):end])

    return Clade(l, r)
end

"""
	balance(root::Clade)

Rotate a binary tree to a balanced configuration.
Preserves the topological ordering of the leaves.
"""
function balance(root::Clade)
    tips = collect(leaves(root))
    return bisect(tips)
end

function Base.length(node::Clade)
    isleaf(node) && return 1
    return 1 + Base.length(node.left) + Base.length(node.right)
end

# ------------------------------------------------------------------------
# ordering functions

# TODO: assumes the input graphs are singletons! generalize
"""
	ordering(compare, Gs...)

Return a guide tree based upon distances computed from a collection of graphs `Gs`, using method `compare`.
The signature of `compare` is expected to be `compare(G::Graphs....) -> distance, names`.
Graphs `Gs...` are expected to be singleton graphs.
"""
function ordering(compare, Gs...)
    distance, names = compare(Gs...)
    return Clade(distance, names; algo=:nj)
end

# ------------------------------------------------------------------------
# align functions

function preprocess(hits, skip, energy, blocks!)
    hits = [
        let
            qry₀, ref₀ = blocks!(hit)

            hit.qry.seq = qry₀.sequence
            hit.ref.seq = ref₀.sequence

            qry, strand =
            if !hit.orientation
                reverse_complement!(hit)
                reverse_complement(qry₀; keepid=true), false
            else
                qry₀, true
            end
            ref = ref₀

            (
                hit=hit,
                qry=qry,
                ref=ref,
                qry₀=qry₀,
                ref₀=ref₀,
                strand=strand,
           )
        end for hit in hits if (energy(hit) < 0 && !skip(hit))
    ]

    return hits
end

function do_align(G₁::Graph, G₂::Graph, energy::Function, minblock::Int, sensitivity::String)
    hits = Minimap.align(pancontigs(G₁), pancontigs(G₂), minblock, sensitivity)
    sort!(hits; by=energy)

    return hits
end

# TODO: make thread safe
function align_kernel(matches, minblock, replace!, verbose)
    blocks   = Array{Array{Block,1},1}(undef, length(matches))
    for (i, match) in enumerate(matches)
        # destructure precomputed hit
        hit,qry,ref,qry₀,ref₀,strand = match

        verbose && log(hit)

        enforce_cutoff!(hit, minblock)
        blks = combine(qry, ref, hit; minblock=minblock)

        qrys = map(b -> b.block, filter(b -> b.kind != :ref, blks))
        refs = map(b -> b.block, filter(b -> b.kind != :qry, blks))

        replace!((qry=qry₀, ref=ref₀), (qry=qrys, ref=refs), strand)

        blocks[i] = map(b->b.block, blks)
    end

    return Dict(blk.uuid => blk for block in blocks for blk in block)
end

"""
	align_self(G₁::Graph, energy::Function, minblock::Int, verify::Function, verbose::Bool; maxiter=100)

Align graph `G₁` to itself by looking for homology between blocks.
This is a low-level function.

`energy` is to be a function that takes an alignment between two blocks and produces a score.
The _lower_ the score, the _better_ the alignment. Only negative energies are considered.

`minblock` is the minimum size block that will be produced from the algorithm.
`maxiter` is maximum number of duplications that will be considered during this alignment.
"""
function align_self(G₁::Graph, energy::Function, minblock::Int, verify::Function, verbose::Bool; maxiter=100, sensitivity="asm10")
    G₀ = G₁

    for niter in 1:maxiter
        # calculate pairwise hits
        hits = do_align(G₀, G₀, energy, minblock, sensitivity)

        skip = (hit) -> (
               (hit.qry.name == hit.ref.name)
            || (hit.length < minblock)
            || (!(hit.qry.name in keys(G₀.block)) || !(hit.ref.name in keys(G₀.block)))
        )
        block = (hit) -> (
            qry = pop!(G₀.block, hit.qry.name),
            ref = pop!(G₀.block, hit.ref.name),
        )
        replace = (old, new, orientation) -> let
            for path in values(G₀.sequence)
                replace!(path, old.qry, new.qry, orientation)
            end

            for path in values(G₀.sequence)
                replace!(path, old.ref, new.ref, true)
            end
        end

        # filter out according to energy & skip conditions
        # only continue if at least one hit survived
        merges = preprocess(hits, skip, energy, block)
        length(merges) > 0 || break

        blocks = align_kernel(merges, minblock, replace, verbose)
        merge!(blocks, G₀.block)

        G₀ = Graph(
            blocks,
            G₀.sequence,
        )
        detransitive!(G₀)
        purge!(G₀)
        prune!(G₀)
    end

    return G₀
end

function compare(old, new, strand)
    for node in keys(old)
        oldseq = sequence(old, node)

        δ         = 1
        newblkseq = UInt8[]
        oldblkseq = UInt8[]
        oldstrand = node.strand

        if oldstrand && strand
            for blk in new
                newblkseq = sequence(blk, node)
                oldblkseq = oldseq[δ:δ+length(newblkseq)-1]
                if !all(oldblkseq .== newblkseq)
                    @goto bad
                end

                δ += length(newblkseq)
            end
        elseif oldstrand && !strand
            for blk in reverse(new)
                newblkseq = reverse_complement(sequence(blk, node))
                oldblkseq = oldseq[δ:δ+length(newblkseq)-1]
                if !all(oldblkseq .== newblkseq)
                    @goto bad
                end

                δ += length(newblkseq)
            end
        elseif !oldstrand && strand
            for blk in reverse(new)
                newblkseq = sequence(blk, node)
                oldblkseq = oldseq[δ:δ+length(newblkseq)-1]
                if !all(oldblkseq .== newblkseq)
                    @goto bad
                end

                δ += length(newblkseq)
            end
        else
            for blk in new
                newblkseq = reverse_complement(sequence(blk, node))
                oldblkseq = oldseq[δ:δ+length(newblkseq)-1]
                if !all(oldblkseq .== newblkseq)
                    @goto bad
                end

                δ += length(newblkseq)
            end
        end

        continue

    @label bad
        badblkloci = Int[]
        for i ∈ 1:min(length(newblkseq),length(oldblkseq))
            if newblkseq[i] != oldblkseq[i]
                push!(badblkloci, i)
            end
        end
        left, right = max(badblkloci[1]-10, 1), min(badblkloci[1]+10, length(newblkseq))

        println("--> length:           ref($(length(oldblkseq))) <=> seq($(length(newblkseq)))")
        println("--> number bad:       $(length(badblkloci))")
        println("--> window:           $(left):$(badblkloci[1]):$(right)")
        println("--> old:              $(String(oldblkseq[left:right]))")
        println("--> new:              $(String(newblkseq[left:right]))")

        @show strand, oldstrand

        error("bad sequence reconstruction")
    end
end

"""
	align_pair(G₁::Graph, G₂::Graph, energy::Function, minblock::Int, verify::Function, verbose::Bool; maxiter=100)

Align graph `G₁` to graph `G₂` by looking for homology between consensus sequences of blocks.
This is a low-level function.

`energy` is to be a function that takes an alignment between two blocks and produces a score.
The _lower_ the score, the _better_ the alignment. Only negative energies are considered.

`minblock` is the minimum size block that will be produced from the algorithm.
`maxiter` is maximum number of duplications that will be considered during this alignment.
"""
function align_pair(G₁::Graph, G₂::Graph, energy::Function, minblock::Int, verify::Function, verbose::Bool; sensitivity="asm10")
    hits = do_align(G₁, G₂, energy, minblock, sensitivity)

    skip = (hit) -> (
            !(hit.qry.name in keys(G₁.block))
         || !(hit.ref.name in keys(G₂.block))
         || (hit.length < minblock)
    )
    block = (hit) -> (
        qry = pop!(G₁.block, hit.qry.name),
        ref = pop!(G₂.block, hit.ref.name),
    )

    replace = (old, new, orientation) -> let
        for path in values(G₁.sequence)
            replace!(path, old.qry, new.qry, orientation)
        end

        for path in values(G₂.sequence)
            replace!(path, old.ref, new.ref, true)
        end
    end

    merges = preprocess(hits, skip, energy, block)

    blocks   = align_kernel(merges, minblock, replace, verbose)
    sequence = merge(G₁.sequence, G₂.sequence)

    # XXX: worry about uuid collision?
    merge!(blocks, G₁.block)
    merge!(blocks, G₂.block)

    G = Graph(
        blocks,
        sequence,
    )

    detransitive!(G)
    purge!(G)
    prune!(G)

    return G
end

# TODO: the associative array is a bit hacky...
#       can we push it directly into the channel?
"""
	align(G::Graph...; compare=Mash.distance, energy=(hit)->(-Inf), minblock=100)

Align graph `G₁` to itself by looking for homology between blocks.
Multithreaded by default.
This is usually the function you want.

`energy` is to be a function that takes an alignment between two blocks and produces a score.
The _lower_ the score, the _better_ the alignment. Only negative energies are considered.

`minblock` is the minimum size block that will be produced from the algorithm.
`maxiter` is maximum number of duplications that will be considered during this alignment.

`compare` is the function to be used to generate pairwise distances that generate the internal guide tree.
"""
function align(Gs::Graph...; compare=Mash.distance, energy=(hit)->(-Inf), minblock=100, reference=nothing, sensitivity="asm10", maxiter=100)
    function verify(graph, msg="")
        if reference !== nothing
            for (name,path) ∈ graph.sequence
                seq = sequence(path)
                ref = reference[name]
                if seq != ref
                    badloci = Int[]
                    for i ∈ 1:min(length(seq),length(ref))
                        if seq[i] != ref[i]
                            push!(badloci, i)
                        end
                    end
                    left, right = max(badloci[1]-10, 1), min(badloci[1]+10, length(seq))
                    cumulative_lengths = cumsum([length(n.block, n) for n in path.node])

                    println("--> length:           ref($(length(ref))) <=> seq($(length(seq)))")
                    println("--> number of nodes:  $(length(path.node))")
                    println("--> path offset:      $(path.offset)")
                    println("--> |badloci|:        $(length(badloci))")
                    println("--> window:           $(left):$(badloci[1]):$(right)")
                    println("--> ref:              $(ref[left:right])")
                    println("--> seq:              $(seq[left:right])")

                    i = 1
                    while i < length(path.node)
                        if cumulative_lengths[i] > badloci[1]
                            break
                        else
                            i += 1
                        end
                    end

                    println("--> i(bad):           $(i)")
                    println("--> mutate:           $(path.node[i].block.mutate[path.node[i]])")
                    println("--> insert:           $(path.node[i].block.insert[path.node[i]])")
                    println("--> delete:           $(path.node[i].block.delete[path.node[i]])")

                    error("--> isolate '$name' incorrectly reconstructed")
                end
            end
        end

        graph
    end

    log("--> ordering")
    tree = ordering(compare, Gs...) |> balance
    log("--> tree: ", tree)

    meter = Progress(length(tree); desc="alignment progress", output=stderr)
    tips  = Dict{String,Graph}(collect(keys(G.sequence))[1] => G for G in Gs)

    log("--> aligning pairs")

    events = Channel{Bool}(10);

    @spawn let
        while isopen(events)
            take!(events)
            next!(meter)
        end
    end

    error_channel = Channel(1)

    G = nothing
    for clade ∈ postorder(tree)
        @spawn try
            if isleaf(clade)
                close(clade.graph)
                put!(clade.parent.graph, tips[clade.name])
                put!(events, true)
            else
                Gₗ = take!(clade.graph)
                Gᵣ = take!(clade.graph)
                close(clade.graph)

                G₀ = align_pair(Gₗ, Gᵣ, energy, minblock, verify, false; sensitivity=sensitivity)
                G₀ = align_self(G₀, energy, minblock, verify, false; sensitivity=sensitivity, maxiter=maxiter)

                put!(events, true)
                if clade.parent !== nothing
                    put!(clade.parent.graph, G₀)
                else
                    G = G₀
                    close(error_channel)
                end
            end
        catch err
            put!(error_channel, [err, catch_backtrace()])
            close(error_channel)
        end
    end

    for (err, backtrace) in error_channel
        @error "In-thread error during graph building:" exception=(err, backtrace)
        error("graph construction failed, see above for stacktrace")
    end

    close(events)

    return G
end

# ------------------------------------------------------------------------
# testing

function test()
    log("> testing mash...")
    distance, names = mash("data/generated/assemblies/isolates.fna.gz")
    tree = Clade(distance, names; algo=:nj)
    log("done!")

end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
@testset "StaticDimer vs BBDimer" begin

heading1("TEST: StaticDimer, BBDimer with and without Preconditioning")

heading2("Muller potential")

V = MullerPotential()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)

dimer = StaticDimer(a_trans=0.002, a_rot=0.002, len=1e-3, verbose=verbose)
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

bbdimer = BBDimer(a0_trans=0.002, a0_rot=0.002, verbose=verbose)
xbb, vbb, bblog = run!(bbdimer, E, dE, x0, v0)
@test res_trans(bblog)[end] <= dimer.tol_trans
@test res_rot(bblog)[end] <= dimer.tol_rot
@test vecnorm(xbb - x, Inf) < 1e-4

heading2("Muller potential and preconditioning")

dimer = StaticDimer( a_trans=0.25, a_rot=0.5, len=1e-3, verbose=verbose,
      precon=eye(2), precon_rot=true, precon_prep! = (P,x) -> hessprecond(V, x) )
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

dimer = BBDimer( a0_trans=0.25, a0_rot=0.5, len=1e-3, verbose=verbose,
      precon=eye(2), precon_rot=true, precon_prep! = (P,x) -> hessprecond(V, x) )
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

heading2("Standard double-well")

V = DoubleWell()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)

dimer = StaticDimer(a_trans=0.66, a_rot=0.4, len=1e-3, verbose=verbose)
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

bbdimer = BBDimer(a0_trans=0.5, a0_rot=0.5, verbose=verbose)
xbb, vbb, bblog = run!(bbdimer, E, dE, x0, v0)
@test res_trans(bblog)[end] <= dimer.tol_trans
@test res_rot(bblog)[end] <= dimer.tol_rot
@test vecnorm(xbb - x, Inf) < 1e-4

heading2("Ill-conditioned double-well")

V = DoubleWell(diagm([1.0, 10.0]))
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)

dimer = StaticDimer(a_trans=0.1, a_rot=0.1, len=1e-3, verbose=verbose)
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

bbdimer = BBDimer(a0_trans=0.1, a0_rot=0.1, verbose=verbose)
xbb, vbb, bblog = run!(bbdimer, E, dE, x0, v0)
@test res_trans(bblog)[end] <= dimer.tol_trans
@test res_rot(bblog)[end] <= dimer.tol_rot
@test vecnorm(xbb - x, Inf) < 1e-4

heading2("Ill-conditioned double-well, but now with preconditioner")

dimer = StaticDimer(a_trans=0.66, a_rot=0.4, len=1e-3, verbose=verbose,
      precon=V.A, precon_rot=true)
x, v, log = run!(dimer, E, dE, x0, v0)
@test res_trans(log)[end] <= dimer.tol_trans
@test res_rot(log)[end] <= dimer.tol_rot

bbdimer = BBDimer(a0_trans=0.66, a0_rot=0.4, precon=V.A,
      precon_rot=true, verbose=verbose)
xbb, vbb, bblog = run!(bbdimer, E, dE, x0, v0)
@test res_trans(bblog)[end] <= dimer.tol_trans
@test res_rot(bblog)[end] <= dimer.tol_rot
@test vecnorm(xbb - x, Inf) < 1e-4

end


text:
module Bisection

using LinearAlgebra
import Base.@kwdef

@kwdef struct Results{T<:Real}
	best::T				# The best x that solves the optimization problem
	val::T				# The best value returned (e.g. val = f(best))
	iter::Int			# Number of iterations used in the method
	f_calls::Int
	info::String		# Converged?
end

@kwdef mutable struct Options
	iter::Int
	f_tol::Float64
	tol::Float64
end

function Options(;iter::Int=1000,f_tol::Float64 = 1e-3,tol::Float64 = 1e-3)

	Options(iter,f_tol,tol)

end

function bisectionSearch(a::Real,b::Real,f::Function,options::Options=Options())

	# Promote the endpoints to the same type
	a,b = promote(a,b)

	# Default (not converged)
	converged = false

	# get initial values for high and low
	fa,fb,f_calls = f(a),f(b),2

	# If they have the same sign then this method will not work
	if sign(fa) == sign(fb)
		@warn("sgn(f(a)) = sgn(f(b)) is not allowed for bisection.")
		return a,fa,b,fb
	end

	# Introduce midpoint and f(midpoint) to be returned to the results struct
	local iterConverged,c, fc

	# Begin bisection
	for iter = 1 : options.iter

		# Compute midpoint
		c = (a + b)/2.0

		# Evaluate f
		fc, f_calls = f(c), f_calls+1

		# Check which sign f(c) has and replace the endpoint with the same sign
		if sign(fa) == sign(fc)
			fa = fc
			a = c
		else
			fb = fc
			b = c
		end

		# Converged
		if abs(fc) < options.f_tol || (b - a) < options.tol
			iterConverged = iter
			converged = true
			break
		end

	end

	if converged
		info = string("Converged in ", iterConverged, " iterations.")
		return Results(c,fc,iterConverged,f_calls,info)
	else
		info = string("Did not converge in ", options.iter, " iterations.")
		return Results(c,fc,options.iter,f_calls,info)
	end

end

end

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
module Autodiff

# package code goes here
include("auto-diff.jl")

end # module


text:
"""
    $SIGNATURES

Return true or false depending on whether child cliques are all up solved.
"""
function areCliqChildrenAllUpSolved(treel::AbstractBayesTree,
                                    prnt::TreeClique)::Bool
  #
  for ch in getChildren(treel, prnt)
    if !isCliqUpSolved(ch)
      return false
    end
  end
  return true
end

"""
    $SIGNATURES

Return `true` if any of the children cliques have status `:needdownmsg`.
"""
function doAnyChildrenNeedDwnMsg(children::Vector{TreeClique})::Bool
  for ch in children
    if getCliqueStatus(ch) == :needdownmsg
      return true
    end
  end
  return false
end

function doAnyChildrenNeedDwnMsg(tree::AbstractBayesTree, cliq::TreeClique)::Bool
  doAnyChildrenNeedDwnMsg( getChildren(tree, cliq) )
end

"""
    $SIGNATURES

Return true if has parent with status `:needdownmsg`.
"""
function isCliqParentNeedDownMsg(tree::AbstractBayesTree, cliq::TreeClique, logger=ConsoleLogger())
  prnt = getParent(tree, cliq)
  if length(prnt) == 0
    return false
  end
  prstat = getCliqueStatus(prnt[1])
  with_logger(logger) do
    @info "$(current_task()) Clique $(cliq.index), isCliqParentNeedDownMsg -- parent status: $(prstat)"
  end
  return prstat == :needdownmsg
end

"""
    $SIGNATURES

Wait here if all siblings and the parent status are `:needdownmsg`.
Return true when parent is `INITIALIZED` after all were `:needdownmsg`

Notes
- used for regulating long need down message chains.
- exit strategy is parent becomes status `INITIALIZED`.
"""
function blockCliqSiblingsParentNeedDown( tree::AbstractBayesTree,
                                          cliq::TreeClique,
                                          prnt_::TreeClique;
                                          logger=ConsoleLogger())
  #

  allneeddwn = true
  prstat = getCliqueStatus(prnt_)
  if prstat == :needdownmsg
    for ch in getChildren(tree, prnt_)
      chst = getCliqueStatus(ch)
      if chst != :needdownmsg
        allneeddwn = false
        break;
      end
    end

    if allneeddwn
      # do actual fetch
      prtmsg = fetchDwnMsgConsolidated(prnt_).status
      if prtmsg == INITIALIZED
        return true
      end
    end
  end
  return false
end