text:
module Entry

import Base: thispatch, nextpatch, nextminor, nextmajor, check_new_version
import ..Git, ..Reqs, ..Read, ..Query, ..Resolve, ..Cache, ..Write, ..GitHub
using ..Types

macro recover(ex)
    quote
        try $(esc(ex))
        catch err
            show(err)
        end
    end
end

function edit(f::Function, pkg::String, args...)
    r = Reqs.read("REQUIRE")
    reqs = Reqs.parse(r)
    avail = Read.available()
    !haskey(avail,pkg) && !haskey(reqs,pkg) && return false
    rʹ = f(r,pkg,args...)
    rʹ == r && return false
    reqsʹ = Reqs.parse(rʹ)
    reqsʹ != reqs && resolve(reqsʹ,avail)
    Reqs.write("REQUIRE",rʹ)
    info("Package database updated")
    return true
end

function edit()
    editor = get(ENV,"VISUAL",get(ENV,"EDITOR",nothing))
    editor != nothing ||
        error("set the EDITOR environment variable to an edit command")
    editor = Base.shell_split(editor)
    reqs = Reqs.parse("REQUIRE")
    run(`$editor REQUIRE`)
    reqsʹ = Reqs.parse("REQUIRE")
    reqs == reqsʹ && return info("Nothing to be done")
    info("Computing changes...")
    resolve(reqsʹ)
end

function add(pkg::String, vers::VersionSet)
    outdated = :maybe
    @sync begin
        @async if !edit(Reqs.add,pkg,vers)
            ispath(pkg) || error("unknown package $pkg")
            info("Nothing to be done")
        end
        branch = Pkg.META_BRANCH
        if Git.branch(dir="METADATA") == branch
            if !Git.success(`diff --quiet origin/$branch`, dir="METADATA")
                outdated = :yes
            else
                try
                    run(Git.cmd(`fetch -q --all`, dir="METADATA") |>DevNull .>DevNull)
                    outdated = Git.success(`diff --quiet origin/$branch`, dir="METADATA") ?
                        (:no) : (:yes)
                end
            end
        else
            outdated = :no # user is doing something funky with METADATA
        end
    end
    if outdated != :no
        is = outdated == :yes ? "is" : "might be"
        info("METADATA $is out-of-date — you may not have the latest version of $pkg")
        info("Use `Pkg.update()` to get the latest versions of your packages")
    end
end
add(pkg::String, vers::VersionNumber...) = add(pkg,VersionSet(vers...))

function rm(pkg::String)
    edit(Reqs.rm,pkg) && return
    ispath(pkg) || return info("Nothing to be done")
    info("Removing $pkg (unregistered)")
    Write.remove(pkg)
end

available() = sort!(ASCIIString[keys(Read.available())...], by=lowercase)

function available(pkg::String)
    avail = Read.available(pkg)
    if !isempty(avail) || Read.isinstalled(pkg)
        return sort!([keys(avail)...])
    end
    error("$pkg is not a package (not registered or installed)")
end

function installed()
    pkgs = Dict{ASCIIString,VersionNumber}()
    for (pkg,(ver,fix)) in Read.installed()
        pkgs[pkg] = ver
    end
    return pkgs
end

function installed(pkg::String)
    avail = Read.available(pkg)
    Read.isinstalled(pkg) && return Read.installed_version(pkg,avail)
    isempty(avail) && error("$pkg is not a package (not registered or installed)")
    return nothing # registered but not installed
end

function status(io::IO)
    reqs = Reqs.parse("REQUIRE")
    instd = Read.installed()
    required = sort!([keys(reqs)...])
    if !isempty(required)
        println(io, "$(length(required)) required packages:")
        for pkg in required
            ver,fix = pop!(instd,pkg)
            status(io,pkg,ver,fix)
        end
    end
    additional = sort!([keys(instd)...])
    if !isempty(additional)
        println(io, "$(length(additional)) additional packages:")
        for pkg in additional
            ver,fix = instd[pkg]
            status(io,pkg,ver,fix)
        end
    end
    if isempty(required) && isempty(additional)
        println(io, "No packages installed")
    end
end
# TODO: status(io::IO, pkg::String)

function status(io::IO, pkg::String, ver::VersionNumber, fix::Bool)
    @printf io " - %-29s " pkg
    fix || return println(io,ver)
    @printf io "%-19s" ver
    if ispath(pkg,".git")
        print(io, Git.attached(dir=pkg) ? Git.branch(dir=pkg) : Git.head(dir=pkg)[1:8])
        attrs = String[]
        isfile("METADATA",pkg,"url") || push!(attrs,"unregistered")
        Git.dirty(dir=pkg) && push!(attrs,"dirty")
        isempty(attrs) || print(io, " (",join(attrs,", "),")")
    else
        print(io, "non-repo (unregistered)")
    end
    println(io)
end

function clone(url::String, pkg::String)
    info("Cloning $pkg from $url")
    ispath(pkg) && error("$pkg already exists")
    try
        Git.run(`clone -q $url $pkg`)
        Git.set_remote_url(url, dir=pkg)
    catch
        run(`rm -rf $pkg`)
        rethrow()
    end
    isempty(Reqs.parse("$pkg/REQUIRE")) && return
    info("Computing changes...")
    resolve()
end

function clone(url_or_pkg::String)
    urlpath = joinpath("METADATA",url_or_pkg,"url")
    if isfile(urlpath)
        pkg = url_or_pkg
        url = readchomp(urlpath)
        # TODO: Cache.prefetch(pkg,url)
    else
        url = url_or_pkg
        m = match(r"(?:^|[/\\])(\w+?)(?:\.jl)?(?:\.git)?$", url)
        m != nothing || error("can't determine package name from URL: $url")
        pkg = m.captures[1]
    end
    clone(url,pkg)
end

function _checkout(pkg::String, what::String, merge::Bool=false, pull::Bool=false)
    Git.transact(dir=pkg) do
        Git.dirty(dir=pkg) && error("$pkg is dirty, bailing")
        Git.run(`checkout -q $what`, dir=pkg)
        merge && Git.run(`merge -q --ff-only $what`, dir=pkg)
        if pull
            info("Pulling $pkg latest $what...")
            Git.run(`pull -q --ff-only`, dir=pkg)
        end
        resolve()
    end
end

function checkout(pkg::String, branch::String, merge::Bool, pull::Bool)
    ispath(pkg,".git") || error("$pkg is not a git repo")
    info("Checking out $pkg $branch...")
    _checkout(pkg,branch,merge,pull)
end

function free(pkg::String)
    ispath(pkg,".git") || error("$pkg is not a git repo")
    Read.isinstalled(pkg) || error("$pkg cannot be freed – not an installed package")
    avail = Read.available(pkg)
    isempty(avail) && error("$pkg cannot be freed – not a registered package")
    Git.dirty(dir=pkg) && error("$pkg cannot be freed – repo is dirty")
    info("Freeing $pkg")
    vers = sort!([keys(avail)...], rev=true)
    while true
        for ver in vers
            sha1 = avail[ver].sha1
            Git.iscommit(sha1, dir=pkg) || continue
            return _checkout(pkg,sha1)
        end
        isempty(Cache.prefetch(pkg, Read.url(pkg), [a.sha1 for (v,a)=avail])) && continue
        error("can't find any registered versions of $pkg to checkout")
    end
end

function pin(pkg::String, head::String)
    ispath(pkg,".git") || error("$pkg is not a git repo")
    branch = "pinned.$(head[1:8]).tmp"
    rslv = (head != Git.head(dir=pkg))
    info("Creating $pkg branch $branch")
    Git.run(`checkout -q -B $branch $head`, dir=pkg)
    rslv ? resolve() : nothing
end
pin(pkg::String) = pin(pkg,Git.head(dir=pkg))

function pin(pkg::String, ver::VersionNumber)
    ispath(pkg,".git") || error("$pkg is not a git repo")
    Read.isinstalled(pkg) || error("$pkg cannot be pinned – not an installed package".tmp)
    avail = Read.available(pkg)
    isempty(avail) && error("$pkg cannot be pinned – not a registered package".tmp)
    haskey(avail,ver) || error("$pkg – $ver is not a registered version")
    pin(pkg,avail[ver].sha1)
end

function update(branch::String)
    info("Updating METADATA...")
    cd("METADATA") do
        if Git.branch() != branch
            Git.dirty() && error("METADATA is dirty and not on $branch, bailing")
            Git.attached() || error("METADATA is detached not on $branch, bailing")
            Git.run(`fetch -q --all`)
            Git.run(`checkout -q HEAD^0`)
            Git.run(`branch -f $branch refs/remotes/origin/$branch`)
            Git.run(`checkout -q $branch`)
        end
        # TODO: handle merge conflicts
        Base.with_env("GIT_MERGE_AUTOEDIT","no") do
            Git.run(`pull -q`, out=DevNull)
        end
    end
    avail = Read.available()
    # this has to happen before computing free/fixed
    @sync for pkg in filter!(Read.isinstalled,[keys(avail)...])
        @async Cache.prefetch(pkg, Read.url(pkg), [a.sha1 for (v,a)=avail[pkg]])
    end
    instd = Read.installed(avail)
    free = Read.free(instd)
    @sync for (pkg,ver) in free
        @async Cache.prefetch(pkg, Read.url(pkg), [a.sha1 for (v,a)=avail[pkg]])
    end
    fixed = Read.fixed(avail,instd)
    @sync for (pkg,ver) in fixed
        ispath(pkg,".git") || continue
        @async begin
            if Git.attached(dir=pkg) && !Git.dirty(dir=pkg)
                info("Updating $pkg...")
                @recover begin
                    Git.run(`fetch -q --all`, dir=pkg)
                    Git.success(`pull -q --ff-only`, dir=pkg) # suppress output
                end
            end
            if haskey(avail,pkg)
                Cache.prefetch(pkg, Read.url(pkg), [a.sha1 for (v,a)=avail[pkg]])
            end
        end
    end
    info("Computing changes...")
    resolve(Reqs.parse("REQUIRE"), avail, instd, fixed, free)
    # Don't use instd here since it may have changed
    updatehook(sort!([keys(installed())...]))
end

function pull_request(dir::String, commit::String="", url::String="")
    commit = isempty(commit) ? Git.head(dir=dir) :
        Git.readchomp(`rev-parse --verify $commit`, dir=dir)
    isempty(url) && (url = Git.readchomp(`config remote.origin.url`, dir=dir))
    m = match(Git.GITHUB_REGEX, url)
    m == nothing && error("not a GitHub repo URL, can't make a pull request: $url")
    owner, repo = m.captures[2:3]
    user = GitHub.user()
    info("Forking $owner/$repo to $user")
    response = GitHub.fork(owner,repo)
    fork = response["ssh_url"]
    branch = "pull-request/$(commit[1:8])"
    info("Pushing changes as branch $branch")
    Git.run(`push -q $fork $commit:refs/heads/$branch`, dir=dir)
    pr_url = "$(response["html_url"])/compare/$branch"
    @osx? run(`open $pr_url`) : info("To create a pull-request open:\n\n  $pr_url\n")
end

function submit(pkg::String, commit::String="")
    urlpath = joinpath("METADATA",pkg,"url")
    url = ispath(urlpath) ? readchomp(urlpath) : ""
    pull_request(pkg, commit, url)
end

function publish(branch::String)
    Git.branch(dir="METADATA") == branch ||
        error("METADATA must be on $branch to publish changes")
    Git.success(`push -q -n origin $branch`, dir="METADATA") ||
        error("METADATA is behind origin/$branch – run `Pkg.update()` before publishing")
    Git.run(`fetch -q`, dir="METADATA")
    info("Validating METADATA")
    check_metadata()
    tags = Dict{ASCIIString,Vector{ASCIIString}}()
    Git.run(`update-index -q --really-refresh`, dir="METADATA")
    cmd = `diff --name-only --diff-filter=AMR origin/$branch HEAD --`
    for line in eachline(Git.cmd(cmd, dir="METADATA"))
        path = chomp(line)
        m = match(r"^(.+?)/versions/([^/]+)/sha1$", path)
        m != nothing && ismatch(Base.VERSION_REGEX, m.captures[2]) || continue
        pkg, ver = m.captures; ver = convert(VersionNumber,ver)
        sha1 = readchomp(joinpath("METADATA",path))
        if Git.success(`cat-file -e origin/$branch:$path`, dir="METADATA")
            old = Git.readchomp(`cat-file blob origin/$branch:$path`, dir="METADATA")
            old == sha1 || error("$pkg v$ver SHA1 changed in METADATA – refusing to publish")
        end
        any(split(Git.readall(`tag --contains $sha1`, dir=pkg))) do tag
            ver == convert(VersionNumber,tag) || return false
            haskey(tags,pkg) || (tags[pkg] = ASCIIString[])
            push!(tags[pkg], tag)
            return true
        end || error("$pkg v$ver is incorrectly tagged – $sha1 expected")
    end
    isempty(tags) && info("No new package versions to publish")
    @sync for pkg in sort!([keys(tags)...])
        @async begin
            forced = ASCIIString[]
            unforced = ASCIIString[]
            for tag in tags[pkg]
                ver = convert(VersionNumber,tag)
                push!(isrewritable(ver) ? forced : unforced, tag)
            end
            if !isempty(forced)
                info("Pushing $pkg temporary tags: ", join(forced,", "))
                refspecs = ["refs/tags/$tag:refs/tags/$tag" for tag in forced]
                Git.run(`push -q --force origin $refspecs`, dir=pkg)
            end
            if !isempty(unforced)
                info("Pushing $pkg permanent tags: ", join(unforced,", "))
                refspecs = ["refs/tags/$tag:refs/tags/$tag" for tag in unforced]
                Git.run(`push -q origin $refspecs`, dir=pkg)
            end
        end
    end
    info("Submitting METADATA changes")
    pull_request("METADATA")
