text:
module Entry

import Base: thispatch, nextpatch, nextminor, nextmajor, check_new_version
import Base.Pkg: Reqs, Read, Query, Resolve, PkgError
import Base.LibGit2
importall Base.LibGit2
using Base.Pkg.Types
import ..PkgDev
import ..PkgDev.GitHub
using ..PkgDev: getrepohttpurl

function pull_request(dir::AbstractString; commit::AbstractString="", url::AbstractString="", branch::AbstractString="")
    with(GitRepo, dir) do repo
        if isempty(commit)
            commit = string(LibGit2.head_oid(repo))
        else
            !LibGit2.iscommit(commit, repo) && throw(PkgError("Cannot find pull commit: $commit"))
        end
        if isempty(url)
            url = LibGit2.getconfig(repo, "remote.origin.url", "")
        end
        force_branch = !isempty(branch)
        if !force_branch
            branch = "pull-request/$(commit[1:8])"
        end

        m = match(LibGit2.GITHUB_REGEX, url)
        m === nothing && throw(PkgError("not a GitHub repo URL, can't make a pull request: $url"))
        owner, owner_repo = m.captures[2:3]
        user = GitHub.user()
        info("Forking $owner/$owner_repo to $user")
        response = GitHub.fork(owner,owner_repo)
        fork = response["clone_url"]
        info("Pushing changes as branch $branch")
        refspecs = ["HEAD:refs/heads/$branch"]  # workaround for $commit:refs/heads/$branch
        LibGit2.push(repo, remoteurl=fork, refspecs=refspecs, force=force_branch, payload=GitHub.credentials())
        pr_url = "$(response["html_url"])/compare/$branch"
        info("To create a pull-request, open:\n\n  $pr_url\n")
    end
end

function submit(pkg::AbstractString, commit::AbstractString="")
    urlpath = Pkg.dir("METADATA",pkg,"url")
    url = ispath(urlpath) ? readchomp(urlpath) : ""
    pull_request(PkgDev.dir(pkg), commit=commit, url=url)
end

function publish(branch::AbstractString, prbranch::AbstractString="")
    tags = Dict{ByteString,Vector{ASCIIString}}()
    metapath = Pkg.dir("METADATA")
    with(GitRepo, metapath) do repo
        LibGit2.branch(repo) == branch ||
            throw(PkgError("METADATA must be on $branch to publish changes"))
        LibGit2.fetch(repo)

        ahead_remote, ahead_local = LibGit2.revcount(repo, "origin/$branch", branch)
        rcount = min(ahead_remote, ahead_local)
        ahead_remote-rcount > 0 && throw(PkgError("METADATA is behind origin/$branch – run `Pkg.update()` before publishing"))
        ahead_local-rcount == 0 && throw(PkgError("There are no METADATA changes to publish"))

        # get changed files
        for path in LibGit2.diff_files(repo, "origin/$branch", LibGit2.Consts.HEAD_FILE)
            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(metapath,path))
            old = LibGit2.cat(repo, LibGit2.GitBlob, "origin/$branch:$path")
            old !== nothing && old != sha1 && throw(PkgError("$pkg v$ver SHA1 changed in METADATA – refusing to publish"))
            with(GitRepo, PkgDev.dir(pkg)) do pkg_repo
                tag_name = "v$ver"
                tag_commit = LibGit2.revparseid(pkg_repo, "$(tag_name)^{commit}")
                LibGit2.iszero(tag_commit) || string(tag_commit) == sha1 || return false
                haskey(tags,pkg) || (tags[pkg] = ASCIIString[])
                push!(tags[pkg], tag_name)
                return true
            end || throw(PkgError("$pkg v$ver is incorrectly tagged – $sha1 expected"))
        end
        isempty(tags) && info("No new package versions to publish")
        info("Validating METADATA")
        check_metadata(Set(keys(tags)))
    end

    for pkg in sort!(collect(keys(tags)))
        with(GitRepo, PkgDev.dir(pkg)) do pkg_repo
            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,", "))
                LibGit2.push(pkg_repo, remote="origin", force=true,
                             refspecs=["refs/tags/$tag:refs/tags/$tag" for tag in forced],
                             payload=GitHub.credentials())
            end
            if !isempty(unforced)
                info("Pushing $pkg permanent tags: ", join(unforced,", "))
                LibGit2.push(pkg_repo, remote="origin",
                             refspecs=["refs/tags/$tag:refs/tags/$tag" for tag in unforced],
                             payload=GitHub.credentials())
            end
        end
    end
    info("Submitting METADATA changes")
    pull_request(metapath, branch=prbranch)
end

function write_tag_metadata(repo::GitRepo, pkg::AbstractString, ver::VersionNumber, commit::AbstractString, force::Bool=false)
    pkgdir = PkgDev.dir(pkg)
    content = with(GitRepo, pkgdir) do pkg_repo
        LibGit2.cat(pkg_repo, LibGit2.GitBlob, "$commit:REQUIRE")
    end
    reqs = content !== nothing ? Reqs.read(split(content, '\n', keep=false)) : 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 ||
                throw(PkgError("$pkg v$ver is already registered as $current, bailing"))
        end
        open(io->println(io,commit), sha1file, "w")
        LibGit2.add!(repo, sha1file)
        reqsfile = joinpath(d,"requires")
        if isempty(reqs)
            ispath(reqsfile) && LibGit2.remove!(repo, reqsfile)
        else
            Reqs.write(reqsfile,reqs)
            LibGit2.add!(repo, reqsfile)
        end
    end
    return nothing
end

function register(pkg::AbstractString, url::AbstractString)
    pkgdir = PkgDev.dir(pkg)
    isempty(pkgdir) && throw(PkgError("$pkg does not exist"))
    metapath = Pkg.dir("METADATA")
    ispath(pkgdir,".git") || throw(PkgError("$pkg is not a git repo"))
    isfile(metapath,pkg,"url") && throw(PkgError("$pkg already registered"))
    LibGit2.transact(GitRepo(metapath)) do repo
        # Get versions from package repo
        versions = with(GitRepo, pkgdir) do pkg_repo
            tags = filter(t->startswith(t,"v"), LibGit2.tag_list(pkg_repo))
            filter!(tag->ismatch(Base.VERSION_REGEX,tag), tags)
            [
                convert(VersionNumber,tag) => string(LibGit2.revparseid(pkg_repo, "$tag^{commit}"))
                for tag in tags
            ]
        end
        # Register package url in METADATA
        cd(metapath) do
            info("Registering $pkg at $url")
            mkdir(pkg)
            path = joinpath(pkg,"url")
            open(io->println(io,url), path, "w")
            LibGit2.add!(repo, path)
        end
        # Register package version in METADATA
        vers = sort!(collect(keys(versions)))
        for ver in vers
            info("Tagging $pkg v$ver")
            write_tag_metadata(repo, pkg,ver,versions[ver])
        end
        # Commit changes in METADATA
        if LibGit2.isdirty(repo)
            info("Committing METADATA for $pkg")
            msg = "Register $pkg [$url]"
            if !isempty(versions)
                msg *= ": $(join(map(v->"v$v", vers),", "))"
            end
            LibGit2.commit(repo, msg)
        else
            info("No METADATA changes to commit")
        end
    end
    return
end

