text:
using PkgEval
using DataFrames
using Feather
using JSON
using Base: UUID
using LibGit2

################################
# Package Selection Validation #
################################

# The package selection is valid if it is simply a single package, "ALL", or an
# array expression that lists several packages. This validation is
# only to prevent server-side evaluation of arbitrary code. No check is
# performed to ensure that the tag predicate is grammatically correct.

function is_valid_pkgsel(pkgsel::AbstractString)
    parsed = Meta.parse(pkgsel)
    if isa(parsed, Expr)
        return is_valid_stringvector(parsed)
    elseif parsed == :ALL
        return true
    else
        return isa(parsed, AbstractString)
    end
end

function is_valid_stringvector(pkgsel::Expr)
    if pkgsel.head != :vect
        return false
    else
        for item in pkgsel.args
            if !isa(item, AbstractString)
                return false
            end
        end
    end
    return true
end

##############
# PkgEvalJob #
##############

mutable struct PkgEvalJob <: AbstractJob
    submission::JobSubmission        # the original submission
    pkgsel::String                   # selection of packages
    against::Union{BuildRef,Nothing} # the comparison build (if available)
    date::Dates.Date                 # the date of the submitted job
    isdaily::Bool                    # is the job a daily job?
    buildflags::Vector{String}       # a list of build flags for Make.user generation
    against_buildflags::Vector{String}
    # FIXME: put build flags in BuildRef? currently created too early for that (when the
    #        GitHub event is parsed, while we get the build flags from the comment)
end

function PkgEvalJob(submission::JobSubmission)
    if haskey(submission.kwargs, :vs)
        againststr = Meta.parse(submission.kwargs[:vs])
        if in(SHA_SEPARATOR, againststr) # e.g. againststr == christopher-dG/julia@e83b7559df94b3050603847dbd6f3674058027e6
            reporef, againstsha = split(againststr, SHA_SEPARATOR)
            againstrepo = isempty(reporef) ? submission.config.trackrepo : reporef
            againstbuild = BuildRef(againstrepo, againstsha)
        elseif in(BRANCH_SEPARATOR, againststr)
            reporef, againstbranch = split(againststr, BRANCH_SEPARATOR)
            againstrepo = isempty(reporef) ? submission.config.trackrepo : reporef
            againstbuild = branchref(submission.config, againstrepo, againstbranch)
        elseif in(TAG_SEPARATOR, againststr)
            reporef, againsttag = split(againststr, TAG_SEPARATOR)
            againstrepo = isempty(reporef) ? submission.config.trackrepo : reporef
            againstbuild = tagref(submission.config, againstrepo, againsttag)
        else
            error("invalid argument to `vs` keyword")
        end
        against = againstbuild
    else
        against = nothing
    end

    if haskey(submission.kwargs, :isdaily)
        isdaily = submission.kwargs[:isdaily] == "true"
    else
        isdaily = false
    end

    if haskey(submission.kwargs, :buildflags)
        expr = Meta.parse(submission.kwargs[:buildflags])
        if !is_valid_stringvector(expr)
            error("invalid argument to `buildflags` keyword")
        end
        buildflags = eval(expr)
    else
        buildflags = String[]
    end

    if haskey(submission.kwargs, :vs_buildflags)
        expr = Meta.parse(submission.kwargs[:vs_buildflags])
        if !is_valid_stringvector(expr)
            error("invalid argument to `vs_buildflags` keyword")
        end
        against_buildflags = eval(expr)
    else
        against_buildflags = String[]
    end

    return PkgEvalJob(submission, first(submission.args), against,
                      Dates.today(), isdaily, buildflags, against_buildflags)
end

function Base.summary(job::PkgEvalJob)
    result = "PkgEvalJob $(summary(submission(job).build))"
    if job.isdaily
        result *= " [daily]"
    elseif job.against !== nothing
        result *= " vs. $(summary(job.against))"
    end
    return result
end

function isvalid(submission::JobSubmission, ::Type{PkgEvalJob})
    allowed_kwargs = (:vs, :isdaily, :buildflags, :vs_buildflags)
    args, kwargs = submission.args, submission.kwargs
    has_valid_args = length(args) == 1 && is_valid_pkgsel(first(args))
    has_valid_kwargs = (all(in(allowed_kwargs), keys(kwargs)) &&
                        (length(kwargs) <= length(allowed_kwargs)))
    return (submission.func == "runtests") && has_valid_args && has_valid_kwargs
end

submission(job::PkgEvalJob) = job.submission

#############
# Utilities #
#############

function jobdirname(job::PkgEvalJob; latest::Bool=false)
    if job.isdaily
        joinpath("by_date", latest ? "latest" : datedirname(job.date))
    else
        primarysha = snipsha(submission(job).build.sha)
        tag = if job.against === nothing
            primarysha
        else
            againstsha = snipsha(job.against.sha)
            string(primarysha, "_vs_", againstsha)
        end
        joinpath("by_hash", tag)
    end
end

reportdir(job::PkgEvalJob; kwargs...) = joinpath(reportdir(submission(job).config), "pkgeval", jobdirname(job; kwargs...))
tmpdir(job::PkgEvalJob) = joinpath(workdir(submission(job).config), "tmpresults")
tmplogdir(job::PkgEvalJob) = joinpath(tmpdir(job), "logs")
tmpdatadir(job::PkgEvalJob) = joinpath(tmpdir(job), "data")

########################
# PkgEvalJob Execution #
########################

# execute the tests of all packages specified by a PkgEvalJob on one or more Julia builds
function execute_tests!(job::PkgEvalJob, builds::Dict, flags::Dict, results::Dict)
    node = myid()
    cfg = submission(job).config

    # determine Julia versions to use
    julia_versions = Dict{String,VersionNumber}()
    for (whichbuild, build) in builds
        # obtain Julia version matching requested BuildRef
        julia = nothing
        if whichbuild == "primary" && submission(job).fromkind == :pr
            # if we're dealing with a PR, try the merge commit
            pr = submission(job).prnumber
            if pr !== nothing
                try
                    # NOTE: the merge head only exists in the upstream Julia repository,
                    #       and not in the repository where the pull request originated.
                    julia =
                        PkgEval.perform_julia_build("pull/$pr/merge", "JuliaLang/julia";
                                                       buildflags=flags[whichbuild])
                    nodelog(cfg, node, "Resolved $whichbuild build to Julia $julia (merge head of PR $pr)")
                catch err
                    isa(err, LibGit2.GitError) || rethrow()
                    # there might not be a merge commit (e.g. in the case of merge conflicts)
                end
            end
        end
        if julia === nothing
            # fall back to the last commit in the PR
            julia = if isempty(flags[whichbuild])
                PkgEval.obtain_julia_build(build.sha, build.repo)
            else
                PkgEval.perform_julia_build(build.sha, build.repo;
                                               buildflags=flags[whichbuild])
            end
            nodelog(cfg, node, "Resolved $whichbuild build to Julia $julia (commit $(build.sha) at $(build.repo))")
        end
        @assert !in(julia, values(julia_versions)) "Cannot compare identical Julia builds"
        julia_versions[whichbuild] = julia

        # get some version info
        mktempdir() do install
            PkgEval.prepare_julia(julia, install)
            try
                out = Pipe()
                PkgEval.run_sandboxed_julia(install, ```-e '
                        VERSION >= v"0.7.0-DEV.3630" && using InteractiveUtils
                        VERSION >= v"0.7.0-DEV.467" ? versioninfo(verbose=true) : versioninfo(true)
                        '
                    ```; stdout=out, stderr=out, stdin=devnull, interactive=false)
                close(out.in)
                build.vinfo = first(split(read(out, String), "Environment"))
            catch err
                build.vinfo = string("retrieving versioninfo() failed: ", sprint(showerror, err))
            end
        end
    end

    # determine packages to test
    pkgsel = Meta.parse(job.pkgsel)
    pkg_names = if pkgsel == :ALL
        String[]
    else
        eval(pkgsel)    # should be safe, it's a :vec of Strings
    end
    pkgs = PkgEval.read_pkgs(pkg_names)

    # run tests
    all_tests = withenv("CI" => true) do
        cpus = mycpus(submission(job).config)
        PkgEval.run(collect(values(julia_versions)), pkgs; ninstances=length(cpus))
    end

    # process the results for each Julia version separately
    for (whichbuild, build) in builds
        tests = all_tests[all_tests[!, :julia] .== julia_versions[whichbuild], :]
        results[whichbuild] = tests

        # write logs
        if cfg.bucket !== nothing
            for test in eachrow(tests)
                if !ismissing(test.log)
                    S3.put_object("$(cfg.bucket)/pkgeval/$(jobdirname(job))",
                                  "$(test.name).$(test.julia).log",
                                  Dict("body"       => test.log,
                                       "x-amz-acl"  => "public-read",
                                       "headers"    => Dict("Content-Type"=>"text/plain")))
                end
            end
        else
            cd(tmplogdir(job)) do
                for test in eachrow(tests)
                    isdir(test.name) || mkdir(test.name)
                    open(joinpath(test.name, "$(test.julia).log"), "w") do io
                        if !ismissing(test.log)
                            write(io, test.log)
                        end
                    end
                end
            end
        end

        # write data
        cd(tmpdatadir(job)) do
            # dataframe with test results
            let tests = copy(tests)
                # Feather can't handle non-primitive types, so stringify them
                for col in (:julia, :version, :status, :reason, :uuid)
                    tests[!, col] = map(repr, tests[!, col])
                end
                Feather.write("$(whichbuild).feather", tests)
            end

            # dict with build properties
            open("$(whichbuild).json", "w") do io
                json = Dict{String,Any}(
                    "build" => Dict(
                        "repo"  => build.repo,
                        "sha"   => build.sha,
                    )
                )
                JSON.print(io, json)
            end
        end
    end
end

function Base.run(job::PkgEvalJob)
    node = myid()
    cfg = submission(job).config

    # make temporary directory for job results
    # Why not create the job's actual report directory now instead? The answer is that
    # the commit SHA that currently describes the job might change if we find out that
    # we should use a merge commit instead. To avoid confusion, we dump all the results
    # to this temporary directory first, then move the data to the correct location
    # in the reporting phase.
    nodelog(cfg, node, "creating temporary directory for benchmark results")
    if isdir(tmpdir(job))
        nodelog(cfg, node, "...removing old temporary directory...")
        rm(tmpdir(job), recursive=true)
    end
    nodelog(cfg, node, "...creating $(tmpdir(job))...")
    mkdir(tmpdir(job))
    if cfg.bucket === nothing
        nodelog(cfg, node, "...creating $(tmplogdir(job))...")
        mkdir(tmplogdir(job))
    end
    nodelog(cfg, node, "...creating $(tmpdatadir(job))...")
    mkdir(tmpdatadir(job))

    # prepare PkgEval
    PkgEval.prepare_registry("General"; update=true)

    # instantiate the dictionary that will hold all of the info needed by `report`
    results = Dict{Any,Any}()

    if job.isdaily
        # get build from previous day
        try
            nodelog(cfg, node, "retrieving results from previous daily build")
            latest_dir = reportdir(job; latest=true)
            latest_db = joinpath(latest_dir, "db.json")
            if isfile(latest_db)
                latest = JSON.parsefile(latest_db)

                # NOTE: we don't actually use the results from the previous day, just the
                #       build properties, since packages upgrades might cause failures too.
                results["against_date"] = parse(Date, latest["date"])
                job.against = BuildRef(latest["build"]["repo"], latest["build"]["sha"])
            else
                nodelog(cfg, node, "didn't find previous daily build data")
            end
        catch err
            rethrow(NanosoldierError("encountered error when retrieving old daily build data", err))
        end
    end

    # refuse to test against an identical build
    if job.against !== nothing && job.against.sha == submission(job).build.sha && job.against_buildflags == job.buildflags
        nodelog(cfg, node, "refusing to compare identical builds, demoting to non-comparing evaluation")
        delete!(results, "against_date")
        job.against = nothing
    end

    # run tests
    builds = Dict("primary" => submission(job).build)
    flags = Dict("primary" => job.buildflags)
    if job.against !== nothing
        builds["against"] = job.against
        flags["against"] = job.against_buildflags
    end
    try
        nodelog(cfg, node, "running tests for $(summary(job))")
        execute_tests!(job, builds, flags, results)
        nodelog(cfg, node, "running tests for $(summary(job))")
    catch err
        results["error"] = NanosoldierError("failed to run tests", err)
        results["backtrace"] = catch_backtrace()
    end

    PkgEval.purge()

    # report results
    nodelog(cfg, node, "reporting results for $(summary(job))")
    report(job, results)
    nodelog(cfg, node, "completed $(summary(job))")

    return
end

########################
# PkgEvalJob Reporting #
########################