end

function resolve(
    reqs  :: Dict = Reqs.parse("REQUIRE"),
    avail :: Dict = Read.available(),
    instd :: Dict = Read.installed(avail),
    fixed :: Dict = Read.fixed(avail,instd),
    have  :: Dict = Read.free(instd),
)
    reqs = Query.requirements(reqs,fixed,avail)
    deps, conflicts = Query.dependencies(avail,fixed)

    for pkg in keys(reqs)
        if !haskey(deps,pkg)
            error("$pkg's requirements can't be satisfied because of the following fixed packages: ",
                   join(conflicts[pkg], ", ", " and "))
        end
    end

    Query.check_requirements(reqs,deps,fixed)

    deps = Query.prune_dependencies(reqs,deps)
    want = Resolve.resolve(reqs,deps)

    # compare what is installed with what should be
    changes = Query.diff(have, want, avail, fixed)
    isempty(changes) && return info("No packages to install, update or remove")

    # prefetch phase isolates network activity, nothing to roll back
    missing = {}
    for (pkg,(ver1,ver2)) in changes
        vers = ASCIIString[]
        ver1 !== nothing && push!(vers,Git.head(dir=pkg))
        ver2 !== nothing && push!(vers,Read.sha1(pkg,ver2))
        append!(missing,
            map(sha1->(pkg,(ver1,ver2),sha1),
                Cache.prefetch(pkg, Read.url(pkg), vers)))
    end
    if !isempty(missing)
        msg = "Missing package versions (possible metadata misconfiguration):"
        for (pkg,ver,sha1) in missing
            msg *= "  $pkg v$ver [$sha1[1:10]]\n"
        end
        error(msg)
    end

    # try applying changes, roll back everything if anything fails
    changed = {}
    try
        for (pkg,(ver1,ver2)) in changes
            if ver1 === nothing
                info("Installing $pkg v$ver2")
                Write.install(pkg, Read.sha1(pkg,ver2))
            elseif ver2 == nothing
                info("Removing $pkg v$ver1")
                Write.remove(pkg)
            else
                up = ver1 <= ver2 ? "Up" : "Down"
                info("$(up)grading $pkg: v$ver1 => v$ver2")
                Write.update(pkg, Read.sha1(pkg,ver2))
            end
            push!(changed,(pkg,(ver1,ver2)))
        end
    catch
        for (pkg,(ver1,ver2)) in reverse!(changed)
            if ver1 == nothing
                info("Rolling back install of $pkg")
                @recover Write.remove(pkg)
            elseif ver2 == nothing
                info("Rolling back deleted $pkg to v$ver1")
                @recover Write.install(pkg, Read.sha1(pkg,ver1))
            else
                info("Rolling back $pkg from v$ver2 to v$ver1")
                @recover Write.update(pkg, Read.sha1(pkg,ver1))
            end
        end
        rethrow()
    end
    # re/build all updated/installed packages
    build(map(x->x[1],filter(x->x[2][2]!=nothing,changes)))
end

function write_tag_metadata(pkg::String, ver::VersionNumber, commit::String, force::Bool=false)
    cmd = Git.cmd(`cat-file blob $commit:REQUIRE`, dir=pkg)
    reqs = success(cmd) ? Reqs.read(cmd) : Reqs.Line[]
    cd("METADATA") do
        d = joinpath(pkg,"versions",string(ver))
        mkpath(d)
        sha1file = joinpath(d,"sha1")
        if !force && ispath(sha1file)
            current = readchomp(sha1file)
            current == commit ||
                error("$pkg v$ver is already registered as $current, bailing")
        end
        open(io->println(io,commit), sha1file, "w")
        Git.run(`add $sha1file`)
        reqsfile = joinpath(d,"requires")
        if isempty(reqs)
            ispath(reqsfile) && Git.run(`rm -f -q $reqsfile`)
        else
            Reqs.write(reqsfile,reqs)
            Git.run(`add $reqsfile`)
        end
    end
    return nothing
end

function register(pkg::String, url::String)
    ispath(pkg,".git") || error("$pkg is not a git repo")
    isfile("METADATA",pkg,"url") && error("$pkg already registered")
    tags = split(Git.readall(`tag -l v*`, dir=pkg))
    filter!(tag->ismatch(Base.VERSION_REGEX,tag), tags)
    versions = [
        convert(VersionNumber,tag) =>
        Git.readchomp(`rev-parse --verify $tag^{commit}`, dir=pkg)
        for tag in tags
    ]
    Git.transact(dir="METADATA") do
        cd("METADATA") do
            info("Registering $pkg at $url")
            mkdir(pkg)
            path = joinpath(pkg,"url")
            open(io->println(io,url), path, "w")
            Git.run(`add $path`)
        end
        vers = sort!([keys(versions)...])
        for ver in vers
            info("Tagging $pkg v$ver")
            write_tag_metadata(pkg,ver,versions[ver])
        end
        if Git.staged(dir="METADATA")
            info("Committing METADATA for $pkg")
            msg = "Register $pkg"
            if !isempty(versions)
                msg *= ": $(join(map(v->"v$v", vers),", "))"
            end
            Git.run(`commit -q -m $msg -- $pkg`, dir="METADATA")
        else
            info("No METADATA changes to commit")
        end
    end
end

function register(pkg::String)
    Git.success(`config remote.origin.url`, dir=pkg) ||
        error("$pkg: no URL configured")
    url = Git.readchomp(`config remote.origin.url`, dir=pkg)
    register(pkg,Git.normalize_url(url))
end

function isrewritable(v::VersionNumber)
    thispatch(v)==v"0" ||
    length(v.prerelease)==1 && isempty(v.prerelease[1]) ||
    length(v.build)==1 && isempty(v.build[1])
end

nextbump(v::VersionNumber) = isrewritable(v) ? v : nextpatch(v)

function tag(pkg::String, ver::Union(Symbol,VersionNumber), force::Bool=false, commit::String="HEAD")
    ispath(pkg,".git") || error("$pkg is not a git repo")
    Git.dirty(dir=pkg) &&
        error("$pkg is dirty – commit or stash changes to tag")
    Git.dirty(pkg, dir="METADATA") &&
        error("METADATA/$pkg is dirty – commit or stash changes to tag")
    commit = Git.readchomp(`rev-parse $commit`, dir=pkg)
    registered = isfile("METADATA",pkg,"url")
    if !force
        if registered
            avail = Read.available(pkg)
            existing = VersionNumber[keys(Read.available(pkg))...]
            ancestors = filter(v->Git.is_ancestor_of(avail[v].sha1,commit,dir=pkg), existing)
        else
            tags = split(Git.readall(`tag -l v*`, dir=pkg))
            filter!(tag->ismatch(Base.VERSION_REGEX,tag), tags)
            existing = VersionNumber[tags...]
            filter!(tags) do tag
                sha1 = Git.readchomp(`rev-parse --verify $tag^{commit}`, dir=pkg)
                Git.is_ancestor_of(sha1,commit,dir=pkg)
            end
            ancestors = VersionNumber[tags...]
        end
        sort!(existing)
        if isa(ver,Symbol)
            prv = isempty(existing) ? v"0" :
                  isempty(ancestors) ? maximum(existing) : maximum(ancestors)
            ver = (ver == :bump ) ? nextbump(prv)  :
                  (ver == :patch) ? nextpatch(prv) :
                  (ver == :minor) ? nextminor(prv) :
                  (ver == :major) ? nextmajor(prv) :
                                    error("invalid version selector: $ver")
        end
        isrewritable(ver) && filter!(v->v!=ver,existing)
        check_new_version(existing,ver)
    end
    # TODO: check that SHA1 isn't the same as another version
    info("Tagging $pkg v$ver")
    opts = ``
    if force || isrewritable(ver)
        opts = `$opts --force`
    end
    if !isrewritable(ver)
        opts = `$opts --annotate --message "$pkg v$ver [$(commit[1:10])]"`
    end
    Git.run(`tag $opts v$ver $commit`, dir=pkg, out=DevNull)
    registered || return
    try
        Git.transact(dir="METADATA") do
            write_tag_metadata(pkg,ver,commit,force)
            if Git.staged(dir="METADATA")
                info("Committing METADATA for $pkg")
                Git.run(`commit -q -m "Tag $pkg v$ver" -- $pkg`, dir="METADATA")
            else
                info("No METADATA changes to commit")
            end
        end
    catch
        Git.run(`tag -d v$ver`, dir=pkg)
        rethrow()
    end
end

function check_metadata()
    avail = Read.available()
    deps, conflicts = Query.dependencies(avail)

    for (dp,dv) in deps, (v,a) in dv, p in keys(a.requires)
        haskey(deps, p) || error("package $dp v$v requires a non-registered package: $p")
    end

    problematic = Resolve.sanity_check(deps)
    if !isempty(problematic)
        msg = "packages with unsatisfiable requirements found:\n"
        for (p, vn, rp) in problematic
            msg *= "    $p v$vn – no valid versions exist for package $rp\n"
        end
        error(msg)
    end
    return
end

function warnbanner(msg...; label="[ WARNING ]", prefix="")
    warn(prefix="", Base.cpad(label,Base.tty_cols(),"="))
    println(STDERR)
    warn(prefix=prefix, msg...)
    println(STDERR)
    warn(prefix="", "="^Base.tty_cols())
end

function build!(pkgs::Vector, errs::Dict, seen::Set=Set())
    for pkg in pkgs
        pkg in seen && continue
        build!(Read.requires_list(pkg),errs,push!(seen,pkg))
        path = abspath(pkg,"deps","build.jl")
        isfile(path) || continue
        info("Building $pkg")
        cd(dirname(path)) do
            try evalfile(path)
            catch err
                warnbanner(err, label="[ ERROR: $pkg ]")
                errs[pkg] = err
            end
        end
    end
end

function build(pkgs::Vector)
    errs = Dict()
    build!(pkgs,errs)
    isempty(errs) && return
    println(STDERR)
    warnbanner(label="[ BUILD ERRORS ]", """
    WARNING: $(join(map(x->x[1],errs),", "," and ")) had build errors.

     - packages with build errors remain installed in $(pwd())
     - build a package and all its dependencies with `Pkg.build(pkg)`
     - build a single package by running its `deps/build.jl` script
    """)
end
build() = build(sort!([keys(installed())...]))

function updatehook!(pkgs::Vector, errs::Dict, seen::Set=Set())
    for pkg in pkgs
        pkg in seen && continue
        updatehook!(Read.requires_list(pkg),errs,push!(seen,pkg))
        path = abspath(pkg,"deps","update.jl")
        isfile(path) || continue
        info("Running update script for $pkg")
        cd(dirname(path)) do
            try evalfile(path)
            catch err
                warnbanner(err, label="[ ERROR: $pkg ]")
                errs[pkg] = err
            end
        end
    end
end

function updatehook(pkgs::Vector)
    errs = Dict()
    updatehook!(pkgs,errs)
    isempty(errs) && return
    println(STDERR)
    warnbanner(label="[ UPDATE ERRORS ]", """
    WARNING: $(join(map(x->x[1],errs),", "," and ")) had update errors.

     - Unrelated packages are unaffected
     - To retry, run Pkg.update() again
    """)
end

const JULIA = joinpath(JULIA_HOME, "julia")

function test!(pkg::String, errs::Vector{String}, notests::Vector{String})
    const reqs_path = abspath(pkg,"test","REQUIRE")
    if isfile(reqs_path)
        const tests_require = Reqs.parse(reqs_path)
        if (!isempty(tests_require))
            info("Computing test dependencies for $pkg...")
            resolve(tests_require)
        end
    end
    const test_path = abspath(pkg,"test","runtests.jl")
    if isfile(test_path)
        info("Testing $pkg")
        cd(dirname(test_path)) do
            try
                run(`$JULIA $test_path`)
                info("$pkg tests passed")
            catch err
                warnbanner(err, label="[ ERROR: $pkg ]")
                push!(errs,pkg)
            end
        end
    else
        push!(notests,pkg)
    end
    resolve()
end

function test(pkgs::Vector{String})
    errs = String[]
    notests = String[]
    for pkg in pkgs
        test!(pkg,errs,notests)
    end
    if !isempty(errs) || !isempty(notests)
        messages = String[]
        isempty(errs) || push!(messages, "$(join(errs,", "," and ")) had test errors")
        isempty(notests) || push!(messages, "$(join(notests,", "," and ")) did not provide a test/runtests.jl file")
        error(join(messages, "and"))
    end
end

test() = test(sort!(String[keys(installed())...]))

end # module

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
import ...BalanceLaws:
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

using Printf