function register(pkg::AbstractString)
    pkgdir = PkgDev.dir(pkg)
    isempty(pkgdir) && throw(PkgError("$pkg does not exist"))
    url = ""
    try
        url = LibGit2.getconfig(pkgdir, "remote.origin.url", "")
    catch err
        throw(PkgError("$pkg: $err"))
    end
    !isempty(url) || throw(PkgError("$pkg: no URL configured"))
    register(pkg, Pkg.Cache.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::AbstractString, ver::Union{Symbol,VersionNumber}, force::Bool=false, commitish::AbstractString="HEAD")
    pkgdir = PkgDev.dir(pkg)
    ispath(pkgdir,".git") || throw(PkgError("$pkg is not a git repo"))
    metapath = Pkg.dir("METADATA")
    with(GitRepo,metapath) do repo
        LibGit2.isdirty(repo, pkg) && throw(PkgError("METADATA/$pkg is dirty – commit or stash changes to tag"))
    end
    with(GitRepo, pkgdir) do repo
        LibGit2.isdirty(repo) && throw(PkgError("$pkg is dirty – commit or stash changes to tag"))
        commit = string(LibGit2.revparseid(repo, commitish))
        urlfile = joinpath(metapath,pkg,"url")
        registered = isfile(urlfile)

        if !force
            if registered
                avail = Pkg.cd(Pkg.Read.available, pkg)
                existing = VersionNumber[keys(avail)...]
                ancestors = filter(v->LibGit2.is_ancestor_of(avail[v].sha1, commit, repo), existing)
            else
                tags = filter(t->startswith(t,"v"), Pkg.LibGit2.tag_list(repo))
                filter!(tag->ismatch(Base.VERSION_REGEX,tag), tags)
                existing = VersionNumber[tags...]
                filter!(tags) do tag
                    sha1 = string(LibGit2.revparseid(repo, "$tag^{commit}"))
                    LibGit2.is_ancestor_of(sha1, commit, repo)
                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) :
                                        throw(PkgError("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")
        LibGit2.tag_create(repo, "v$ver", commit,
                           msg=(!isrewritable(ver) ? "$pkg v$ver [$(commit[1:10])]" : ""),
                           force=(force || isrewritable(ver)) )
        registered || return
        try
            LibGit2.transact(GitRepo(metapath)) do repo
                write_tag_metadata(repo, pkg, ver, commit, force)
                if LibGit2.isdirty(repo)
                    info("Committing METADATA for $pkg")
                    # Convert repo url into proper http url
                    repourl = getrepohttpurl(readchomp(urlfile))
                    LibGit2.commit(repo, "Tag $pkg v$ver [$repourl]")
                else
                    info("No METADATA changes to commit")
                end
            end
        catch
            LibGit2.tag_delete(repo, "v$ver")
            rethrow()
        end
    end
    return
end

function check_metadata(pkgs::Set{ByteString} = Set{ByteString}())
    avail = Pkg.cd(Read.available)
    deps, conflicts = Query.dependencies(avail)

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

    problematic = Resolve.sanity_check(deps, pkgs)
    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
        throw(PkgError(msg))
    end
    return
end

end


text:
using Unitful
using UnitfulEquivalences
using Test

struct Equiv1 <: Equivalence end
UnitfulEquivalences.edconvert(::typeof(dimension(u"m")), x::Unitful.Time,   ::Equiv1) = x * 10u"m/s"
UnitfulEquivalences.edconvert(::typeof(dimension(u"s")), x::Unitful.Length, ::Equiv1) = x * (1//10)u"s/m"
UnitfulEquivalences.edconvert(::typeof(dimension(u"g")), x::Unitful.Temperature, ::Equiv1) = x * 2.0u"g/K"
UnitfulEquivalences.edconvert(::typeof(dimension(u"K")), x::Unitful.Mass,        ::Equiv1) = x * 0.5u"K/g"
UnitfulEquivalences.edconvert(::typeof(NoDims),          x::Unitful.Energy,        ::Equiv1) = x / u"eV"
UnitfulEquivalences.edconvert(::typeof(dimension(u"J")), x::DimensionlessQuantity, ::Equiv1) = x * u"eV"

struct Equiv2 <: Equivalence end
UnitfulEquivalences.edconvert(::typeof(dimension(u"m")), x::Unitful.Time,   ::Equiv2) = ustrip(u"s", x)^3 * u"m"
UnitfulEquivalences.edconvert(::typeof(dimension(u"s")), x::Unitful.Length, ::Equiv2) = cbrt(ustrip(u"m", x)) * u"s"

struct Equiv3 <: Equivalence; val::Int; end
UnitfulEquivalences.edconvert(::typeof(dimension(u"J")),   x::Unitful.Area,   e::Equiv3) = e.val*u"J" / ustrip(u"m^2", x)
UnitfulEquivalences.edconvert(::typeof(dimension(u"m^2")), x::Unitful.Energy, e::Equiv3) = e.val*u"m^2" / ustrip(u"J", x)

struct NoEquiv end
UnitfulEquivalences.edconvert(::typeof(dimension(u"m")), x::Unitful.Time,   ::NoEquiv) = x * 10u"m/s"
UnitfulEquivalences.edconvert(::typeof(dimension(u"s")), x::Unitful.Length, ::NoEquiv) = x * (1//10)u"s/m"

@testset "Conversion" begin
    # Equiv1
    @test uconvert(u"ms", 1u"inch", Equiv1()) === (254//100)u"ms"
    @test uconvert(u"km", 1.0u"d", Equiv1()) === 864.0u"km"
    @test uconvert(u"kg", 10u"K", Equiv1()) === 0.02u"kg"
    @test uconvert(u"K", 10u"kg", Equiv1()) === 5000.0u"K"
    @test uconvert(u"eV", 10u"mm/m", Equiv1()) === (1//100)u"eV"
    @test uconvert(u"eV", 10, Equiv1()) === 10u"eV"
    @test uconvert(NoUnits, 10u"eV", Equiv1()) === 10
    @test ustrip(u"ms", 1u"inch", Equiv1()) === (254//100)
    @test ustrip(u"km", 1.0u"d", Equiv1()) === 864.0
    @test ustrip(u"kg", 10u"K", Equiv1()) === 0.02
    @test ustrip(u"K", 10u"kg", Equiv1()) === 5000.0
    @test ustrip(u"eV", 10u"mm/m", Equiv1()) === 1//100
    @test ustrip(u"eV", 10, Equiv1()) === 10
    @test ustrip(NoUnits, 10u"eV", Equiv1()) === 10
    @test ustrip(Float64, u"ms", 1u"inch", Equiv1()) === 2.54
    @test ustrip(Rational{Int}, u"km", 1.0u"d", Equiv1()) === 864//1
    @test ustrip(Float32, u"kg", 10u"K", Equiv1()) === 0.02f0
    @test ustrip(Int, u"K", 10u"kg", Equiv1()) === 5000
    @test ustrip(Float64, u"eV", 10, Equiv1()) === 10.0
    @test_throws ArgumentError uconvert(u"s", 10u"s", Equiv1())
    @test_throws ArgumentError uconvert(u"kg", 1u"s", Equiv1())
    @test_throws ArgumentError uconvert(u"m", 1, Equiv1())
    @test_throws MethodError   uconvert(u"km", 1u"s", Equiv1) # need instance, not type
    @test_throws ArgumentError ustrip(u"s", 10u"s", Equiv1())
    @test_throws ArgumentError ustrip(u"kg", 1u"s", Equiv1())
    @test_throws ArgumentError ustrip(u"m", 1, Equiv1())
    @test_throws MethodError   ustrip(u"km", 1u"s", Equiv1) # need instance, not type
    @test_throws ArgumentError ustrip(Float64, u"s", 10u"s", Equiv1())
    @test_throws ArgumentError ustrip(Float64, u"kg", 1u"s", Equiv1())
    @test_throws ArgumentError ustrip(Float64, u"m", 1, Equiv1())
    @test_throws MethodError   ustrip(Float64, u"km", 1u"s", Equiv1) # need instance, not type

    # Equiv2
    @test uconvert(u"cm", 1u"minute", Equiv2()) === 21_600_000u"cm"
    @test uconvert(u"km", 2u"minute", Equiv2()) === (1728//1)u"km"
    @test uconvert(u"s", -8u"m", Equiv2()) === -2.0u"s"
    @test uconvert(u"s", 27u"m", Equiv2()) === 3.0u"s"

    # Equiv3
    @test uconvert(u"J", 20u"ha", Equiv3(1)) === 0.000_005u"J"
    @test uconvert(u"J", 20u"ha", Equiv3(-4)) === -0.000_02u"J"
    @test uconvert(u"cm^2", 1u"J", Equiv3(2)) === 20_000.0u"cm^2"
    @test uconvert(u"cm^2", 1u"J", Equiv3(5)) === 50_000.0u"cm^2"

    # NoEquiv
    @test_throws MethodError uconvert(u"km", 1u"minute", NoEquiv()) # !(NoEquiv <: Equivalence)
    @test_throws MethodError ustrip(u"km", 1u"minute", NoEquiv()) # !(NoEquiv <: Equivalence)
    @test_throws MethodError ustrip(Float64, u"km", 1u"minute", NoEquiv()) # !(NoEquiv <: Equivalence)

    @testset "Broadcasting" begin
        @test uconvert.(u"ms", [5u"cm", 10u"cm"], Equiv1()) == [5u"ms", 10u"ms"]
        @test uconvert.(u"km", Quantity[1u"hr", -1u"s"], Equiv1()) == [36u"km", (-1//100)u"km"]
        @test ustrip.(u"ms", [5u"cm", 10u"cm"], Equiv1()) == [5, 10]
        @test ustrip.(u"km", Quantity[1u"hr", -1u"s"], Equiv1()) == [36, -1//100]
        @test ustrip.(Float64, u"ms", [5u"cm", 10u"cm"], Equiv1()) == [5.0, 10.0]
        @test ustrip.(Float64, u"km", Quantity[1u"hr", -1u"s"], Equiv1()) == [36.0, -0.01]
    end
end

struct Equiv4 <: Equivalence end
@eqrelation Equiv4 Unitful.Velocity/Unitful.Voltage = 1.5*u"m/(V*s)"
@eqrelation Equiv4 Unitful.Force*Unitful.Volume = -1u"N*m^3"
@eqrelation Equiv4 Unitful.Energy/DimensionlessQuantity = u"eV"

@testset "@eqrelation" begin
    @test uconvert(u"m/s", 10u"V", Equiv4()) === 15.0u"m/s"
    @test uconvert(u"V", 15u"m/s", Equiv4()) === 10.0u"V"
    @test uconvert(u"N", 10u"m^3", Equiv4()) === -0.1u"N"
    @test uconvert(u"m^3", 5u"N", Equiv4()) === -0.2u"m^3"
    @test uconvert(NoUnits, 5u"eV", Equiv4()) === 5
    @test uconvert(u"keV", 5u"km/m", Equiv4()) === 5u"keV"
    @test uconvert(u"MeV", 5, Equiv4()) === (1//200_000)u"MeV"
    @test_throws ArgumentError uconvert(u"m^3", 1u"V", Equiv4())
    @test_throws LoadError @macroexpand @eqrelation Equiv4 Unitful.Energy = Unitful.Mass * Unitful.c0^2
    @test_throws LoadError @macroexpand @eqrelation Equiv4 Unitful.Energy + Unitful.Mass = Unitful.c0^2
end

≈ᵤ(x, y; kwargs...) = unit(x) === unit(y) && ≈(x, y; kwargs...)

@testset "MassEnergy" begin
    @test uconvert(u"keV", 1u"me", MassEnergy()) ≈ᵤ 510.998_95u"keV"
    @test uconvert(u"kg", 938.272_088_16u"MeV", MassEnergy()) ≈ᵤ 1u"mp"
end

@testset "PhotonEnergy" begin
    @test PhotonEnergy() === PhotonEnergy(frequency=:linear, wavelength=:linear, wavenumber=:linear)
    @test PhotonEnergy(frequency=:angular) === PhotonEnergy(frequency=:angular, wavelength=:linear, wavenumber=:linear)
    @test sprint(show, PhotonEnergy()) === "PhotonEnergy(frequency=:linear, wavelength=:linear, wavenumber=:linear)"

    # Energy ↔ frequency
    for L = (:linear, :angular), N = (:linear, :angular)
        @test uconvert(u"eV", 1u"MHz", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=N)) ≈ᵤ 4.135_667_697e-09u"eV" # E = h*f
        @test uconvert(u"J", 1u"ns^-1", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=N)) ≈ᵤ 1.054_571_818e-25u"J" # E = ħ*ω
        @test uconvert(u"Hz", 1u"eV", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=N)) ≈ᵤ 2.417_989_242e+14u"Hz" # f = E/h
        @test uconvert(u"fs^-1", 1u"J", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=N)) ≈ᵤ 9.482_521_562e+18u"fs^-1" # ω = E/ħ
    end

    # Energy ↔ wavelength
    for F = (:linear, :angular), N = (:linear, :angular)
        @test uconvert(u"eV", 1u"nm", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=N)) ≈ᵤ 1239.841_984u"eV" # E = h*c/λ
        @test uconvert(u"J", 1u"Å", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=N)) ≈ᵤ 3.161_526_773e-16u"J" # E = ħ*c/ƛ
        @test uconvert(u"nm", 1u"keV", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=N)) ≈ᵤ 1.239_841_984u"nm" # λ = h*c/E
        @test uconvert(u"μm", 1u"J", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=N)) ≈ᵤ 3.161_526_773e-20u"μm" # ƛ = ħ*c/E
    end

    # Energy ↔ wavenumber
    for F = (:linear, :angular), L = (:linear, :angular)
        @test uconvert(u"eV", 1u"cm^-1", PhotonEnergy(frequency=F, wavelength=L, wavenumber=:linear)) ≈ᵤ 0.000_123_984_1984u"eV" # E = h*c*ν̃
        @test uconvert(u"J", 1u"m^-1", PhotonEnergy(frequency=F, wavelength=L, wavenumber=:angular)) ≈ᵤ 3.161_526_773e-26u"J" # E = ħ*c*k
        @test uconvert(u"cm^-1", 1u"J", PhotonEnergy(frequency=F, wavelength=L, wavenumber=:linear)) ≈ᵤ 5.034_116_568e+22u"cm^-1" # ν̃ = E/(h*c)
        @test uconvert(u"m^-1", 1u"eV", PhotonEnergy(frequency=F, wavelength=L, wavenumber=:angular)) ≈ᵤ 5.067_730_716e+06u"m^-1" # k = E/(ħ*c)
    end

    # Frequency ↔ wavelength
    for N = (:linear, :angular)
        @test uconvert(u"MHz", 1u"nm", PhotonEnergy(frequency=:linear, wavelength=:linear, wavenumber=N)) ≈ᵤ 299_792_458_000u"MHz" # f = c/λ
        @test uconvert(u"ns^-1", 1u"Å", PhotonEnergy(frequency=:linear, wavelength=:angular, wavenumber=N)) ≈ᵤ 4.771_345_159e+08u"ns^-1" # f = c/(2π*ƛ)
        @test uconvert(u"fs^-1", 1u"nm", PhotonEnergy(frequency=:angular, wavelength=:linear, wavenumber=N)) ≈ᵤ 1883.651_567u"fs^-1" # ω = 2π*c/λ
        @test uconvert(u"Hz", 1u"m", PhotonEnergy(frequency=:angular, wavelength=:angular, wavenumber=N)) ≈ᵤ 299_792_458u"Hz" # ω = c/ƛ

        @test uconvert(u"nm", 1u"MHz", PhotonEnergy(frequency=:linear, wavelength=:linear, wavenumber=N)) ≈ᵤ 299_792_458_000u"nm" # λ = c/f
        @test uconvert(u"Å", 1u"ns^-1", PhotonEnergy(frequency=:linear, wavelength=:angular, wavenumber=N)) ≈ᵤ 4.771_345_159e+08u"Å" # ƛ = c/(2π*f)
        @test uconvert(u"nm", 1u"fs^-1", PhotonEnergy(frequency=:angular, wavelength=:linear, wavenumber=N)) ≈ᵤ 1883.651_567u"nm" # λ = 2π*c/ω
        @test uconvert(u"m", 1u"Hz", PhotonEnergy(frequency=:angular, wavelength=:angular, wavenumber=N)) ≈ᵤ 299_792_458u"m" # ƛ = c/ω
    end

    # Frequency ↔ wavenumber
    for L = (:linear, :angular)
        @test uconvert(u"MHz", 1u"cm^-1", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=:linear)) ≈ᵤ 29_979.2458u"MHz" # f = c*ν̃
        @test uconvert(u"ns^-1", 1u"m^-1", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=:angular)) ≈ᵤ 0.047_713_451_59u"ns^-1" # f = c*k/2π
        @test uconvert(u"fs^-1", 1u"cm^-1", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=:linear)) ≈ᵤ 0.000_188_365_1567u"fs^-1" # ω = 2π*c*ν̃
        @test uconvert(u"Hz", 1u"m^-1", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=:angular)) ≈ᵤ 299_792_458u"Hz" # ω = c*k

        @test uconvert(u"cm^-1", 1u"MHz", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=:linear)) ≈ᵤ 3.335_640_952e-05u"cm^-1" # ν̃ = f/c
        @test uconvert(u"m^-1", 1u"ns^-1", PhotonEnergy(frequency=:linear, wavelength=L, wavenumber=:angular)) ≈ᵤ 20.958_450_22u"m^-1" # k = 2π*f/c
        @test uconvert(u"cm^-1", 1u"fs^-1", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=:linear)) ≈ᵤ 5308.837_459u"cm^-1" # ν̃ = ω/(2π*c)
        @test uconvert(u"m^-1", 1u"Hz", PhotonEnergy(frequency=:angular, wavelength=L, wavenumber=:angular)) ≈ᵤ 3.335_640_952e-09u"m^-1" # k = ω/c
    end

    # Wavelength ↔ wavenumber
    for F = (:linear, :angular)
        @test uconvert(u"nm", 1u"cm^-1", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=:linear)) ≈ᵤ 10_000_000u"nm" # λ = 1/ν̃
        @test uconvert(u"Å", 1u"m^-1", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=:angular)) ≈ᵤ 6.283_185_307e+10u"Å" # λ = 2π/k
        @test uconvert(u"nm", 1u"cm^-1", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=:linear)) ≈ᵤ 1.591_549_431e+06u"nm" # ƛ = 1/(2π*ν̃)
        @test uconvert(u"m", 1u"m^-1", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=:angular)) ≈ᵤ 1u"m" # ƛ = 1/k

        @test uconvert(u"cm^-1", 1u"nm", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=:linear)) ≈ᵤ 10_000_000u"cm^-1" # ν̃ = 1/λ
        @test uconvert(u"m^-1", 1u"Å", PhotonEnergy(frequency=F, wavelength=:linear, wavenumber=:angular)) ≈ᵤ 6.283_185_307e+10u"m^-1" # k = 2π/λ
        @test uconvert(u"cm^-1", 1u"nm", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=:linear)) ≈ᵤ 1.591_549_431e+06u"cm^-1" # ν̃ = 1/(2π*ƛ)
        @test uconvert(u"m^-1", 1u"m", PhotonEnergy(frequency=F, wavelength=:angular, wavenumber=:angular)) ≈ᵤ 1u"m^-1" # k = 1/ƛ
    end
end

ext: jl