# report job results back to GitHub
function report(job::PkgEvalJob, results)
    node = myid()
    cfg = submission(job).config
    if haskey(results, "primary") && isempty(results["primary"])
        reply_status(job, "error", "no tests were executed")
        reply_comment(job, "[Your package evaluation job]($(submission(job).url)) has completed, " *
                      "but no tests were actually executed. Perhaps your package selection " *
                      "contains misspelled names? cc @$(cfg.admin)")
    else
        #  prepare report + data and push it to report repo
        target_url = ""
        try
            nodelog(cfg, node, "...generating report...")
            reportname = "report.md"
            open(joinpath(tmpdir(job), reportname), "w") do file
                printreport(file, job, results)
            end
            if job.isdaily && !haskey(results, "error")
                nodelog(cfg, node, "...generating database...")
                dbname = "db.json"
                open(joinpath(tmpdir(job), dbname), "w") do file
                    printdb(file, job, results)
                end
            end
            nodelog(cfg, node, "...tarring data...")
            cd(tmpdir(job)) do
                run(`tar -cf data.tar data`)
                run(`xz --compress -9 --extreme data.tar`)
                rm(tmpdatadir(job), recursive=true)
            end
            nodelog(cfg, node, "...moving $(tmpdir(job)) to $(reportdir(job))...")
            mkpath(reportdir(job))
            mv(tmpdir(job), reportdir(job); force=true)
            if job.isdaily
                latest = reportdir(job; latest=true)
                islink(latest) && rm(latest)
                symlink(datedirname(job.date), latest)
            end
            nodelog(cfg, node, "...pushing $(reportdir(job)) to GitHub...")
            target_url = upload_report_repo!(job, joinpath("pkgeval", jobdirname(job), reportname),
                                             "upload report for $(summary(job))")
        catch err
            rethrow(NanosoldierError("error when preparing/pushing to report repo", err))
        end

        if haskey(results, "error")
            # TODO: throw with backtrace?
            if haskey(results, "backtrace")
                @error("An exception occurred during job execution",
                       exception=(results["error"], results["backtrace"]))
            else
                @error("An exception occurred during job execution",
                       exception=results["error"])
            end
            err = results["error"]
            err.url = target_url
            throw(err)
        else
            # determine the job's final status
            state = results["has_issues"] ? "failure" : "success"
            if job.against !== nothing
                status = results["has_issues"] ? "possible new issues were detected" :
                                                 "no new issues were detected"
            else
                status = results["has_issues"] ? "possible issues were detected" :
                                                 "no issues were detected"
            end
            # reply with the job's final status
            reply_status(job, state, status, target_url)
            if isempty(target_url)
                comment = "[Your package evaluation job]($(submission(job).url)) has completed, but " *
                          "something went wrong when trying to upload the result data. cc @$(cfg.admin)"
            else
                comment = "[Your package evaluation job]($(submission(job).url)) has completed - " *
                          "$(status). A full report can be found [here]($(target_url)). cc @$(cfg.admin)"
            end
            reply_comment(job, comment)
        end
    end
end

# Markdown Report Generation #
#----------------------------#

function printreport(io::IO, job::PkgEvalJob, results)
    cfg = submission(job).config
    build = submission(job).build
    buildname = string(build.repo, SHA_SEPARATOR, build.sha)
    buildlink = "https://github.com/$(build.repo)/commit/$(build.sha)"
    joblink = "[$(buildname)]($(buildlink))"
    hasagainstbuild = job.against !== nothing

    # in contrast to BenchmarkJob, comparison jobs always have an against build, even daily
    # ones (so we don't need `iscomparisonjob`). in the case of a daily comparison job,
    # `results["against_date"]` is guaranteed to be set (so we don't need `hasprevdate`).

    if hasagainstbuild
        againstbuild = job.against
        againstname = string(againstbuild.repo, SHA_SEPARATOR, againstbuild.sha)
        againstlink = "https://github.com/$(againstbuild.repo)/commit/$(againstbuild.sha)"
        joblink = "$(joblink) vs [$(againstname)]($(againstlink))"
    end

    # print report preface + job properties #
    #---------------------------------------#

    println(io, """
                # Package Evaluation Report

                ## Job Properties

                *Commit(s):* $(joblink)

                *Triggered By:* [link]($(submission(job).url))

                *Package Selection:* `$(job.pkgsel)`
                """)

    if job.isdaily
        if hasagainstbuild
            dailystr = string(job.date, " vs ", results["against_date"])
        else
            dailystr = string(job.date)
        end
        println(io, """
                    *Daily Job:* $(dailystr)
                    """)
    end

    # if errors are found, end the report now #
    #-----------------------------------------#

    if haskey(results, "error")
        println(io, """
                    ## Error

                    The build could not finish due to an error:

                    ```""")

        Base.showerror(io, results["error"])
        if haskey(results, "backtrace")
            Base.show_backtrace(io, results["backtrace"])
        end
        println(io)

        println(io, """
                    ```

                    Check the logs folder in this directory for more detailed output.
                    """)
        return nothing
    end

    # print summary of tested packages #
    #----------------------------------#

    # we don't care about the distinction between failed and killed tests,
    # so lump them together
    for key in ("primary", "against", "previous")
        if haskey(results, key)
            df = results[key]
            df[df[!, :status] .== :kill, :status] .= :fail
        end
    end

    o = count(==(:ok),      results["primary"].status)
    s = count(==(:skip),    results["primary"].status)
    f = count(==(:fail),    results["primary"].status)
    x = nrow(results["primary"])

    println(io, """
                In total, $x packages were tested, out of which $o succeeded, $f failed and $s were skipped.
                """)

    println(io)

    # print result list #
    #-------------------#

    if hasagainstbuild
        package_results = join(results["primary"], results["against"],
                               on=:uuid, kind=:left, makeunique=true, indicator=:source)
    else
        package_results = results["primary"]
        package_results[!, :source] .= "left_only" # fake a left join
    end

    results["has_issues"] = false

    # report test results in groups based on the test status
    for (status, (verb, emoji)) in (:fail   => ("failed tests", ":heavy_multiplication_x:"),
                                    :ok     => ("passed tests", ":heavy_check_mark:"),
                                    :skip   => ("were skipped", ":heavy_minus_sign:"))
        # NOTE: no `groupby(package_results, :status)` because we can't impose ordering
        group = package_results[package_results[!, :status] .== status, :]
        sort!(group, :name)

        if !isempty(group)
            println(io, "## $emoji Packages that $verb\n")

            # report on a single test
            function reportrow(test)
                verstr(version) = ismissing(version) ? "" : " v$(version)"

                primary_log = if cfg.bucket !== nothing
                    "https://s3.amazonaws.com/$(cfg.bucket)/pkgeval/$(jobdirname(job))/$(test.name).$(test.julia).log"
                else
                    "logs/$(test.name)/$(test.julia).log"
                end
                print(io, "- [$(test.name)$(verstr(test.version))]($primary_log)")

                # "against" entries are suffixed with `_1` because of the join
                if test.source == "both"
                    against_log = if cfg.bucket !== nothing
                        "https://s3.amazonaws.com/$(cfg.bucket)/pkgeval/$(jobdirname(job))/$(test.name_1).$(test.julia_1).log"
                    else
                        "logs/$(test.name_1)/$(test.julia_1).log"
                    end
                    print(io, " vs. [$(test.name_1)$(verstr(test.version_1))]($against_log)")

                    print(io, " ($(PkgEval.statusses[test.status_1])")
                    if !ismissing(test.reason_1)
                        print(io, ", $(PkgEval.reasons[test.reason_1])")
                    end
                    print(io, ")")
                end

                println(io)
            end

            # report on a group of tests, prefixed with the reason
            function reportgroup(group)
                subgroups = groupby(group, :reason; skipmissing=true)
                for subgroup in subgroups
                    println(io, uppercasefirst(PkgEval.reasons[first(subgroup).reason]), ":")
                    println(io)
                    foreach(reportrow, eachrow(subgroup))
                    println(io)
                end

                # print tests without a reason separately, at the end
                subgroup = group[group[!, :reason] .=== missing, :]
                if !isempty(subgroup)
                    if length(subgroups) > 0
                        println(io, "Other:")
                        println(io)
                    end
                    foreach(reportrow, eachrow(subgroup))
                    println(io)
                end
            end

            if hasagainstbuild
                # first report on tests that changed status
                changed_tests = filter(test->test.source == "both" &&
                                             test.status != test.status_1, group)
                if !isempty(changed_tests)
                    println(io, "**$(nrow(changed_tests)) packages $verb only on the current version.**")
                    println(io)
                    reportgroup(changed_tests)

                    if status == :fail
                        results["has_issues"] |= true

                        # if this was an explicit "vs" build (i.e., not a daily comparison
                        # against a previous day), give the syntax to re-test failures.
                        if haskey(submission(job).kwargs, :vs)
                            vs = submission(job).kwargs[:vs]
                            println(io,  """
                                <details><summary>Click here for the Nanosoldier invocation to re-run these tests.</summary>
                                <p>

                                ```
                                @nanosoldier `runtests($(repr(changed_tests.name)), vs = $vs)`
                                ```

                                </p>
                                </details>
                                """)

                            println(io)
                        end
                    end
                end

                # now report the other ones
                unchanged_tests = filter(test->test.source == "left_only" ||
                                               test.status == test.status_1, group)
                if !isempty(unchanged_tests)
                    println(io, """
                        <details><summary><strong>$(nrow(unchanged_tests)) packages $verb on the previous version too.</strong></summary>
                        <p>
                        """)
                    unchanged_tests = copy(unchanged_tests)     # only report the
                    unchanged_tests[!, :source] .= "left_only"  # primary result
                    reportgroup(unchanged_tests)
                    println(io, """
                        </p>
                        </details>
                        """)
                end
            else
                # just report on all tests
                println(io, "$(nrow(group)) packages $verb.")
                println(io)
                reportgroup(group)

                if status == :fail
                    results["has_issues"] |= true
                end
            end

            println(io)
        end
    end

    # print build version info #
    #--------------------------#

    print(io, """
              ## Version Info

              #### Primary Build

              ```
              $(build.vinfo)
              ```
              """)

    if !isempty(job.buildflags)
        println(io, "Build flags: ", join(map(flag->"`$flag`", job.buildflags), ", "))
    end

    if hasagainstbuild
        println(io)
        print(io, """
                  #### Comparison Build

                  ```
                  $(job.against.vinfo)
                  ```
                  """)

        if !isempty(job.against_buildflags)
            println(io, "Build flags: ", join(map(flag->"`$flag`", job.against_buildflags), ", "))
        end
    end

    println(io, "<!-- Generated on $(now()) -->")

    return nothing
end

# JSON Database Generation #
#--------------------------#

function printdb(io::IO, job::PkgEvalJob, results)
    build = submission(job).build

    # build information
    json = Dict{String,Any}(
        "build" => Dict(
            "repo"  => build.repo,
            "sha"   => build.sha,
        ),
        "date" => job.date,
    )

    # test results
    tests = Dict()
    for test in eachrow(results["primary"])
        tests[test.uuid] = Dict(
            "julia"         => test.julia,
            "name"          => test.name,
            "version"       => test.version,
            "status"        => test.status,
            "reason"        => test.reason,
            "duration"      => test.duration,
        )
    end
    json["tests"] = tests

    JSON.print(io, json)

    return
end


text:
module ORBIT

using ..ParticleSwarm
using ..Sequencing # Includes Clifford, Waveforms, and Qubit
using ..Waveforms

export tuneGate

# Fitness calls here are of the form
# "How well does <waveform> implement <pulse> on <qubit>?"

# The overarching "ORBIT" call will include a qubit and a gate, and define the
# fitness anonymously through the below.
#function fitness(params::Vector{UInt16}, q::Qubit, p::Pulse)
#  pulselength = div(length(params), 3)
#  q[p] = ExactWaveform(params[1:pulselength], params[pulselength:2*pulselength],
#                       params[2:pulselength:3*pulselength], true)
#  pulseIndex = findfirst(SQPulse, p)
#
#  # Every "1" in the readout detracts from the fitness score.
#  1 - mean([sendSequenceAndReadout(q, benchmark1Qubit(
#                                   pulseIndex, orbitSequenceLength))
#        for _ in 1:orbitAveraging])
#end

# Same but for FloatWaveforms, where we are optimizing vectors of Float64
#function fitness(params::Vector{Float64}, q::Qubit, p::Pulse)
#  q[p] = params
#  pulseIndex = findfirst(SQPulse, p)
#  1 - mean([sendSequenceAndReadout(q, benchmark1Qubit(
#                                   pulseIndex, orbitSequenceLength))
#        for _ in 1:orbitAveraging])
#end

const numClifs = 80
const numSeq = 6
import InstrumentControl: IQSoftwareResponse, Stimulus, Response, Instrument, source, measure, sweep, status, Done
type rbStim <: Stimulus
    q::Qubit
    r::Readout
    p::Pulse
    ind::Integer
end

import ..Clifford: SQPulseNoZ, Ypi
function rbStim(q::Qubit, r::Readout, p::Pulse)
    ind = findfirst(SQPulseNoZ, p)
    rbStim(q, r, p, ind)
end

function source(x::rbStim, g::ExactWaveform)
    gprime = ExactWaveform(g.XYQ, 16382-g.XYI, UInt16[], false, true)
    x.q[x.p] = g
    x.q[Ypi] = gprime
    sendNormalSequence(x.q, benchmark1Qubit(numClifs, x.ind, false), x.r)