@inline function set_fast_for_stepping!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgFast,
    Qfast,
    S_fast,
    slow_dt,
    rkC,
    rkW,
    s,
    nStages,
    fast_time_rec,
    fast_steps,
)
    FT = typeof(slow_dt)

    #- inverse ratio of additional fast time steps (for weighted average)
    #  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
    add = slow.add_fast_substeps

    #- set time-step
    #  Warning: only make sense for LS3NRK33Heuns
    #  where 12 is lowest common mutiple (LCM) of all RK-Coeff inverse
    fast_dt = fast_time_rec[1]
    steps = fast_dt > 0 ? ceil(Int, slow_dt / fast_dt / FT(12)) : 1
    ntsFull = 12 * steps
    fast_dt = slow_dt / ntsFull
    add = add > 0 ? floor(Int, ntsFull / add) : 0

    #- time to start fast time-stepping (fast_time_rec[3]) for this stage:
    if s == nStages
        #  Warning: only works with few RK-scheme such as LS3NRK33Heuns
        fast_time_rec[3] = rkW[1]
        fract_dt = (1 - fast_time_rec[3])
        fast_time_rec[3] *= slow_dt
        fast_steps[2] = 1
    else
        fast_time_rec[3] = 0.0
        fract_dt = rkC[s + 1] - fast_time_rec[3]
        fast_steps[2] = 0
    end

    #- set number of sub-steps we need
    #  will time-average fast over: fast_steps[1] , fast_steps[3]
    #  centered on fract_dt*slow_dt which corresponds to advance in time of slow
    steps = ceil(Int, fract_dt * slow_dt / fast_dt)
    add = min(add, steps - 1)
    fast_steps[1] = steps - add
    fast_steps[3] = steps + add
    fast_time_rec[1] = fract_dt * slow_dt / steps

    #- select which fast time-step (fast_steps[2]) solution to save for next time-step
    #  Warning: only works with few RK-scheme such as LS3NRK33Heuns
    fast_steps[2] *= steps
    if s == 1
        fast_steps[2] = round(Int, ntsFull * rkW[1])
    end
    # @printf("Update @ s= %i : frac_dt = %.6f , dt_fast = %.1f , steps= %i , add= %i\n",
    #          s, fract_dt, fast_time_rec[1], steps, add)
    # println(" fast_time_rec = ",fast_time_rec)
    # println(" fast_steps = ",fast_steps)

    # set starting point for fast-state solution
    #  Warning: only works with few RK-scheme such as LS3NRK33Heuns
    if s == 1
        S_fast.η .= Qfast.η
        S_fast.U .= Qfast.U
    elseif s == nStages
        Qfast.η .= dgFast.state_auxiliary.η_s
        Qfast.U .= dgFast.state_auxiliary.U_s
    else
        Qfast.η .= S_fast.η
        Qfast.U .= S_fast.U
    end

    # initialise cumulative arrays
    fast_time_rec[2] = 0.0
    dgFast.state_auxiliary.η_c .= -0
    dgFast.state_auxiliary.U_c .= (@SVector [-0, -0])'

    return nothing
end

@inline function initialize_fast_state!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    slow_dt,
    fast_time_rec,
    fast_steps;
    firstStage = false,
)

    #- inverse ratio of additional fast time steps (for weighted average)
    #  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
    add = slow.add_fast_substeps

    #- set time-step and number of sub-steps we need
    #  will time-average fast over: fast_steps[1] , fast_steps[3]
    #  centered on fast_steps[2] which corresponds to advance in time of slow
    fast_dt = fast_time_rec[1]
    if add == 0
        steps = fast_dt > 0 ? ceil(Int, slow_dt / fast_dt) : 1
        fast_steps[1:3] = [1 1 1] * steps
    else
        steps = fast_dt > 0 ? ceil(Int, slow_dt / fast_dt / add) : 1
        fast_steps[2] = add * steps
        fast_steps[1] = (add - 1) * steps
        fast_steps[3] = (add + 1) * steps
    end
    fast_time_rec[1] = slow_dt / fast_steps[2]
    fast_time_rec[2] = 0.0
    # @printf("Update: frac_dt = %.1f , dt_fast = %.1f , nsubsteps= %i\n",
    #          slow_dt,fast_time_rec[1],fast_steps[3])
    # println(" fast_steps = ",fast_steps)

    dgFast.state_auxiliary.η_c .= -0
    dgFast.state_auxiliary.U_c .= (@SVector [-0, -0])'

    # set fast-state to previously stored value
    if !firstStage
        Qfast.η .= dgFast.state_auxiliary.η_s
        Qfast.U .= dgFast.state_auxiliary.U_s
    end

    return nothing
end

@inline function initialize_adjustment!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
)
    ## reset tendency adjustment before calling Baroclinic Model
    dgSlow.state_auxiliary.ΔGu .= 0

    return nothing
end

@inline function tendency_from_slow_to_fast!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    dQslow2fast,
)
    FT = eltype(Qslow)

    # integrate the tendency
    tendency_dg = dgSlow.modeldata.tendency_dg
    tend = tendency_dg.balance_law
    grid = dgSlow.grid
    elems = grid.topology.elems
    update_auxiliary_state!(tendency_dg, tend, dQslow2fast, 0, elems)

    info = basic_grid_info(dgSlow)
    Nqh, Nqk = info.Nqh, info.Nqk
    nelemv, nelemh = info.nvertelem, info.nhorzelem
    nrealelemh = info.nhorzrealelem

    ## get top value (=integral over full depth) of ∫du
    nb_aux_tnd = number_states(tend, Auxiliary())
    data_tnd = reshape(
        tendency_dg.state_auxiliary.data,
        Nqh,
        Nqk,
        nb_aux_tnd,
        nelemv,
        nelemh,
    )
    index_∫du = varsindex(vars_state(tend, Auxiliary(), FT), :∫du)
    flat_∫du = @view data_tnd[:, end, index_∫du, end, 1:nrealelemh]

    ## copy into Gᵁ of dgFast
    nb_aux_fst = number_states(fast, Auxiliary())
    data_fst = reshape(dgFast.state_auxiliary.data, Nqh, nb_aux_fst, nelemh)
    index_Gᵁ = varsindex(vars_state(fast, Auxiliary(), FT), :Gᵁ)
    boxy_Gᵁ = @view data_fst[:, index_Gᵁ, 1:nrealelemh]
    boxy_Gᵁ .= flat_∫du

    ## scale by -1/H and copy back to ΔGu
    # note: since tendency_dg.state_auxiliary.∫du is not used after this, could be
    #   re-used to store a 3-D copy of "-Gu"
    nb_aux_slw = number_states(slow, Auxiliary())
    data_slw = reshape(
        dgSlow.state_auxiliary.data,
        Nqh,
        Nqk,
        nb_aux_slw,
        nelemv,
        nelemh,
    )
    index_ΔGu = varsindex(vars_state(slow, Auxiliary(), FT), :ΔGu)
    boxy_ΔGu = @view data_slw[:, :, index_ΔGu, :, 1:nrealelemh]
    boxy_∫du = @view data_tnd[:, end:end, index_∫du, end:end, 1:nrealelemh]
    boxy_ΔGu .= -boxy_∫du / slow.problem.H

    return nothing
end

@inline function cummulate_fast_solution!(
    fast::BarotropicModel,
    dgFast,
    Qfast,
    fast_time,
    fast_dt,
    substep,
    fast_steps,
    fast_time_rec,
)
    #- might want to use some of the weighting factors: weights_η & weights_U
    #- should account for case where fast_dt < fast.param.dt

    # cumulate Fast solution:
    if substep >= fast_steps[1]
        dgFast.state_auxiliary.U_c .+= Qfast.U
        dgFast.state_auxiliary.η_c .+= Qfast.η
        fast_time_rec[2] += 1.0
    end

    # save mid-point solution to start from the next time-step
    if substep == fast_steps[2]
        dgFast.state_auxiliary.U_s .= Qfast.U
        dgFast.state_auxiliary.η_s .= Qfast.η
    end

    return nothing
end

@inline function reconcile_from_fast_to_slow!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    fast_time_rec;
    lastStage = false,
)
    FT = eltype(Qslow)
    info = basic_grid_info(dgSlow)
    Nqh, Nqk = info.Nqh, info.Nqk
    nelemv, nelemh = info.nvertelem, info.nhorzelem
    nrealelemh = info.nhorzrealelem
    grid = dgSlow.grid
    elems = grid.topology.elems

    # need to calculate int_u using integral kernels
    # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow)

    ## get time weighted averaged out of cumulative arrays
    dgFast.state_auxiliary.U_c .*= 1 / fast_time_rec[2]
    dgFast.state_auxiliary.η_c .*= 1 / fast_time_rec[2]
    # @printf(" reconcile_from_fast_to_slow! @ s= %i : time_Count = %6.3f\n",
    #          0, fast_time_rec[2])
    #   #      s, fast_time_rec[2])

    ## Compute: \int_{-H}^{0} u_slow

    # integrate vertically horizontal velocity
    flowintegral_dg = dgSlow.modeldata.flowintegral_dg
    flowint = flowintegral_dg.balance_law
    update_auxiliary_state!(flowintegral_dg, flowint, Qslow, 0, elems)

    # get top value (=integral over full depth)
    nb_aux_flw = number_states(flowint, Auxiliary())
    data_flw = reshape(
        flowintegral_dg.state_auxiliary.data,
        Nqh,
        Nqk,
        nb_aux_flw,
        nelemv,
        nelemh,
    )
    index_∫u = varsindex(vars_state(flowint, Auxiliary(), FT), :∫u)
    flat_∫u = @view data_flw[:, end, index_∫u, end, 1:nrealelemh]

    ## substract ∫u from U and divide by H

    # Δu is a place holder for 1/H * (Ū - ∫u)
    Δu = dgFast.state_auxiliary.Δu
    Δu .= dgFast.state_auxiliary.U_c

    nb_aux_fst = number_states(fast, Auxiliary())
    data_fst = reshape(dgFast.state_auxiliary.data, Nqh, nb_aux_fst, nelemh)
    index_Δu = varsindex(vars_state(fast, Auxiliary(), FT), :Δu)
    boxy_Δu = @view data_fst[:, index_Δu, 1:nrealelemh]
    boxy_Δu .-= flat_∫u
    boxy_Δu ./= slow.problem.H

    ## apply the 2D correction to the 3D solution
    nb_cons_slw = number_states(slow, Prognostic())
    data_slw = reshape(Qslow.data, Nqh, Nqk, nb_cons_slw, nelemv, nelemh)
    index_u = varsindex(vars_state(slow, Prognostic(), FT), :u)
    boxy_u = @view data_slw[:, :, index_u, :, 1:nrealelemh]
    boxy_u .+= reshape(boxy_Δu, Nqh, 1, 2, 1, nrealelemh)

    ## save Eta from 3D model into η_diag (aux var of 2D model)
    ## and store difference between η from Barotropic Model and η_diag
    ## Note: since 3D η is not used (just in output), only do this at last stage
    ##       (save computation and get Δη diagnose over full time-step)
    if lastStage
        index_η = varsindex(vars_state(slow, Prognostic(), FT), :η)
        boxy_η_3D = @view data_slw[:, :, index_η, :, 1:nrealelemh]
        flat_η = @view data_slw[:, end, index_η, end, 1:nrealelemh]
        index_η_diag = varsindex(vars_state(fast, Auxiliary(), FT), :η_diag)
        boxy_η_diag = @view data_fst[:, index_η_diag, 1:nrealelemh]
        boxy_η_diag .= flat_η
        dgFast.state_auxiliary.Δη .=
            dgFast.state_auxiliary.η_c - dgFast.state_auxiliary.η_diag

        ## copy 2D model Eta over to 3D model
        index_η_c = varsindex(vars_state(fast, Auxiliary(), FT), :η_c)
        boxy_η_2D = @view data_fst[:, index_η_c, 1:nrealelemh]
        boxy_η_3D .= reshape(boxy_η_2D, Nqh, 1, 1, 1, nrealelemh)

        # reset fast-state to end of time-step value
        Qfast.η .= dgFast.state_auxiliary.η_s
        Qfast.U .= dgFast.state_auxiliary.U_s
    end

    return nothing
end

ext: jl


text:
#!/usr/bin/env julia

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "..", "matrix"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "core"))

using MatrixBase
using EigenBase
using MatrixMarket
using LinearAlgebra

using QJuliaSolvers
using QJuliaEnums
using MPI
using Random

#initialize MPI
MPI.Init()

Random.seed!(2018)

#matrix_path = "/home/astrel/data/nasa2146/nasa2146.mtx"
matrix_path = "/home/astrel/data/bcsstk03/bcsstk03.mtx"
#matrix_path = "/home/alex/data/nasa2910/nasa2910.mtx"
#matrix_path = "/home/astrel/data/nasa4704/nasa4704.mtx"
#matrix_path = "/home/astrel/data/smt/smt.mtx"
#matrix_path = "/home/astrel/data/nos4/nos4.mtx"

prec         = Float64
prec_sloppy  = Float64
prec_precond = Float64

csrM = MatrixBase.CSRMat{prec}(matrix_path)
csrMpre = MatrixBase.CSRMat{prec_precond}(matrix_path)

#ilu0csrM = MatrixBase.ilu0(csrM)
#MatrixBase.print_CSRMat_info(ilu0csrM)

data_type = typeof(csrM.csrVals[1])
pre_data_type = typeof(csrMpre.csrVals[1])

b = Vector{data_type}(undef, csrM.Dims[2])
x = Vector{data_type}(undef, csrM.Dims[2])
z = Vector{data_type}(undef, csrM.Dims[2])

x .= 0.0
z .= 0.0
b .= rand(data_type, csrM.Dims[2])

M(b,x) = MatrixBase.csrmv(b, csrM, x)
#or just:
#M=MatrixMarket.mmread(matrix_path)
Mpre(b,z) = MatrixBase.csrmv(b, csrMpre, z)

#setup preconditioner
pre_solv_param = QJuliaSolvers.QJuliaSolverParam_qj()