text:
include("../examples/rijke.jl")
include("../src/adjoint_lss.jl")
using JLD
using OrdinaryDiffEq
function Rijke_adjoint_sensitivity(n_spe)
    s = [6.9, 0.2]
    n = 2000
    d = N

    J_trj = ones(n)
    dJ_trj = zeros(d,n)
    f_trj = zeros(d,n)
    X_trj = zeros(2,d,n)
    du_trj = zeros(d, d, n)

    dJds = zeros(2)
    #vsh = zeros(d, n)
    nRunup = 1000000
    u = Rijke_ODE(rand(d),s,nRunup)

    for i = 1:n
        du_trj[:,:, i] = dRijke(u, s, 1.e-6)
        X_trj[1,:,i] = perturbation(u, s)
        X_trj[2,:,i] = tau_perturbation(u, s)

        #this last parameter is t to be compatible with
        # ODE problem
        f!(view(f_trj, :, i), u, s, 1.)

        vel_p_i = view(u, 1:2*Ng)
        J_trj[i] = 0.25*sum(x->x*x,vel_p_i)
        dJ_trj[1:2*Ng,i] = 0.5*vel_p_i

        u = Rijke_ODE(u, s, 1)

    end
    du_trj = reverse(du_trj, dims=3)
    du_trj = permutedims(du_trj, [2, 1, 3])

    u = Rijke_ODE(u, s, 1)
    X_trj = reverse(X_trj, dims=3)
    #X1_end = perturbation(u, s)
    #X2_end = tau_perturbation(u, s)
    #X_trj[1,:,:] = [X1_end X_trj[1,:,1:end-1]]
    #X_trj[2,:,:] = [X2_end X_trj[2,:,1:end-1]]

    dJ_trj = reverse(dJ_trj, dims=2)
    dJ_trj /= n

    f_trj = reverse(f_trj, dims=2)
    f_end = zeros(d)
    f!(f_end, u, s, 1.)
    f_trj = [f_end f_trj[:,1:end-1]]

    y, xi = lss(du_trj, dJ_trj,
    zeros(d,n), s, 2, f_trj)
    dJds = compute_sens(y, zeros(n), X_trj, zeros(d,n))
    #vsh[:,:] = y
	@show dJds
    save(string("../data/rijke_adjoint_sensitivity/",
        "vsh_and_dJds_", string(n_spe), ".jld"),
             "dJds", dJds)
end
function collect_adjoint_sensitivities()
    pmap(Rijke_adjoint_sensitivity, 9001:10000)
end

#=
    n_samples = 1000
    dJacds_ens, dJrayds_ens = zeros(2,n_samples),
                    zeros(2,n_samples)
    up = zeros(d)
    um = zeros(d)
    eps = 1.e-3
    n = 100 # produces a growth of ≈ e
    for i = 1:n_samples
        @time dJacds_ens[:,i] = Zygote.gradient(s ->Jac(u, s, n), s)[1]
    dJrayds_ens[:,i] = Zygote.gradient(s ->Jray(u, s, n),
                               s)[1]
    u .= Rijke(u, [7.0, 0.2], n)
    end
=#

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
function nnet_layers(nin::Int, nout::Int;
              nh = Int[],
              inner_activation = :tanh,
              final_activation = :identity,
              layernorm = false,
              inputnorm = false,
              kw...)
    ns = vcat(nin, nh, nout)
    num_affine = length(ns) - 1
    layers = Transformation[]

    for i=1:num_affine
        # push!(layers, (layernorm ? LayerNorm : Affine)(ns[i], ns[i+1]))
        if inputnorm
            push!(layers, InputNorm(ns[i]; kw...))
        end

        if layernorm && !(i==num_affine && final_activation == :identity)
            push!(layers, Linear(ns[i], ns[i+1]))
            push!(layers, LayerNorm(ns[i+1]; kw...))
        else
            push!(layers, Affine(ns[i], ns[i+1]))
        end

        if inner_activation != :identity && i < num_affine
            push!(layers, Activation(inner_activation, ns[i+1]))
        end
    end
    if final_activation != :identity
        push!(layers, Activation(final_activation, nout))
    end
    layers
end

function nnet(nin::Int, nout::Int, nh::AbstractVector{Int},
              iact = :tanh, fact = :identity;
              prep = NoPreprocessing(),
              grad_calc = :backprop,
              kw...)
    Chain(Float64, nnet_layers(
        nin, nout;
        nh=nh,
        inner_activation=iact,
        final_activation=fact,
        kw...
    ), prep=prep, grad_calc=grad_calc)
end

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

"""
A ResNet layer: y = F(x) + x

Kaiming He et al (2015): "Deep Residual Learning for Image Recognition"
    https://arxiv.org/pdf/1512.03385v1.pdf

Notes:
    - x and y have the same dimension
    - F can be any transformation mapping n-->n
"""
type ResidualLayer{T,F<:Transformation,P} <: Learnable
    n::Int
    f::F
    input::SumNode{T,1}
    output::OutputNode{T,1}
    params::P
end
ResidualLayer(n::Int; kw...) = ResidualLayer(Float64, n; kw...)
ResidualLayer(n::Int, f::Transformation) = ResidualLayer(Float64, n, f)
ResidualLayer{T}(::Type{T}, n::Int; nh = Int[], kw...) = ResidualLayer(T, n, nnet(n, n, nh; kw...))

params(t::ResidualLayer) = params(t.f)

function ResidualLayer{T}(::Type{T}, n::Int, f::Transformation)
    @assert input_length(f) == n
    @assert output_length(f) == n
    ResidualLayer(
        n,
        f,
        InputNode(T,n),
        OutputNode(T,n),
        f.params
    )
end

function transform!(t::ResidualLayer)
    x = input_value(t)
    y = output_value(t)

    yᶠ = transform!(t.f, x)
    for i=1:t.n
        y[i] = yᶠ[i] + x[i]
    end
    y
end

function grad!{T}(t::ResidualLayer{T})
    ∇y = output_grad(t)
    ∇x = input_grad(t)
    ∇xᶠ = input_grad(t.f)

    grad!(t.f, ∇y)
    for i=1:t.n
        ∇x[i] = ∇xᶠ[i] + ∇y[i]
    end
    ∇x
end

function reset_params!{T}(t::ResidualLayer{T}, θ::AbstractVector, ∇::AbstractVector)
    t.params = reset_params!(t.f, θ, ∇)
end

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

"Construct nblocks layers of stacked ResidualLayer blocks, with a final affine mapping to nout"
function resnet(nin::Int, nout::Int, nblocks::Int; kw...)
    layers = Transformation[ResidualLayer(nin; kw...) for i=1:nblocks]
    push!(layers, Affine(nin, nout))
    Chain(Float64, layers)
end


text:
using StochasticDiffEq, DiffEqDevTools, DiffEqBase, Test
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear
Random.seed!(100)
prob = deepcopy(prob_sde_linear)
prob2 = SDEProblem(prob.f,prob.g,prob.u0,(1.0,0.0))