end

#const binary1 = 28.26*cis(-0.786)
#const binary0 = 27.62*cis(-0.513)
# 1 microsecond readout has different numbers
const binary0 = 49.4 - 26.1im
const binary1 = 41.3 - 40.2im
binaryOut(IQdata) = abs(IQdata-binary1) < abs(IQdata-binary0) ? 1 : 0
type zeroProb <: Response
    r::IQSoftwareResponse
end

function measure(x::zeroProb)
    m = measure(x.r)::Array{Complex{Float32},1}
    1-mean(binaryOut.(m))
end

function fitnessMap(candidates::Vector{ExactWaveform}, stim::rbStim, out::zeroProb)
    pool = length(candidates)
    theJob = sweep(out, (stim, vcat(fill(candidates, numSeq)...)))
    while status(theJob) != Done
        sleep(1)
    end
    bare = theJob.sweep.result.data
    return [mean(bare[a:pool:(numSeq*pool)]) for a in 1:pool]
end

toParams(x::ExactWaveform) = Int64[x.XYI; x.XYQ]
function fromParams(x::Vector{Int64})
    numpoints = div(length(x), 2)
    y = UInt16.(x)
    ExactWaveform(y[1:numpoints], y[(numpoints+1):(2*numpoints)], UInt16[], false, true)
end

# Generalize later.  For now, let's just get it working with AWG and
# a lot of assumptions.
function tuneGate(q::Qubit, p::Pulse, r::Readout, iq::IQSoftwareResponse)
  init = q[p]
  numpoints = length(init.XYI)

  rb = rbStim(q, r, p)
  zp = zeroProb(iq)
  fitness(pool) = fitnessMap(fromParams.(pool), rb, zp)
  lowBound = map(x -> max(x - 800, 0), toParams(init))
  hiBound = map(x -> min(x + 800, 16382), toParams(init))
  winner, score = PSO(fitness, lowBound, hiBound, candidate=toParams(init))
  q[p] = fromParams(winner)
  return score

end

end # End module ORBIT

ext: jl

max_stars_repo_name: PainterQubits/PulseController

max_issues_repo_name: PainterQubits/PulseController

max_forks_repo_name: PainterQubits/PulseController


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

@testset "graph.clique" begin
    nullg = plaingraph()
    @test issetequal(maximumclique(maximalcliques(nullg)), [])

    noedges = plaingraph(5, Tuple{Int,Int}[])
    @test length(maximumclique(maximalcliques(noedges))) == 1

    g1 = cyclegraph(5)
    @test length(maximumclique(maximalcliques(g1))) == 2

    g2 = plaingraph(5, [(1, 2), (2, 3), (3, 1), (3, 4), (4, 5)])
    @test issetequal(maximumclique(maximalcliques(g2)), [1, 2, 3])

    g3 = plaingraph(
        5, [(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4), (4, 5)])
    @test issetequal(maximumclique(maximalcliques(g3)), [1, 2, 3, 4])
end # graph.clique


text:
(AxisVector{T, A, SVector{1, T}} where {T})(
    a::Real,
    ::LocalGeometry,
) where {A} = AxisVector(A.instance, SVector(a))

# standard conversions
ContravariantVector(u::ContravariantVector, ::LocalGeometry) = u
CovariantVector(u::CovariantVector, ::LocalGeometry) = u
LocalVector(u::LocalVector, ::LocalGeometry) = u

ContravariantVector(
    u::LocalVector{T, I},
    local_geometry::LocalGeometry{I},
) where {T, I} = local_geometry.∂ξ∂x * u
LocalVector(
    u::ContravariantVector{T, I},
    local_geometry::LocalGeometry{I},
) where {T, I} = local_geometry.∂x∂ξ * u
CovariantVector(
    u::LocalVector{T, I},
    local_geometry::LocalGeometry{I},
) where {T, I} = local_geometry.∂x∂ξ' * u
LocalVector(
    u::CovariantVector{T, I},
    local_geometry::LocalGeometry{I},
) where {T, I} = local_geometry.∂ξ∂x' * u

# Converting to specific dimension types
(::Type{<:ContravariantVector{<:Any, I}})(
    u::ContravariantVector{<:Any, I},
    ::LocalGeometry{I},
) where {I} = u

(::Type{<:ContravariantVector{<:Any, I}})(
    u::ContravariantVector,
    ::LocalGeometry,
) where {I} = transform(ContravariantAxis{I}(), u)

(::Type{<:ContravariantVector{<:Any, I}})(
    u::AxisVector,
    local_geometry::LocalGeometry,
) where {I} =
    transform(ContravariantAxis{I}(), ContravariantVector(u, local_geometry))

(::Type{<:CovariantVector{<:Any, I}})(
    u::CovariantVector{<:Any, I},
    ::LocalGeometry{I},
) where {I} = u

(::Type{<:CovariantVector{<:Any, I}})(
    u::CovariantVector,
    ::LocalGeometry,
) where {I} = transform(CovariantAxis{I}(), u)

(::Type{<:CovariantVector{<:Any, I}})(
    u::AxisVector,
    local_geometry::LocalGeometry,
) where {I} = transform(CovariantAxis{I}(), CovariantVector(u, local_geometry))

(::Type{<:LocalVector{<:Any, I}})(
    u::LocalVector{<:Any, I},
    ::LocalGeometry{I},
) where {I} = u

(::Type{<:LocalVector{<:Any, I}})(u::LocalVector, ::LocalGeometry) where {I} =
    transform(LocalAxis{I}(), u)

(::Type{<:LocalVector{<:Any, I}})(
    u::AxisVector,
    local_geometry::LocalGeometry,
) where {I} = transform(LocalAxis{I}(), LocalVector(u, local_geometry))

# Generic N-axis conversion functions,
# Convert to specific local geometry dimension then convert vector type
LocalVector(u::CovariantVector, local_geometry::LocalGeometry{I}) where {I} =
    transform(LocalAxis{I}(), transform(CovariantAxis{I}(), u), local_geometry)

LocalVector(
    u::ContravariantVector,
    local_geometry::LocalGeometry{I},
) where {I} = transform(
    LocalAxis{I}(),
    transform(ContravariantAxis{I}(), u),
    local_geometry,
)

CovariantVector(u::LocalVector, local_geometry::LocalGeometry{I}) where {I} =
    transform(CovariantAxis{I}(), transform(LocalAxis{I}(), u), local_geometry)

CovariantVector(
    u::ContravariantVector,
    local_geometry::LocalGeometry{I},
) where {I} = transform(
    CovariantAxis{I}(),
    transform(ContravariantAxis{I}(), u),
    local_geometry,
)

ContravariantVector(
    u::LocalVector,
    local_geometry::LocalGeometry{I},
) where {I} = transform(
    ContravariantAxis{I}(),
    transform(LocalAxis{I}(), u),
    local_geometry,
)

ContravariantVector(
    u::CovariantVector,
    local_geometry::LocalGeometry{I},
) where {I} = transform(
    ContravariantAxis{I}(),
    transform(CovariantAxis{I}(), u),
    local_geometry,
)