pre_solv_param.inv_type  = QJuliaEnums.QJULIA_LANMR_INVERTER
pre_solv_param.dtype     = pre_data_type
pre_solv_param.tol       = 1e-2
#
pre_solv_param.maxiter   = 10
pre_solv_param.Nsteps    = 1
pre_solv_param.global_reduction = false

do_ilu = false

if do_ilu == false

K(pPre, rPre) = QJuliaSolvers.solve(pPre, rPre, Mpre, Mpre, pre_solv_param)

else

eigpc = EigenBase.PrecondDescr()
eigpc.pctype = EigenBase.EIGENBASE_PRECOND_ICC
eigpc.rows   = csrM.Dims[1]
eigpc.cols   = csrM.Dims[2]
eigpc.nnz    = length(csrMpre.csrVals)
eigpc.precision = typeof(csrMpre.csrVals[1]) == Float64 ? EigenBase.EIGENBASE_FLOAT_DOUBLE : EigenBase.EIGENBASE_FLOAT_SINGLE
# We need to reformat arrays for the external c/c++ library:

csrRows=Vector{Int32}(undef, length(csrMpre.csrRows))
csrCols=Vector{Int32}(undef, length(csrMpre.csrCols))

csrRows .=@. csrMpre.csrRows - 1
csrCols .=@. csrMpre.csrCols - 1

EigenBase.CreateEigenPreconditioner_qj(csrMpre.csrVals, csrCols, csrRows, eigpc)
K(pPre, rPre) = EigenBase.ApplyEigenPreconditioner_qj(pPre, rPre, eigpc)

end

solv_param = QJuliaSolvers.QJuliaSolverParam_qj()
# Set up parameters
#solv_param.inv_type               = QJuliaEnums.QJULIA_PIPEPCG_INVERTER
solv_param.inv_type               = QJuliaEnums.QJULIA_FCG_INVERTER
solv_param.dtype                  = data_type
solv_param.inv_type_precondition  = pre_solv_param.inv_type
#solv_param.inv_type_precondition  = QJuliaEnums.QJULIA_INVALID_INVERTER
solv_param.dtype_precondition     = pre_solv_param.dtype
solv_param.tol                    = 1e-10
solv_param.delta                  = 1e-1
#
solv_param.maxiter                = 1000
solv_param.Nsteps                 = 2
#nKrylov = 5 is the best choice for the NOTAY option
solv_param.nKrylov                = 5

do_unit_source = true
if do_unit_source == true
  z .=@. 1.0 / sqrt(csrM.Dims[2])
else
  z .= rand(data_type, csrM.Dims[2])
end

M(b,z)

println("Initial residual:  ", norm(b))

if (solv_param.inv_type_precondition == QJuliaEnums.QJULIA_INVALID_INVERTER)
@time QJuliaSolvers.solve(x, b, M, M, solv_param)
else
@time QJuliaSolvers.solve(x, b, M, M, solv_param, K)
end

#compute error:
z .=@. z - x
println("Solution error:  ", norm(z))

#compute residual:
M(z, x)
z .= b - z

println("Solution norm: ", norm(x), ", absolute residual: ", norm(z), " relative ", norm(z) / norm(b))

MPI.Finalize()

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# Dependencies
using BenchmarkTools
include("mccall_utilities.jl")

# Benchmark the QuantEcon default solution
@btime default_quantecon()

# Benchmark the QuantEcon iteration scheme, with update_bellman! replaced with the stacked bellman_operator!
@btime stacked_quantecon()

# Benchmark the version using fixedpoint() with m = 0 (simple iteration).
@btime fixedpoint_iteration()

# Benchmark the version using fixedpoint() and m = 1.
@btime fixedpoint_anderson_m(1)

# Benchmark the version using fixedpoint() and m = 2.
@btime fixedpoint_anderson_m(2)

# Benchmark the version using fixedpoint() and m = 7.
@btime fixedpoint_anderson_m(7)

# Benchmark the newton version.
@btime fixedpoint_newton()

# Get solutions.
default_quantecon_results = default_quantecon()
stacked_quantecon_results = stacked_quantecon()
default_fixedpoint_results = fixedpoint_iteration()
fixedpoint_results_1 = fixedpoint_anderson_m(1)
fixedpoint_results_2 = fixedpoint_anderson_m(2)
fixedpoint_results_7 = fixedpoint_anderson_m(7)

# Test for accuracy.

# Differences between Anderson results.
@test default_quantecon_results == (stacked_quantecon_results[2:end], stacked_quantecon_results[1])
@show norm(default_fixedpoint_results.zero - fixedpoint_results_1.zero, Inf)
@show norm(fixedpoint_results_1.zero - fixedpoint_results_2.zero, Inf)
@show norm(fixedpoint_results_2.zero - fixedpoint_results_7.zero, Inf)

# Differences between Anderson results and QuantEcon.
@show norm(default_fixedpoint_results.zero - cat(1, default_quantecon_results[2], default_quantecon_results[1]), Inf)
@show norm(fixedpoint_results_1.zero - cat(1, default_quantecon_results[2], default_quantecon_results[1]), Inf)
@show norm(fixedpoint_results_2.zero - cat(1, default_quantecon_results[2], default_quantecon_results[1]), Inf)
@show norm(fixedpoint_results_7.zero - cat(1, default_quantecon_results[2], default_quantecon_results[1]), Inf)

# Test fixed-point-ness of each solution
    # Utility function.
    function u(c::Real, σ::Real)
        if c > 0
            return (c^(1 - σ) - 1) / (1 - σ)
        else
            return -10e6
        end
    end
    # Instantiate model object.
    mcm = McCallModel()
    # QuantEcon Test
        # Bellman function from QuantEcon
        function update_bellman!(mcm::McCallModel,
                                V::AbstractVector,
                                V_new::AbstractVector,
                                U::Real)
            # Simplify notation
            α, β, σ, c, γ = mcm.α, mcm.β, mcm.σ, mcm.c, mcm.γ

            for (w_idx, w) in enumerate(mcm.w_vec)
                # w_idx indexes the vector of possible wages
                V_new[w_idx] = u(w, σ) + β * ((1 - α) * V[w_idx] + α * U)
            end

            U_new = u(c, σ) + β * (1 - γ) * U +
                            β * γ * dot(max.(U, V), mcm.p_vec)

            return U_new
        end
        # Difference
        old = cat(1, default_quantecon_results[2], default_quantecon_results[1])
        V_new = similar(default_quantecon_results[1])
        U_new = update_bellman!(mcm, default_quantecon_results[1], V_new, default_quantecon_results[2]);
        @show norm(old - cat(1, U_new, V_new), Inf)
    # Other test
        # New Bellman Operator Function
        function bellman_operator!(newVec, oldVec, model = mcm)
            # Unpack parameters.
            α, β, σ, c, γ = mcm.α, mcm.β, mcm.σ, mcm.c, mcm.γ
            # Add new V(w) values to newVec.
            for (w_idx, w) in enumerate(mcm.w_vec)
                # w_idx indexes the vector of possible wages
                newVec[w_idx + 1] = u(w, σ) + β * ((1 - α) * oldVec[1 + w_idx] + α * oldVec[1])
            end
            # Add new U value to newVec.
            newVec[1] = u(c, σ) + β * (1 - γ) * oldVec[1] + β * γ * dot(max.(oldVec[1], oldVec[2:end]), mcm.p_vec)
        end
        newVec = similar(default_fixedpoint_results.zero)
        bellman_operator!(newVec, default_fixedpoint_results.zero)
        @show norm(default_fixedpoint_results.zero - newVec, Inf)
        bellman_operator!(newVec, fixedpoint_results_1.zero)
        @show norm(fixedpoint_results_1.zero - newVec, Inf)
        bellman_operator!(newVec, fixedpoint_results_2.zero)
        @show norm(fixedpoint_results_2.zero - newVec, Inf)
        bellman_operator!(newVec, fixedpoint_results_7.zero)
@show norm(fixedpoint_results_7.zero - newVec, Inf)


text:
using NASAPrecipitation
using Test

@testset "NASAPrecipitation.jl" begin
    # Write your tests here.
end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

function cholesky_L(m::AbstractSparseMatrix)
    sparse(cholesky(m).L)
end

function cholesky_L(m::AbstractMatrix)
    cholesky(m).L
end

function cholesky_L(m::Diagonal)
    Diagonal(sqrt.(m.diag))
end

struct PDLinMapWithChol{
    T,
    PMap<:LinearMaps.LinearMap{T},
    CMap<:LinearMaps.LinearMap{T}
} <: LinearMaps.LinearMap{T}
    parent::PMap
    chol_L::CMap
end

function PDLinMapWithChol(A::AbstractMatrix, chol_A::AbstractMatrix)
    chol_L = LinearMap(
        chol_A,
        issymmetric = false, ishermitian = false, isposdef = false
    )

    wrapped_A = LinearMap(
        A,
        issymmetric = true, ishermitian = true, isposdef = true
    )

    PDLinMapWithChol(wrapped_A, chol_L)
end

PDLinMapWithChol(A::AbstractMatrix) = PDLinMapWithChol(A, cholesky_L(A))

Base.parent(A::PDLinMapWithChol) = A.parent
cholesky_L(A::PDLinMapWithChol) = A.chol_L

Base.Matrix(A::PDLinMapWithChol) = Matrix(parent(A))

Base.convert(::Type{AbstractMatrix}, A::PDLinMapWithChol) =
    convert(AbstractMatrix, parent(A))
Base.convert(::Type{Matrix}, A::PDLinMapWithChol) =
    convert(Matrix, parent(A))
Base.convert(::Type{SparseMatrixCSC}, A::PDLinMapWithChol) =
    convert(SparseMatrixCSC, parent(A))
SparseArrays.sparse(A::PDLinMapWithChol) =
    sparse(parent(A))

Base.size(A::PDLinMapWithChol) = size(parent(A))

LinearAlgebra.issymmetric(A::PDLinMapWithChol) = true
LinearAlgebra.ishermitian(A::PDLinMapWithChol) = true
LinearAlgebra.isposdef(A::PDLinMapWithChol) = true

LinearAlgebra.adjoint(A::PDLinMapWithChol) = A

LinearAlgebra.transpose(A::PDLinMapWithChol{<:Real}) = A
LinearAlgebra.transpose(A::PDLinMapWithChol) = PDLinMapWithChol(transpose(A.parent), A.chol_L)

LinearMaps.MulStyle(A::PDLinMapWithChol) = LinearMaps.MulStyle(parent(A))

Base.@propagate_inbounds LinearAlgebra.mul!(
    y::AbstractVector, A::PDLinMapWithChol, x::AbstractVector, α::Number, β::Number
) = mul!(y, parent(A), x, α, β)

Base.@propagate_inbounds LinearAlgebra.mul!(
    y::AbstractMatrix, A::PDLinMapWithChol, x::AbstractMatrix, α::Number, β::Number
) = mul!(y, parent(A), x, α, β)

Base.:(*)(A::LinearMap, B::PDLinMapWithChol) = A * parent(B)
Base.:(*)(A::PDLinMapWithChol, B::LinearMap) = parent(A) * B
Base.:(*)(A::PDLinMapWithChol, B::PDLinMapWithChol) = parent(A) * parent(B)

function _blockdiag(As::Tuple{Vararg{<:PDLinMapWithChol}})
    PDLinMapWithChol(blockdiag(map(parent, As)...), blockdiag(map(A -> A.chol_L, As)...))
end

const DiagLinearMap{T} = LinearMaps.WrappedMap{T,<:Diagonal}
const DiagPDLinMapWithChol = PDLinMapWithChol{T,<:DiagLinearMap,<:DiagLinearMap} where T

function _blockdiag(As::Tuple{Vararg{<:DiagPDLinMapWithChol}}) where T
    PDLinMapWithChol(
        LinearMap(Diagonal(vcat(map(A -> A.parent.lmap.diag, As)...)), issymmetric = true, ishermitian = true, isposdef = true),
        LinearMap(Diagonal(vcat(map(A -> A.chol_L.lmap.diag, As)...)), issymmetric = true, ishermitian = true, isposdef = true)
    )
end

blockdiag(As::PDLinMapWithChol...) = _blockdiag(As)

function _blockdiag(A::AbstractVector{<:DiagPDLinMapWithChol})
    d = reduce(vcat, map(x -> parent(x).lmap.diag, A))
    PDLinMapWithChol(Diagonal(d))
end


text:
abstract type AbstractDecisionRule{S<:Grid,T<:Grid,nx} end

function Base.show(io::IO, dr::AbstractDecisionRule)
    println(io, typeof(dr))
end

outdim(dr::AbstractDecisionRule{<:Grid,<:Grid,nx}) where nx = nx

# ---------------------- #
# Constant decision rule #
# ---------------------- #

mutable struct ConstantDecisionRule{nx} <: AbstractDecisionRule{EmptyGrid,EmptyGrid,nx}
    constants::SVector{nx,Float64}
end

function ConstantDecisionRule(constants::Vector{Float64})
    nx = length(constants)
    ConstantDecisionRule{nx}(SVector{nx,Float64}(constants...))
end

function ConstantDecisionRule(model::Model)
    ConstantDecisionRule(model.calibration[:controls])
end

(dr::ConstantDecisionRule)(x::Point) = dr.constants
(dr::ConstantDecisionRule)(x::Vector{Point{d}}) where d = [dr.constants for n=1:length(x)]
(dr::ConstantDecisionRule)(i::Int, x::Union{Point{d},Vector{Point{d}}}) where d = dr(x)
(dr::ConstantDecisionRule)(i::Int, j::Int, x::Union{Point{d},Vector{Point{d}}}) where d = dr(x)