sol = solve(prob2,EM(),dt=-1//2^(4),tstops = [0.33])

@test sol.t[15] > sol.t[16]

sol = solve(prob2,SRIW1(),tstops = [0.33])

@test sol.t[15] > sol.t[16]

prob = deepcopy(prob_sde_2Dlinear)
prob2 = SDEProblem(prob.f,prob.g,prob.u0,(1.0,0.0))

sol = solve(prob2,EM(),dt=-1//2^(4),tstops = [0.33])

@test sol.t[15] > sol.t[16]

sol = solve(prob2,SRIW1(),tstops = [0.33])

@test sol.t[15] > sol.t[16]


text:
struct RegularLatLonGrid{TX, TY, TZ, R, A} <: AbstractGrid{TX, TY, TZ}
    # corrdinates at cell centers, unit: degree
    xC::R
    yC::R
    # corrdinates at cell centers, unit: meter
    zC::R
    # corrdinates at cell faces, unit: degree
    xF::R
    yF::R
    # corrdinates at cell faces, unit: meter
    zF::R
    # grid spacing, unit: degree
    Δx::Float64
    Δy::Float64
    # grid spacing, unit: meter
    Δz::Float64
    # grid spacing from center to center, unit: meter
    dxC::A
    dyC::A
    # grid spacing from face to face, unit: meter
    dxF::A
    dyF::A
    # areas and volume, unit: m² or m³
    Ax::A
    Ay::A
    Az::A
    Vol::A
    # number of grid points
    Nx::Int
    Ny::Int
    Nz::Int
    # number of halo points
    Hx::Int
    Hy::Int
    Hz::Int
end

"""
    RegularLatLonGrid(;size, lat, lon, z,
                       radius = 6370.0e3,
                       halo = (2, 2, 2))
Creats a `RegularLatLonGrid` struct with `size = (Nx, Ny, Nz)` grid points.

Keyword Arguments (Required)
============================
- `size` : A tuple prescribing the number of grid points.
                `size` is a 3-tuple no matter for 3D, 2D, or 1D model.
- `lat` : A 2-tuple specifying the startind and ending points in latitudinal direction.
                Possible values are from -80 (80S) to 80 (80N).
- `lon` : A 2-tuple specifying the startind and ending points in longitudinal direction.
                Possible values are from -180 (180W) to 180 (180E).
- `z` : A 2-tuple specifying the startind and ending points in vertical direction.
                Vertical indexing starts from surface and use negative numbers for depth.

Keyword Arguments (Optional)
============================
- `radius` : Specify the radius of the Earth used in the model, 6370.0e3 meters by default.
- `halo` : A tuple of integers that specifies the size of the halo region of cells
                surrounding the physical interior for each direction.
                `halo` is a 3-tuple no matter for 3D, 2D, or 1D model.
                At least 2 halo points are needed for DST3FL advection scheme.
"""
function RegularLatLonGrid(;size, lat, lon, z,
                            radius = 6370.0e3,
                            halo = (2, 2, 2))
    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo
    lat₁, lat₂ = lat
    lon₁, lon₂ = lon
    z₁, z₂ = z

    @assert -180 <= lon₁ < lon₂ <= 180
    @assert -80 <= lat₁ < lat₂ <= 80
    @assert z₁ > z₂

    TX = lon₁ == -180 && lon₂ == 180 ? Periodic : Bounded
    TY = Bounded
    TZ = Bounded

    Δx = (lon₂ - lon₁) / Nx
    Δy = (lat₂ - lat₁) / Ny
    Δz = -(z₂ - z₁) / Nz

    xF = range(lon₁ - Hx * Δx, lon₁ + (Nx + Hx - 1) * Δx, length = Nx + 2 * Hx)
    yF = range(lat₁ - Hy * Δy, lat₁ + (Ny + Hy - 1) * Δy, length = Ny + 2 * Hy)
    zF = -range(z₁ - Hz * Δz, z₁ + (Nz + Hz - 1) * Δz, length = Nz + 2 * Hz)

    xC = range(lon₁ + (0.5 - Hx) * Δx, lon₁ + (Nx + Hx - 0.5) * Δx, length = Nx + 2 * Hx)
    yC = range(lat₁ + (0.5 - Hy) * Δy, lat₁ + (Ny + Hy - 0.5) * Δy, length = Ny + 2 * Hy)
    zC = -range(z₁ + (0.5 - Hz) * Δz, z₂ + (Nz + Hz - 0.5) * Δz, length = Nz + 2 * Hz)

    # inclue halo points
    dxC = zeros(Nx+2*Hx, Ny+2*Hy)
    dyC = zeros(Nx+2*Hx, Ny+2*Hy)
    dxF = zeros(Nx+2*Hx, Ny+2*Hy)
    dyF = zeros(Nx+2*Hx, Ny+2*Hy)
    Ax = zeros(Nx+2*Hx, Ny+2*Hy)
    Ay = zeros(Nx+2*Hx, Ny+2*Hy)
    Az = zeros(Nx+2*Hx, Ny+2*Hy)
    Vol = zeros(Nx+2*Hx, Ny+2*Hy)

    for i in 1:Nx+2*Hx
        for j in 1:Ny+2*Hy
            dxC[i, j] = radius * cos(yC[j]*π/180) * deg2rad(Δx)
            dyC[i, j] = radius * deg2rad(Δy)
            dxF[i, j] = radius * cos(yC[j]*π/180) * deg2rad(Δx)
            dyF[i, j] = radius * deg2rad(Δy)
        end
    end
    for i in 1:Nx+2*Hx-1
        for j in 1:Ny+2*Hy-1
            Ax[i,j]   = dyF[i,j] * Δz
            Ay[i,j]   = dxF[i,j] * Δz
            Az[i,j]   = radius^2 * deg2rad(Δx) * (sin(yF[j+1]*π/180) - sin(yF[j]*π/180))
        end
    end
    for i in 1:Nx+2*Hx-1
        for j in 1:Ny+2*Hy-1
            Vol[i,j]   = Az[i,j] * Δz
        end
    end

    return RegularLatLonGrid{TX, TY, TZ, typeof(xF), typeof(dxC)}(
        xC, yC, zC, xF, yF, zF, Δx, Δy, Δz, dxC, dyC, dxF, dyF, Ax, Ay, Az, Vol, Nx, Ny, Nz, Hx, Hy, Hz)
end

function show(io::IO, g::RegularLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ}
    xL, xR = g.xF[g.Hx+1], g.xF[g.Hx+1+g.Nx]
    yL, yR = g.yF[g.Hy+1], g.yF[g.Hy+1+g.Ny]
    zL, zR = g.zF[g.Hz+1], g.zF[g.Hz+1+g.Nz]
    print(io, "RegularLatLonGrid{$TX, $TY, $TZ}\n",
              "domain: x ∈ [$xL, $xR], y ∈ [$yL, $yR], z ∈ [$zL, $zR]\n",
              "topology (Tx, Ty, Tz):     ", (TX, TY, TZ), '\n',
              "resolution (Nx, Ny, Nz):   ", (g.Nx, g.Ny, g.Nz), '\n',
              "halo size (Hx, Hy, Hz):    ", (g.Hx, g.Hy, g.Hz), '\n',
              "grid spacing (Δx, Δy, Δz): ", (g.Δx, g.Δy, g.Δz))
end

function short_show(grid::RegularLatLonGrid{TX, TY, TZ}) where {TX, TY, TZ}
    return "RegularLatLonGrid{$TX, $TY, $TZ}(Nx=$(grid.Nx), Ny=$(grid.Ny), Nz=$(grid.Nz))"
end


text:
using MPI
using CLIMA
using Logging
using Test
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LinearSolvers
using CLIMA.GeneralizedMinimalResidualSolver
using CLIMA.ColumnwiseLUSolver: SingleColumnLU, ManyColumnLU, banded_matrix,
                                banded_matrix_vector_product!
using CLIMA.AdditiveRungeKuttaMethod
using LinearAlgebra
using Printf
using Dates
using CLIMA.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using CLIMA.ODESolvers: solve!, gettime
using CLIMA.VTK: writevtk, writepvtu
using CLIMA.Mesh.Grids: EveryDirection, HorizontalDirection, VerticalDirection

const ArrayType = CLIMA.array_type()

if !@isdefined integration_testing
  if length(ARGS) > 0
    const integration_testing = parse(Bool, ARGS[1])
  else
    const integration_testing =
      parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
  end
end

const output = parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_OUTPUT","false")))

include("advection_diffusion_model.jl")

struct Pseudo1D{n, α, β, μ, δ} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(::Pseudo1D{n, α, β}, aux::Vars,
                                  geom::LocalGeometry) where {n, α, β}
  # Direction of flow is n with magnitude α
  aux.u = α * n

  # diffusion of strength β in the n direction
  aux.D = β * n * n'
end

function initial_condition!(::Pseudo1D{n, α, β, μ, δ}, state, aux, x,
                            t) where {n, α, β, μ, δ}
  ξn = dot(n, x)
  # ξT = SVector(x) - ξn * n
  state.ρ = exp(-(ξn - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
end

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model, testname)
  ## name of the file that this MPI rank will write
  filename = @sprintf("%s/%s_mpirank%04d_step%04d",
                      vtkdir, testname, MPI.Comm_rank(mpicomm), vtkstep)

  statenames = flattenednames(vars_state(model, eltype(Q)))
  exactnames = statenames .* "_exact"

  writevtk(filename, Q, dg, statenames, Qe, exactnames)

  ## Generate the pvtu file for these vtk files
  if MPI.Comm_rank(mpicomm) == 0
    ## name of the pvtu file
    pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

    ## name of each of the ranks vtk files
    prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
      @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
    end

    writepvtu(pvtuprefix, prefixes, (statenames..., exactnames...))

    @info "Done writing VTK: $pvtuprefix"
  end
end

function run(mpicomm, ArrayType, dim, topl, N, timeend, FT, dt,
             n, α, β, μ, δ, vtkdir, outputtime, linearsolvertype)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )
  model = AdvectionDiffusion{dim}(Pseudo1D{n, α, β, μ, δ}())
  dg = DGModel(model,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty(),
               direction=EveryDirection())

  vdg = DGModel(model,
                grid,
                Rusanov(),
                CentralNumericalFluxDiffusive(),
                CentralGradPenalty(),
                auxstate=dg.auxstate,
                direction=VerticalDirection())

  Q = init_ode_state(dg, FT(0))

  ode_solver = ARK548L2SA2KennedyCarpenter(dg, vdg, linearsolvertype(), Q;
                                           dt = dt, t0 = 0,
                                           split_nonlinear_linear=false)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", gettime(ode_solver),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end
  callbacks = (cbinfo,)
  if ~isnothing(vtkdir)
    # create vtk dir
    mkpath(vtkdir)

    vtkstep = 0
    # output initial step
    do_output(mpicomm, vtkdir, vtkstep, dg, Q, Q, model, "advection_diffusion")

    # setup the output callback
    cbvtk = EveryXSimulationSteps(floor(outputtime/dt)) do
      vtkstep += 1
      Qe = init_ode_state(dg, gettime(ode_solver))
      do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model,
                "advection_diffusion")
    end
    callbacks = (callbacks..., cbvtk)
  end

  numberofsteps = convert(Int64, cld(timeend, dt))
  dt = timeend / numberofsteps

  @info "time step" dt numberofsteps dt*numberofsteps timeend

  solve!(Q, ode_solver; numberofsteps=numberofsteps, callbacks=callbacks,
         adjustfinalstep=false)

  # Print some end of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, FT(timeend))

  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
  errf
end

let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 4
  base_num_elem = 4

  expected_result = Dict()
  expected_result[2, 1, Float64] = 7.2741067278166746e-02
  expected_result[2, 2, Float64] = 6.8111330473941490e-03
  expected_result[2, 3, Float64] = 1.4428968208587861e-04
  expected_result[2, 4, Float64] = 2.4253713555277215e-06
  expected_result[3, 1, Float64] = 1.0446155345232400e-01
  expected_result[3, 2, Float64] = 1.0267824273059278e-02
  expected_result[3, 3, Float64] = 2.0613600120241151e-04
  expected_result[3, 4, Float64] = 3.3441210801005078e-06
  expected_result[2, 1, Float32] = 7.2741046547889709e-02
  expected_result[2, 2, Float32] = 6.8110809661448002e-03
  expected_result[2, 3, Float32] = 1.4426559209823608e-04
  # This is near roundoff so we will not check it
  # expected_result[2, 4, Float32] = 2.7985299766442040e-06
  expected_result[3, 1, Float32] = 1.0446154326200485e-01
  expected_result[3, 2, Float32] = 1.0267823934555054e-02
  expected_result[3, 3, Float32] = 2.0618981216102839e-04
  expected_result[3, 4, Float32] = 1.9863322449964471e-05

  numlevels = integration_testing ? 4 : 1

  @testset "$(@__FILE__)" begin
    for FT in (Float64, Float32)
      result = zeros(FT, numlevels)
      for dim = 2:3
        for linearsolvertype in (SingleColumnLU, ManyColumnLU)
          d = dim == 2 ? FT[1, 10, 0] : FT[1, 1, 10]
          n = SVector{3, FT}(d ./ norm(d))

          α = FT(1)
          β = FT(1 // 100)
          μ = FT(-1 // 2)
          δ = FT(1 // 10)
          for l = 1:numlevels
            Ne = 2^(l-1) * base_num_elem
            brickrange = (ntuple(j->range(FT(-1); length=Ne+1, stop=1), dim-1)...,
                          range(FT(-5); length=5Ne+1, stop=5))

            periodicity = ntuple(j->false, dim)
            topl = StackedBrickTopology(mpicomm, brickrange;
                                        periodicity = periodicity,
                                        boundary = (ntuple(j->(1,1), dim-1)...,
                                                    (3,3)))
            dt = (α/4) / (Ne * polynomialorder^2)

            outputtime = 0.01
            timeend = 0.5

            @info (ArrayType, FT, dim, linearsolvertype, l)
            vtkdir = output ? "vtk_advection" *
                              "_poly$(polynomialorder)" *
                              "_dim$(dim)_$(ArrayType)_$(FT)" *
                              "_$(linearsolvertype)_level$(l)" : nothing
            result[l] = run(mpicomm, ArrayType, dim, topl, polynomialorder,
                            timeend, FT, dt, n, α, β, μ, δ, vtkdir,
                            outputtime, linearsolvertype)
            # test the errors significantly larger than floating point epsilon
            if !(dim == 2 && l == 4 && FT == Float32)
              @test result[l] ≈ FT(expected_result[dim, l, FT])
            end
          end
          @info begin
            msg = ""
            for l = 1:numlevels-1
              rate = log2(result[l]) - log2(result[l+1])
              msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
            end
            msg
          end
        end
      end
    end
  end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# Implementation of ideas found in F.F. Wu's "Two stage" paper.
using Test
using ExaPF
using FiniteDiff
using ForwardDiff
using LinearAlgebra
using LinearOperators
using Printf
using KernelAbstractions
using UnicodePlots

import ExaPF: PowerSystem

# This function computes the Davidon cubic interpolation as seen in:
# http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.693.272&rep=rep1&type=pdf
# Notes:
#
#  - In some cases v^2 - F0p*FMp < 0. What to do if we need to keep alpha_m??
function davidon_ls(polar, xk, uk, p, delta_x, delta_u, alpha_m)

    F0 = ExaPF.cost_production(polar, xk, uk, p)
    FM = ExaPF.cost_production(polar, xk + alpha_m*delta_x, uk + alpha_m*delta_u, p)

    function cost_a(a)
        ud = uk + a*delta_u
        xkk, convergence = ExaPF.powerflow(polar, xk, ud, p)
        return ExaPF.cost_production(polar, xkk, ud, p)
    end

    F0p = FiniteDiff.finite_difference_derivative(cost_a, 0.0)
    FMp = FiniteDiff.finite_difference_derivative(cost_a, alpha_m)

    v = (3.0/alpha_m)*(F0 - FM) + F0p + FMp
    w = sqrt(v^2 - F0p*FMp)

    scale = (FMp + w - v)/(FMp - F0p + 2*w)
    alpha = alpha_m - scale*alpha_m

    return alpha
end

# sample along descent line and find minimum.
function sample_ls(nlp, uk, delta_x, delta_u, alpha_m; sample_max=30)

    xk = nlp.x
    alpha = 0.0

    function cost_a(a)
        ud = uk + a*delta_u
        ExaPF.update!(nlp, ud)
        return ExaPF.objective(nlp, ud)
    end

    alpha_vec = collect(range(0, stop=alpha_m, length=sample_max))
    f_vec = zeros(sample_max)

    for i=1:sample_max
        a = alpha_vec[i]
        f_vec[i] = cost_a(a)
    end

    (val, ind) = findmin(f_vec)

    return alpha_vec[ind]
end

# This function computes an approximation of delta_x
# given delta_u
function deltax_approx(delta_u, dGdx, dGdu)
    b = -dGdu*delta_u
    delta_x = dGdx\b
    return delta_x
end

# This function determines the descent direction given gradient
# and limits on control variables.
function descent_direction(nlp, rk, u; damping=false, scale=2.0)
    pf = nlp.model.network
    u_min, u_max = nlp.u_min, nlp.u_max

    dim = length(u)
    delta_u = zeros(dim)
    nref = length(pf.ref)
    npv = length(pf.pv)
    npq = length(pf.pq)

    for i=1:dim
        if u[i] < u_max[i] && u[i] > u_min[i]
            delta_u[i] = -rk[i]
        elseif isapprox(u[i], u_max[i]) && rk[i] > 0.0
            delta_u[i] = -rk[i]
        elseif isapprox(u[i], u_min[i]) && rk[i] < 0.0
            delta_u[i] = -rk[i]
        end
    end

    # u = [VMAG^{REF}, P^{PV}, VMAG^{PV}]

    # scale ratio
    for i=1:npv
        delta_u[nref + i] = scale*delta_u[nref + i]
    end

    # damping factor
    if damping
        for i=1:npv
            idx_u = nref + npv + i
            if u[idx_u] < u_max[idx_u] && u[idx_u] > u_min[idx_u]
                if rk[idx_u] < 0.0
                    damp = min((u_max[idx_u] - u[idx_u]), 1.0)
                    delta_u[idx_u] = damp*delta_u[idx_u]
                elseif rk[idx_u] > 0.0
                    damp = min((u[idx_u] - u_min[idx_u]), 1.0)
                    delta_u[idx_u] = damp*delta_u[idx_u]
                end
            end
        end
    end

    return delta_u
end

# This function determines if the algorithm has converged
function check_convergence(rk, u, u_min, u_max; eps=1e-4)
    dim = length(rk)

    for i=1:dim
        if u[i] < u_max[i] && u[i] > u_min[i] && abs(rk[i]) > eps
            return false
        elseif isapprox(u[i], u_max[i]) && rk[i] > eps
            return false
        elseif isapprox(u[i], u_min[i]) && rk[i] < eps
            return false
        end
    end

    return true
end

# obtain the maximum alpha along the descent line that satisfies all the
# constraints.
function alpha_max(nlp, delta_x, uk, delta_u)
    xk = nlp.x
    x_min, x_max = nlp.x_max, nlp.x_min
    u_min, u_max = nlp.u_max, nlp.u_min

    x_dim = length(delta_x)
    u_dim = length(delta_u)

    alpha_x = 1e10
    alpha_u = 1e10

    for i=1:u_dim
        if abs(delta_u[i]) > 0.0
            am = (u_max[i] - uk[i])/delta_u[i]
            al = (u_min[i] - uk[i])/delta_u[i]
            #@printf("i: %d. delta_u: %f. uk= %f, u_min=%f, u_max=%f, am=%f, al=%f\n",
            #        i, delta_u[i], uk[i], u_min[i], u_max[i], am, al)
            # alpha needs to be positive
            a_prop = max(am, al)
            # need to find alpha that satisfies all constraints
            alpha_u = min(alpha_u, a_prop)
        end
    end

    for i=1:x_dim
        if abs(delta_x[i]) > 0.0 && (x_max[i] > x_min[i])
            am = (x_max[i] - xk[i])/delta_x[i]
            al = (x_min[i] - xk[i])/delta_x[i]
            # alpha needs to be positive
            a_prop = max(am, al)
            # need to find alpha that satisfies all constraints
            alpha_x = min(alpha_x, a_prop)
        end
    end

    # in some cases there seems to be no feasible step with the computed
    # delta_x. maybe this is a problem with using linearized approximation.
    if alpha_x < 0.0
        return alpha_u
    end

    # TODO: we will only calculate alpha_max for u for now.
    # Need to take a closer look to delta_x approximation.
    return alpha_u
end

# given limit alpha, compute costs along a direction.
function cost_direction(nlp, u, delta_u, alpha_max, alpha_dav; points=10)

    alphas = zeros(points)
    costs = zeros(points)
    k = 1

    for a=range(0.0, stop=alpha_max, length=points)
        u_prop = u + a*delta_u
        ExaPF.update!(nlp, u_prop)
        c = ExaPF.objective(nlp, u_prop)

        alphas[k] = a
        costs[k] = c
        k += 1
    end

    #UNCOMMENT THESE LINES
    plt = lineplot(alphas, costs, title = "Cost along alpha", width=80);
    alpha_dav_vert = alpha_dav*ones(points)
    scatterplot!(plt, alpha_dav_vert, costs)
    println(plt)

end

function run_ffwu(datafile; bfgs=false)
    pf = PowerSystem.PowerNetwork(datafile, 1)
    polar = PolarForm(pf, CPU())

    xk = ExaPF.initial(polar, State())
    uk = ExaPF.initial(polar, Control())
    p = ExaPF.initial(polar, Parameters())

    buffer = ExaPF.get(polar, ExaPF.PhysicalState())
    constraints = Function[ExaPF.state_constraint, ExaPF.power_constraints]
    nlp = @time ExaPF.ReducedSpaceEvaluator(polar, xk, uk, p; constraints=constraints)
    jx, ju = ExaPF.init_ad_factory(polar, buffer)

    # reduced gradient method
    iterations = 0
    iter_max = 30
    step = 0.0001
    converged = false
    norm_tol = 1e-5

    # initialize arrays
    grad = similar(uk)
    fill!(grad, 0)
    grad_prev = copy(grad)

    cost_history = zeros(iter_max)
    grad_history = zeros(iter_max)

    if bfgs
        H = InverseLBFGSOperator(Float64, length(uk), 50, scaling=true)
    end

    ExaPF.gradient!(nlp, grad_prev, uk)
    grad = copy(grad_prev)

    iter = 1
    @printf("%6s %12s %6s %6s\n", "iter", "objective", "αₘₐₓ", "αₗₛ")
    while !converged && iter <= iter_max
        # solve power flow and compute gradients
        ExaPF.update!(nlp, uk)

        # evaluate cost
        c = ExaPF.objective(nlp, uk)
        cost_history[iter] = c

        # check convergence
        converged = check_convergence(grad_prev, uk, nlp.u_min, nlp.u_max)

        # compute descent direction
        if bfgs
            delta_u = -H * grad_prev
        else
            ExaPF.gradient!(nlp, grad, uk)
            delta_u = descent_direction(nlp, grad, uk, damping=false, scale=2.0)
        end

        # compute gradients of G
        dGdx = ExaPF.jacobian(polar, jx, nlp.buffer)
        dGdu = ExaPF.jacobian(polar, ju, nlp.buffer)

        # Line search
        # 1st - compute approximate delta_x
        delta_x = deltax_approx(delta_u, dGdx, dGdu)
        # 2nd - compute a_m such that x + a_m*delta_x, u + a_m*delta_u
        # does not violate constraints
        a_m = alpha_max(nlp, delta_x, uk, delta_u)

        # 3rd - fit cubic and find its minimum.
        #a_dav = davidon_ls(pf, nlp.x, uk, p, delta_x, delta_u, a_m)
        a_ls = sample_ls(nlp, uk, delta_x, delta_u, a_m, sample_max=20)
        @printf("%6d %.8e %.6e %.6e %.4e\n", iter, c, a_m, a_ls, norm(delta_u))

        # Uncomment function below to plot the cost function along the descent direction
        # and the calculated alpha.
        # cost_direction(nlp, uk, delta_u, a_m, a_ls; points=100)

        # compute control step
        uk .= uk .+ a_ls*delta_u
        grad_history[iter] = norm(grad)

        # compute gradient
        if bfgs
            ExaPF.gradient!(nlp, grad, uk)
            bfgs && push!(H, a_ls * delta_u, grad - grad_prev)
            grad_prev .= grad
        end
        iter += 1
    end
    # uncomment to plot cost evolution
    # plt = lineplot(cost_history[1:iter - 1], title = "Cost history", width=80);
    # plt = lineplot(log10.(grad_history[1:iter - 1]), title = "Cost history", width=80);
    # println(plt)

    return
end
datafile = joinpath(dirname(@__FILE__), "..", "test", "data", "case9.m")
run_ffwu(datafile; bfgs=false)


text:
#load after ssobol.jl
import RandomNumbers: AbstractRNG
import Random: rand!
import AQFED.Math: norminv
export DigitalSobolSeq, ModuloSobolSeq

mutable struct DigitalSobolSeq{N} <: AbstractSeq{N}
    delegate::ScrambledSobolSeq{N, NoScrambling}
    rng::AbstractRNG{UInt32}
    upoint::Array{UInt32, 1}
end

function DigitalSobolSeq(dimension::Int, n::Int, rng::AbstractRNG{UInt32})
    delegate = ScrambledSobolSeq(dimension, n, NoScrambling())
    upoint = Vector{UInt32}(undef, dimension)
    rand!(rng, upoint)
    return DigitalSobolSeq(delegate, rng, upoint)
end

@inline normalize(s::DigitalSobolSeq, x::UInt32) = ldexp(Float64(x)+0.5, Int32(-32))

#next vector at counter containing all dimensions
@inline function next!(s::DigitalSobolSeq, points::AbstractVector{<:AbstractFloat})
    next(s.delegate, UInt32)
    sx = s.delegate.x
    @inbounds for j = 1:ndims(s)
        points[j] = normalize(s, (sx[j]<< (32-s.delegate.l)) ⊻ s.upoint[j])
    end
    return points
end

@inline function nextn!(s::DigitalSobolSeq, points::AbstractVector{<:AbstractFloat})
    next!(s, points)
    @. points = norminv(points)
    points
end

@inline function nextn!(s::DigitalSobolSeq, dim::Int, points::AbstractVector{<:AbstractFloat})
    next!(s, dim, points)
    @. points = norminv(points)
    points
end
#next vector for a given dimension (vertical) from counter to counter + length(points)
@inline function next!(s::DigitalSobolSeq, dim::Int, points::AbstractVector{<:AbstractFloat})
    j = dim
    sx = s.delegate.x
    sv = s.delegate.v
    @inbounds for i = 1:length(points)
        if s.delegate.counter != 0
            c = ffz(s.delegate.counter)
            sx[j] ⊻= sv[j, c]
        end
        points[i] = normalize(s, (sx[j]<< (32-s.delegate.l)) ⊻ s.upoint[j])
        s.delegate.counter += one(s.delegate.counter)
    end
    return points
end

@inline skipTo(s::DigitalSobolSeq, n::Int) = skipTo(s.delegate, n)
@inline skipTo(s::DigitalSobolSeq, dim::Int, n::Int) = skipTo(s.delegate, dim, n)

mutable struct ModuloSobolSeq{N} <: AbstractSeq{N}
    delegate::ScrambledSobolSeq{N, NoScrambling}
    rng::AbstractRNG
    upoint::Array{Float64, 1}
end

function ModuloSobolSeq(dimension::Int, n::Int, rng::AbstractRNG)
    delegate = ScrambledSobolSeq(dimension, n, NoScrambling())
    upoint = Vector{Float64}(undef, dimension)
    rand!(rng, upoint)
    return ModuloSobolSeq(delegate, rng, upoint)
end

#next vector at counter containing all dimensions
@inline function next!(s::ModuloSobolSeq, points::AbstractVector{<:AbstractFloat})
    next!(s.delegate, points)
    @. points = (points + s.upoint) % 1
    return points
end

#next vector for a given dimension (vertical) from counter to counter + length(points)
@inline function next!(s::ModuloSobolSeq, dim::Int, points::AbstractVector{<:AbstractFloat})
    next!(s.delegate, dim, points)
    uj = s.upoint[j]
    @. points = (points + uj) % 1
    return points
end

@inline skipTo(s::ModuloSobolSeq, n::Int) = skipTo(s.delegate, n)
@inline skipTo(s::ModuloSobolSeq, dim::Int, n::Int) = skipTo(s.delegate, dim, n)


text:
########################################################################### GLP

function getmean!{T<:AbstractFloat}(c::Vector{T}, p::GLP,
                                    features::Vector{T})
   W, b, _ = split_theta!(p)
   A_mul_B!(c, W, features)
   c .+= b
end

# noise is in c vector already
function getaction!{T<:AbstractFloat}(c::Vector{T}, p::GLP,
                                      features::Vector{T})
   W, b, ls = split_theta!(p)
   @. c = c * exp(ls) + b
   # c += A*features
   DT = eltype(p.theta)
   if DT == T
      LinAlg.gemv!('N', DT(1.0), W, features, DT(1.0), c)
   else
      c += W*features
   end
end

function loglikelihood!{T<:AbstractFloat}(ll,
                                          theta::Vector{T}, p::GLP,
                                          features::Matrix{T},
                                          actions::Matrix{T})
   W = reshape( theta[p.ranges[:W]],  p.m, p.n)
   b = reshape( theta[p.ranges[:b]],  p.m)
   ls= reshape( theta[p.ranges[:ls]], p.m)

   dt = eltype(theta)
   terms = dt(-0.5*p.m*log(2pi)) - sum(ls)

   mu = broadcast(+, W*features, b) # W*features + b

   zs = ((actions-mu) ./ exp.(ls)).^2
   ll .= sum(zs, 1)
   ll .= ll.*dt(-0.5) .+ terms
end

function loglikelihood(theta, p::GLP, features, actions)
   W = reshape( theta[p.ranges[:W]],  p.m, p.n)
   b = reshape( theta[p.ranges[:b]],  p.m)
   ls= reshape( theta[p.ranges[:ls]], p.m)

   dt = eltype(theta)
   #terms = dt(-0.5*p.m*log(2pi)) - sum(ls)
   terms = -0.5*p.m*log(2pi) - sum(ls)

   mu = broadcast(+, W*features, b) # W*features + b

   zs = ((actions-mu) ./ exp.(ls)).^2
   ll = sum(zs, 1)
   #ll .= ll.*dt(-0.5) .+ terms
   ll .= ll.*-0.5 .+ terms

   return ll
end

function gradLL!{T<:AbstractFloat}(gradll::Matrix{T},
                                   p::GLP,
                                   features::Matrix{T},
                                   actions::Matrix{T})
   #gradloglikelihood(theta) = loglikelihood(theta, p, features, actions)
   #f_jc = ForwardDiff.JacobianConfig(gradloglikelihood, p.theta)
   #@time ForwardDiff.jacobian!(gradll, gradloglikelihood, p.theta, f_jc)

   #function ll_vec{GT<:PolType}(theta::GT, ff::GT, aa::GT)
   function ll_vec(theta, ff, aa)
      loglikelihood(theta, p, ff, aa)
   end

   const n = p.n
   const m = p.m
   const N = size(features, 2) # N = T*numT

   BLAS.set_num_threads(1)
   nthread = min(N, Threads.nthreads())
   Threads.@threads for tid=1:nthread
   #   for tid=1:nthread
      thread_range = Distributed.splitrange(N, nthread)[tid]

      const ff = zeros(T, n)
      const aa = zeros(T, m)

      const ctp = ReverseDiff.GradientTape(ll_vec, (p.theta, ff, aa))

      gr = map(similar, (p.theta, ff, aa))
      for i=thread_range
         ff .= features[:,i] # implicit conversion and memory saving
         aa .= actions[:,i]

         ReverseDiff.gradient!(gr, ctp, (p.theta, ff, aa)) # taped
         gradll[:,i] .= gr[1][:]
      end
   end
   BLAS.set_num_threads(Threads.nthreads())
end

function split_theta!(p::GLP)
   w_r = p.ranges[:W]
   b_r = p.ranges[:b]
   l_r = p.ranges[:ls]

   W   = reshape( view(p.theta, w_r), p.m, p.n)
   b   = reshape( view(p.theta, b_r), p.m)
   ls  = reshape( view(p.theta, l_r), p.m)

   return W, b, ls
end

function split_theta(p::GLP)
   w_r = p.ranges[:W]
   b_r = p.ranges[:b]
   l_r = p.ranges[:ls]

   W   = reshape( p.theta[w_r], p.m, p.n)
   b   = reshape( p.theta[b_r], p.m)
   ls  = reshape( p.theta[l_r], p.m)

   return W, b, ls
end


text:
using DiffEqJump, DiffEqBase
using Test
using LinearAlgebra
using LightGraphs

function get_mean_sol(jump_prob, Nsims, saveat)
    sol = solve(jump_prob, SSAStepper(), saveat = saveat).u
    for i in 1:Nsims-1
        sol += solve(jump_prob, SSAStepper(), saveat = saveat).u
    end
    sol/Nsims
end

# assume sites are labeled from 1 to num_sites(spatial_system)
function discrete_laplacian_from_spatial_system(spatial_system, hopping_rate)
    sites = 1:num_sites(spatial_system)
    laplacian = zeros(length(sites), length(sites))
    for site in sites
        laplacian[site,site] = -outdegree(spatial_system, site)
        for nb in neighbors(spatial_system, site)
            laplacian[site, nb] = 1
        end
    end
    laplacian .*= hopping_rate
    laplacian
end

# problem setup
majumps = JumpSet(nothing)
tf = 0.5
u0 = [100]
num_species = 1

domain_size = 1.0 #μ-meter
linear_size = 5
diffusivity = 0.1
dim = 2
dims = Tuple([linear_size for i in 1:dim])
num_nodes = prod(dims)

# Starting state setup
starting_state = zeros(Int, length(u0), num_nodes)
center_node = trunc(Int,(num_nodes+1)/2)
starting_state[:,center_node] = copy(u0)
tspan = (0.0, tf)
prob = DiscreteProblem(starting_state,tspan, [])

hopping_rate = diffusivity * (linear_size/domain_size)^2
hopping_constants = [hopping_rate for i in starting_state]

# analytic solution
lap = discrete_laplacian_from_spatial_system(LightGraphs.grid(dims), hopping_rate)
evals, B = eigen(lap) # lap == B*diagm(evals)*B'
Bt = B'
analytic_solution(t) = B*diagm(ℯ.^(t*evals))*Bt * reshape(prob.u0, num_nodes, 1)

num_time_points = 10
Nsims = 10000
rel_tol = 0.01
times = 0.0:tf/num_time_points:tf

algs = [NSM(), DirectCRDirect()]
grids = [CartesianGridRej(dims), CartesianGridIter(dims), LightGraphs.grid(dims)]
jump_problems = JumpProblem[JumpProblem(prob, algs[2], majumps, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for grid in grids]
sizehint!(jump_problems, length(algs) * 3 + length(jump_problems))

# flattenned jump prob
push!(jump_problems, JumpProblem(prob, NRM(), majumps, hopping_constants=hopping_constants, spatial_system = grids[1], save_positions=(false,false)))

# hop rates of form L_{s,i,j}
hop_constants = Matrix{Vector{Float64}}(undef, size(hopping_constants))
for ci in CartesianIndices(hop_constants)
    (species, site) = Tuple(ci)
    hop_constants[ci] = repeat([hopping_constants[species, site]], outdegree(grids[1], site))
end
for alg in algs
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=hop_constants, spatial_system=grids[1], save_positions=(false,false)))
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=hop_constants, spatial_system=grids[3], save_positions=(false,false)))
end

# hop rates of form D_s * L_{i,j}
species_hop_constants = [hopping_rate]
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = repeat([1.0], DiffEqJump.outdegree(grids[1], site))
end
for alg in algs
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=Pair(species_hop_constants, site_hop_constants), spatial_system=grids[1], save_positions=(false,false)))
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=Pair(species_hop_constants, site_hop_constants), spatial_system=grids[3], save_positions=(false,false)))
end