# In order to make curls and cross products work in 2D, we define the 3rd
# dimension to be orthogonal to the exisiting dimensions, and have unit length
# (so covariant, contravariant and local axes are equal)
ContravariantVector(u::LocalVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(Contravariant3Axis(), components(u))
ContravariantVector(u::CovariantVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(Contravariant3Axis(), components(u))

CovariantVector(u::LocalVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(Covariant3Axis(), components(u))
CovariantVector(u::ContravariantVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(Covariant3Axis(), components(u))

LocalVector(u::CovariantVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(WAxis(), components(u))
LocalVector(u::ContravariantVector{<:Any, (3,)}, ::LocalGeometry{(1, 2)}) =
    AxisVector(WAxis(), components(u))

covariant1(u::AxisVector, local_geometry::LocalGeometry) =
    CovariantVector(u, local_geometry).u₁
covariant2(u::AxisVector, local_geometry::LocalGeometry) =
    CovariantVector(u, local_geometry).u₂
covariant3(u::AxisVector, local_geometry::LocalGeometry) =
    CovariantVector(u, local_geometry).u₃

contravariant1(u::AxisVector, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry).u¹
contravariant2(u::AxisVector, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry).u²
contravariant3(u::AxisVector, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry).u³

contravariant1(u::Axis2Tensor, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry)[1, :]
contravariant2(u::Axis2Tensor, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry)[2, :]
contravariant3(u::Axis2Tensor, local_geometry::LocalGeometry) =
    transform(Contravariant123Axis(), u, local_geometry)[3, :]

Jcontravariant3(u::AxisTensor, local_geometry::LocalGeometry) =
    local_geometry.J * contravariant3(u, local_geometry)

# required for curl-curl
covariant3(u::Contravariant3Vector, local_geometry::LocalGeometry{(1, 2)}) =
    contravariant3(u, local_geometry)

"""
    transform(axis, V[, local_geometry])

Transform the first axis of the vector or tensor `V` to `axis`.
"""
function transform end

"""
    project(axis, V[, local_geometry])

Project the first axis component of the vector or tensor `V` to `axis`
"""
function project end

for op in (:transform, :project)
    @eval begin
        # Covariant <-> Cartesian
        $op(
            ax::CartesianAxis,
            v::CovariantTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂ξ∂x' * $op(dual(axes(local_geometry.∂ξ∂x, 1)), v),
        )
        $op(
            ax::CovariantAxis,
            v::CartesianTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂x∂ξ' * $op(dual(axes(local_geometry.∂x∂ξ, 1)), v),
        )
        $op(ax::LocalAxis, v::CovariantTensor, local_geometry::LocalGeometry) =
            $op(
                ax,
                local_geometry.∂ξ∂x' *
                $op(dual(axes(local_geometry.∂ξ∂x, 1)), v),
            )
        $op(ax::CovariantAxis, v::LocalTensor, local_geometry::LocalGeometry) =
            $op(
                ax,
                local_geometry.∂x∂ξ' *
                $op(dual(axes(local_geometry.∂x∂ξ, 1)), v),
            )

        # Contravariant <-> Cartesian
        $op(
            ax::ContravariantAxis,
            v::CartesianTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂ξ∂x * $op(dual(axes(local_geometry.∂ξ∂x, 2)), v),
        )
        $op(
            ax::CartesianAxis,
            v::ContravariantTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂x∂ξ * $op(dual(axes(local_geometry.∂x∂ξ, 2)), v),
        )
        $op(
            ax::ContravariantAxis,
            v::LocalTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂ξ∂x * $op(dual(axes(local_geometry.∂ξ∂x, 2)), v),
        )

        $op(
            ax::LocalAxis,
            v::ContravariantTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂x∂ξ * $op(dual(axes(local_geometry.∂x∂ξ, 2)), v),
        )

        # Covariant <-> Contravariant
        $op(
            ax::ContravariantAxis,
            v::CovariantTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂ξ∂x *
            local_geometry.∂ξ∂x' *
            $op(dual(axes(local_geometry.∂ξ∂x, 1)), v),
        )
        $op(
            ax::CovariantAxis,
            v::ContravariantTensor,
            local_geometry::LocalGeometry,
        ) = $op(
            ax,
            local_geometry.∂x∂ξ' *
            local_geometry.∂x∂ξ *
            $op(dual(axes(local_geometry.∂x∂ξ, 2)), v),
        )

        $op(ato::CovariantAxis, v::CovariantTensor, ::LocalGeometry) =
            $op(ato, v)
        $op(ato::ContravariantAxis, v::ContravariantTensor, ::LocalGeometry) =
            $op(ato, v)
        $op(ato::CartesianAxis, v::CartesianTensor, ::LocalGeometry) =
            $op(ato, v)
        $op(ato::LocalAxis, v::LocalTensor, ::LocalGeometry) = $op(ato, v)
    end
end

"""
    divergence_result_type(V)

The return type when taking the divergence of a field of `V`.

Required for statically infering the result type of the divergence operation for `AxisVector` subtypes.
"""
@inline divergence_result_type(::Type{V}) where {V <: AxisVector} = eltype(V)

# this isn't quite right: it only is true when the Christoffel symbols are zero
@inline divergence_result_type(
    ::Type{Axis2Tensor{FT, Tuple{A1, A2}, S}},
) where {FT, A1, A2 <: LocalAxis, S <: StaticMatrix{S1, S2}} where {S1, S2} =
    AxisVector{FT, A2, SVector{S2, FT}}

"""
    gradient_result_type(Val(I), V)

The return type when taking the gradient along dimension `I` of a field `V`.

Required for statically infering the result type of the gradient operation for `AxisVector` subtypes.
    """
@inline function gradient_result_type(
    ::Val{I},
    ::Type{V},
) where {I, V <: Number}
    N = length(I)
    AxisVector{V, CovariantAxis{I}, SVector{N, V}}
end
@inline function gradient_result_type(
    ::Val{I},
    ::Type{V},
) where {I, V <: Geometry.AxisVector{T, A, SVector{N, T}}} where {T, A, N}
    M = length(I)
    Axis2Tensor{T, Tuple{CovariantAxis{I}, A}, SMatrix{M, N, T, M * N}}
end

"""
    curl_result_type(Val(I), Val(L), V)

The return type when taking the curl along dimensions `I` of a field of eltype `V`, defined on dimensions `L`

Required for statically infering the result type of the curl operation for `AxisVector` subtypes.
"""
@inline curl_result_type(
    ::Val{(1, 2)},
    ::Type{Covariant12Vector{FT}},
) where {FT} = Contravariant3Vector{FT}
@inline curl_result_type(
    ::Val{(1, 2)},
    ::Type{Covariant3Vector{FT}},
) where {FT} = Contravariant12Vector{FT}
@inline curl_result_type(::Val{(1,)}, ::Type{Covariant3Vector{FT}}) where {FT} =
    Contravariant2Vector{FT}
@inline curl_result_type(
    ::Val{(3,)},
    ::Type{Covariant12Vector{FT}},
) where {FT} = Contravariant12Vector{FT}
@inline curl_result_type(::Val{(3,)}, ::Type{Covariant1Vector{FT}}) where {FT} =
    Contravariant2Vector{FT}
@inline curl_result_type(::Val{(3,)}, ::Type{Covariant2Vector{FT}}) where {FT} =
    Contravariant1Vector{FT}

# these are only true in 2D
@inline curl_result_type(
    ::Val{(1, 2)},
    ::Type{Contravariant12Vector{FT}},
) where {FT} = Contravariant3Vector{FT}
@inline curl_result_type(::Val{(1, 2)}, ::Type{UVVector{FT}}) where {FT} =
    Contravariant3Vector{FT}
@inline curl_result_type(
    ::Val{(1, 2)},
    ::Type{Contravariant3Vector{FT}},
) where {FT} = Contravariant12Vector{FT}
@inline curl_result_type(::Val{(1, 2)}, ::Type{WVector{FT}}) where {FT} =
    Contravariant12Vector{FT}

_norm_sqr(x, local_geometry::LocalGeometry) =
    sum(x -> _norm_sqr(x, local_geometry), x)
_norm_sqr(x::Number, ::LocalGeometry) = LinearAlgebra.norm_sqr(x)
_norm_sqr(x::AbstractArray, ::LocalGeometry) = LinearAlgebra.norm_sqr(x)

function _norm_sqr(uᵢ::CovariantVector, local_geometry::LocalGeometry)
    LinearAlgebra.norm_sqr(LocalVector(uᵢ, local_geometry))
end

function _norm_sqr(uᵢ::ContravariantVector, local_geometry::LocalGeometry)
    LinearAlgebra.norm_sqr(LocalVector(uᵢ, local_geometry))
end

_norm_sqr(u::Contravariant2Vector, ::LocalGeometry{(1,)}) =
    LinearAlgebra.norm_sqr(u.u²)
_norm_sqr(u::Contravariant2Vector, ::LocalGeometry{(3,)}) =
    LinearAlgebra.norm_sqr(u.u²)
_norm_sqr(u::Contravariant2Vector, ::LocalGeometry{(1, 3)}) =
    LinearAlgebra.norm_sqr(u.u²)

_norm_sqr(u::Contravariant3Vector, ::LocalGeometry{(1,)}) =
    LinearAlgebra.norm_sqr(u.u³)
_norm_sqr(u::Contravariant3Vector, ::LocalGeometry{(1, 2)}) =
    LinearAlgebra.norm_sqr(u.u³)

_norm(u::AxisVector, local_geometry::LocalGeometry) =
    sqrt(_norm_sqr(u, local_geometry))

_cross(u::AxisVector, v::AxisVector, local_geometry::LocalGeometry) = _cross(
    ContravariantVector(u, local_geometry),
    ContravariantVector(v, local_geometry),
    local_geometry,
)

_cross(
    x::ContravariantVector,
    y::ContravariantVector,
    local_geometry::LocalGeometry,
) =
    local_geometry.J * Covariant123Vector(
        x.u² * y.u³ - x.u³ * y.u²,
        x.u³ * y.u¹ - x.u¹ * y.u³,
        x.u¹ * y.u² - x.u² * y.u¹,
    )
_cross(
    x::Contravariant12Vector,
    y::Contravariant12Vector,
    local_geometry::LocalGeometry,
) = local_geometry.J * Covariant3Vector(x.u¹ * y.u² - x.u² * y.u¹)
_cross(
    x::Contravariant2Vector,
    y::Contravariant1Vector,
    local_geometry::LocalGeometry,
) = local_geometry.J * Covariant3Vector(-x.u² * y.u¹)
_cross(
    x::Contravariant12Vector,
    y::Contravariant3Vector,
    local_geometry::LocalGeometry,
) = local_geometry.J * Covariant12Vector(x.u² * y.u³, -x.u¹ * y.u³)
_cross(
    x::Contravariant3Vector,
    y::Contravariant12Vector,
    local_geometry::LocalGeometry,
) = local_geometry.J * Covariant12Vector(-x.u³ * y.u², x.u³ * y.u¹)

_cross(
    x::Contravariant2Vector,
    y::Contravariant3Vector,
    local_geometry::LocalGeometry,
) = local_geometry.J * Covariant1Vector(x.u² * y.u³)

_cross(u::CartesianVector, v::CartesianVector, ::LocalGeometry) =
    LinearAlgebra.cross(u, v)
_cross(u::LocalVector, v::LocalVector, ::LocalGeometry) =
    LinearAlgebra.cross(u, v)


text:
immutable LocationIterator{w,T}
    range::UnitRange{Int}
    index::FMIndex{w,T}
end

length(iter::LocationIterator) = length(iter.range)
@inline start(iter::LocationIterator) = 1
@inline done(iter::LocationIterator, i) = i > length(iter)
@inline next(iter::LocationIterator, i) = sa_value(iter.range[i], iter.index) + 1, i + 1


text:
# This example plots msf plots of Lorenz system
import LyapunovExponents as LE
using LinearAlgebra
using PyPlot
pygui(true)
close("all")

# Lyapunov exponents calculation
ds = LE.Lorenz(ρ = 28)
P = diagm([1.0, 0.0, 0.0])
nsteps = Int(3e4)       # Number of steps
ntrsteps = Int(1e4)     # Number of transient steps
dt = 0.001              # Step size of the integration
ηr = collect(range(0, 20, length = 21))
Λr = map(ηi -> LE.msf(ds, ηi, P, nsteps = nsteps, ntrsteps = ntrsteps, dt = dt), ηr)

# Findin the threshold
ηt = ηr[findfirst(Λr .< 0)]

# Plots
figure()
plot(ηr, Λr, marker = "*")
axhline(0, color = "red")
xlabel(L"\eta")
ylabel(L"\Lambda")
tight_layout()


text:
module Projeto2Solvers

# standard lib
using LinearAlgebra
using ForwardDiff

# JSO packages
using NLPModels
using SolverTools
using JuMP
using Ipopt

include("uncsolver.jl")
include("newtoncombusca.jl")
include("l_bfgs_rcst.jl")
include("basicos.jl")
include("StrongWolfe.jl")
include("Newton_rc_bissec.jl")
include("Lbfgs_StrongWolfe.jl")

end # module


text:
struct graph_coloring
    graph::SimpleGraph
    num_cols::UInt8
end

Base.eltype(::Type{graph_coloring}) = Vector{UInt8}
Base.length(GC::graph_coloring) = GC.num_cols*((GC.num_cols-1)^(nv(GC.graph)-1))

function Base.iterate(GC::graph_coloring, prev::Vector{UInt8} = UInt8[])::Union{Nothing, Tuple{Vector{UInt8}, Vector{UInt8}}}

    if isempty(prev)
        prev = ones(UInt8, nv(GC.graph))
    end

    while true

        index = findfirst(x -> prev[x] < GC.num_cols, 1:length(prev))

        index === nothing && return nothing

        prev[index] += UInt8(1)
        for i in 1:(index-1)
            @inbounds prev[i] = UInt8(1)
        end

        #is_coloring(GC.graph, prev) && break
        all(e -> prev[src(e)] != prev[dst(e)], edges(GC.graph)) && break
    end

    return prev, prev #use copy(prev), prev in order to collect the elements
end

########

const Read = Dict("1"=>1,"2"=>2,"3"=>3,"4"=>4,"5"=>5,"6"=>6,"7"=>7,"8"=>8,"9"=>9,"a"=>10,
            "b"=>11,"c"=>12,"d"=>13,"e"=>14,"f"=>15,"g"=>16,"h"=>17,"i"=>18,"j"=>19,"k"=>20,
            "l"=>21,"m"=>22,"n"=>23,"o"=>24,"p"=>25,"q"=>26,"r"=>27,"s"=>28,"t"=>29,"u"=>30)

struct graph_coloring_from_file
    color_file_gz::CodecZlib.TranscodingStreams.TranscodingStream{GzipDecompressor, IOStream}
end

Base.eltype(::Type{graph_coloring_from_file}) = Tuple{Vector{UInt8}, Int64}
# function Base.length(GC::graph_coloring_from_file)
#     res = countlines(GC.color_file_gz) - 1
#     seekstart(GC.color_file_gz)
#     return res
# end
function Base.iterate(GC::graph_coloring_from_file, ::Int64 = 0)::Union{Nothing, Tuple{Tuple{Vector{UInt8}, Int64}, Int64}}

    st = readline(GC.color_file_gz)
    if st == "STOP"
        close(GC.color_file_gz)
        return nothing
    end

    st_s = split(st, ",")

    current_color = [UInt8(Read[x]) for x in split(st_s[1], "")]
    current_aut = parse(Int64, st_s[2])

    return (current_color, current_aut), 0
end

#=
"""
Compute the graph with a specific Prufer sequence.
"""=#
function PruferToGraph(prufer::Vector{UInt8})::SimpleGraph

    num_v = length(prufer)+2 #number of vertices
    g = Graph(num_v) #initialize a graph with specific number of vertices

    degree::Vector{Int64} = zeros(Int64,num_v)

    for i in 1:num_v-2
        degree[prufer[i]] += 1
    end

    for i in 1:(num_v-2)
        for j in 1:num_v
            if degree[j] == 0 #If j is not present in prufer set
                add_edge!(g, j, prufer[i]) #add edge
                degree[j] = -1 #Remove from Prufer
                degree[prufer[i]] -= 1 #low the degree of the vertex i
                break
            end
        end
    end

    # For the last element, we look at those vertices still with degree
    i = indexin(0, degree)[1]
    j = num_v + 1 - indexin(0, reverse(degree))[1]
    add_edge!(g, i, j)

    return g
end
#=
"""
Extract a Prufer sequence and a number from a string. The number is the number of automorphisms of the graph with that Prufer sequence.
"""=#
# function get_graph(str::String)::Tuple{SimpleGraph{Int64}, Int64}

#     s = split(str, ',')
#     g = PruferToGraph(UInt8[Read[string(s[1][i])] for i in 1:length(s[1])]) #get the graph
#     a = parse(Int64,s[2])

#     return (g, a)
# end
#=
"""
Return all the arrays of length `l` such that the sum of all elements of the array is `d`.
"""=#
function get_weights(l::Int64, d::Int64)::Vector{Vector{Int64}}

    return vcat(unique.(permutations.(partitions(d,l)))...)
end

struct graph_gen
    list_g::IOStream
    tot_number::Int64
    in_dir::Union{Vector{String},Nothing}
end

Base.eltype(::Type{graph_gen}) = Tuple{String, Int64, Bool}
Base.length(Gg::graph_gen) = Gg.tot_number

function Base.iterate(Gg::graph_gen, this_num::Int64 = 1)::Union{Nothing, Tuple{Tuple{String, Int64, Bool}, Int64}}

    if Gg.tot_number < this_num
        close(Gg.list_g)
        return nothing
    end

    this_num += 1

    st = split(readline(Gg.list_g), ",")

    has_file::Bool = true

    if Gg.in_dir !== nothing
        has_file = st[1]*"0.gz" in Gg.in_dir
    end

    top_aut::Int64 = parse(Int64, st[2])

    return (string(st[1]), top_aut, has_file), this_num
end


text:
@doc raw"""
    SlackEvaluator{Evaluator<:AbstractNLPEvaluator, T, VT} <: AbstractNLPEvaluator

Reformulate a problem with inequality constraints as an
equality constrained problem, by introducing a set of slack
variables.

### Description
A `SlackEvaluator` takes as input an original [`AbstractNLPEvaluator`](@ref),
subject to inequality constraints
```math
\begin{aligned}
       \min_{u \in \mathbb{R}^n} \quad & f(u)\\
\mathrm{s.t.} \quad & h^♭ ≤ h(u) ≤ h^♯,\\
                    & u^♭ ≤  u   ≤ u^♯.
\end{aligned}
```
The `SlackEvaluator` instance rewrites this problem with inequalities
as a new problem comprising only *equality constraints*, by introducing
$m$ slack variables $s_1, ⋯, s_m$. The new problem writes out
```math
\begin{aligned}
       \min_{u \in \mathbb{R}^n, s \in \mathbb{R}^m} \quad & f(u)\\
    \mathrm{s.t.} \quad & h(u) - s = 0 \\
                    & u^♭ ≤  u   ≤ u^♯, \\
                    & h^♭ ≤  s   ≤ h^♯.
\end{aligned}
```

### Attributes

- `inner::Evaluator`: original evaluator
- `s_min::VT`: stores lower bounds for slack variables
- `s_max::VT`: stores upper bounds for slack variables
- `nv::Int`: number of original variables
- `ns::Int`: number of slack variables

"""
mutable struct SlackEvaluator{Evaluator<:AbstractNLPEvaluator, T, VT} <: AbstractNLPEvaluator
    inner::Evaluator
    s_min::VT
    s_max::VT
    nv::Int
    ns::Int
end
function SlackEvaluator(nlp::AbstractNLPEvaluator)
    if !is_constrained(nlp)
        error("Input problem must have inequality constraints")
    end
    nv, ns = n_variables(nlp), n_constraints(nlp)
    s_min, s_max = bounds(nlp, Constraints())
    return SlackEvaluator{typeof(nlp), eltype(s_min), typeof(s_min)}(nlp, s_min, s_max, nv, ns)
end
function SlackEvaluator(
    datafile::String;
    device=KA.CPU(), options...
)
    nlp = ReducedSpaceEvaluator(datafile; device=device)
    return SlackEvaluator(nlp)
end

n_variables(nlp::SlackEvaluator) = nlp.nv + nlp.ns
n_constraints(nlp::SlackEvaluator) = n_constraints(nlp.inner)

constraints_type(::SlackEvaluator) = :equality
has_hessian(nlp::SlackEvaluator) = has_hessian(nlp.inner)

# Getters
get(nlp::SlackEvaluator, attr::AbstractNLPAttribute) = get(nlp.inner, attr)
get(nlp::SlackEvaluator, attr::AbstractVariable) = get(nlp.inner, attr)
get(nlp::SlackEvaluator, attr::PS.AbstractNetworkAttribute) = get(nlp.inner, attr)

# Setters
function setvalues!(nlp::SlackEvaluator, attr::PS.AbstractNetworkValues, values)
    setvalues!(nlp.inner, attr, values)
end

# Bounds
function bounds(nlp::SlackEvaluator, ::Variables)
    u♭, u♯ = bounds(nlp.inner, Variables())
    return [u♭; nlp.s_min], [u♯; nlp.s_max]
end
bounds(nlp::SlackEvaluator, ::Constraints) = (zeros(nlp.ns), zeros(nlp.ns))

function initial(nlp::SlackEvaluator)
    u0 = initial(nlp.inner)
    update!(nlp.inner, u0)
    cons = constraint(nlp.inner, u0)
    return [u0; -cons]
end

function update!(nlp::SlackEvaluator, w)
    u = @view w[1:nlp.nv]
    return update!(nlp.inner, u)
end

# f(x) = f₀(u)   , with x = (u, s)
function objective(nlp::SlackEvaluator, w)
    u = @view w[1:nlp.nv]
    return objective(nlp.inner, u)
end

# h(x) = h₀(u) - s
function constraint!(nlp::SlackEvaluator, cons, w)
    u = @view w[1:nlp.nv]
    s = @view w[nlp.nv+1:end]
    constraint!(nlp.inner, cons, u)
    cons .-= s
    return
end

# Gradient
# ∇f = [ ∇f₀ ; 0 ]
function gradient!(nlp::SlackEvaluator, grad, w)
    # w.r.t. u
    u = @view w[1:nlp.nv]
    gu = @view grad[1:nlp.nv]
    gradient!(nlp.inner, gu, u)
    # w.r.t. s
    gs = @view grad[nlp.nv+1:end]
    fill!(gs, 0.0)
    return
end

## Transpose Jacobian-vector product
# N.B.: constraints are specified as h(u) - s = 0
# J = [J₀  -I], so
# J' = [ J₀' ]
#      [ -I  ]
function jtprod!(nlp::SlackEvaluator, jv, w, v)
    # w.r.t. u
    u = @view w[1:nlp.nv]
    jvu = @view jv[1:nlp.nv]
    jtprod!(nlp.inner, jvu, u, v)
    # w.r.t. s
    jvs = @view jv[nlp.nv+1:end]
    jvs .= -v
    return
end

function jacobian!(nlp::SlackEvaluator, jac, w)
    fill!(jac, 0)
    u = @view w[1:nlp.nv]
    # w.r.t. u
    Jᵤ = @view jac[:, 1:nlp.nv]
    jacobian!(nlp.inner, Jᵤ, u)
    # w.r.t. s
    for i in 1:nlp.ns
        jac[i, i + nlp.nv] = -1.0
    end
end

function ojtprod!(nlp::SlackEvaluator, jv, w, σ, v)
    # w.r.t. u
    u = @view w[1:nlp.nv]
    jvu = @view jv[1:nlp.nv]
    ojtprod!(nlp.inner, jvu, u, σ, v)
    # w.r.t. s
    jvs = @view jv[nlp.nv+1:end]
    jvs .= -v
    return
end

# H = [ H₀   0 ]
#     [ 0    0 ]
function hessprod!(nlp::SlackEvaluator, hessvec, w, v)
    # w.r.t. u
    @views hessprod!(nlp.inner, hessvec[1:nlp.nv], w[1:nlp.nv], v[1:nlp.nv])
    # w.r.t. s
    hus = @view hessvec[nlp.nv+1:end]
    hus .= 0.0
    return
end

# J = [Jᵤ -I] , hence
# J' * J = [ Jᵤ' * Jᵤ    - Jᵤ']
#          [ - Jᵤ           I ]
function hessian_lagrangian_penalty_prod!(
    nlp::SlackEvaluator, hessvec, x, y, σ, v, ρ,
)
    @views begin
        u   = x[1:nlp.nv]
        vᵤ  = v[1:nlp.nv]
        vₛ  = v[nlp.nv+1:end]
        hvu = hessvec[1:nlp.nv]
        hvs = hessvec[nlp.nv+1:end]
    end
    fill!(hessvec, 0)
    u_buf = similar(u) ; fill!(u_buf, 0)
    y_buf = similar(y) ; fill!(y_buf, 0)
    # w.r.t. uu
    # ∇²L + ρ Jᵤ' * Jᵤ
    hessian_lagrangian_penalty_prod!(nlp.inner, hvu, u, y, σ, vᵤ, ρ)

    if !iszero(ρ)
        # w.r.t. us
        y_buf .= ρ .* vₛ
        # - Jᵤ' * vₛ
        jtprod!(nlp.inner, u_buf, u, y_buf)
        hvu .-=  u_buf
        # w.r.t. su
        jprod!(nlp.inner, y_buf, u, vᵤ)
        hvs .-= ρ .* y_buf
        # w.r.t. ss
        hvs .+= ρ .* vₛ
    end
    return
end

function hessian_lagrangian_penalty!(
    nlp::SlackEvaluator, H, x, y, σ, w,
)

    n = n_variables(nlp)
    @views begin
        u   = x[1:nlp.nv]
        Hᵤᵤ = H[1:nlp.nv, 1:nlp.nv]
        Hᵤᵥ = H[1:nlp.nv, 1+nlp.nv:n]
        Hᵥᵤ = H[1+nlp.nv:n, 1:nlp.nv]
        Hᵥᵥ = H[1+nlp.nv:n, 1+nlp.nv:n]
    end
    # w.r.t. uu
    # ∇²L + ρ Jᵤ' * Jᵤ
    hessian_lagrangian_penalty!(nlp.inner, Hᵤᵤ, u, y, σ, w)

    if !iszero(w)
        D = Diagonal(w)
        Jᵤ = Hᵥᵤ
        jacobian!(nlp.inner, Jᵤ, u)
        mul!(Hᵤᵥ, Jᵤ', -D)
        mul!(Hᵥᵤ, - D, Jᵤ)
        fill!(Hᵥᵥ, 0)
        @inbounds for i in 1:nlp.ns
            Hᵥᵥ[i, i] = w[i]
        end
    else
        fill!(Hᵤᵥ, 0)
        fill!(Hᵥᵤ, 0)
        fill!(Hᵥᵥ, 0)
    end
end

# TODO: return sparse sparsity pattern for bottom-left block
function hessian_structure(nlp::SlackEvaluator)
    n = n_variables(nlp)
    rows = Int[r for r in 1:n for c in 1:r]
    cols = Int[c for r in 1:n for c in 1:r]
    return rows, cols
end

function reset!(nlp::SlackEvaluator)
    reset!(nlp.inner)
end

# Utils function
function primal_infeasibility!(nlp::SlackEvaluator, cons, u)
    constraint!(nlp, cons, u) # Evaluate constraints
    return norm(cons, Inf)
end
function primal_infeasibility(nlp::SlackEvaluator, u)
    cons = similar(nlp.s_min) ; fill!(cons, 0)
    return primal_infeasibility!(nlp, cons, u)
end


text:
@testset "DependencyGraph" begin
  model1 = kendall()
  model1_rxnrxn = [
    [1, 2], # X --> X + X
    [1, 2], # X --> 0
    [1, 2]  # 0 --> X
  ]
  model1_spcrxn = [
    [1, 2]
  ]

  model2 = neg_autoreg() # gene P2_gene RNA P P2
  model2_rxnrxn = [
    [1, 2, 3, 6], # 1. gene + P2 --> P2_gene
    [1, 2, 3, 6], # 2. P2_gene --> gene + P2
    [4, 7],       # 3. gene --> gene + RNA  / gene does not change number
    [5, 8],       # 4. RNA --> RNA + P      / RNA does not change number
    [1, 5, 6, 8], # 5. P + P --> P2
    [1, 5, 6, 8], # 6. P2 --> P + P
    [4, 7],       # 7. RNA --> 0
    [5, 8]        # 8. P --> 0
  ]
  model2_spcrxn = [
    [1, 3],
    [2],
    [4, 7],
    [5, 8],
    [1, 6]
  ]

  @testset "DGVector{$(T)}" for T in [DGView, DGLazy]
    @testset "Birth-Death-Immigration Process" begin
      @testset "reaction-reaction" begin
        dg = rxnrxn_depgraph(T(), model1)

        for j in eachindex(model1_rxnrxn)
          @test collect(dependents(dg, j)) == model1_rxnrxn[j]
        end
      end

      @testset "species-reaction" begin
        dg = spcrxn_depgraph(T(), model1)

        for i in eachindex(model1_spcrxn)
          @test collect(dependents(dg, i)) == model1_spcrxn[i]
        end
      end
    end

    @testset "Negative Gene Autoregulation" begin
      @testset "reaction-reaction" begin
        dg = rxnrxn_depgraph(T(), model2)

        for j in eachindex(model2_rxnrxn)
          @test collect(dependents(dg, j)) == model2_rxnrxn[j]
        end
      end

      @testset "species-reaction" begin
        dg = spcrxn_depgraph(T(), model2)

        for i in eachindex(model2_spcrxn)
          @test collect(dependents(dg, i)) == model2_spcrxn[i]
        end
      end
    end
  end
end


text:
@testset "Solution objects" begin
    @testset "Graph 1" begin
        g = MetaDiGraph()
        add_vertex!(g, :name, "A")
        add_vertex!(g, :name, "B")
        add_vertex!(g, :name, "C")
        add_vertex!(g, :name, "D")
        add_vertex!(g, :name, "E")
        add_vertex!(g, :name, "F")
        add_edge!(g, 1, 3, :capacity, 2.0)
        add_edge!(g, 2, 3, :capacity, 2.0)
        add_edge!(g, 2, 6, :capacity, 2.0)
        add_edge!(g, 3, 4, :capacity, 20.0) # High capacity, shared.
        add_edge!(g, 4, 5, :capacity, 2.0)
        add_edge!(g, 1, 5, :capacity, 2.0) # Direct links.
        add_edge!(g, 4, 6, :capacity, 2.0)

        k = SimpleDiGraph(6)
        add_edge!(k, 1, 5)
        add_edge!(k, 2, 6)
        d1 = Edge(1, 5)
        d2 = Edge(2, 6)

        paths_d1 = [[Edge(1, 5)], [Edge(1, 3), Edge(3, 4), Edge(4, 5)]]
        paths_d2 = [[Edge(2, 6)], [Edge(2, 3), Edge(3, 4), Edge(4, 6)]]
        paths_d1_idx = [1, 2]
        paths_d2_idx = [3, 4]

        @testset "Flow solution" begin
            # Taken from an actual solution, with some numerical difficulties removed.
            # This is already a poor-quality solution, expect deviations like 0.001.
            sm = JuMP.Containers.DenseAxisArray{Float64}(undef, edges(k), edges(g))
            sm[d1, Edge(1, 3)] = 0.5003495146197712
            sm[d1, Edge(2, 3)] = 1.1217083702874072e-11
            sm[d1, Edge(2, 6)] = 1.1217083702874072e-11
            sm[d1, Edge(3, 4)] = 0.5003495146197712
            sm[d1, Edge(4, 5)] = 0.5003495146197712
            sm[d1, Edge(1, 5)] = 0.49965048538022855
            sm[d1, Edge(4, 6)] = 1.1217083702874072e-11
            sm[d2, Edge(1, 3)] = -1.1163257993389714e-11
            sm[d2, Edge(2, 3)] = 0.5003495146197711
            sm[d2, Edge(2, 6)] = 0.49965048538022877
            sm[d2, Edge(3, 4)] = 0.5003495146197711
            sm[d2, Edge(4, 5)] = -1.1163257993389714e-11
            sm[d2, Edge(1, 5)] = -1.1163257993389714e-11
            sm[d2, Edge(4, 6)] = 0.5003495146197711

            mt = ModelType(Load(), MinimumMaximum(), FlowFormulation(), false, Automatic(), NoUncertaintyHandling(), NoUncertainty())
            rd = RoutingData(g, k, nothing, mt)

            @testset "flow_routing_to_path" begin
                @test_throws ErrorException path_routing_to_flow(rd, sm)
                path_flows, edge_flows, paths = flow_routing_to_path(rd, sm)

                @test length(path_flows) == 2
                @test length(edge_flows) == 2
                @test length(paths) == 4

                @test [Edge(1, 5)] in paths
                @test [Edge(2, 6)] in paths
                @test [Edge(1, 3), Edge(3, 4), Edge(4, 5)] in paths
                @test [Edge(2, 3), Edge(3, 4), Edge(4, 6)] in paths

                for d in [d1, d2]
                    @test length(path_flows[d]) == 2
                    @test length(edge_flows[d]) == 4
                end

                # Paths with almost zero flow (below numerical accuracy) should not exist.
                for i in paths_d1_idx
                    @test path_flows[d1][i] ≈ 0.5 atol=1.e-3
                    @test_throws KeyError path_flows[d2][i]
                end
                for i in paths_d2_idx
                    @test_throws KeyError path_flows[d1][i]
                    @test path_flows[d2][i] ≈ 0.5 atol=1.e-3
                end
            end

            @testset "routing_to_matrix" begin
                r = Routing(rd, sm)

                mat = routing_to_matrix(r)
                @test mat[1, 1] ≈ 0.5 atol=1.e-3
                @test mat[1, 2] ≈ 0.5 atol=1.e-3
                @test mat[1, 3] ≈ 0.0 atol=1.e-3
                @test mat[1, 4] ≈ 0.0 atol=1.e-3
                @test mat[2, 1] ≈ 0.0 atol=1.e-3
                @test mat[2, 2] ≈ 0.0 atol=1.e-3
                @test mat[2, 3] ≈ 0.5 atol=1.e-3
                @test mat[2, 4] ≈ 0.5 atol=1.e-3
            end
        end

        @testset "Path solution" begin
            # Very similar solution to the flow one.
            sm = JuMP.Containers.SparseAxisArray{Float64, 2, Tuple{Edge{Int}, Int}}(Dict())
            for p in paths_d1_idx
                sm[d1, p] = 0.5003495146197712
                sm[d2, p] = 1.1217083702874072e-11
            end
            for p in paths_d2_idx
                sm[d1, p] = -1.1163257993389714e-11
                sm[d2, p] = 0.49965048538022855
            end

            mt = ModelType(Load(), MinimumMaximum(), PathFormulation(), false, Automatic(), NoUncertaintyHandling(), NoUncertainty())
            rd = RoutingData(g, k, nothing, mt)

            @testset "path_routing_to_flow" begin
                @test_throws ErrorException flow_routing_to_path(rd, sm)
                path_flows, edge_flows = path_routing_to_flow(rd, sm)

                @test length(path_flows) == 2
                @test length(edge_flows) == 2

                for d in [d1, d2]
                    @test length(path_flows[d]) == 2
                    @test length(edge_flows[d]) == 4
                end

                # Paths with almost zero flow (below numerical accuracy) should not exist.
                for i in paths_d1_idx
                    @test path_flows[d1][i] ≈ 0.5 atol=1.e-3
                    @test_throws KeyError path_flows[d2][i]
                end
                for i in paths_d2_idx
                    @test_throws KeyError path_flows[d1][i]
                    @test path_flows[d2][i] ≈ 0.5 atol=1.e-3
                end
            end

            @testset "routing_to_matrix" begin
                r = Routing(rd, sm)

                mat = routing_to_matrix(r)
                @test mat[1, 1] ≈ 0.5 atol=1.e-3
                @test mat[1, 2] ≈ 0.5 atol=1.e-3
                @test mat[1, 3] ≈ 0.0 atol=1.e-3
                @test mat[1, 4] ≈ 0.0 atol=1.e-3
                @test mat[2, 1] ≈ 0.0 atol=1.e-3
                @test mat[2, 2] ≈ 0.0 atol=1.e-3
                @test mat[2, 3] ≈ 0.5 atol=1.e-3
                @test mat[2, 4] ≈ 0.5 atol=1.e-3
            end
        end
    end

    @testset "Graph 2" begin # Same graph (and solution) as in compute.jl, for α-fairness.
        g = MetaDiGraph()
        add_vertex!(g, :name, "A")
        add_vertex!(g, :name, "B")
        add_vertex!(g, :name, "C")
        add_vertex!(g, :name, "D")
        add_edge!(g, 1, 2, :capacity, 2.0)
        add_edge!(g, 1, 3, :capacity, 2.0)
        add_edge!(g, 2, 4, :capacity, 2.0)
        add_edge!(g, 3, 4, :capacity, 2.0)
        add_edge!(g, 1, 4, :capacity, 3.0)

        k = SimpleDiGraph(4)
        add_edge!(k, 1, 4)
        d = Edge(1, 4)
        dm = Dict{Edge{Int}, Float64}(d => 4.0)

        paths_d = [[Edge(1, 2), Edge(2, 4)], [Edge(1, 3), Edge(3, 4)], [Edge(1, 4)]]

        @testset "Flow solution" begin
            # Taken from an actual solution, with some numerical difficulties removed.
            # This is already a poor-quality solution, expect deviations like 0.001.
            sm = JuMP.Containers.DenseAxisArray{Float64}(undef, edges(k), edges(g))
            sm[d, Edge(1, 4)] = 0.30413786305612406
            sm[d, Edge(1, 2)] = 0.347931068471938
            sm[d, Edge(1, 3)] = 0.347931068471938
            sm[d, Edge(2, 4)] = 0.347931068471938
            sm[d, Edge(3, 4)] = 0.347931068471938

            mt = ModelType(Load(), MinimumMaximum(), FlowFormulation(), false, Automatic(), NoUncertaintyHandling(), NoUncertainty())
            rd = RoutingData(g, k, nothing, mt)

            @testset "flow_routing_to_path" begin
                @test_throws ErrorException path_routing_to_flow(rd, sm)
                path_flows, edge_flows, paths = flow_routing_to_path(rd, sm)

                @test length(path_flows) == 1
                @test length(edge_flows) == 1
                @test length(paths) == 3

                for p in paths_d
                    @test p in paths
                end

                @test length(path_flows[d]) == 3
                @test length(edge_flows[d]) == 5

                @test path_flows[d][1] ≈ 0.30413786305612406
                @test path_flows[d][2] ≈ 0.347931068471938
                @test path_flows[d][3] ≈ 0.347931068471938
            end

            @testset "routing_to_matrix" begin
                r = Routing(rd, sm)

                mat = routing_to_matrix(r)
                @test mat[1, 1] ≈ 0.30413786305612406
                @test mat[1, 2] ≈ 0.347931068471938
                @test mat[1, 3] ≈ 0.347931068471938
            end
        end

        @testset "Path solution" begin
            # Very similar solution to the flow one.
            sm = JuMP.Containers.SparseAxisArray{Float64, 2, Tuple{Edge{Int}, Int}}(Dict())
            sm[d, 1] = 0.30413786305612406
            sm[d, 2] = 0.347931068471938
            sm[d, 3] = 0.347931068471938

            mt = ModelType(Load(), MinimumMaximum(), PathFormulation(), false, Automatic(), NoUncertaintyHandling(), NoUncertainty())
            rd = RoutingData(g, k, nothing, mt)

            @testset "path_routing_to_flow" begin
                @test_throws ErrorException flow_routing_to_path(rd, sm)
                path_flows, edge_flows = path_routing_to_flow(rd, sm)

                @test length(path_flows) == 1
                @test length(edge_flows) == 1

                @test length(path_flows[d]) == 3
                @test length(edge_flows[d]) == 5

                @test path_flows[d][1] ≈ 0.30413786305612406
                @test path_flows[d][2] ≈ 0.347931068471938
                @test path_flows[d][3] ≈ 0.347931068471938
            end

            @testset "routing_to_matrix" begin
                r = Routing(rd, sm)

                mat = routing_to_matrix(r)
                @test mat[1, 1] ≈ 0.30413786305612406
                @test mat[1, 2] ≈ 0.347931068471938
                @test mat[1, 3] ≈ 0.347931068471938
            end
        end
    end

    @testset "Graph 3" begin # Solution from the "basic" graph in the examples, the solution has a loop.
        g = MetaDiGraph()

        add_vertex!(g, :name, "A")
        add_vertex!(g, :name, "B")
        add_vertex!(g, :name, "C")
        add_vertex!(g, :name, "D")
        add_vertex!(g, :name, "E")
        add_vertex!(g, :name, "F")

        add_edge!(g, 1, 2, :capacity, 1.0)
        add_edge!(g, 2, 3, :capacity, 0.9)
        add_edge!(g, 2, 6, :capacity, 1.0)
        add_edge!(g, 3, 4, :capacity, 1.0)
        add_edge!(g, 3, 6, :capacity, 1.0)
        add_edge!(g, 5, 6, :capacity, 1.0)

        add_edge!(g, 2, 1, :capacity, 1.0) # Duplicate edges.
        add_edge!(g, 3, 2, :capacity, 0.9)
        add_edge!(g, 6, 2, :capacity, 1.0)
        add_edge!(g, 4, 3, :capacity, 1.0)
        add_edge!(g, 6, 3, :capacity, 1.0)
        add_edge!(g, 6, 5, :capacity, 1.0)

        paths_all = [
            [Edge(1, 2), Edge(2, 3)],
            [Edge(1, 2), Edge(2, 3), Edge(3, 6)],
            [Edge(1, 2), Edge(2, 6)],
            [Edge(4, 3), Edge(3, 2), Edge(2, 1)],
            [Edge(4, 3), Edge(3, 2), Edge(2, 6), Edge(6, 5)],
            [Edge(4, 3), Edge(3, 6), Edge(6, 5)],
            [Edge(5, 6), Edge(6, 3), Edge(3, 2)],
            [Edge(5, 6), Edge(6, 2)],
            [Edge(6, 2), Edge(2, 3), Edge(3, 4)],
            [Edge(6, 3), Edge(3, 4)]
        ]

        k = SimpleDiGraph(6)
        add_edge!(k, 1, 3)
        add_edge!(k, 1, 6)
        add_edge!(k, 4, 1)
        add_edge!(k, 4, 5)
        add_edge!(k, 5, 2)
        add_edge!(k, 6, 4)
        d1 = Edge(1, 3)
        d2 = Edge(1, 6)
        d3 = Edge(4, 1)
        d4 = Edge(4, 5)
        d5 = Edge(5, 2)
        d6 = Edge(6, 4)
        dm = Dict{Edge{Int}, Float64}(d1 => 0.5, d2 => 0.5, d3 => 0.5, d4 => 0.5, d5 => 0.5, d6 => 0.5)

        paths_d1 = [[Edge(1, 2), Edge(2, 3)]]
        paths_d2 = [[Edge(1, 2), Edge(2, 3), Edge(3, 6)], [Edge(1, 2), Edge(2, 6)]]
        paths_d3 = [[Edge(4, 3), Edge(3, 2), Edge(2, 1)]]
        paths_d4 = [[Edge(4, 3), Edge(3, 2), Edge(2, 6), Edge(6, 5)], [Edge(4, 3), Edge(3, 6), Edge(6, 5)]]
        paths_d5 = [[Edge(5, 6), Edge(6, 3), Edge(3, 2)], [Edge(5, 6), Edge(6, 2)]]
        paths_d6 = [[Edge(6, 2), Edge(2, 3), Edge(3, 4)], [Edge(6, 3), Edge(3, 4)]]

        @testset "Flow solution" begin
            # Taken from an actual solution, numerical difficulties not removed.
            sm = JuMP.Containers.DenseAxisArray{Float64}(undef, edges(k), edges(g))

            sm[d1, Edge(1, 2)] = 1.0
            sm[d1, Edge(2, 1)] = 0.0
            sm[d1, Edge(2, 3)] = 0.9999988912829403
            sm[d1, Edge(2, 6)] = 2.120611460099956e-6
            sm[d1, Edge(3, 2)] = 0.0
            sm[d1, Edge(3, 4)] = 0.0
            sm[d1, Edge(3, 6)] = 0.0
            sm[d1, Edge(4, 3)] = 0.0
            sm[d1, Edge(5, 6)] = 1.1094507307655902e-8
            sm[d1, Edge(6, 2)] = 1.0122885971469296e-6
            sm[d1, Edge(6, 3)] = 1.1087170597439012e-6
            sm[d1, Edge(6, 5)] = 1.1094507307655902e-8

            sm[d2, Edge(1, 2)] = 1.0
            sm[d2, Edge(2, 1)] = 0.0
            sm[d2, Edge(2, 3)] = 0.3177048579857707
            sm[d2, Edge(2, 6)] = 0.7726791282105232
            sm[d2, Edge(3, 2)] = 0.09038398659069616
            sm[d2, Edge(3, 4)] = 2.476831093836052e-7
            sm[d2, Edge(3, 6)] = 0.22732087178947685
            sm[d2, Edge(4, 3)] = 2.476831093836052e-7
            sm[d2, Edge(5, 6)] = 0.0
            sm[d2, Edge(6, 2)] = 0.0
            sm[d2, Edge(6, 3)] = 0.0
            sm[d2, Edge(6, 5)] = 0.0

            sm[d3, Edge(1, 2)] = 0.0
            sm[d3, Edge(2, 1)] = 1.0
            sm[d3, Edge(2, 3)] = 1.2410080677635113e-6
            sm[d3, Edge(2, 6)] = 8.844169471686035e-7
            sm[d3, Edge(3, 2)] = 0.9999987504657342
            sm[d3, Edge(3, 4)] = 0.0
            sm[d3, Edge(3, 6)] = 3.3707368767873062e-6
            sm[d3, Edge(4, 3)] = 1.0
            sm[d3, Edge(5, 6)] = 0.0
            sm[d3, Edge(6, 2)] = 3.374629109674882e-6
            sm[d3, Edge(6, 3)] = 8.805247142810213e-7
            sm[d3, Edge(6, 5)] = 0.0

            sm[d4, Edge(1, 2)] = 3.087618250488776e-7
            sm[d4, Edge(2, 1)] = 3.087618250488776e-7
            sm[d4, Edge(2, 3)] = 0.08955509114375726
            sm[d4, Edge(2, 6)] = 0.2343737049053851
            sm[d4, Edge(3, 2)] = 0.3239280973532021
            sm[d4, Edge(3, 4)] = 0.0
            sm[d4, Edge(3, 6)] = 0.76562769985842
            sm[d4, Edge(4, 3)] = 1.0
            sm[d4, Edge(5, 6)] = 0.0
            sm[d4, Edge(6, 2)] = 6.986957449624771e-7
            sm[d4, Edge(6, 3)] = 7.063964009965341e-7
            sm[d4, Edge(6, 5)] = 1.0

            sm[d5, Edge(1, 2)] = 0.0
            sm[d5, Edge(2, 1)] = 0.0
            sm[d5, Edge(2, 3)] = 0.0
            sm[d5, Edge(2, 6)] = 0.0
            sm[d5, Edge(3, 2)] = 0.2126523775599113
            sm[d5, Edge(3, 4)] = 0.0
            sm[d5, Edge(3, 6)] = 6.99368822994739e-7
            sm[d5, Edge(4, 3)] = 0.0
            sm[d5, Edge(5, 6)] = 1.0
            sm[d5, Edge(6, 2)] = 0.7873476224400887
            sm[d5, Edge(6, 3)] = 0.2126530769949393
            sm[d5, Edge(6, 5)] = 0.0

            sm[d6, Edge(1, 2)] = 0.0
            sm[d6, Edge(2, 1)] = 0.0
            sm[d6, Edge(2, 3)] = 0.3037925494376403
            sm[d6, Edge(2, 6)] = 0.0
            sm[d6, Edge(3, 2)] = 0.09114078583113316
            sm[d6, Edge(3, 4)] = 1.0
            sm[d6, Edge(3, 6)] = 0.0
            sm[d6, Edge(4, 3)] = 0.0
            sm[d6, Edge(5, 6)] = 0.0
            sm[d6, Edge(6, 2)] = 0.21265176354084359
            sm[d6, Edge(6, 3)] = 0.7873482364591564
            sm[d6, Edge(6, 5)] = 0.0

            mt = ModelType(Load(), MinimumMaximum(), FlowFormulation(), false, Automatic(), NoUncertaintyHandling(), NoUncertainty())
            rd = RoutingData(g, k, nothing, mt)

            @testset "flow_routing_to_path" begin
                @test_throws ErrorException path_routing_to_flow(rd, sm)
                path_flows, edge_flows, paths = flow_routing_to_path(rd, sm)

                @test length(path_flows) == 6
                @test length(edge_flows) == 6
                @test length(paths) == length(paths_all)

                for p in paths_all
                    @test p in paths
                end

                @test length(path_flows[d1]) == 1
                @test length(edge_flows[d1]) == 2
                @test length(path_flows[d2]) == 2
                @test length(edge_flows[d2]) == 4
                @test length(path_flows[d3]) == 1
                @test length(edge_flows[d3]) == 3
                @test length(path_flows[d4]) == 2
                @test length(edge_flows[d4]) == 5
                @test length(path_flows[d5]) == 2
                @test length(edge_flows[d5]) == 4
                @test length(path_flows[d6]) == 2
                @test length(edge_flows[d6]) == 4

                for p in keys(path_flows[d1])
                    @test paths[p] in paths_d1
                end
                for p in keys(path_flows[d2])
                    @test paths[p] in paths_d2
                end
                for p in keys(path_flows[d3])
                    @test paths[p] in paths_d3
                end
                for p in keys(path_flows[d4])
                    @test paths[p] in paths_d4
                end
                for p in keys(path_flows[d5])
                    @test paths[p] in paths_d5
                end
                for p in keys(path_flows[d6])
                    @test paths[p] in paths_d6
                end
            end

            @testset "routing_to_matrix" begin
                r = Routing(rd, sm)

                mat = routing_to_matrix(r)
                for i in 1:6
                    @test sum(mat[i, :]) ≈ 1.0 atol=1.0e-5
                end
            end
        end
    end
end


text:
using StaticArrays
using PMDPs
using Test

@testset "product.jl" begin

    p1 = PMDPs.Product(SA[true, false], 6)
    @test typeof(p1) <: PMDPs.Product

    p2 = PMDPs.Product(SA[true, true], 3)
    @test p1+p2 == [2,1]

    products = [PMDPs.Product(SA[true, false], 6),
                PMDPs.Product(SA[false, true], 8),
                PMDPs.Product(SA[true, true], 6)]
    @test PMDPs.are_unique(products)

    products_not_unique = [PMDPs.Product(SA[true, false], 6),
                PMDPs.Product(SA[false, true], 8),
                PMDPs.Product(SA[false, true], 6)]
    @test !PMDPs.are_unique(products_not_unique)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "ImageMagick"
version = v"6.9.10-12"
# Collection of sources required to build imagemagick
sources = [
    "https://github.com/ImageMagick/ImageMagick6.git" =>
    "d9d6f94ba8a40ca50d3b2fb748c1ec3014e335ef",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ImageMagick6/
./configure --prefix=$prefix --host=$target --disable-openmp --disable-install --disable-dependency-tracking --without-frozenpaths --without-perl
make -j${ncore}
make install

if [[ ${target} == "x86_64-apple-darwin14" ]]; then
    echo "-- Modifying link references for ImageMagick libraries"
    opts=""
    # for some reason libtiff and libpng don't need help?
    for XLIB in libMagick++-6.Q16.8 libMagickCore-6.Q16.6 libMagickWand-6.Q16.6 libjpeg.9 libz.1
    do
  	  opts="${opts} -change ${WORKSPACE}/destdir/lib/${XLIB}.dylib @rpath/${XLIB}.dylib"
    done
    for YLIB in libMagickWand-6.Q16.6 libMagickCore-6.Q16.6 libMagick++-6.Q16.8
    do
	  install_name_tool ${opts} ${prefix}/lib/${YLIB}.dylib
    done
fi
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libMagickWand", :libwand),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl",
    "https://github.com/SimonDanisch/LibpngBuilder/releases/download/v1.0.1/build_libpng.v1.6.31.jl",
    "https://github.com/SimonDanisch/LibJPEGBuilder/releases/download/v9b/build_libjpeg.v9.0.0-b.jl",
    "https://github.com/SimonDanisch/LibTIFFBuilder/releases/download/v1.0.0/build_libtiff.v4.0.9.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using DataFrames, CSV, Statistics

data = CSV.read("../data/temperatures.csv", copycols=true)
brisT = data.Brisbane
gcT = data.GoldCoast

sigB = std(brisT)
sigG = std(gcT)
covBG = cov(brisT, gcT)

meanVect = [mean(brisT) , mean(gcT)]
covMat = [sigB^2  covBG
          covBG   sigG^2]

outfile = open("../data/mvParams.jl","w")
write(outfile,"meanVect = $meanVect \ncovMat = $covMat")
close(outfile)
print(read("../data/mvParams.jl", String))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using XRayQuadrics

"""
    in_out(p::Particle, q::Quadric)

Computes entry and exit times for a `Particle` passing through a `Quadric` surface. Returns an Array of length 2, with the entry and exit times in increasing order, unless there is only one intersection with the surface, in which case the time is the first element. If there is no interaction, returns two zeros.
"""
function in_out(p::Particle, q::Quadric)
    # Project to homogeneous coordinate
    rh = [p.r0;1]   # AN ISSUE WITH THIS: if r = r0 + vΔt, and r0 and v have 1 as their     homogenous coord, get (1 + Δt) as the homogeneous coord for r. FIX
    vh = [p.v;0]    # REMOVE TIME TERM IN HOMEOGENEOUS COORDS?

    Q = q.Q
    Qr = Q[1:3,1:3]
    qd = Q[1:3,end]
    q0 = Q[end]

    # check for coplanar/axial alignment
    denom = vh'*Q*vh

    if all(Qr .== 0)                # plane
        if qd'*p.v == 0             # ray parallel to plane
            return [0, 0]
        else                        # ray intersects plane
            a = 2*qd
            return [(-q0 - a'*p.r0)/(a'*p.v), 0]
        end
    else                            # nonreducible quadric
        if denom == 0               # ray parallel to axis, not hyperboloid
            if qd'*p.v != 0         # paraboloid
                return [(p.r0'*Qr*p.r0 + 2*qd'*p.r0 + q0)/(2*qd'*p.v), 0]
            else                    # cylinder, but ray parallel to axis
                return [0, 0]
            end
        else                        # ray oblique to axis, or hyperboloid
            post = (sqrt(Complex((vh'*Q*rh)^2 - (vh'*Q*vh)*(rh'*Q*rh))))/denom
            pred = (-vh'*Q*rh)/denom
            if imag(post) == 0      # ray hits quadric
                post = real(post)
                return [min(pred - post, pred + post), max(pred - post, pred + post)]
            else                    # imaginary solution ⟹ no intersection with quadric
                return [0, 0]
            end
        end
    end

    # if denom == 0.0
    #     if all(Q[1:3, 1:3] .== 0)
    #         # a plane
    #         Qd = 2*Q[1:3,end]
    #         if Qd'*p.v == 0
    #             return (0,0)
    #         else
    #             return (0, (-Qd'*p.r0 - Q[end])/(Qd'*p.v))
    #         end
    #     else
    #         # colinear with nondegenerate quadric axis
    #         return (0,0)
    #     end
    # end

    # post = (sqrt(Complex((vh'*Q*rh)^2 - (vh'*Q*vh)*(rh'*Q*rh))))/denom
    # pred = (-vh'*Q*rh)/denom

    # if imag(post) == 0
    #     post = real(post)
    #     # sort roots in ascending order
    #     return extrema([pred - post, pred + post])
    # else
    #     # imaginary component implies no intersection
    #     return (0,0)
    # end
end

"""
    interacts(p::Particle, tq::TruncatedQuadric)

!!! warning
    FINISH IMPLEMENTING: should return a copy of `p` and a probability (of reflection or transmission, depending on )
"""
function interacts(p::Particle, tq::TruncatedQuadric)
    quadric = tq.q
    Qr = quadric.Q[1:3,1:3]
    qd = quadric.Q[1:3,end]
    q0 = quadric.Q[end]

    plane1 = tq.p[1]
    plane2 = tq.p[2]

    qtimes = in_out(p,quadric)
    p1time = in_out(p,plane1)[1]
    p2time = in_out(p,plane2)[1]

    nzeros = length(findall(qtimes .== 0))

    if nzeros == 2                  # parallel to quadric axis
        # project into both planes, check if within quadric at those positions
        # PROBLEM: HOW TO TELL IF IT'S *WITHIN* QUADRIC?

    elseif nzeros == 1              # quadric must be paraboloid
        # badorders:
        #   spp (anti-axis)
        #   pps (with axis)
        #   ppss    CANNOT APPEAR HERE
        #   sspp    CANNOT APPEAR HERE
        if tq.caps
            dir = sign(qd'*p.v)
            if dir == 1             # path parallel to axis

            elseif dir == 0
                error("Particle was supposed to travel parallel to Paraboloid axis.")
            else                    # path antiparallel to axis

            end
        else
            # psp is only passing case
        end
    else                            # hits walls of quadric twice

    end
end

"""
    cylinderentryexit(particle, cylinder)

Compute times at which a ray passes through a cylinder. The ray has initial position `particle.r0` and velocity `particle.v`. The cylinder has radius `cylinder.R`, height `cylinder.h`, a point one one of its caps `cylinder.c`, and a unit-length axis pointing from `cylinder.c` towards its other cap `cylinder.a`.

Two scalars are returned: the time at which the ray enters the cylinder and the time at which it leaves the cylinder. If the ray misses the cylinder, zeros are returned for both entry and exit time.
"""
function cylinderentryexit(particle, cylinder)
    r0 = particle.r0
    v = particle.v
    c = cylinder.c
    a = cylinder.a
    R = cylinder.R
    h = cylinder.h

    ttop = (c'*a - r0'*a)/(v'*a)
    tbottom = (c'*a + h - r0'*a)/(v'*a)

    topcapradius = norm(r0 + v*ttop - c)
    bottomcapradius = norm(r0 + v*tbottom - c - a*h)

    qa = v'*v - (v'*a)^2
    qb = 2*((r0'*v - c'*v) - ((v'*a)*(r0'*a - c'*a)))
    qc = norm(r0 - c)^2 - (r0'*a - c'*a)^2 - R^2

    # particle coaxial with cylinder:
    if qa == 0
        # check if it hits caps, if so return caps times
        if topcapradius ≤ R
            if bottomcapradius ≤ R
                # rectify to zero if particle starts within cylinder
                ttop = max(ttop, 0)
                tbottom = max(tbottom, 0)
                return (min(ttop, tbottom), max(ttop, tbottom))
            else
                error("particle supposedly parallel to cylinder axis")
            end
        else
            # outside cap radius, miss
            if bottomcapradius ≤ R
                error("particle supposedly parallel to cylinder axis")
            end

            return (0, 0)
        end

    # photon oblique to cylinder axis:
    else
        root1 = (-qb + sqrt(Complex(qb^2 - 4*qa*qc)))/(2*qa)
        root2 = (-qb - sqrt(Complex(qb^2 - 4*qa*qc)))/(2*qa)

        if imag(root1) == 0
            root1 = real(root1)
            root2 = real(root2)
            # hits infinite cylinder somewhere
            t1 = min(root1, root2)
            t2 = max(root1, root2)

            # list of all intersection times with cap planes and cylinder
            dt = [tbottom, ttop, t1, t2]
            # order intersection times
            I = sortperm(dt)

            # check if first two intersection times are with cylinder (t1, t2)
            if length(intersect(I[1:2], [3, 4])) == 2
                # miss
                return (0, 0)
            end
            # check if last two intersection times are with cylinder (t1, t2)
            if length(intersect(I[3:4], [3, 4])) == 2
                # miss
                return (0, 0)
            end

            # if you've made it this far, cylinder is hit somewhere
            tin = dt[I[2]]
            tout = dt[I[3]]

            tin = max(tin, 0)
            tout = max(tout, 0)
            return (tin, tout)

        # miss altogether (complex roots for intersection with cylinder)
        else
            return (0, 0)
        end
    end
end

"""
    lengthsinattenuator(photon, attenuator)

Get distances a `photon::Particle` travels within an `attenuator::PixelatedAttenuator` which has cylinders cut out of it.
"""
function lengthsinattenuator(photon, attenuator)
    topcylinderstimes = Array{Float64}(undef, 0, 2)
    # filter cylinders along path
    hits = falses(length(attenuator.holes[:]))
    for c = 1:length(attenuator.holes[:])
        if norm(cross(photon.v,attenuator.holes[c].a)) == 0
            # if cylinder parallel to ray
            hits[c] = norm(attenuator.holes[c].c - photon.r0 - photon.v*((photon.v'*(attenuator.holes[c].c - photon.r0))./(photon.v'*photon.v))) ≤ attenuator.holes[c].R
        else
            # if cylinder oblique to ray
            hits[c] = cross(photon.v, attenuator.holes[c].a)'*(photon.r0 - attenuator.holes[c].c)/norm(cross(photon.v, attenuator.holes[c].a)) ≤ attenuator.holes[c].R
        end

        # if this photon flies within this cylinder's radius, check for collision
        if hits[c]
            (tin, tout) = cylinderentryexit(photon, attenuator.holes[c])

            # make sure it actually hits cylinder
            if tout - tin > 0
                topcylinderstimes = cat(topcylinderstimes, [tin tout], dims=1)
            end
        end
    end

    # sort cylinders visited in order of visit time (increasing distance from r0 along v)
    sortedcylinderstimes = zeros(size(topcylinderstimes))
    if length(topcylinderstimes[:,1]) > 1
        topcylinderorder = sortperm(topcylinderstimes[:,1])
        sortedcylinderstimes = topcylinderstimes[topcylinderorder,:]
    else
        sortedcylinderstimes = topcylinderstimes
    end

    # times at which photon hits attenuator top surface and bottom surface
    tslabbottom = (attenuator.bottompoint'*attenuator.normal - photon.r0'*attenuator.normal)/(photon.v'*attenuator.normal)
    tslabtop = (attenuator.toppoint'*attenuator.normal - photon.r0'*attenuator.normal)/(photon.v'*attenuator.normal)

    tslabin = min(tslabbottom, tslabtop)
    tslabout = max(tslabbottom, tslabtop)

    if length(sortedcylinderstimes) > 0
        # times spent within attenuator material (complement of time spent in cylinders)
        timediffs = [sortedcylinderstimes[:,1]; tslabout] - [tslabin; sortedcylinderstimes[:,2]]
    else
        timediffs = tslabout - tslabin
    end

    # path lengths in material
    lengths = norm(photon.v).*timediffs

    return lengths
end

"""
    absorption_probability(photon, attenuator)

Return the probability of a photon::Particle being absorbed in attenuator::PixelatedAttenuator.
"""
function absorption_probability(photon, attenuator)
    lengths = lengthsinattenuator(photon, attenuator)
    if length(lengths) != 0
        transmissionlikelihoods = zeros(BigFloat, length(lengths))
        for i = 1:length(lengths)
            # interpolate attenuation coefficients from table data:
            attenuationcoeff = interpolateattenuation(photon.E, attenuator.massattenuation)

            # probability of passing through material:
            transmissionlikelihoods[i] = exp(BigFloat(-lengths[i]*attenuationcoeff))
        end

        # complement of total transmission likelihood is absorption likelihood
        absorptionlikelihood = 1 - prod(transmissionlikelihoods)

        return absorptionlikelihood
    else
        return 0.0
    end
end

"""
    transmission_probability(photon, attenuator)

Return the probability of a photon::Particle being transmitted through attenuator::PixelatedAttenuator.
"""
function transmission_probability(photon, attenuator)
    lengths = lengthsinattenuator(photon, attenuator)
    if length(lengths) != 0
        transmissionlikelihoods = zeros(BigFloat, length(lengths))
        for i = 1:length(lengths)
            # interpolate attenuation coefficients from table data:
            attenuationcoeff = interpolateattenuation(photon.E, attenuator.massattenuation)

            # probability of passing through material:
            transmissionlikelihoods[i] = exp(BigFloat(-lengths[i]*attenuationcoeff))
        end

        # get total transmission likelihood
        return prod(transmissionlikelihoods)

    else
        return 1.0
    end
end

"""
    batch_photons(photons, attenuator)

Compute absorption likelihood for a set of `photons` passing through `attenuator`. If multiple processes are available, computation will be parallelized. To execute in parallel, use the `Distributed` module in the calling context, add processes using `addprocs(n)`, and include this source with the `@everywhere` macro before `using` the module:

    @everywhere include("Attenuator3D.jl")
    using .Attenuator3D
"""
function batch_photons(photons, attenuator)
    # check if multiple processes available:
    if nprocs() > 1
        # need to factor jobsize by nworkers()
        spares = nworkers() - rem(length(photons), nworkers())
        # populate extra "dummy" photons to make array correct size
        appendix = Array{Particle}(undef, spares, 1)
        fill!(appendix, photons[1])
        inphotons = [photons[:]; appendix]

        # size of batch to run on each worker
        batchsize = length(inphotons[:]) ÷ nworkers()

        if rem(length(inphotons[:]), nworkers()) != 0
            error("Cannot distribute job to number of workers")
        end

        # shape array so each worker gets one column
        inphotons = reshape(inphotons, batchsize, nworkers())
        # make distributed array for result on all workers
        transmitlikelihood = dzeros(BigFloat, size(inphotons), workers(), [1, nworkers()])

        @sync @distributed for j in 1:nworkers()
            localtransmitlikelihood = localpart(transmitlikelihood)
            locali = DistributedArrays.localindices(transmitlikelihood)

            for i = 1:length(localtransmitlikelihood)
                localtransmitlikelihood[i] = transmission_probability(inphotons[i,locali[2][1]], attenuator)
            end
        end

        return Array(transmitlikelihood[1:(end - spares)])

    else
        transmitlikelihood = zeros(size(photons))

        @inbounds for i = 1:length(photons)
            transmitlikelihood[i] = transmission_probability(photons[i], attenuator)
        end
        return transmitlikelihood
    end
end


text:
@testset "MeshData" begin
  etable = (a=rand(100), b=rand(100))
  d = meshdata(CartesianGrid{T}(10,10), Dict(2 => etable))
  @test domain(d) == CartesianGrid{T}(10,10)
  @test values(d, 0) === nothing
  @test values(d, 1) === nothing
  @test values(d, 2) == etable

  vertices = P2[(0,0),(1,0),(1,1),(0,1)]
  elements = connect.([(1,2,3),(3,4,1)])
  d = meshdata(vertices, elements,
    Dict(
      0 => (temperature=[1.0,2.0,3.0,4.0], pressure=[4.0,3.0,2.0,1.0]),
      2 => (quality=["A","B"], state=[true,false])
    )
  )
  @test values(d, 0) == (temperature=[1.0,2.0,3.0,4.0], pressure=[4.0,3.0,2.0,1.0])
  @test values(d, 1) === nothing
  @test values(d, 2) == (quality=["A","B"], state=[true,false])

  etable = (a=rand(100), b=rand(100))
  d = meshdata(CartesianGrid{T}(10,10), etable = etable)
  @test domain(d) == CartesianGrid{T}(10,10)
  @test values(d, 0) === nothing
  @test values(d, 1) === nothing
  @test values(d, 2) == etable

  vertices = P2[(0,0),(1,0),(1,1),(0,1)]
  elements = connect.([(1,2,3),(3,4,1)])
  d = meshdata(vertices, elements,
    vtable = (temperature=[1.0,2.0,3.0,4.0], pressure=[4.0,3.0,2.0,1.0]),
    etable = (quality=["A","B"], state=[true,false])
  )
  @test domain(d) == SimpleMesh(vertices, elements)
  @test values(d, 0) == (temperature=[1.0,2.0,3.0,4.0], pressure=[4.0,3.0,2.0,1.0])
  @test values(d, 1) === nothing
  @test values(d, 2) == (quality=["A","B"], state=[true,false])
end


text:
module TestBooleans
    using Base.Test
    using DataArrays

    @assert NA | true == true
    @assert isna(NA | false)
    @assert isna(NA | NA)
    @assert true | NA == true
    @assert isna(false | NA)

    @assert isna(NA & true)
    @assert NA & false == false
    @assert isna(NA & NA)
    @assert isna(true & NA)
    @assert false & NA == false

    @assert any((@data [1, 2, NA]) .== 1) == true
    @assert any((@data [NA, 1, 2]) .== 1) == true
    @assert isna(any((@data [1, 2, NA]) .== 3))
    @assert any((@data [1, 2, 3] ).== 4) == false

    @assert isna(all((@data [1, 1, NA]) .== 1))
    @assert isna(all((@data [NA, 1, 1]) .== 1))
    @assert all((@data [1, 1, 1]) .== 1) == true
    @assert all((@data [1, 2, 1]) .== 1) == false
end


text:
function estimate_weights(profile::ValueProfileAnalysis, data::OneDimensionalTimeSeries, relative_weights)
    mean_variances = [mean_variance(data, field) for field in data.relevant_fields]
    weights = [1/σ for σ in mean_variances]

    if !isnothing(relative_weights)
        weights .*= relative_weights
    end

    return weights
end

function estimate_weights(profile::GradientProfileAnalysis, data::OneDimensionalTimeSeries, relative_weights)
    gradient_weight = profile.gradient_weight
    value_weight = profile.value_weight

    @warn "Dividing the gradient weight of profile by data.grid.Lz = $(data.grid.Lz))"
    gradient_weight = profile.gradient_weight = gradient_weight / data.grid.Lz

    fields = data.relevant_fields

    max_variances = [max_variance(data, field) for field in fields]
    #max_gradient_variances = [max_gradient_variance(data, field) for field in fields]

    weights = zeros(length(fields))
    for i in 1:length(fields)
        σ = max_variances[i]
        #ς = max_gradient_variances[i]
        weights[i] = 1/σ * (value_weight + gradient_weight / data.grid.Lz)
    end

    if relative_weights != nothing
        weights .*= relative_weights
    end

    max_variances = [max_variance(data, field) for field in fields]
    weights = [1/σ for σ in max_variances]

    if !isnothing(relative_weights)
        weights .*= relative_weights
    end

    return weights
end