# (dr::ConstantDecisionRule)(x::AbstractVector, y::AbstractVector) = dr.constants
# (dr::ConstantDecisionRule)(x::AbstractVector, y::AbstractMatrix) = repeat(dr.constants', size(y, 1), 1)
# (dr::ConstantDecisionRule)(x::AbstractMatrix, y::AbstractVector) = repeat(dr.constants', size(x, 1), 1)
# (dr::ConstantDecisionRule)(x::AbstractMatrix, y::AbstractMatrix) = repeat(dr.constants', size(x, 1), 1)
# (dr::ConstantDecisionRule)(i::Int, x::Union{AbstractVector,AbstractMatrix}) = dr(x)
# (dr::ConstantDecisionRule)(i::Int, j::Int, x::Union{AbstractVector,AbstractMatrix}) = dr(x)

# ------------------------------ #
# 2-dimensional Taylor Expansion #
# ------------------------------ #

mutable struct BiTaylorExpansion{nx} <: AbstractDecisionRule{EmptyGrid,EmptyGrid,nx}
    m0::Vector{Float64}
    s0::Vector{Float64}
    x0::Vector{Float64}
    x_m::Matrix{Float64}
    x_s::Matrix{Float64}
end

(dr::BiTaylorExpansion)(m::Point, s::Point) = SVector( (dr.x0 + dr.x_m*(m-dr.m0) + dr.x_s*(s-dr.s0))... )

(dr::BiTaylorExpansion)(m::AbstractVector{<:Point}, s::AbstractVector{<:Point}) = [dr(m[i], s[i]) for i=1:length(s)]