# hop rates of form D_{s,i} * L_{i,j}
species_hop_constants = hopping_rate * ones(1, num_nodes)
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = repeat([1.0], DiffEqJump.outdegree(grids[1], site))
end
for alg in algs
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=Pair(species_hop_constants, site_hop_constants), spatial_system=grids[1], save_positions=(false,false)))
    push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants=Pair(species_hop_constants, site_hop_constants), spatial_system=grids[3], save_positions=(false,false)))
end

# testing
for (j,spatial_jump_prob) in enumerate(jump_problems)
    mean_sol = get_mean_sol(spatial_jump_prob, Nsims, tf/num_time_points)
    for (i,t) in enumerate(times)
        local diff = analytic_solution(t) - reshape(mean_sol[i], num_nodes, 1)
        @test abs(sum(diff[1:center_node])/sum(analytic_solution(t)[1:center_node])) < rel_tol
    end
end

# testing non-uniform hopping rates
dims = (2,2)
num_nodes = prod(dims)
grid = LightGraphs.grid(dims)
hopping_constants = Matrix{Vector{Float64}}(undef, 1, num_nodes)
for ci in CartesianIndices(hopping_constants)
    (species, site) = Tuple(ci)
    hopping_constants[species, site] = zeros(outdegree(grid, site))
    for (n, nb) in enumerate(neighbors(grid, site))
        if nb < site
            hopping_constants[species, site][n] = 1.0
        end
    end