(dr::BiTaylorExpansion)(m::AbstractVector, s::AbstractVector) = Array(dr(SVector(m...), SVector(s...)))
(dr::BiTaylorExpansion)(m::AbstractMatrix, s::AbstractVector) = vcat([(dr(m[i, :], s))' for i=1:size(m, 1) ]...)
(dr::BiTaylorExpansion)(m::AbstractVector, s::AbstractMatrix) = vcat([(dr(m, s[i, :]))' for i=1:size(s, 1) ]...)
(dr::BiTaylorExpansion)(m::AbstractMatrix, s::AbstractMatrix) = vcat([(dr(m[i, :], s[i, :]))' for i=1:size(m, 1) ]...)

# User defined functions
struct CFunDR{S,T,nx} <: AbstractDecisionRule{S,T,nx}
    fun::Function
    dprocess::AbstractDiscretizedProcess
end
CFunDR(fun::Function, dprocess::AbstractDiscretizedProcess, n_x::Int) = CFunDR{typeof(dprocess.grid), EmptyGrid, n_x}(fun, dprocess)
(cfdr::CFunDR)(i::Int, x::Point{d}) where d = cfdr.fun(node(Point, cfdr.dprocess,i),x)
(cfdr::CFunDR)(i::Int, x::Vector{Point{d}}) where d = [cfdr(i,e) for e in x]


text:
## Simulation of SMLM Data using Nmer Pattern

using Revise
using SMLMSim
using SMLMData
using PlotlyJS

# Simulation paramters use physical units
# smld structures are in units of pixels and frames

smld_true, smld_model, smld_noisy=SMLMSim.sim(;
ρ=1.0,
σ_PSF=.13, #micron
minphotons=50,
ndatasets=10,
nframes=1000,
framerate=50.0, # 1/s
pattern=SMLMSim.Nmer2D(),
molecule=SMLMSim.GenericFluor(;q=[0 50; 1e-2 0]), #1/s
camera=SMLMSim.IdealCamera(;xpixels=256,ypixels=256,pixelsize=0.1) #pixelsize is microns
)

plt=PlotlyJS.plot(scattergl(x=smld_noisy.x, y=smld_noisy.y, mode="markers"))
display(plt)


text:
### Compute mean recursively given previous mean and new sample

recursive_mean(lastmean::Real, k::Integer, x::Real) = ((k-1)*lastmean+x)/k

recursive_mean!(m::RealVector, lastmean::RealVector, k::Integer, x::RealVector) = (m[:] = ((k-1)*lastmean+x)/k)

mean(s::VariableNState{Univariate}) = mean(s.value)

mean(s::VariableNState{Multivariate}, i::Integer) = mean(s.value[i, :])

mean(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size) = eltype(s)[mean(s, i) for i in r]


text:
using BenchmarkTools, ConstrainedSystems

Ns, Nc = 1000, 100
A,B2,B1t,C,rhs1v,rhs2v,solex = ConstrainedSystems.basic_linalg_problem(Ns=Ns,Nc=Nc);

Δt = 1.0e-2
prob1, xexact, yexact = ConstrainedSystems.basic_constrained_problem(tmax=3Δt)
prob2, xexact, yexact = ConstrainedSystems.cartesian_pendulum_problem(tmax=3Δt)

BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

SUITE = BenchmarkGroup()

SUITE["saddle setup"] = BenchmarkGroup()
SUITE["saddle solve"] = BenchmarkGroup()

rhs = SaddleVector(rhs1v,rhs2v)
sol = similar(rhs)

SUITE["saddle setup"]["basic problem"] = @benchmarkable ConstrainedSystems.SaddleSystem(A,B2,B1t,C,rhs)

As = ConstrainedSystems.SaddleSystem(A,B2,B1t,C,rhs)

SUITE["saddle solve"]["basic problem"] = @benchmarkable sol .= As\rhs

SUITE["timemarch setup"] = BenchmarkGroup()
SUITE["timemarch basic problem"] = BenchmarkGroup()
SUITE["timemarch pendulum problem"] = BenchmarkGroup()

SUITE["timemarch setup"]["IFHEEuler"] = @benchmarkable ConstrainedSystems.init(prob2, IFHEEuler(),dt=Δt)
SUITE["timemarch setup"]["LiskaIFHERK"] = @benchmarkable ConstrainedSystems.init(prob2, LiskaIFHERK(),dt=Δt)

integrator11 = ConstrainedSystems.init(prob1, IFHEEuler(),dt=Δt)
integrator12 = ConstrainedSystems.init(prob1, LiskaIFHERK(),dt=Δt)
integrator21 = ConstrainedSystems.init(prob2, IFHEEuler(),dt=Δt)
integrator22 = ConstrainedSystems.init(prob2, LiskaIFHERK(),dt=Δt)

#SUITE["basic problem"]["LiskaIFHERK"] = @benchmarkable solve(prob1, LiskaIFHERK(),dt=Δt)
#SUITE["basic problem"]["IFHEEuler"] = @benchmarkable solve(prob1, IFHEEuler(),dt=Δt)
#SUITE["pendulum problem"]["LiskaIFHERK"] = @benchmarkable solve(prob2, LiskaIFHERK(),dt=Δt)
#SUITE["pendulum problem"]["IFHEEuler"] = @benchmarkable solve(prob2, IFHEEuler(),dt=Δt)

SUITE["timemarch basic problem"]["IFHEEuler"] = @benchmarkable step!(integrator11,Δt)
SUITE["timemarch basic problem"]["LiskaIFHERK"] = @benchmarkable step!(integrator12,Δt)
SUITE["timemarch pendulum problem"]["IFHEEuler"] = @benchmarkable step!(integrator21,Δt)
SUITE["timemarch pendulum problem"]["LiskaIFHERK"] = @benchmarkable step!(integrator22,Δt)

run(SUITE)


text:
using Flux
using WindMixing
using JLD2
using FileIO
using OceanParameterizations
using OrdinaryDiffEq
using Random
using GalacticOptim
using LinearAlgebra
using WindMixing: solve_NDE_nonmutating
using WindMixing: solve_NDE_nonmutating_backprop
using WindMixing: calculate_profile_gradient
using BenchmarkTools

train_files = ["-1e-3"]
# train_files = [
#                "wind_-1e-3_cooling_4e-8",
#             #    "wind_-2e-4_cooling_1e-8",
#             #    "wind_-1e-3_cooling_2e-8",
#                "wind_-2e-4_cooling_5e-8",
#                "wind_-5e-4_cooling_3e-8"
#                ]

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

PATH = pwd()
FILE_PATH_NN = joinpath(PATH, "extracted_training_output",
                        "NDE_training_modified_pacanowski_philander_1sim_-1e-3_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_1e-2_rate_1e-4_extracted.jld2")

@assert isfile(FILE_PATH_NN)
file = jldopen(FILE_PATH_NN, "r")

uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]
close(file)

# N_inputs = 96
# hidden_units = 400
# N_outputs = 31

# # weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs)))
# weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs)))

# uw_NN = re(weights ./ 1f5)
# vw_NN = re(weights ./ 1f5)
# wT_NN = re(weights ./ 1f5)

tsteps = 1:1:10
timestepper = ROCK4()
n_simulations = length(train_files)
solve_NDE_nonmutating(uw_NN, vw_NN, wT_NN, 𝒟train, tsteps, timestepper;
                                n_simulations=n_simulations, ν₀=1f-4, ν₋=1f-1, ΔRi=1f0, Riᶜ=0.25, Pr=1f0, κ=10f0, f=1f-4, α=1.67f-4, g=9.81f0)

optimizer = ADAM(2e-4)
maxiters = 3
solve_NDE_nonmutating_backprop(uw_NN, vw_NN, wT_NN, 𝒟train, tsteps, timestepper, optimizer,
                                maxiters=maxiters, n_simulations=n_simulations, gradient_scaling=1f-2, ν₀=1f-4, ν₋=1f-1, ΔRi=1f0, Riᶜ=0.25, Pr=1f0, κ=10f0, f=1f-4, α=1.67f-4, g=9.81f0)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    GenericDensity{F<:Function} <: AbstractDensity

*BAT-internal, not part of stable public API.*

Constructors:

    GenericDensity(log_f)

Turns the logarithmic density function `log_f` into a BAT-compatible
`AbstractDensity`. `log_f` must support

    `log_f(params::Any)::Real`

It must be safe to execute `log_f` in parallel on multiple threads and
processes.
"""
struct GenericDensity{F<:Function} <: AbstractDensity
    log_f::F
end

Base.convert(::Type{GenericDensity}, log_f::Function) = GenericDensity(log_f)
Base.convert(::Type{AbstractDensity}, log_f::Function) = GenericDensity(log_f)

Base.parent(density::GenericDensity) = density.log_f

function density_logval(
    density::GenericDensity,
    params::Any
)
    density.log_f(params)
end


text:
module forceintensitytest1
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   vector = [10.0]
   fi = ForceIntensity(vector)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest1
forceintensitytest1.test()

module forceintensitytest2
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      v .= [10.0]
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest2
forceintensitytest2.test()

module forceintensitytest3
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
       return (time < 5.0 ?  v .= [10.0] : v .= [0.0])
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!, 0.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
   settime!(fi, 6.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [0.0]
end
end
using .forceintensitytest3
forceintensitytest3.test()

module forceintensitytest4
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
      v .= [10.0]
      return v
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
   settime!(fi, 6.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest4
forceintensitytest4.test()


text:
using PyCall
using Flux
using BetaML.Nn
using BenchmarkTools
using Distributions
import LinearAlgebra.BLAS

BLAS.set_num_threads(1)

torch = pyimport("torch")
torch.set_num_threads(1)

NN = torch.nn.Sequential(
    torch.nn.Linear(8, 64),
    torch.nn.Tanh(),
    torch.nn.Linear(64, 32),
    torch.nn.Tanh(),
    torch.nn.Linear(32, 2),
    torch.nn.Tanh()
)

torch_nn(in) = NN(in)

Flux_nn = Chain(Dense(8,64,tanh),
                Dense(64,32,tanh),
                Dense(32,2,tanh))

BetaML_nn = buildNetwork([
               DenseLayer(8,64,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(8+64),sqrt(6)/sqrt(8+64)),64,8),wb=rand(Uniform(-sqrt(6)/sqrt(8+64),sqrt(6)/sqrt(8+64)),64)),
               DenseLayer(64,32,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(64+32),sqrt(6)/sqrt(64+32)),32,64),wb=rand(Uniform(-sqrt(6)/sqrt(64+32),sqrt(6)/sqrt(64+32)),32)),
               DenseLayer(32,2,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(32+2),sqrt(6)/sqrt(32+2)),2,32),wb=rand(Uniform(-sqrt(6)/sqrt(32+2),sqrt(6)/sqrt(32+2)),2))],
               squaredCost,name="Bike sharing regression model")

BetaML_nn = buildNetwork([
              DenseLayer(8,64,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(8+64),sqrt(6)/sqrt(8+64)),64,8),wb=zeros(64)),
              DenseLayer(64,32,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(64+32),sqrt(6)/sqrt(64+32)),32,64),wb=zeros(32)),
              DenseLayer(32,2,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(32+2),sqrt(6)/sqrt(32+2)),2,32),wb=zeros(2))],
              squaredCost,name="Bike sharing regression model")

BetaML_nn2 = buildNetwork([
            DenseLayer(8,64,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(8+64),sqrt(6)/sqrt(8+64)),64,8),wb=zeros(64), df=dtanh),
            DenseLayer(64,32,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(64+32),sqrt(6)/sqrt(64+32)),32,64),wb=zeros(32),df=dtanh),
            DenseLayer(32,2,f=tanh,w=rand(Uniform(-sqrt(6)/sqrt(32+2),sqrt(6)/sqrt(32+2)),2,32),wb=zeros(2),df=dtanh)],
            squaredCost,name="Bike sharing regression model",dcf=dSquaredCost)

for i in [1, 10, 100, 1000]
    println("Batch size: $i")
    torch_in = torch.rand(i,8)
    flux_in = rand(Float32,8,i)
    betaml_in = rand(Float32,i,8)
    print("pytorch:")
    @btime torch_nn($torch_in)
    print("flux   :")
    @btime Flux_nn($flux_in)
    print("betaml   :")
    @btime predict($BetaML_nn,$betaml_in)
    print("betaml2   :")
    @btime predict($BetaML_nn2,$betaml_in)
end

#= Output:
Batch size: 1
pytorch:  89.920 μs (6 allocations: 192 bytes)
flux   :  3.426 μs (6 allocations: 1.25 KiB)
betaml   :  3.046 μs (19 allocations: 3.55 KiB)
betaml2   :  3.089 μs (19 allocations: 3.55 KiB)
Batch size: 10
pytorch:  100.737 μs (6 allocations: 192 bytes)
flux   :  19.743 μs (6 allocations: 8.22 KiB)
betaml   :  33.137 μs (181 allocations: 34.77 KiB)
betaml2   :  32.259 μs (181 allocations: 34.77 KiB)
Batch size: 100
pytorch:  132.689 μs (6 allocations: 192 bytes)
flux   :  184.807 μs (8 allocations: 77.16 KiB)
betaml   :  306.326 μs (1801 allocations: 347.08 KiB)
betaml2   :  310.554 μs (1801 allocations: 347.08 KiB)
Batch size: 1000
pytorch:  392.295 μs (6 allocations: 192 bytes)
flux   :  1.838 ms (10 allocations: 766.19 KiB)
betaml   :  3.172 ms (18490 allocations: 3.40 MiB)
betaml2   :  3.116 ms (18490 allocations: 3.40 MiB)
=#


text:
@testset "Interval" begin
    @test 0..1 isa Interval
    interval = 0..1

    @test 0 ∈ interval
    @test 1 ∈ interval
    @test 2 ∉ interval
    @test -1 ∉ interval
    @test 0.5 ∈ interval
    @test 1.1 ∉ interval
    @test -0.1 ∉ interval

    @test 0 ∈ 0 .. 0
    @test 0 ∉ 0 .. -1

    @test interval == interval
    @test interval != -1..1
    @test interval != 0..2
    @test interval != 0..0
    @test interval != 1..1
    @test interval == 0.0 .. 1.0
    @test 0 .. 0 == 0 .. 0
    @test 0 .. 0 != 1 .. 0
    @test 0 .. -1 == 1 .. 0

    @test isequal(interval, interval)
    @test !isequal(interval, -1..1)
    @test !isequal(interval, 0..2)
    @test !isequal(interval, 0..0)
    @test !isequal(interval, 1..1)
    @test isequal(interval, 0.0 .. 1.0)
    @test isequal(0 .. 0, 0 .. 0)
    @test !isequal(0 .. 0, 1 .. 0)
    @test isequal(0 .. -1, 1 .. 0)

    @test isless(0..1, 2..3)
    @test !isless(2..3, 0..1)
    @test isless(0..1, 0..2)
    @test !isless(0..2, 0..1)

    @test !isless(0..0, 0..0)
    @test isless(0..0, 1..1)
    @test isless(0..0, 1..0)
    @test isless(0..0, 1 .. -1)
    @test isless(0..0, 0..1)
    @test isless(0..0, 0 .. -1)
    @test !isless(1..1, 0..0)
    @test !isless(1..0, 0..0)
    @test !isless(1 .. -1, 0..0)
    @test !isless(0..1, 0..0)
    @test !(isless(0 .. -1, 0..0))
end

@testset "lessthan and greaterthan" begin
    interval = greaterthan(0) .. lessthan(1)

    @test -1 ∉ interval
    @test 0 ∉ interval
    @test 0.5 ∈ interval
    @test 1 ∉ interval
    @test 2 ∉ interval
end


text:
function check_module_node!(cur_package::Module, cur_info::FileInfo, cur_shard::Expr)
  cur_shard = _clean_module_shard!(cur_package, deepcopy(cur_shard))

  cur_undef = _check_for_bad_imports(cur_package, cur_info, cur_shard)
  println(123)
  println(cur_undef)

  if cur_undef == nothing
    cur_undef = _check_for_bad_types(cur_package, cur_info, cur_shard)
  end
  println(1234)
  println(cur_undef)

  if cur_undef == nothing
    cur_undef = _check_for_bad_methods(cur_package, cur_info, cur_shard)
  end
  println(1235)
  println(cur_undef)

  if cur_undef == nothing
    has_unloaded_file = false

    cur_folder = dirname(cur_info.name)

    for cur_sub_shard in cur_shard.args[3].args
      ( cur_sub_shard.head == :call ) || continue
      ( first(cur_sub_shard.args) == :include ) || continue

      full_include_path = joinpath(cur_folder, cur_sub_shard.args[2])

      if in(full_include_path, map(tmp_info -> tmp_info.name, cur_info.parent.file_infos))
        has_unloaded_file = true
        break
      end

      cur_sub_shard.args[2] = full_include_path
    end

    has_unloaded_file || return true
  end

  # put modules in the middle of the pack

  if iszero(length(cur_info.parent.file_infos))
    cur_info.time = 0.0
  else
    cur_info.time = median(map(tmp_info -> tmp_info.time, cur_info.parent.file_infos))
  end

  cur_info.undef = cur_undef

  return false

end

function _check_for_bad_methods(cur_package::Module, cur_info::FileInfo, cur_shard::Expr)
  method_list = filter(
    cur_arg -> isdefined(cur_arg, :head) && (
      cur_arg.head == :function ||
      ( cur_arg.head == :(=) && isdefined(cur_arg.args[1], :head) && cur_arg.args[1].head == :call )
    ),
    cur_shard.args[3].args
  )

  method_list = map(cur_method -> first(cur_method.args), method_list)

  isempty(method_list) && return

  base_excused_vars = Symbol[]

  type_list = filter(
    cur_arg -> isdefined(cur_arg, :head) && ( cur_arg.head == :type || cur_arg.head == :abstract ),
    cur_shard.args[3].args
  )

  for cur_type in type_list

    cur_type_symbol = nothing

    ( cur_type.head == :abstract ) &&
      ( cur_type_symbol = cur_type.args[1] )

    ( cur_type.head == :type ) &&
      ( cur_type_symbol = cur_type.args[2] )

    ( cur_type_symbol == nothing ) &&
      error("Unable to determine what to do with cur_type: $cur_type")

    if isa(cur_type_symbol, Expr) && cur_type_symbol.head == :curly
      println("rookie")
      main_shard, cur_type_symbol = cur_type_symbol.args
      push!(base_excused_vars, main_shard)
    end

    if isa(cur_type_symbol, Expr) && cur_type_symbol.head == :(<:)
      println("sup")
      cur_type_symbol = first(cur_type_symbol.args)
    end

    isa(cur_type_symbol, Symbol) ||
      error("Unable to determine what to do with cur_type: $cur_type_symbol")

    push!(base_excused_vars, cur_type_symbol)
  end

  for cur_method in method_list
    ( cur_method.head == :call || cur_method.head == :where ) ||
      error("Unable to determine what to do with non-call function shard: $(cur_method.head)")

    cur_excused_vars = deepcopy(base_excused_vars)

    if cur_method.head == :where
      ( length(cur_method.args) == 2 ) ||
        error("Currently invalid type shard in module load [3/4]")

      cur_method, cur_type_definition = cur_method.args

      ( cur_type_definition.head == :(<:) ) ||
        error("Invalid where clause in module load")

      ( length(cur_type_definition.args) == 2 ) ||
        error("Currently invalid type shard in module load [4/4]")

      cur_param_type, cur_type_association = cur_type_definition.args

      deep_isdefined(cur_package, cur_info, cur_type_association, cur_excused_vars) ||
        return cur_type_association

      push!(cur_excused_vars, cur_param_type)
    end

    cur_call = first(cur_method.args)
    cur_args = cur_method.args[2:end]

    if isa(cur_call, Expr) && cur_call.head != :(.)
      ( cur_call.head == :curly ) ||
        error("Unable to check type with non-curly head type of: $(cur_call.head) [3/4]")

      new_excused_vars = cur_call.args[2:end]

      for (cur_index, cur_sub_shard) in enumerate(new_excused_vars)
        isa(cur_sub_shard, Expr) || continue

        cur_param_type, cur_type_association = cur_sub_shard.args

        deep_isdefined(cur_package, cur_info, cur_type_association, cur_excused_vars) ||
          return cur_type_association

        new_excused_vars[cur_index] = cur_param_type
      end

      append!(cur_excused_vars, new_excused_vars)
    end

    for cur_arg in cur_args
      isa(cur_arg, Expr) || continue
      cur_arg.head == :(::) || continue

      ( length(cur_arg.args) < 3 ) || error("Invalid :(::) args size in module call")

      cur_type_association = last(cur_arg.args)

      deep_isdefined(cur_package, cur_info, cur_type_association, cur_excused_vars) ||
        return cur_type_association
    end
  end

  return
end

function _get_import_list(cur_shard::Expr)
  import_list = filter(
    cur_arg -> isdefined(cur_arg, :head) && ( cur_arg.head == :import || cur_arg.head == :using ),
    cur_shard.args[3].args
  )

  cur_import_exprs = filter(
    cur_arg -> isdefined(cur_arg, :head) && cur_arg.head == :toplevel,
    cur_shard.args[3].args
  )

  for cur_expr in cur_import_exprs
    for cur_sub_expr in cur_expr.args
      ( cur_sub_expr.head == :import || cur_sub_expr.head == :using ) &&
        push!(import_list, cur_sub_expr)
    end
  end

  import_list
end

function _check_for_bad_imports(cur_package::Module, cur_info::FileInfo, cur_shard::Expr)
  import_list = _get_import_list(cur_shard)

  for cur_import in import_list

    dot_count = 0
    for cur_arg in cur_import.args
      ( cur_arg == :(.) ) || break
      dot_count += 1
    end

    if !iszero(dot_count)
      ( dot_count == 2 ) ||
        error("Invalid module import")

      shift!(cur_import.args)
      cur_import.args[1] = Symbol(cur_package)
    end

    if cur_import.args[1] != Symbol(cur_package)
      isdefined(cur_package, last(cur_import.args)) && continue

      cur_info.parent.boot_module.include_string(string(cur_import))

      continue
    end

    work_method = cur_import.args[end]

    work_module = nothing
    bad_expression = nothing

    tmp_modules = [ cur_package , cur_info.parent.boot_module , cur_package ]

    for tmp_module in tmp_modules
      work_module = tmp_module

      for cur_sub_sub_expr in cur_import.args[2+(dot_count-1):end-1]
        if !isdefined(work_module, cur_sub_sub_expr)
          bad_expression = cur_sub_sub_expr
          work_module = nothing
          break
        end

        work_module = getfield(work_module, cur_sub_sub_expr)
      end

      ( work_module == nothing ) && continue

      isdefined(work_module, work_method) && break
    end

    ( work_module == nothing ) && return bad_expression

    isdefined(work_module, work_method) && continue

    return work_method
  end

  return
end

function _check_for_bad_types(cur_package::Module, cur_info::FileInfo, cur_shard::Expr)
  type_list = filter(
    cur_arg -> isdefined(cur_arg, :head) && ( cur_arg.head == :type || cur_arg.head == :abstract ),
    cur_shard.args[3].args
  )

  excused_vars = Symbol[]

  for cur_type in type_list
    println("asdf ", cur_type)
    cur_type_definition = nothing

    if cur_type.head == :abstract
      cur_type_definition = cur_type.args[1]
      push!(excused_vars, cur_type_definition)
      continue
    end

    cur_type_definition = cur_type.args[2]

    if isa(cur_type_definition, Symbol)
      push!(excused_vars, cur_type_definition)
    end

    if isa(cur_type_definition, Expr)
      if cur_type_definition.head == :(<:)
        ( length(cur_type_definition.args) == 2 ) ||
          error("Currently invalid type shard in module load [1/4]")

        cur_param_type, cur_type_association = cur_type_definition.args

        if isa(cur_param_type, Symbol)
          push!(excused_vars, cur_param_type)
        elseif cur_param_type.head == :curly
          cur_type_definition = cur_param_type
        else
          error("Unable to determine what to do with param type of: $(cur_param_type.head)")
        end

        deep_isdefined(cur_package, cur_info, cur_type_association, excused_vars) ||
          ( println("axd") ; return cur_type_association )
      end

      if cur_type_definition.head != :(<:)
        ( cur_type_definition.head == :curly ) ||
          error("Unable to check type with non-curly head type of: $(cur_type_definition.head) [2/4]")

        for cur_sub_shard in cur_type_definition.args[2:end]
          if isa(cur_sub_shard, Symbol)
            push!(excused_vars, cur_sub_shard)
            continue
          end

          ( cur_sub_shard.head == :(<:) ) ||
            error("Unable to check type with a sub-head type of: $(cur_type_definition.head)")

          ( length(cur_sub_shard.args) == 2 ) ||
            error("Currently invalid type shard in module load [2/4]")

          cur_param_type, cur_type_association = cur_sub_shard.args

          push!(excused_vars, cur_param_type)

          deep_isdefined(cur_package, cur_info, cur_type_association, excused_vars) ||
            ( println("nimubs") ; return cur_type_association )
        end
      end
    end

    for cur_field in cur_type.args[3].args
      isa(cur_field, Symbol) && continue

      ( cur_field.head == :function ) && continue
      ( cur_field.head == :line ) && continue
      ( cur_field.head == :(=) ) && continue

      ( cur_field.head == :(::) ) ||
        error("Invalid type definition for loading")

      ( length(cur_field.args) == 2 ) ||
        error("Invalid field definition for loading")

      cur_undef = _check_type_field(cur_package, cur_info, cur_field.args[2], excused_vars)

      ( cur_undef == nothing ) || return cur_undef
    end
  end

  return
end

function _check_type_field(cur_package::Module, cur_info::FileInfo, cur_shard::Expr, excused_vars::Array{Symbol})
  ( cur_shard.head == :(.) ) &&
    return _check_dot_field(cur_package, cur_info, cur_shard)

  for cur_sub_shard in cur_shard.args
    undef_type = _check_type_field(cur_package, cur_info, cur_sub_shard, excused_vars)
    ( undef_type == nothing ) || return undef_type
  end

  return
end

function _check_type_field(cur_package::Module, cur_info::FileInfo, cur_shard::Symbol, excused_vars::Array{Symbol})
  in(cur_shard, excused_vars) && return

  deep_isdefined(cur_package, cur_info, cur_shard, excused_vars) && return

  cur_shard
end

function _check_dot_field(cur_package::Module, cur_info::FileInfo, cur_shard::Expr)
  cur_error = nothing

  try
    cur_package.eval(cur_shard)
  catch tmp_error
    cur_error = tmp_error
  end

  ( cur_error == nothing ) && return

  isa(cur_error, UndefVarError) || rethrow(cur_error)

  return cur_error.var
end

function deep_isdefined(cur_package::Module, cur_info::FileInfo, cur_shard::Expr, excused_vars::Vector{Symbol})
  if cur_shard.head == :(.)
    cur_dot_reference = _expand_dot_var(cur_package, cur_shard)

    tmp_packages = [ cur_package , cur_info.parent.boot_module ]

    cur_field = nothing

    first_reference = first(cur_dot_reference)

    for tmp_package in tmp_packages
      isdefined(tmp_package, first_reference) || continue
      cur_field = getfield(tmp_package, first_reference)
      break
    end

    ( cur_field == nothing ) && return false

    for cur_reference in cur_dot_reference[2:end]
      isdefined(cur_field, cur_reference) || return false
      cur_field = getfield(cur_field, cur_reference)
    end

    return true
  end

  ( cur_shard.head == :curly ) ||
    error("Unable to check type with non-curly head type of: $(cur_shard.head) [1/4]")

  for cur_sub_shard in cur_shard.args
    isa(cur_sub_shard, Symbol) || continue

    deep_isdefined(cur_package, cur_info, cur_sub_shard, excused_vars) ||
      return false
  end

  return true
end

function deep_isdefined(cur_package::Module, cur_info::FileInfo, cur_shard::Symbol, excused_vars::Vector{Symbol})
  tmp_packages = [ cur_package , cur_info.parent.boot_module ]

  in(cur_shard, excused_vars) && return true

  for tmp_package in tmp_packages
    isdefined(tmp_package, cur_shard) && return true
  end

  return false
end

_clean_module_shard!(cur_package::Module, cur_shard::Any) = cur_shard

function _clean_module_shard!(cur_package::Module, cur_shard::Expr)
  if cur_shard.head == :macrocall && string(cur_shard.args[1]) == "Core.@doc"
    cur_shard = cur_shard.args[3]
    isa(cur_shard, Symbol) && return cur_shard
  end

  expanded_shard = cur_package.macroexpand(cur_shard)

  if expanded_shard.head != :error
    cur_shard = expanded_shard
  end

  for (cur_index, cur_sub_shard) in enumerate(cur_shard.args)
    cur_shard.args[cur_index] = _clean_module_shard!(cur_package, cur_sub_shard)
  end

  cur_shard
end

function _expand_dot_var(cur_package::Module, cur_shard::Expr, cur_list::Vector{Symbol}=Symbol[])
  ( length(cur_shard.args) == 2 ) ||
    error("Currently invalid dot shard in module load")

  first_shard, second_shard = cur_shard.args

  if isa(first_shard, Symbol)
    push!(cur_list, first_shard)
  elseif isa(first_shard, Expr)
    _expand_dot_var(cur_package, first_shard, cur_list)
  else
    error("Unable to determine what to do with dot shard: $cur_shard [1/2]")
  end

  if isa(second_shard, Symbol)
    push!(cur_list, second_shard)
  elseif isa(second_shard, QuoteNode)
    push!(cur_list, cur_package.eval(second_shard))
  else
    error("Unable to determine what to do with dot shard: $cur_shard [2/2]")
  end

  cur_list
end


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpReduce.jl")
TabLevel=TabLevel*"    "

##### Reducing expressions

function OpReduce(MajSp_in::MajoranaSp;verbose=false)
    verbose && println(" ........................  ")
    verbose && println("Reduce MajoranaSp")
    verbose && println("MajSp_in=")
    verbose && OpPrintln(MajSp_in)
    ###loop over all the terms and then reduce them
    MajSp=OpCopy(MajSp_in)
    for indx in 1:length(MajSp.SpList)
        #println("indx=",indx)
        #println("MajSp.SpList[$indx]=")
        #OpPrintln(MajSp.SpList[indx])
        MajSp.SpList[indx]=OpReduce(MajSp.SpList[indx])
        #println("MajSp.SpList[$indx]=")
        #OpPrintln(MajSp.SpList[indx])
    end
    verbose && println("Reduced the terms")
    ###Loop over the terms again and check that they are correctly ordered
    NumElems=length(MajSp.SpList)
    CheckStep=NumElems-1
    while CheckStep>0
        verbose && println("CheckStep=",CheckStep," of ",NumElems)
        if CheckStep==NumElems
            verbose && println("moved to faar to right")
            ###Move to the left and continue
            CheckStep-=1
            continue
        end
        ElemA=MajSp.SpList[CheckStep]
        ElemB=MajSp.SpList[CheckStep+1]
        verbose && println("ElemA:",ElemA," ElemB:",ElemB)
        if OpProp(ElemA,ElemB) ###If the elements are proportional merge them
            deleteat!(MajSp.SpList,CheckStep)
            ###Remove on of the elements
            MajSp.SpList[CheckStep]=ElemA+ElemB
            NumElems-=1
            verbose && println("Proprtional - merge pair")
        elseif OpOrdered(ElemA,ElemB)
            ###Ordered move to the left and continute
            verbose && println("Correctly ordered")
            CheckStep-=1
            continue
        else ###Reverse order
            MajSp.SpList[CheckStep]=ElemB
            MajSp.SpList[CheckStep+1]=ElemA
            verbose && println("Reverse order")
            CheckStep+=1
            continue
        end
        CheckStep-=1
    end
    MajSp.Reduced=true
    #println("MajoranaSp-reduced")
    return MajSp
end

function OpReduce(MajOp_in::MajoranaOp;verbose=false)
    MajOp=OpCopy(MajOp_in)
    scale=MajOp.scale
    OpList=MajOp.OpList
    NumElems=length(OpList)
    verbose && print("----------------\nReduce the expresison:")
    verbose && OpPrintln(MajOp)
    CheckStep=NumElems
    while CheckStep>0
        verbose && println("CheckStep=",CheckStep," of ",NumElems)
        ElemA=OpList[CheckStep]
        if CheckStep==NumElems
            verbose && println("movement to the left not possible")
            if MajOp.KetState
                verbose && println("There is a |0> at the end. Can the operator $ElemA be transformed?")
                OpType=ElemA.Type
                if OpType == "by"
                    verbose && println("Transforming $ElemA as by")
                    scale*=-im
                    OpList[CheckStep]=MajElem("bx",ElemA.Index)
                elseif OpType == "bz"
                    verbose && println("Transforming $ElemA as bz")
                    scale*=-im
                    OpList[CheckStep]=MajElem("c",ElemA.Index)
                elseif OpType == "bx" || OpType == "c"
                    verbose && println("Not Transforming $ElemA")
                else
                    error("Unknown operators $Op")
                end
            else
            end
            ###Move to the left and continue
            CheckStep-=1
            continue
        end
        ElemB=OpList[CheckStep+1]
        verbose && println("ElemA:",ElemA," ElemB:",ElemB)
        if ElemA==ElemB ###If the elements are the same remove them
            deleteat!(OpList,CheckStep)
            deleteat!(OpList,CheckStep)
            NumElems-=2
            verbose && OpPrintln(MajoranaOp(scale,OpList))
            verbose && println("remove pair")
        elseif ElemA < ElemB
            ###Ordered move to the left and continute
            verbose && println("Correctly ordered")
            CheckStep-=1
            continue
        else ###Reverse order
            OpList[CheckStep]=ElemB
            OpList[CheckStep+1]=ElemA
            scale*=-1
            verbose && println("Reverse order")
            verbose && OpPrintln(MajoranaOp(scale,OpList))
            CheckStep+=1
            continue
        end
        CheckStep-=1
    end
    MajOp.scale=scale
    MajOp.OpList=OpList
    return MajOp
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpReduce.jl")


text:
include("testexamples.jl")
include("mathprogtest.jl")


text:
module TestNewick

using Phylo
using Test

@testset "newick and nexus parsing" begin
    @testset "For $TreeType" for TreeType in [NamedBinaryTree]
        @test length(nodeiter(parsenewick("((,),(,));"))) == 7
        @test ["where", "when it's good", "Not mine", "MyLeaf"] ⊆
            nodenameiter(parsenewick("""((MyLeaf,"when it's good"),
                                             ('Not mine',where));""",
                                     TreeType))
        tree =
            parsenewick("((MyLeaf[&Real=23,'Not real'={5,4}]:4.0,)Parent,(,));",
                        TreeType)
        branches = branchfilter(tree) do tree, branch
            return getnodename(tree, src(tree, branch)) == "Parent" &&
                getnodename(tree, dst(tree, branch)) == "MyLeaf"
        end
        @test length(branches) == 1
        @test getlength(tree, first(branches)) ≈ 4.0
        @test getnodedata(tree, "MyLeaf")["Real"] == 23
        @test 5 ∈ getnodedata(tree, "MyLeaf")["Not real"]
        @test 4 ∈ getnodedata(tree, "MyLeaf")["Not real"]
        @test_throws Exception parsenewick("((,),(,)));", TreeType)
        @test_throws Exception parsenewick("((,),(,))", TreeType)
        @test_throws Exception parsenewick("((,),(,);", TreeType)
        @test_throws Exception parsenewick("((MyLeaf:-4.0,)Parent,(,));",
                                           TreeType)
        tree = open(f -> parsenewick(f, TreeType), Phylo.path("H1N1.newick"))
        @test nleaves(tree) == 507
        @test ntrees(tree) == 1
        treex = open(f -> parsenexus(f, TreeType), Phylo.path("H1N1.trees"))
        @test nleaves(tree) == nleaves(treex) ==
            nleaves(treex["TREE1"]) == nleaves(treex["TREE2"])
        @test ntrees(treex) == 2
    end

    @testset "For $TreeType" for TreeType in
        [NamedTree, RootedTree, ManyRootTree]
        @test nnodes(parsenewick("((,),(,,));", TreeType)) == 8
        @test ["where", "when it's good", "Not mine", "MyLeaf", "next"] ⊆
            nodenameiter(parsenewick("""((MyLeaf,"when it's good",next),
                                             ('Not mine',where));""",
                                     TreeType))
        tree = parsenewick("((MyLeaf:4.0,)Parent,());", TreeType)
        branches = branchfilter(tree) do tree, branch
            return getnodename(tree, src(tree, branch)) == "Parent" &&
                getnodename(tree, dst(tree, branch)) == "MyLeaf"
        end
        @test length(branches) == 1
        @test getlength(tree, first(branches)) ≈ 4.0
        @test_throws Exception parsenewick("((,),(,)));", TreeType)
        @test_throws Exception parsenewick("((,),(,))", TreeType)
        @test_throws Exception parsenewick("((,),(,);", TreeType)
        tree = open(f -> parsenewick(f, TreeType), Phylo.path("H1N1.newick"))
        @test nleaves(tree) == 507
        @test ntrees(tree) == 1
        treex = open(f -> parsenexus(f, TreeType), Phylo.path("H1N1.trees"))
        @test nleaves(tree) == nleaves(treex) ==
            nleaves(treex["TREE1"]) == nleaves(treex["TREE2"])
        @test ntrees(treex) == 2
    end
end

end


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

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

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

baremodule API

using Base
#import Base: call

if VERSION >= v"0.6.0-pre"
using SpecialFunctions
end

eval(x) = Core.eval(API, x)

@noinline function setindex!{T}(A::DenseArray{T}, args...)
  Base.setindex!(A, args...)
end

@inline function setindex!(A, args...)
  Base.setindex!(A, args...)
end

@noinline function getindex{T}(A::DenseArray{T}, args...)
  Base.getindex(A, args...)
end

@inline function getindex{T}(A::Type{T}, X...)
  T[ X[i] for i=1:length(X) ]
end

@inline function getindex(A, args...)
  Base.getindex(A, args...)
end

@noinline function reshape{T}(A::DenseArray{T}, args)
  Base.reshape(A, args...)
end

# inline reshape to version with dims as tuple for easier handling in compiler
@inline function reshape{T}(A::DenseArray{T}, args...)
  reshape(A, args)
end

# Unary operators/functions
const unary_map_operators_normal = Symbol[
    :-, :+, :acos, :acosh, :angle, :asin, :asinh, :atan, :atanh, :cbrt,
    :cis, :cos, :cosh, :exp10, :exp2, :exp, :expm1, :lgamma,
    :log10, :log1p, :log2, :log, :sin, :sinh, :sqrt, :tan, :tanh,
    :abs, :isnan]

const special_functions = Symbol[:erf]

const unary_map_operators = vcat(unary_map_operators_normal, special_functions)

const reduce_operators = Symbol[:sum, :prod, :minimum, :maximum, :any, :all]

const unary_operators = vcat(unary_map_operators, reduce_operators, Symbol[:copy])

# Binary operators/functions
const comparison_map_operators = Symbol[ :.>, :.<, :.<=, :.>=, :.== ]

const binary_map_operators = vcat(comparison_map_operators, Symbol[ :*, :/,
    :-, :+, :.+, :.-, :.*, :./, :.\, :.%, :.<<, :.>>, :.^,
    :div, :mod, :rem, :&, :|, :$, :min, :max])

const binary_operators = binary_map_operators

const elem_ops = Symbol[ :.-, :.+, :.*, :./, :.>, :.<, :.<=, :.>=, :.==, :.\, :.%, :.<<, :.>>, :.^ ]
const elem_bare_ops = Symbol[ :-, :+, :*, :/, :>, :<, :<=, :>=, :(==), :\, :%, :<<, :>>, :^ ]
const elem_map = Dict{Symbol,Symbol}(zip(elem_ops,elem_bare_ops))

const rename_op_from = Symbol[ :-, :+, :*, :/, :.-, :.+, :.*, :./, :.>, :.<, :.<=, :.>=, :.==, :.\, :.%, :.<<, :.>>, :.^, :&, :|, :$ ]
const rename_op_to = Symbol[ :pa_api_sub, :pa_api_add, :pa_api_mul, :pa_api_div, :pa_api_elem_sub, :pa_api_elem_add, :pa_api_elem_mul, :pa_api_elem_div, :pa_api_elem_gt, :pa_api_elem_lt, :pa_api_elem_lte, :pa_api_elem_gte, :pa_api_elem_equal, :pa_api_elem_back, :pa_api_elem_mod, :pa_api_elem_lshift, :pa_api_elem_rshift, :pa_api_elem_carat, :pa_api_ampersand, :pa_api_pipe, :pa_api_dollar ]
const rename_forward = Dict{Symbol,Symbol}(zip(rename_op_from,rename_op_to))
const rename_back = Dict{Symbol,Symbol}(zip(rename_op_to,rename_op_from))

@inline sum(A::DenseArray{Bool}) = sum(pa_api_elem_mul(1,A))
@inline sum(A::DenseArray{Bool}, x::Int) = sum(pa_api_elem_mul(1,A), x)

function rename_if_needed(f)
    if haskey(rename_forward, f)
        return rename_forward[f]
    else
        return f
    end
end

function rename_back_if_needed(f)
    if haskey(rename_back, f)
        return rename_back[f]
    else
        return f
    end
end

# reduction across a dimension
for f in reduce_operators
    nf = rename_if_needed(f)
    @eval begin
        @noinline function ($nf){T<:Number}(A::DenseArray{T}, x::Int)
            (Base.$f)(A, x)
        end
    end
end

for f in unary_operators
    nf = rename_if_needed(f)
    @eval begin
        @noinline function ($nf){T<:Number}(A::DenseArray{T})
            (Base.$f)(A)
        end
    end
    if !(in(f, binary_operators))
        @eval begin
            @inline function ($nf)(A...)
                (Base.$f)(A...)
            end
        end
    end
end

for f in unary_map_operators_normal
    nf = rename_if_needed(f)
    @eval begin
        @noinline function ($nf){T<:Number}(A::T)
            (Base.$f)(A)
        end
    end
end

for f in special_functions
    nf = rename_if_needed(f)
    if VERSION >= v"0.6.0-pre"
        @eval begin
            @noinline function ($nf){T<:Number}(A::T)
                (SpecialFunctions.$f)(A)
            end
        end
    else
        @eval begin
            @noinline function ($nf){T<:Number}(A::T)
                (Base.$f)(A)
            end
        end
    end
end

for f in binary_operators
    nf = rename_if_needed(f)
    bare_op = f
    if haskey(elem_map, bare_op)
        bare_op = elem_map[bare_op]
    end

    @eval begin
        @noinline function ($nf){T1<:Number,T2<:Number}(A::T1, B::DenseArray{T2})
            (Base.$f)(A, B)
        end
    end

    if bare_op == f
		@eval begin
			@noinline function ($nf){T1<:Number,T2<:Number}(B::DenseArray{T1}, A::T2)
				(Base.$f)(B, A)
			end
		end
    else
		@eval begin
			@noinline function ($nf){T1<:Number,T2<:Number}(B::DenseArray{T1}, A::T2)
				broadcast($bare_op, B, A)
			end
		end
    end

    if f != :* && f != :/
        if bare_op == f
			@eval begin
				@noinline function ($nf){T1<:Number,T2<:Number}(A::DenseArray{T1}, B::DenseArray{T2})
					(Base.$f)(A, B)
				end
			end
        else
			@eval begin
				@noinline function ($nf){T1<:Number,T2<:Number}(A::DenseArray{T1}, B::DenseArray{T2})
					broadcast($bare_op, A, B)
				end
			end
        end
    end

    @eval begin
        @inline function ($nf)(args...)
            (Base.$f)(args...)
        end
    end
end

@noinline function cartesianmapreduce(body, ndims, reductions...)
  # reductions are ignored in sequential semantics
  for (redFunc, redVar) in reductions
      # redVar holds initial value that has to be neutral
      vs = deepcopy(redVar)
      @assert (vs == redFunc(redVar)) ("initial value " * string(vs) * " not idempotent with respect to the reduction function")
  end
  for I in CartesianRange(ndims)
    body(I.I...)
  end
end

if VERSION >= v"0.5.0-dev+5306"
@inline function to_tuple_type(t)
  if t <: Tuple
    t
  else
    Tuple{t}
  end
end

@inline function cartesianarray(body, L::Int)
   if 1<0
     x=body(L)
   end
   cartesianarray(body, to_tuple_type(typeof(x)), (L,))
end

@inline function cartesianarray(body, L::Int, M::Int)
   if 1<0
     x=body(L, M)
   end
   cartesianarray(body, to_tuple_type(typeof(x)), (L,M,))
end

@inline function cartesianarray(body, L::Int, M::Int, N::Int)
   if 1<0
     x=body(L, M, N)
   end
   cartesianarray(body, to_tuple_type(typeof(x)), (L,M,N))
end

end

function cartesianarray{T}(body, ::Type{Tuple{T}}, ndims)
  a = Array{T}(ndims...)
  for I in CartesianRange(ndims)
    a[I.I...] = body(I.I...)
  end
  a
end

function cartesianarray{T1,T2}(body, ::Type{Tuple{T1, T2}}, ndims)
  a = Array{T1}(ndims...)
  b = Array{T2}(ndims...)
  for I in CartesianRange(ndims)
    (u, v) = body(I.I...)
    a[I.I...] = u
    b[I.I...] = v
  end
  return a, b
end

function cartesianarray{T1,T2,T3}(body, ::Type{Tuple{T1, T2, T3}}, ndims)
  a = Array{T1}(ndims...)
  b = Array{T2}(ndims...)
  c = Array{T3}(ndims...)
  for I in CartesianRange(ndims)
    (u, v, w) = body(I.I...)
    a[I.I...] = u
    b[I.I...] = v
    c[I.I...] = w
  end
  return a, b, c
end

function parallel_for(loopvar, range, body)
  throw("Not Implemented")
end

@noinline function map{T<:Number}(f::Function, a::DenseArray{T})
    Base.map(f, a)
end

@noinline function map{T<:Number}(f::Function, a::DenseArray{T}, b...)
    Base.map(f, a, b...)
end

@noinline function map{T<:Array}(f::Function, a::Array{T})
    Base.map(f, a)
end

@noinline function map{T<:Array}(f::Function, a::Array{T}, b...)
    Base.map(f, a, b...)
end

@noinline function map(f, a::Range)
    Base.map(f, a)
end

@noinline function map(f, a::Range, b...)
    Base.map(f, a, b...)
end

@inline function map(f, a...)
    Base.map(f, a...)
end

@noinline function map!{T<:Number}(f::Function, a::DenseArray{T})
    Base.map!(f, a)
end

@noinline function map!{T<:Number}(f::Function, a::DenseArray{T}, b...)
    Base.map!(f, a, b...)
end

@inline function map!(f, a...)
    Base.map!(f, a...)
end

@noinline function reduce{T<:Number}(f::Function, v::T, a::DenseArray{T})
    Base.reduce(f, v, a)
end

@noinline function reduce{T<:Number}(f::Function, v::T, a::DenseArray{T}, b...)
    Base.reduce(f, v, a, b...)
end

@inline function reduce(f, a...)
    Base.reduce(f, a...)
end

@noinline function broadcast{T<:Number}(f::Function, a::DenseArray{T}, b...)
    Base.broadcast(f, a, b...)
end

@noinline function broadcast{T<:Number}(f::Function, a::T, b...)
    Base.broadcast(f, a, b...)
end

@inline function broadcast(f, a...)
    Base.broadcast(f, a...)
end

operators = Set(vcat(unary_operators, binary_operators,
    Symbol[:map, :map!, :reduce, :broadcast, :setindex!, :getindex, :reshape]))

for opr in operators
    @eval export $opr
end

include("api-stencil.jl")

@noinline function pointer(args...)
    Base.pointer(args...)
end

macro par(args...)
    na = length(args)
    @assert (na > 0) "Expect a for loop as argument to @par"
    redVars = Array{Symbol}(na-1)
    redOps = Array{Symbol}(na-1)
    loop = args[end]
    if !isa(loop,Expr) || !(loop.head === :for)
        error("malformed @par loop")
    end
    if !isa(loop.args[1], Expr)
        error("maltformed for loop")
    end
    for i = 1:na-1
        @assert (isa(args[i], Expr) && args[i].head == :call) "Expect @par reduction in the form of var(op), but got " * string(args[i])
        v = args[i].args[1]
        op = args[i].args[2]
        #println("got reduction variable ", v, " with function ", op)
        @assert (isa(v, Symbol)) "Expect reduction variable to be symbols, but got " * string(v)
        @assert (isa(op, Symbol)) "Expect reduction operator to be symbols, but got " * string(op)
        redVars[i] = v
        redOps[i] = op
    end
    if loop.args[1].head == :block
        loopheads = loop.args[1].args
    else
        loopheads = Any[loop.args[1]]
    end
    # Here is a no-op, real implementation is in Capture module, and will run as an OptFramework pass
    esc(loop)
end

import .Stencil.runStencil
export @par, cartesianmapreduce, cartesianarray, parallel_for, runStencil, pointer

include("api-lib.jl")
importall .Lib

function enableLib()
  global operators
  oprs = delete!(Set(names(Lib)), :Lib)
  union!(operators, oprs)
  for opr in oprs
    @eval export $opr
  end
end

function disableLib()
  global operators
  setdiff!(operators, delete!(Set(names(Lib)), :Lib))
end

enableLib()

include("api-capture.jl")

end


text:
export OrderedPiecewiseFunction

struct OrderedPiecewiseFunction{N, T<:Real, F<:Function} <: AbstractPiecewise{N}
    functions::Vector{<:Function}
    breakpoints::Vector{T}
    evalfunc::F
    expressions::Vector{String}
    function OrderedPiecewiseFunction{N}(functions::Vector{Function}, breakpoints::Vector{T}, evalfunc::F, expressions::Vector{String}) where {N, T<:Real, F<:Function}
        @assert length(functions) ==  N + 1
        @assert length(breakpoints) == N
        @assert issorted(breakpoints)
        new{N, T, F}(functions, breakpoints, evalfunc, expressions)
    end
end

function OrderedPiecewiseFunction{N}(functions::Vector{Function}, breakpoints::Vector{T}) where {N, T<:Real}
    evalfunc = x -> evalpiecewise(Val(N), functions, breakpoints, x)
    expressions = [quote $f end for f in functions] .|> string
    OrderedPiecewiseFunction{N}(functions, breakpoints, evalfunc, expressions)
end

(p::OrderedPiecewiseFunction{N, T, F})(x) where {N, T, F} = p.evalfunc(x)

function Base.show(io::IO, p::OrderedPiecewiseFunction{N, T, F}) where {N, T, F}
    println(io, polyname(p) * ": ")
    printpiecewise(io, p.expressions, p.breakpoints)
end


text:
using Documenter
using AnyMOD

makedocs(sitename="AnyMOD.jl",
    authors = "Leonard Goeke",
    pages = [
        "Introduction" => "index.md",
        "Model Elements" => Any["Model object" => "model_object.md",
                                "Sets and Mappings" => "sets.md",
                                "Parts" => "parts.md",
                                "Parameter" => Any["Parameter overview" => "parameter_overview.md","Parameter list" => "parameter_list.md"],
                                "Variables" => "variables.md",
                                "Constraints" => "constraints.md"],
        "Tools" => Any["Reporting" => Any["Error handling" => "error.md","Data files" => "data.md","Plots" => "plots.md",],
                        "Performance and stability" => "performance.md"],
        "Annex" => Any["Tips" => "tips.md",
                        "API" => "api.md",
                        "Related material" => "related.md"]
        ],
    )

deploydocs(
    repo = "github.com/leonardgoeke/AnyMOD.jl.git",
    devbranch = "master"
)

max_stars_repo_name: leonardgoeke/anyMOD

max_issues_repo_name: leonardgoeke/anyMOD

max_forks_repo_name: leonardgoeke/anyMOD