end
starting_state = 25*ones(Int, length(u0), num_nodes)
tspan = (0.0, 10.0)
prob = DiscreteProblem(starting_state,tspan, [])

jp=JumpProblem(prob, NSM(), majumps, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false))
sol = solve(jp, SSAStepper())

@test sol.u[end][1,1] == sum(sol.u[end])


text:
include("hydrogenic_orbitals.jl")

using SymPy
using Polynomials

function Base.:(/)(p::Polynomials.Poly, q::Polynomials.Poly)
    all(q.a .>= 0) || sum(q.a) == 1 ||
        throw(ArgumentError("True polynomial division not implemented"))
    i = findfirst(isequal(1), q.a)
    all(p.a[1:(i-1)] .== 0) ||
        throw(ArgumentError("Polynomial division by $q would not result in a polynomial"))
    Polynomials.Poly(p.a[i:end], p.var)
end

function Ykcomp(n,ℓ,n′,ℓ′,k,Z=1//1)
    N,ρ,P=hydredwfn(Tuple,n,ℓ,Z)
    N′,ρ′,P′=hydredwfn(Tuple,n′,ℓ′,Z)

    a = (polyder(ρ)*polyder(ρ′)).a[1]
    r = Polynomials.Poly([0,1], P.var)
    f = P*P′
    g = k > 0 ? f*r^k : f/r^(abs(k))
    I = Polynomials.Poly([0], f.var)

    i = 1
    while true
        I += (isodd(i) ? 1 : -1)*g/a^i

        g = polyder(g)
        i += 1
        iszero(g) && break
    end
    I,r,a
end

function Y(n,ℓ,n′,ℓ′,k,Z=1//1)
    N,ρ,P=hydredwfn(Tuple,n,ℓ,Z)
    N′,ρ′,P′=hydredwfn(Tuple,n′,ℓ′,Z)

    r = symbols("r", real=true, positive=true)

    ee = exp(-(polyval(ρ, r)+polyval(ρ′, r))/2)
    PP = polyval(P*P′, r)

    Ysmaller = integrate(r^k*PP*ee, (r,0,r))/r^k
    Ygreater = integrate(PP*ee/(r^(k+1)), (r,r,oo))*r^(k+1)
    Yᵏ = N*N′*expand(Ysmaller + Ygreater)
end

function triangle_range(a,b)
    kmin = abs(a-b)
    if !iseven(kmin + a + b)
        kmin += 1
    end
    kmin:2:(a+b)
end

function generate_Yᵏ_references()
    nℓs = [(1,0),
           (2,0), (2,1),
           (3,0), (3,1), (3,2),
           (4,0), (4,1), (4,2), (4,3)]
    ks = 0:3

    open("exact_Yks.jl", "w") do file
        write(file, "exact_Yᵏs = Dict(")
        for ((n,ℓ), (n′,ℓ′), k) in Iterators.product(nℓs, nℓs, ks)
            k ∉ triangle_range(ℓ,ℓ′) && continue
            YY = Y(n,ℓ,n′,ℓ′,k)
            write(file, "(($n,$ℓ), ($(n′),$(ℓ′)), $k) => r -> $YY,\n")
        end
        write(file, ")")
    end
end


text:
using JulieTest

describe ("myTest") do
  it("basic") do
    @is 1 => 1
  end
end


text:
# Script to generate Figure 15.
using Plots, SparseArrays, JLD2, LinearAlgebra, Wavelets, MultiscaleGraphSignalTransforms
include("auxilaries.jl")
include("../../../src/unbalanced_haar_image.jl")

# Set up the resolution and display size
gr(dpi=200, size=(800,600))

# Load the Barbara image data
JLD2.@load "barbara.jld2"

# smaller face region of size 100x100
row_zoom = 71:170
col_zoom = 341:440
display(heatmap(barbara[row_zoom, col_zoom],ratio=1, yaxis =:flip,
    showaxis = false, ticks = false, color = :grays, clims = (0,1)))

# Extract that portion of the image
matrix = deepcopy(barbara[row_zoom, col_zoom])

#
# Perform the image partition using the penalized total variation and then
# compute the expansion coefficients
#
p = 3.0
maxlev = 9 # need this to reach the single node leaves
GPcols = PartitionTreeMatrix_unbalanced_haar_binarytree(matrix, maxlev, p)
maxlev = 8 # need this to reach the single node leaves
GProws = PartitionTreeMatrix_unbalanced_haar_binarytree(copy(matrix'), maxlev, p)
# copy() is necessary to switch the Adjoint type to a regular Matrix{Float64}.
dmatrix = ghwt_analysis_2d(matrix, GProws, GPcols)

#
# Compute the graph Haar coefficients from the previous PTV partition tree
#
BS_haar_rows = bs_haar(GProws)
BS_haar_cols = bs_haar(GPcols)
dvec_haar, loc_haar = BS2loc(dmatrix, GProws, GPcols, BS_haar_rows, BS_haar_cols)

#
# Compute the eGHWT best basis
#
dvec_eghwt, loc_eghwt = eghwt_bestbasis_2d(matrix, GProws, GPcols)

#
# Classical Haar transform via Wavelets.jl with direct input
#
dvec_classichaar = dwt(matrix, wavelet(WT.haar))

#
# Classical Haar transform via Wavelets.jl via putting in the dyadic image
#
matrix2 = zeros(128,128)
matrix2[15:114,15:114] = deepcopy(matrix)
dvec_classichaar2 = dwt(matrix2, wavelet(WT.haar))

#
# Classical Haar transform via Wavelets.jl via putting in the dyadic image
#
matrix3 = zeros(128,128)
matrix3[1:100,1:100] = deepcopy(matrix)
matrix3[101:end,1:100] = deepcopy(matrix[end:-1:end-27,1:100])
matrix3[:,101:end] = matrix3[:,100:-1:73]
dvec_classichaar3 = dwt(matrix3, wavelet(WT.haar))

# Figure 15 (Approximation error plot up to the top 5000 coefficients)
DVEC = [ dvec_classichaar[:], dvec_classichaar2[:], dvec_classichaar3[:],
    dvec_haar[:], dvec_eghwt[:] ]
T = [ "Classic Haar (direct)", "Classic Haar (zero pad)",
    "Classic Haar (even refl)", "Graph Haar (PTV cost)", "eGHWT (PTV cost)" ]
L = [ (:dashdot, :orange), (:dashdot, :red), (:dashdot, :green),
      (:solid, :blue), (:solid, :black) ]
approx_error3(DVEC, T, L, 5000)
savefig("barbara_face100_nondyadic_haar_error.pdf")
savefig("barbara_face100_nondyadic_haar_error.png")


text:
# Autogenerated wrapper script for FastHadamardStructuredTransforms_jll for x86_64-unknown-freebsd
export libfasttransforms

JLLWrappers.@generate_wrapper_header("FastHadamardStructuredTransforms")
JLLWrappers.@declare_library_product(libfasttransforms, "libfasttransforms.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfasttransforms,
        "lib/libfasttransforms.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
Module to enclose [Zernike polynomial](https://en.wikipedia.org/wiki/Zernike_polynomials) specific functionality.
"""
module Zernike

"""
    OSAtoNM(j::Int) -> Tuple{Int, Int}

Convert OSA zernike index `j` to `(N,M)` form according to formula `J = N * (N + 2) + M`.
"""
function OSAtoNM(j::Int)::Tuple{Int, Int}
    n = Int(ceil((-3 + sqrt(9 + 8j)) / 2))
    m = 2j - n * (n + 2)
    return (Int(n), Int(m))
end

"""
    NolltoNM(j::Int) -> Tuple{Int, Int}

Convert Noll zernike index `j` to `(N,M)` form.
"""
function NolltoNM(j::Int)
    n = Int(ceil((-3 + sqrt(1 + 8j)) / 2))
    jr = j - Int(n * (n + 1) / 2)
    if mod(n, 4) ∈ (0, 1)
        m1 = jr
        m2 = -(jr - 1)
        if iseven(n - m1)
            m = m1
        else
            m = m2
        end
    else # mod(n,4) ∈ (2,3)
        m1 = jr - 1
        m2 = -(jr)
        if iseven(n - m1)
            m = m1
        else
            m = m2
        end
    end
    return (Int(n), Int(m))
end

"""
    normalisation(::Type{T}, N::Int, M::Int) -> T

Normalisation coefficient for Zernike polynomial term ``Z_{n}^{m}``.
"""
@inline function normalisation(::Type{T}, N::Int, M::Int) where {T <: Real}
    return T(sqrt((2 * (N + 1)) / (1 + (M == 0 ? 1 : 0))))
end

"""
    ζ(N::Int, M::Int, ρ::T, ϕ::T) -> Tuple{T,T}

Evaluate Zernike polynomial term ``Z_{n}^{m}(\\rho, \\phi)``.
"""
@inline function ζ(N::Int, M::Int, ρ::T, ϕ::T)::T where {T <: Real}
    aM = abs(M)
    if M < 0
        return normalisation(T, N, M) * R(N, aM, ρ) * sin(aM * ϕ)
    else
        return normalisation(T, N, M) * R(N, aM, ρ) * cos(aM * ϕ)
    end
end

"""
    δζ(N::Int, M::Int, ρ::T, ϕ::T) -> Tuple{T,T}

Evaluate partial derivatives of Zernike polynomial term ``Z_{n}^{m}(\\rho, \\phi)``.
"""
@inline function δζ(N::Int, M::Int, ρ::T, ϕ::T)::Tuple{T, T} where {T <: Real}
    n = normalisation(T, N, M)
    aM = abs(M)
    RNM = R(N, aM, ρ)
    ρ2 = ρ^2
    δRNM =
        (
            (2 * N * aM * (ρ2 - 1) + (N - aM) * (aM + N * (2 * ρ2 - 1))) * RNM -
            (N + aM) * (N - aM) * R(N - 2, aM, ρ)
        ) / (2 * N * ρ * (ρ2 - 1))
    Mϕ = aM * ϕ
    if M < 0
        δρ = n * δRNM * sin(Mϕ)
        δϕ = n * RNM * aM * cos(Mϕ)
    else
        δρ = n * δRNM * cos(Mϕ)
        δϕ = n * RNM * aM * -sin(Mϕ)
    end
    return δρ, δϕ
end

"""
    R(N::Int, M::Int, ρ::T) -> T

Evaluate radial polynomial ``R_{n}^{m}(\\rho)``.
"""
@inline function R(N::Int, M::Int, ρ::T)::T where {T <: Real}
    if (N - M) % 2 === 1
        return zero(T)
    end
    total = zero(T)
    @simd for k in 0:((N - M) ÷ 2)
        total +=
            ((-1)^k * factorial(N - k)) / (
                factorial(k) *
                factorial((N + M) ÷ 2 - k) *
                factorial((N - M) ÷ 2 - k)
            ) * ρ^(N - 2 * k)
    end
    return total
end

end # module Zernike

#########################################################################################################

"""
Either `ZernikeIndexingOSA` or `ZernikeIndexingNoll`, see [Zernike polynomials wikipedia entry](https://en.wikipedia.org/wiki/Zernike_polynomials) for details.
"""
@enum ZernikeIndexType ZernikeIndexingOSA ZernikeIndexingNoll

"""
    ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N}

Surface incorporating the Zernike polynomials - radius, conic and aspherics are defined relative to absolute semi-diameter, Zernike terms are normalized according to the `normradius` parameter.
`T` is the datatype, `N` is the dimensionality, `P` is the number of Zernike terms, `Q` is the number of aspheric terms and `M` is the Aspheric Type.

The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space.
Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned.

For convenience the input `zcoeff` can be indexed using either OSA or Noll convention, indicated using the `indexing` argument as either `ZernikeIndexingOSA` or `ZernikeIndexingNoll`.

```julia
ZernikeSurface(semidiameter, radius = Inf, conic = 0, zcoeff = nothing, aspherics = nothing, normradius = semidiameter, indexing = ZernikeIndexingOSA)
```

`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term
for the corresponding `indexing`, or the polynomial power of the aspheric term (may be even or odd) and
`v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively..  `M` will be determined from the terms entered to optimize
the evaluation of the aspheric polynomial.

The sag is defined by the equation

```math
z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir^{2i} + \\sum_{i}^PA_iZ_i(\\rho, \\phi)
```

where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial.
"""
struct ZernikeSurface{T, N, P, Q, M} <: ParametricSurface{T, N}
    asp::AsphericSurface{T, N, Q, M}
    coeffs::SVector{P, Tuple{Int, Int, T}}
    boundingcylinder::Cylinder{T, N}

    function ZernikeSurface(
        semidiameter::T;
        radius::T = typemax(T),
        conic::T = zero(T),
        zcoeff::Union{Nothing, Vector{Tuple{Int, T}}} = nothing,
        aspherics::Union{Nothing, Vector{Tuple{Int, T}}} = nothing,
        normradius::T = semidiameter,
        indexing::ZernikeIndexType = ZernikeIndexingOSA,
    ) where {T <: Real}
        asp =
            AsphericSurface(semidiameter; radius, conic, aspherics, normradius)
        Q = length(asp.aspherics) #this is not the same as the aspherics variable passed to the function!

        zcs = []
        if zcoeff !== nothing
            for (i, k) in zcoeff
                if abs(k) > zero(T)
                    if indexing === ZernikeIndexingOSA
                        R, S = Zernike.OSAtoNM(i)
                    else
                        R, S = Zernike.NolltoNM(i)
                    end
                    push!(zcs, (R, S, k))
                end
            end
        end
        P = length(zcs)
        M = asphericType(asp)
        return new{T, 3, P, Q, M}(
            asp::AsphericSurface{T, 3, Q, M},
            SVector{P, Tuple{Int, Int, T}}(zcs),
            Cylinder(semidiameter, interface = opaqueinterface(T)),
        ) # TODO!! incorrect interface on cylinder
    end
end

function uvrange(
    ::Type{ZernikeSurface{T, N, P, Q, M}},
) where {T <: Real, N, P, Q, M}
    return ((zero(T), one(T)), (-T(π), T(π)))
end # ρ and ϕ

semidiameter(z::ZernikeSurface{T}) where {T <: Real} = z.asp.semidiameter
halfsizeu(z::ZernikeSurface{T}) where {T <: Real} = semidiameter(z)
halfsizev(z::ZernikeSurface{T}) where {T <: Real} = semidiameter(z)

boundingobj(z::ZernikeSurface{T}) where {T <: Real} = z.boundingcylinder

function point(
    z::ZernikeSurface{T, 3, P, Q, M},
    ρ::T,
    ϕ::T,
)::SVector{3, T} where {T <: Real, P, Q, M}
    pnt = point(z.asp, ρ, ϕ)

    # sum zernike
    rad = semidiameter(z.asp)
    r = ρ * rad
    u = r / z.asp.normradius
    h = zero(T)
    @inbounds @simd for m in 1:P
        (R, S, k) = z.coeffs[m]
        h += k * Zernike.ζ(R, S, u, ϕ)
    end
    return SVector{3, T}(pnt[1], pnt[2], pnt[3] + h)
end

function partials(
    z::ZernikeSurface{T, 3, P, Q, M},
    ρ::T,
    ϕ::T,
)::Tuple{SVector{3, T}, SVector{3, T}} where {T <: Real, P, Q, M}
    pρ, pϕ = partials(z.asp, ρ, ϕ)
    # sum zernike partials
    rad = z.asp.semidiameter
    n = rad / z.asp.normradius
    u = ρ * n
    dhdρ = zero(T)
    dhdϕ = zero(T)
    @inbounds @simd for m in 1:P
        (R, S, k) = z.coeffs[m]
        du, dϕ = Zernike.δζ(R, S, u, ϕ)
        dhdρ += k * du * n # want the derivative wrt ρ, not u
        dhdϕ += k * dϕ
    end
    return SVector{3, T}(pρ[1], pρ[2], pρ[3] + dhdρ),
    SVector{3, T}(pϕ[1], pϕ[2], pϕ[3] + dhdϕ)
end

function normal(
    z::ZernikeSurface{T, 3, P, Q, M},
    ρ::T,
    ϕ::T,
)::SVector{3, T} where {T <: Real, P, Q, M}
    du, dv = partials(z, ρ, ϕ)
    if ρ == zero(T) && norm(dv) == zero(T)
        # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric)
        # then we get some big problems, hardcoding this case solves the problems
        return SVector{3, T}(0, 0, 1)
    end
    return normalize(cross(du, dv))
end

function uv(
    z::ZernikeSurface{T, 3, P, Q, M},
    p::SVector{3, T},
) where {T <: Real, P, Q, M}
    # avoid divide by zero for ForwardDiff
    ϕ = NaNsafeatan(p[2], p[1])
    if p[1] == zero(T) && p[2] == zero(T)
        ρ = zero(T)
    else
        ρ = sqrt(p[1]^2 + p[2]^2) / semidiameter(z)
    end
    return SVector{2, T}(ρ, ϕ)
end

function onsurface(
    surf::ZernikeSurface{T, 3, P, Q, M},
    p::SVector{3, T},
) where {T <: Real, P, Q, M}
    ρ, ϕ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, ϕ)
        return samepoint(p[3], surfpoint[3])
    end
end

function inside(
    surf::ZernikeSurface{T, 3, P, Q, M},
    p::SVector{3, T},
) where {T <: Real, P, Q, M}
    ρ, ϕ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, ϕ)
        return p[3] < surfpoint[3]
    end
end

#########################################################################################################

# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(
    surf::AcceleratedParametricSurface{T, 3, ZernikeSurface{T, 3, P, Q, M}},
    r::AbstractRay{T, 3},
) where {T <: Real, P, Q, M}
    cylint = surfaceintersection(surf.surface.boundingcylinder, r)
    if cylint isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        if doesintersect(surf.triangles_bbox, r) ||
           inside(surf.triangles_bbox, origin(r))
            surfint = triangulatedintersection(surf, r)
            if !(surfint isa EmptyInterval{T})
                return intervalintersection(cylint, surfint)
            end
        end
        # hasn't hit the surface
        if lower(cylint) isa RayOrigin{T} && upper(cylint) isa Infinity{T}
            if inside(surf.surface, origin(r))
                return Interval(RayOrigin(T), Infinity(T))
            else
                return EmptyInterval(T)
            end
            # otherwise check that the intersection is underneath the surface
        else
            p = point(closestintersection(cylint, false))
            ρ, ϕ = uv(surf, p)
            surfpoint = point(surf.surface, ρ, ϕ)
            if p[3] < surfpoint[3]
                return cylint # TODO!! UV (and interface) issues?
            else
                return EmptyInterval(T)
            end
        end
    end
end

function AcceleratedParametricSurface(
    surf::T,
    numsamples::Int = 17;
    interface::NullOrFresnel{S} = NullInterface(S),
) where {S <: Real, N, T <: ZernikeSurface{S, N}}
    # Zernike uses ρ, ϕ uv space so need to modify extension of triangulation
    a = AcceleratedParametricSurface(
        surf,
        triangulate(surf, numsamples, true, false, true, false),
        interface = interface,
    )
    emptytrianglepool!(S)
    return a
end

function BoundingBox(
    surf::ZernikeSurface{T, 3, P, Q, M},
) where {T <: Real, P, Q, M}
    bb = BoundingBox(surf.asp)
    # zernike terms have condition than |Zᵢ| <= 1
    # so this gives us a (loose) bounding box
    ak = if P > 0
        sum(abs.(Zernike.normalisation(T, n, m) * k for (n, m, k) in surf.coeffs))
    else
        zero(T)
    end
    bb.zmin -= ak
    bb.zmax += ak
    return BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymax, bb.zmin, bb.zmax) #could just return bb, but this way is safer
end


text:
#==============================================================================#
# OpenSCADIcosahedron
# December 2014
#==============================================================================#

# This module uses OpenSCADPolyhedron to build an icosahedron.
# All points, transforms, and symmetries were derived in Mathematica.

#-------------------------------------------------------------------------------
module OpenSCADIcosahedron
using OpenSCADPolyhedron
export icosahedronPoints, icosahedronFace, icosahedronBase, icosahedronTransforms, icosahedronSymmetries

#-------------------------------------------------------------------------------
# We start with the top face (our icosahedron will have its top and bottom faces
# parallel to the ground, to facilitate 3D printing).
icosahedronPoints = [
  Point3{Float64}([-1., -0.57735, 1.51152]),
  Point3{Float64}([ 0.,  1.1547,  1.51152]),
  Point3{Float64}([ 1., -0.57735, 1.51152]),
]
icosahedronBase = Polyhedron{Float64}(icosahedronPoints)
icosahedronFace = Int[1, 2, 3]
append!(icosahedronBase, icosahedronFace)

icosahedronTransforms = Function[
  #Group 1 ("top 5")
  p::Point3{Float64} -> Point3{Float64}([ 1.        0.        0.;        0.        1.        0.;        0.        0.        1.      ] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.5      -0.645497 -0.57735;   0.866025  0.372678  0.333333;  0.       -0.666667  0.745356] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.309017 -0.178411 -0.934172;  0.755761 -0.64235  -0.127322; -0.57735  -0.745356  0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.309017  0.755761 -0.57735;  -0.178411 -0.64235  -0.745356; -0.934172 -0.127322  0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.5       0.866025  0.;       -0.645497  0.372678 -0.666667; -0.57735   0.333333  0.745356] * p.x),
  #Group 2 ("bottom 5")
  p::Point3{Float64} -> Point3{Float64}([ 1.        0.        0.;        0.       -1.        0.;        0.        0.       -1.      ] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.5       0.645497  0.57735;   0.866025 -0.372678 -0.333333;  0.        0.666667 -0.745356] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.309017  0.178411  0.934172;  0.755761  0.64235   0.127322; -0.57735   0.745356 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.309017 -0.755761  0.57735;  -0.178411  0.64235   0.745356; -0.934172  0.127322 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.5      -0.866025  0.;       -0.645497 -0.372678  0.666667; -0.57735  -0.333333 -0.745356] * p.x),
  #Group 4 ("downward pointing faces")
  p::Point3{Float64} -> Point3{Float64}([ 0.309017 -0.755761  0.57735;   0.755761  0.563661  0.333333; -0.57735   0.333333  0.745356] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.       -0.934172 -0.356822;  0.356822 -0.333333  0.872678; -0.934172 -0.127322  0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.309017 -0.178411 -0.934172; -0.178411 -0.975684  0.127322; -0.934172  0.127322 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.809017  0.467086 -0.356822; -0.110264 -0.475684 -0.872678; -0.57735   0.745356 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.809017  0.110264  0.57735;   0.467086  0.475684 -0.745356; -0.356822  0.872678  0.333333] * p.x),
  #Group 4 ("upward pointing faces")
  p::Point3{Float64} -> Point3{Float64}([-0.809017 -0.467086  0.356822;  0.467086 -0.14235   0.872678; -0.356822  0.872678  0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-0.5      -0.645497 -0.57735;  -0.645497 -0.166667  0.745356; -0.57735   0.745356 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.5      -0.645497 -0.57735;  -0.866025 -0.372678 -0.333333;  0.        0.666667 -0.745356] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.809017 -0.467086  0.356822;  0.110264 -0.475684 -0.872678;  0.57735   0.745356 -0.333333] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.       -0.356822  0.934172;  0.934172 -0.333333 -0.127322;  0.356822  0.872678  0.333333] * p.x),
]

icosahedronSymmetries = [
  #Group 1
  (2,  1) => (1, 1),
  (2,  3) => (1, 2),
  (3,  1) => (1, 1),
  (3,  3) => (2, 2),
  (4,  1) => (1, 1),
  (4,  3) => (3, 2),
  (5,  1) => (1, 1),
  (5,  3) => (4, 2),
  (5,  2) => (1, 3),
  # Group 2
  (7,  1) => (6, 2),
  (7,  3) => (6, 3),
  (8,  1) => (7, 2),
  (8,  3) => (6, 3),
  (9,  1) => (8, 2),
  (9,  3) => (6, 3),
  (10, 1) => (9, 2),
  (10, 3) => (6, 3),
  (10, 2) => (6, 1),
  # Group 3
  (16, 1) => (8, 2),
  (16, 2) => (1, 2),
  (16, 3) => (9, 2),
  (17, 1) => (9, 2),
  (17, 2) => (2, 2),
  (17, 3) => (6, 1),
  (18, 1) => (6, 1),
  (18, 2) => (3, 2),
  (18, 3) => (6, 2),
  (19, 1) => (6, 2),
  (19, 2) => (4, 2),
  (19, 3) => (7, 2),
  (20, 1) => (7, 2),
  (20, 2) => (1, 3),
  (20, 3) => (8, 2),
  # Group 4
  (11, 1) => (1, 3),
  (11, 2) => (1, 2),
  (11, 3) => (8, 2),
  (12, 1) => (1, 2),
  (12, 2) => (2, 2),
  (12, 3) => (9, 2),
  (13, 1) => (2, 2),
  (13, 2) => (3, 2),
  (13, 3) => (6, 1),
  (14, 1) => (3, 2),
  (14, 2) => (4, 2),
  (14, 3) => (6, 2),
  (15, 1) => (4, 2),
  (15, 2) => (1, 3),
  (15, 3) => (7, 2),
]

#-------------------------------------------------------------------------------
# Build the icosahedron polygon.
icosahedron = generatePolyhedron(icosahedronBase, icosahedronTransforms, icosahedronSymmetries)
#println(stringify(icosahedron))

#-------------------------------------------------------------------------------
# Bonus: this base polyhedron and set of transforms can be used to generate an icosahedron
# of a different orientation. The transforms are slightly cleaner, but the resulting model
# only has one edge that touches the ground so cannot stand on its own without support.

icosahedronPoints2 = [
  Point3{Float64}([-1., 0.,      1.61803]),
  Point3{Float64}([ 0., 1.61803, 1.     ]),
  Point3{Float64}([ 1., 0.,      1.61803])
]
icosahedronBase2 = Polyhedron{Float64}(icosahedronPoints)
append!(icosahedronBase2, Int[1, 2, 3])

v1 = 0.309017
v2 = 0.809017
icosahedronTransforms2 = Function[
  # Group 1
  p::Point3{Float64} -> p,
  p::Point3{Float64} -> Point3{Float64}([ .5 -v2 -v1;  v2  v1  .5; -v1 -.5  v2] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-v1 -.5 -v2;  .5 -v2  v1; -v2 -v1  .5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-v1 0.5 -v2; -.5 -v2 -v1; -v2  v1  .5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ .5  v2 -v1; -v2  v1 -.5; -v1  .5  v2] * p.x),
  # Group 2
  p::Point3{Float64} -> Point3{Float64}([ 1.  0.  0.;  0. -1.  0.;  0.  0. -1.] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ .5  v2  v1;  v2 -v1 -.5; -v1  .5 -v2] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-v1  .5  v2;  .5  v2 -v1; -v2  v1 -.5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-v1 -.5  v2; -.5  v2  v1; -v2 -v1 -.5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ .5 -v2  v1; -v2 -v1  .5; -v1 -.5 -v2] * p.x),
  # Group 3
  p::Point3{Float64} -> Point3{Float64}([ v1 -.5  v2;  .5  v2  v1; -v2  v1  .5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0. -1.  0.;  0.  0.  1.; -1.  0.  0.] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ v1 -.5 -v2; -.5 -v2  v1; -v2  v1 -.5] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ v2  v1 -.5; -v1 -.5 -v2; -.5  v2 -v1] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ v2  v1  .5;  v1  .5 -v2; -.5  v2  v1] * p.x),
  # Group 4
  p::Point3{Float64} -> Point3{Float64}([-v2 -v1  .5;  v1  .5  v2; -.5  v2 -v1] * p.x),
  p::Point3{Float64} -> Point3{Float64}([-.5 -v2 -v1; -v2  v1  .5; -v1  .5 -v2] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ .5 -v2 -v1; -v2 -v1 -.5;  v1  .5 -v2] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ v2 -v1  .5;  v1 -.5 -v2;  .5  v2 -v1] * p.x),
  p::Point3{Float64} -> Point3{Float64}([ 0.  0.  1.;  1.  0.  0.;  0.  1.  0.] * p.x)
]

end

max_stars_repo_name: erikstrand/OpenSCADPolyhedron

max_issues_repo_name: erikstrand/OpenSCADPolyhedron

max_forks_repo_name: erikstrand/OpenSCADPolyhedron


text:
using Test
using GoogleCloud

creds = JSONCredentials("/secrets/google-secret.json")

session = GoogleSession(creds, ["devstorage.full_control"])

set_session!(storage, session)

bkts = storage(:Bucket, :list)


text:
using Statistics
using StatsBase:sample
using Revise
using EvoTrees
using MemoryConstrainedTreeBoosting

nrounds = 100

# EvoTrees params
params_evo = EvoTreeRegressor(T=Float32,
        loss=:logistic, metric=:logloss,
        nrounds=nrounds,
        λ=0.5, γ=0.0, η=0.05,
        max_depth=6, min_weight=1.0,
        rowsample=1.0, colsample=0.5, nbins=64)

# MemoryConstrainedTreeBoosting params
params_mctb = (
        weights                 = nothing,
        bin_count               = 128,
        iteration_count         = nrounds,
        min_data_weight_in_leaf = 1.0,
        l2_regularization       = 0.0,
        max_leaves              = 64,
        max_depth               = 6,
        max_delta_score         = 1.0e10, # Before shrinkage.
        learning_rate           = 0.05,
        feature_fraction        = 0.5, # Per tree.
        bagging_temperature     = 0.1,
      )

nobs = Int(1e6)
num_feat = Int(100)
@info "testing with: $nobs observations | $num_feat features."
X = rand(Float32, nobs, num_feat)
Y = Float32.(rand(Bool, size(X, 1)))

@info "evotrees train CPU:"
params_evo.device = "cpu"
@time m, cache = EvoTrees.init_evotree(params_evo, X, Y);
@time EvoTrees.grow_evotree!(m, cache);
@time m, cache = EvoTrees.init_evotree(params_evo, X, Y);
@time EvoTrees.grow_evotree!(m, cache);
@time m_evo = fit_evotree(params_evo, X, Y);
@time fit_evotree(params_evo, X, Y);
@info "evotrees predict CPU:"
@time pred_evo = EvoTrees.predict(m_evo, X);
@time EvoTrees.predict(m_evo, X);

@info "evotrees train GPU:"
params_evo.device = "gpu"
@time m_evo = fit_evotree(params_evo, X, Y);
@time fit_evotree(params_evo, X, Y);
@info "evotrees predict GPU:"
@time pred_evo = EvoTrees.predict(m_evo, X);
@time EvoTrees.predict(m_evo, X);

@info "MemoryConstrainedTreeBoosting train CPU:"
@time bin_splits, trees = MemoryConstrainedTreeBoosting.train(X, Y; params_mctb...);
@time MemoryConstrainedTreeBoosting.train(X, Y; params_mctb...);
@info "MemoryConstrainedTreeBoosting predict CPU, JITed:"
save_path = tempname()
MemoryConstrainedTreeBoosting.save(save_path, bin_splits, trees)
unbinned_predict = MemoryConstrainedTreeBoosting.load_unbinned_predictor(save_path)
@time pred_mctb = unbinned_predict(X)
@time unbinned_predict(X)


text:
pn = TimePetriNet(
    [:red_1,:yellow_1,:green_1, :red_2, :yellow_2, :green_2, :free],
    [:ry1,:yg1,:gr1, :ry2,:yg2,:gr2],
    sparse([
        1;0;0 ; 0;0;0 ; 1 ;; # traffic light 1 red->yellow
        0;1;0 ; 0;0;0 ; 0 ;; # traffic light 1 yellow->green
        0;0;1 ; 0;0;0 ; 0 ;; # traffic light 1 green->red
        0;0;0 ; 1;0;0 ; 1 ;; # traffic light 2 red->yellow
        0;0;0 ; 0;1;0 ; 0 ;; # traffic light 2 yellow->green
        0;0;0 ; 0;0;1 ; 0    # traffic light 2 green->red
    ]),
    sparse([
        -1; 1; 0 ;  0; 0; 0 ; -1 ;; # traffic light 1 red->yellow
         0;-1; 1 ;  0; 0; 0 ;  0 ;; # traffic light 1 yellow->green
         1; 0;-1 ;  0; 0; 0 ;  1 ;; # traffic light 1 green->red
         0; 0; 0 ; -1; 1; 0 ; -1 ;; # traffic light 2 red->yellow
         0; 0; 0 ;  0;-1; 1 ;  0 ;; # traffic light 2 yellow->green
         0; 0; 0 ;  1; 0;-1 ;  1    # traffic light 2 green->red
    ]),
    ones(6),
    ones(6),
    [5;0;0;0;1;0;0]
)

f,ax,p = graphplot(pn); f

#=
trange = 1:60*length(states)
record(f, "traffic_light.mp4", trange; framerate=30) do tt
    mm = states[1+div(tt-1,60)]
    animate_step!(p, pn, mm)
    update!(f.scene)
end
=#


text:
if basename(pwd()) == "aoc"
    cd("2017/10")
end

struct CircularBuffer
    data
end

CircularBuffer(n::Int) = CircularBuffer(collect(0:n-1))

function Base.getindex(buffer::CircularBuffer, index::Integer)
    buffer.data[index%length(buffer.data)+1]
end

function Base.getindex(buffer::CircularBuffer, indices::AbstractRange{<:Integer})
    collect(buffer[index] for index in indices)
end

function Base.setindex!(buffer::CircularBuffer, value, index::Integer)
    buffer.data[index%length(buffer.data)+1] = value
end

function Base.reverse!(buffer::CircularBuffer, start = 0, stop = length(buffer.data) - 1)
    for (i, v) in zip(start:stop, reverse(buffer[start:stop]))
        buffer[i] = v
    end
    buffer
end

function loadlengths(filename::AbstractString)
    parse.(Int, split(readline(filename), ","))
end

function knothash!(buffer::CircularBuffer, lengths)
    position = 0
    for (skip, len) in zip(Iterators.countfrom(0), lengths)
        reverse!(buffer, position, position + len - 1)
        position += skip + len
    end
    buffer
end

function part1(filename::AbstractString, n = 256)
    prod(knothash!(CircularBuffer(n), loadlengths(filename))[0:1])
end

@assert knothash!(CircularBuffer(5), [3, 4, 1, 5])[0:1] == [3, 4]
# part1("input.txt")

function sparsehash!(buffer, input::AbstractString)
    lengths = [(Int(c) for c in input)..., 17, 31, 73, 47, 23]
    position = 0
    for (skip, len) in zip(Iterators.countfrom(0), repeat(lengths, 64))
        reverse!(buffer, position, position + len - 1)
        position += skip + len
    end
    buffer
end
sparsehash(input::AbstractString) = sparsehash!(CircularBuffer(256), input)

function densehash(sparsehash::CircularBuffer)
    join(map(Iterators.partition(sparsehash.data, 16)) do section
        string(reduce(xor, section), base = 16, pad = 2)
    end)
end

densehash(input::AbstractString) = densehash(sparsehash(input))
part2(filename::AbstractString) = densehash(readline(filename))

@assert densehash("") == "a2582a3a0e66e6e86e3812dcb672a272"
@assert densehash("AoC 2017") == "33efeb34ea91902bb2f59c9920caa6cd"
@assert densehash("1,2,3") == "3efbe78a8d82f29979031a4aa0b16a9d"
@assert densehash("1,2,4") == "63960835bcdc130f0b66d7ff4f6a5a8e"

# part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code