text:
import Random

TextToArray(filename::String) = open(filename) do file
    texts = [line for line in eachline(file)]
end

TagBeginning(x::AbstractArray; tag = "<s>") = map(tokens -> vcat([tag], tokens), x)
TagEnding(x::AbstractArray; tag = "</s>") = map(tokens -> vcat(tokens, [tag]), x)

function TokenFrequencies(x::AbstractArray{String, 1}, sep = r"\s+")
    Lookup = Dict{String, Int32}()
    for i in 1:size(x, 1)
        for token in split(x[i], sep)
            if haskey(Lookup, token)
                Lookup[token] += Int32(1)
            else
                Lookup[token] = Int32(1)
            end
        end
    end
    return Lookup
end

GetMaxSequenceLength(x::AbstractArray{String, 1}, sep = r"\s+") = mapreduce(s -> length(split(s, sep)) |> Int32, max, x)
GetTopNTokens(TokenCounts::Dict{String, Int32}, N::Integer) = sort([[kv...] for kv in TokenCounts], lt=(x, y) -> isless(x[2], y[2]), rev = true)[1:N]

function TopNLookup(TokenCounts::Dict{String, Int32}, N::Integer)
    Lookup = Dict{String, Int32}([
        [item[2][1], item[1]]
        for item in enumerate(GetTopNTokens(TokenCounts, N))
    ])
    MaxVal = values(Lookup) |> maximum
    Lookup["<unk>"] = MaxVal + 1
    Lookup["<s>"] = MaxVal + 2
    Lookup["</s>"] = MaxVal + 3
    return Lookup
end

Tokenize(x::AbstractArray{String, 1}, sep = r"\s+") = map(s -> split(s, sep), x)
function TokensHash(s::String, TokenLookup::Dict{String, Int32}, sep = r"\s+")
    stringTokens = split(s, sep)
    out = similar(stringTokens, Int32)

    for item in enumerate(stringTokens)
        pos, token = item
        out[pos] = haskey(TokenLookup, token) ? TokenLookup[token] : TokenLookup["<unk>"]
    end
    return out
end
TokenizeLookup(x::AbstractArray{String, 1}, TokenLookup::Dict{String, Int32}, sep = r"\s+") = map(s -> TokensHash(s, TokenLookup, sep), x)
function DatasetSplit(x::AbstractArray{String, 1}; cvPercent = 0.1)
    if cvPercent >= 0.0 && cvPercent < 1.0
        x = Random.shuffle(x)
        cvSize = Int32(100 * cvPercent)
        return x[1:cvSize], x[cvSize:end]
    end
    return [], x

end


text:
module JuLearn

using JLD

function include_path(path,folder)
	files=readdir(joinpath(path,folder))
	map(file->include(joinpath(path,folder,file)),files)
end

export Layer

abstract Layer

if ENV["JuLearnBackend"]=="native"
	path=Pkg.dir("JuLearn","src","native")
elseif ENV["JuLearnBackend"]=="BLAS"
	path=Pkg.dir("JuLearn","src","BLAS")
elseif ENV["JuLearnBackend"]=="ArrayFire"
	using ArrayFire
	path=Pkg.dir("JuLearn","src","ArrayFire")
end

include_path(path,"layers")
include_path(path,"essentials")

include_path(path,"utils")
include_path(path,"activation")
include_path(path,"chain")
include_path(path,"loss")
include_path(path,"rnn")

end # module

ext: jl


text:
module NiGraphEmbedding

using NiLang, NiLang.AD
import ForwardDiff
using ForwardDiff: Dual
using Optim

export embedding_loss

include("native.jl")
include("reversible.jl")
include("hessian.jl")
include("Adam.jl")
include("train.jl")

end # module


text:
testset = [
    [
        "Umfpack",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPUmfpack,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "LapackCPU-BUNCHKAUFMAN",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPLapackCPU,
            lapackcpu_algorithm=MadNLPLapackCPU.BUNCHKAUFMAN,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "LapackCPU-LU",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPLapackCPU,
            lapackcpu_algorithm=MadNLPLapackCPU.LU,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "LapackCPU-QR",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPLapackCPU,
            lapackcpu_algorithm=MadNLPLapackCPU.QR,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "LapackCPU-CHOLESKY",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPLapackCPU,
            lapackcpu_algorithm=MadNLPLapackCPU.CHOLESKY,
            print_level=MadNLP.ERROR),
        ["infeasible", "lootsma", "eigmina"]
    ],
    [
        "Option: RELAX_BOUND",
        ()->MadNLP.Optimizer(
            fixed_variable_treatment=MadNLP.RELAX_BOUND,
            print_level=MadNLP.ERROR),
        [],
        true
    ],
    [
        "Option: AUGMENTED KKT SYSTEM",
        ()->MadNLP.Optimizer(
            kkt_system=MadNLP.SPARSE_UNREDUCED_KKT_SYSTEM,
            print_level=MadNLP.ERROR),
        ["infeasible","eigmina"] # numerical errors
    ],
    [
        "Option: INERTIA_FREE & AUGMENTED KKT SYSTEM",
        ()->MadNLP.Optimizer(
            inertia_correction_method=MadNLP.INERTIA_FREE,
            kkt_system=MadNLP.SPARSE_UNREDUCED_KKT_SYSTEM,
            print_level=MadNLP.ERROR),
        ["infeasible","eigmina"] # numerical errors
    ],
    [
        "Option: INERTIA_FREE",
        ()->MadNLP.Optimizer(
            inertia_correction_method=MadNLP.INERTIA_FREE,
            print_level=MadNLP.ERROR),
        []
    ],
]

@isdefined(MadNLPPardisoMKL) && push!(
    testset,
    [
        "PardisoMKL",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPPardisoMKL,
            print_level=MadNLP.ERROR),
        ["eigmina"]
    ]
)

for (name,optimizer_constructor,exclude) in testset
    test_madnlp(name,optimizer_constructor,exclude)
end


text:
module ArtifactWrappers

using Downloads
using Pkg.Artifacts
using DocStringExtensions: FIELDS

export ArtifactWrapper, ArtifactFile, get_data_folder

"""
    ArtifactFile

A single data file to be downloaded, containing both the url and the
name to use locally.

# Fields
$(FIELDS)
"""
Base.@kwdef struct ArtifactFile
    "URL pointing to data to be downloaded"
    url::AbstractString = ""
    "Local name used for downloaded online data"
    filename::AbstractString = ""
end

"""
    ArtifactWrapper

A set of data files to be downloaded, grouped by `data_name`. Example:

```
dataset = ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "MyDataSet",
    ArtifactFile[
        ArtifactFile(
            url="https://..../SomeNetCDF1.nc",
            filename="experiment1.nc",
        ),
        ArtifactFile(
            url="https://..../SomeNetCDF2.nc",
            filename="experiment2.nc",
        ),
    ]
)
```

# Fields
$(FIELDS)
"""
struct ArtifactWrapper
    "Directory to store artifact / data"
    artifact_dir::AbstractString
    "Locally running (not using CI)"
    local_run::Bool
    "Path to the used Artifacts.toml"
    artifact_toml::AbstractString
    "Unique name of dataset"
    data_name::AbstractString
    "Array of `ArtifactFile`'s, grouped by this dataset"
    artifact_files::Vector{ArtifactFile}
end
function ArtifactWrapper(artifact_dir, local_run, data_name, artifact_files)
    if !local_run
        artifact_dir = mktempdir(artifact_dir; prefix = "artifact_")
    end
    artifact_toml = joinpath(artifact_dir, "Artifacts.toml")
    return ArtifactWrapper(
        artifact_dir,
        local_run,
        artifact_toml,
        data_name,
        artifact_files,
    )
end

"""
    get_data_folder(art_wrap::ArtifactWrapper)

Get local folder of dataset defined in `art_wrap`.

Example:

```
dataset_path = get_data_folder(dataset)
```
"""
function get_data_folder(art_wrap::ArtifactWrapper)
    if !art_wrap.local_run
        # When running multiple jobs, create_artifact
        # has a race condition when creating/moving
        # files. So, when using CI, just download
        # the data files:
        filenames = [af.filename for af in art_wrap.artifact_files]
        urls = [af.url for af in art_wrap.artifact_files]
        for (url, filename) in zip(urls, filenames)
            Downloads.download(
                "$(url)",
                joinpath(art_wrap.artifact_dir, filename),
            )
        end
        return art_wrap.artifact_dir
    else
        # Query the `Artifacts.toml` file for the hash bound to the name
        # data_name (returns `nothing` if no such binding exists)
        data_hash = artifact_hash(art_wrap.data_name, art_wrap.artifact_toml)

        # If the name was not bound, or the hash it was bound to does not
        # exist, create it!
        if data_hash == nothing || !artifact_exists(data_hash)
            # create_artifact() returns the content-hash of the artifact
            # directory once we're finished creating it
            data_hash = create_artifact() do artifact_dir
                # We create the artifact by simply downloading a few files
                # into the new artifact directory
                filenames = [af.filename for af in art_wrap.artifact_files]
                urls = [af.url for af in art_wrap.artifact_files]
                for (url, filename) in zip(urls, filenames)
                    Downloads.download("$(url)", joinpath(artifact_dir, filename))
                end
            end

            # Now bind that hash within our `Artifacts.toml`. `force = true`
            # means that if it already exists, just overwrite with the new
            # content-hash.  Unless the source files change, we do not expect
            # the content hash to change, so this should not cause
            # unnecessary version control churn.
            bind_artifact!(
                art_wrap.artifact_toml,
                art_wrap.data_name,
                data_hash,
                force = true,
            )
        end

        # Get the path of the dataset, either newly created or
        # previously generated. This should be something like:
        # `~/.julia/artifacts/dbd04e28be047a54fbe9bf67e934be5b5e0d357a`
        dataset_path = artifact_path(data_hash)
        return dataset_path
    end
end

end # module


text:
using Boltzmann
using Distributions
using Test
using MLDatasets.MNIST

function _get_dataset(X)
    X = round.(X ./ (maximum(X) - minimum(X)))
    X_c = copy(X)   # the corrupted images

    # Flip 10% of bits in X_c
    flip_count = round(Int, size(X_c, 1) * 0.1)
    for i in 1:size(X_c, 2)
        flips = map(i -> rand(1:size(X_c, 1)), 1:flip_count)
        for j in flips
            item = Bool(X_c[j,i])
            X_c[j,i] = item ? 0.0 : 1.0
        end
    end

    return vcat(X, X_c)
end

if :integration in TEST_GROUPS
    @testset "Conditional RBM Integration" begin
        model = ConditionalRBM(Bernoulli, Bernoulli, 250, 150; steps=3)
        Boltzmann.test(model; input_size=1000, debug=true)
    end
end

if :acceptance in TEST_GROUPS
    @testset "Conditional RBM Acceptance" begin
        X = MNIST.convert2features(traindata()[1])[:, 1:10000]
        train_X = _get_dataset(X)
	    input_size = round(Int, size(train_X, 1) / 2)

        model = ConditionalRBM(Bernoulli, Bernoulli, input_size, 500; steps=1)
        fit(model, train_X; n_epochs=10, n_gibbs=5)

        test_X = _get_dataset(X)
        corrupt_start = input_size + 1
        forecast = predict(model, test_X[corrupt_start:end,:]; n_gibbs=20)

        mae = mean(test_X[1:corrupt_start-1,:] .!= forecast)

        # At a bare minimum the prediction should have an mae
        # less than the amount of noise added to the condition
        # values. In reality, this should be around ~0.05
        # and the likelihood should be move from >-60 to <-31
        println("MAE = $mae")
        @test mae < 0.1
    end
end

if :benchmark in TEST_GROUPS
    @testset "Conditional RBM Benchmark" begin
        n_vis = 250
        n_cond = 750
        n_hid = 200

        suite = BenchmarkGroup()

        for T in [Float32, Float64]
            model = ConditionalRBM(T, Bernoulli, Bernoulli, n_vis, n_hid, n_cond)
            suite = benchmark!(model, suite; input_size=1000, debug=true)
        end

        tune!(suite)
        results = run(suite, verbose=true, seconds=10)
    end
end


text:
using GraphMakie
using Documenter
using Literate
using CairoMakie

# preload the deps from the examples to supress precompilation output in docs
using JSServe
using NetworkDynamics
using LayeredLayouts
using Graphs
using PkgDeps

DocMeta.setdocmeta!(GraphMakie, :DocTestSetup, :(using GraphMakie); recursive=true)

# generate examples
examples = [
    joinpath(@__DIR__, "examples", "plots.jl"),
    joinpath(@__DIR__, "examples", "interactions.jl"),
    joinpath(@__DIR__, "examples", "depgraph.jl"),
    joinpath(@__DIR__, "examples", "truss.jl"),
    joinpath(@__DIR__, "examples", "syntaxtree.jl"),
]
OUTPUT = joinpath(@__DIR__, "src", "generated")
isdir(OUTPUT) && rm(OUTPUT, recursive=true)
mkpath(OUTPUT)

for ex in examples
    Literate.markdown(ex, OUTPUT)
end

makedocs(; modules=[GraphMakie], authors="Simon Danisch, Hans Würfel",
         repo="https://github.com/JuliaPlots/GraphMakie.jl/blob/{commit}{path}#{line}",
         sitename="GraphMakie.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://JuliaPlots.github.io/GraphMakie.jl", assets=String[]),
         pages=["Home" => "index.md",
                "Examples" => [
                    "Feature Walkthrough" => "generated/plots.md",
                    "Interaction Examples" => "generated/interactions.md",
                    "Dependency Graph" => "generated/depgraph.md",
                    "Stress on Truss" => "generated/truss.md",
                    "Julia AST" => "generated/syntaxtree.md",
                ]
                ])

# if gh_pages branch gets to big, check out
# https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#gh-pages-Branch
deploydocs(;repo="github.com/JuliaPlots/GraphMakie.jl",
           push_preview=true)


text:
include("./utils.jl")

using .Utils
using Knet: progress!, minibatch, adam, relu

EPOCH = 1
BATCHSIZE = 64
SENTENCELEN = 300
DROPOUT = 0.5
LR = 0.001
TRUNCATE = true

EMBEDDINGDIM = 100
HIDDENDIM = 20
OUTPUTDIM = 2

CNN_W1 = 5
CNN_W2 = 5
CNN_CX = 1
CNN_CY = 20

MODELTYPE = 1 # 0: MLP, 1: RNN, 2: CNN

if MODELTYPE == 0
    println("Training MLP...")
elseif MODELTYPE == 1
    println("Training RNN...")
else
    println("Training CNN...")
end

TRAIN_POS_PATH = joinpath(@__DIR__, "aclImdb/train/pos/")
TRAIN_NEG_PATH = joinpath(@__DIR__, "aclImdb/train/neg/")
TEST_POS_PATH = joinpath(@__DIR__, "aclImdb/test/pos/")
TEST_NEG_PATH = joinpath(@__DIR__, "aclImdb/test/neg/")

println("Read Training Data...")

xtrn_pos, ytrn_pos = readfolder(TRAIN_POS_PATH; class=1)
xtrn_neg, ytrn_neg = readfolder(TRAIN_NEG_PATH; class=2)

xtrn = vcat(xtrn_pos, xtrn_neg)
ytrn = vcat(ytrn_pos, ytrn_neg)

println("Read Test Data...")

xtst_pos, ytst_pos = readfolder(TEST_POS_PATH; class=1)
xtst_neg, ytst_neg = readfolder(TEST_NEG_PATH; class=2)

xtst = vcat(xtst_pos, xtst_neg)
ytst = vcat(ytst_pos, ytst_neg)

println("Shuffle...")

xtrn, ytrn = shuffledata(xtrn, ytrn)

println("Build vocabulary...")

voc = vocabulary(frequencies(skipwords(map(doc->map(lowercase, doc), xtrn))), 10)
len_voc = length(voc)

println("Doc2ids and add padding...")

if MODELTYPE == 0
    pos = "post"
    TRUNCATE = false
else
    pos = "pre"
end

xtrn_ids = addpadding!(doc2ids(xtrn, voc), pos=pos)
xtst_ids = addpadding!(doc2ids(xtst, voc), pos=pos)

if TRUNCATE
    println("Truncate sentences...")

    xtrn_ids = truncate!(xtrn_ids, SENTENCELEN, pos=pos)
    xtst_ids = truncate!(xtst_ids, SENTENCELEN, pos=pos)
end

dtrn = minibatch(xtrn_ids, ytrn, BATCHSIZE; shuffle=true)
dtst = minibatch(xtst_ids, ytst, BATCHSIZE; shuffle=true)

println("Build model...")

if MODELTYPE == 0
    model = MLPClassifier(len_voc, EMBEDDINGDIM, HIDDENDIM, OUTPUTDIM, relu; pdrop=0.5)
elseif MODELTYPE == 1
    model = RNNClassifier(len_voc, EMBEDDINGDIM, HIDDENDIM, OUTPUTDIM; pdrop=DROPOUT)
elseif MODELTYPE == 2
    dtrn = minibatch(xtrn_ids, ytrn, BATCHSIZE; shuffle=true)
    dtst = minibatch(xtst_ids, ytst, BATCHSIZE; shuffle=true)

    HIDDENDIM = Int((EMBEDDINGDIM - 4) / 2 * (SENTENCELEN - 4) / 2 * CNN_CY)
    model = CNNClassifier(len_voc, EMBEDDINGDIM, CNN_W1, CNN_W2,
                          CNN_CX, CNN_CY, HIDDENDIM, OUTPUTDIM,
                          relu; pdrop=DROPOUT)
else
    model = MLPClassifier(len_voc, EMBEDDINGDIM, HIDDENDIM, OUTPUTDIM, relu; pdrop=DROPOUT)
end

println("Training...")
progress!(adam(model, repeat(dtrn, EPOCH), lr=LR))

println("Train Accuracy:")
println(accuracy(model, dtrn))

println("Test Accuracy:")
println(accuracy(model, dtst))

println("Finish!")

max_stars_repo_name: alperkesen/imdb-sentiment-analysis

max_issues_repo_name: alperkesen/imdb-sentiment-analysis

max_forks_repo_name: alperkesen/imdb-sentiment-analysis


text:
# Core abstractions for the package

@unit ppm "ppm" PartsPerMillion 1//1000000 false

Unitful.register(MagneticResonanceSignals)

"""
4.7 is a chosen nominal chemical shift of water relative to TMS.  Note that
this is only approximate because the true shift can vary based on pH,
temperature, etc and you should calibrate your spectral axes against known
metabolites.
"""
const _water_tms_offset = 4.7u"ppm"

struct MRMetadata
    protocol_name    # Name of sequence parameter set
    sequence_name    # Name of MR sequence code
    software_version # Version of spectrometer software
    device_id        # Device identifier: (manufacturer, model_name, serial)
    ref_epoch        # Epoch of localizer image
    frequency        # Spectrometer frequency (TODO: or frequencies for MNO?)
end

"""
Extract standard MR metadata about an experiment
"""
function standard_metadata
end

"""
    mr_load(data)

High level function for loading MR `data` and recognizing which experiment was
run. The returned object describes the experimental conditions and acquisition
schedule which can be used during data processing.  `data` may be a file name,
`IO` stream, or a data structure describing raw acquisition data, for example,
Siemens Twix.

If the experiment is not recognized, throw an exception.
"""
mr_load(filename::String, args...; kws...) = open(io->mr_load(io, args...; kws...), filename)

# TODO: A magic number system to figure out what to load...
mr_load(io::IO, args...; kws...) = mr_load(load_twix(io, args...; kws...))

"""
    sampledata(expt, index; downsample=1)

Return the acquired data from `expt` at a given acqusition `index`.  If
`downsample>1`, the data will be subsampled by the given rate by truncating the
tails of the signal in the Fourier spectral domain. This has the effect of
removing noise by filtering away irrelevant high and low frequency components.
"""
function sampledata
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
module NeuralVerification

using JuMP

using GLPK, SCS # SCS only needed for Certify
using PicoSAT # needed for Planet
using LazySets, LazySets.Approximations
using Polyhedra, CDDLib

using LinearAlgebra
using Parameters
using Interpolations # only for PiecewiseLinear

import LazySets: dim, HalfSpace # necessary to avoid conflict with Polyhedra

# only for priority optimization
#import DataStructures: PriorityQueue, Queue, enqueue!, dequeue!

using Requires

abstract type Solver end

# For optimization methods:
import JuMP.MOI.OPTIMAL, JuMP.MOI.INFEASIBLE
JuMP.Model(solver::Solver) = Model(solver.optimizer)
# define a `value` function that recurses so that value(vector) and
# value(VecOfVec) works cleanly. This is only so the code looks nice.
value(var::JuMP.AbstractJuMPScalar) = JuMP.value(var)
value(vars::Vector) = value.(vars)
value(val) = val

include("utils/activation.jl")
include("utils/network.jl")
include("utils/problem.jl")
include("utils/util.jl")

function __init__()
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("utils/flux.jl")
end

export
    Solver,
    Network,
    AbstractActivation,
    # NOTE: not sure if exporting these is a good idea as far as namespace conflicts go:
    # ReLU,
    # Max,
    # Id,
    # Sigmoid,
    # Tanh,
    GeneralAct,
    PiecewiseLinear,
    Problem,
    Result,
    BasicResult,
    CounterExampleResult,
    AdversarialResult,
    ReachabilityResult,
    read_nnet,
    write_nnet,
    solve,
    forward_network,
    check_inclusion

solve(m::Model; kwargs...) = JuMP.solve(m; kwargs...)
export solve

# TODO: consider creating sub-modules for each of these.
include("optimization/utils/constraints.jl")
include("optimization/utils/objectives.jl")
include("optimization/utils/variables.jl")
include("optimization/nsVerify.jl")
include("optimization/convDual.jl")
include("optimization/duality.jl")
include("optimization/certify.jl")
include("optimization/iLP.jl")
include("optimization/mipVerify.jl")
include("optimization/bab.jl")
include("optimization/sherlock.jl")
include("optimization/reluplex.jl")
include("optimization/planet.jl")
export NSVerify, ConvDual, Duality, Certify, ILP, MIPVerify,
       BaB, Sherlock, Reluplex, Planet

include("reachability/utils/reachability.jl")
include("reachability/exactReach.jl")
include("reachability/maxSens.jl")
include("reachability/ai2.jl")
include("reachability/reluVal.jl")
include("reachability/neurify.jl")
include("reachability/fastLin.jl")
include("reachability/fastLip.jl")
include("reachability/dlv.jl")

# added by me
include("reachability/DeepPolyBounds/utils.jl")
include("reachability/deep_poly.jl")
include("reachability/DeepPolyBounds/network_neg_pos_idx.jl")
include("reachability/DeepPolyBounds/symbolic_interval_bounds.jl")
include("reachability/DeepPolyBounds/symbolic_interval_fresh_vars.jl")
include("reachability/DeepPolyBounds/deep_poly_bounds.jl")
include("reachability/DeepPolyBounds/deep_poly_fresh_vars.jl")
include("reachability/DeepPolyBounds/symbolic_interval_heur.jl")
include("reachability/DeepPolyBounds/deep_poly_heuristic.jl")
include("reachability/CROWN/crown.jl")
include("reachability/DeepPolyBounds/fresh_var_heuristic.jl")
include("reachability/AsymESIP/asym_esip.jl")
include("reachability/AsymESIP/esip_solver.jl")
include("reachability/DeepPolyBounds/symbolic_interval_heur_fv.jl")
include("reachability/DeepPolyBounds/dp_neurify_fv.jl")
include("reachability/DeepPolyBounds/dp_neurify_zono.jl")

export ExactReach, MaxSens, Ai2, Ai2h, Ai2z, Box,
       ReluVal, Neurify, FastLin, FastLip, DLV,
       DeepPoly, DeepPolyBounds, DeepPolyFreshVars,
       DeepPolyHeuristic, CROWN, backward_network, calc_bounds, ESIPSolver,
       DPNeurifyFV, DPNeurifyZono

const TOL = Ref(sqrt(eps()))
set_tolerance(x::Real) = (TOL[] = x)
export set_tolerance

end


text:
module TestInference

include("preamble.jl")

madd = ManyAdd()

@test_inferred bcapp(madd, 1:5...)

end  # module


text:
using Documenter, SimpleCache

makedocs(;
    modules=[SimpleCache],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ssfrr/SimpleCache.jl/blob/{commit}{path}#L{line}",
    sitename="SimpleCache.jl",
    authors="Spencer Russell",
    assets=String[],
)

deploydocs(;
    repo="github.com/ssfrr/SimpleCache.jl",
)


text:
using BHAtp, Statistics

ProjDir = @__DIR__
!isdir(joinpath(ProjDir, "plots")) && mkdir(joinpath(ProjDir, "plots"))

actual = DataFrame(
  incls = [40.0, 42.3, 44.0, 44.65, 45.7, 47.25, 49.5],
  atp = [2.55, 2.65, 0.4, 1.15, 1.18, 2.15, 2.35],
)

rpm = DataFrame(
  incls = [35, 43.5,43.6,44.2,44.3,44.9, 45, 48,48.1, 55],
  rpm = [ 100,  100,  60,  60,  70,  70, 80, 80, 120,120]
)

name = "figure_03"

figure_03 = @pgf TikzPicture(
  Axis({
      height = "6cm", width = " 8cm",
      xmin = 30, xmax = 90,
      ymin = 0, ymax = 6.2,
      axis_y_line = "left",
      xlabel = "Inclination [°]",
      ylabel = "TP [°/100ft]",
      title = "Actual vs. Isotropic TP",
      legend_pos = "south east",
    },
    Plot({ color="blue", mark = "." },
      Table(tpf.incl[1:7], tpf.tp[1:7])),
    LegendEntry("TP wob=30"),
    Plot({ color="gray", mark = "." },
      Table(tpf.incl[8:14], tpf.tp[8:14])),
    LegendEntry("TP wob=35"),
    Plot({ color="green", mark = "." },
      Table(tpf.incl[15:21], tpf.tp[15:21])),
    LegendEntry("TP wob=35"),
    Plot({ thick, color="black", mark = "+" },
      Table(actual.incls, actual.atp)),
    LegendEntry("Actual")),
  Axis({
      height = "6cm", width = " 8cm",
      ymin = 0, ymax = 160,
      xmin = 30, xmax = 90,
      hide_x_axis,
      axis_y_line = "right",
      ylabel = "RPM [revs/min]",
      legend_pos = "north east",
    },
    Plot({ color="red", mark = "." },
      Table(rpm.incls, rpm.rpm)),
    LegendEntry("RPM")),
)

savefig(joinpath(ProjDir, "plots", name), figure_03, ProjDir)


text:
"""
    energy_balance(object, meteo, constants = Constants())
    energy_balance!(object, meteo, constants = Constants())

Computes the energy balance of one or several components based on the type of the model it was
parameterized with in `object.energy`, and on one or several meteorology time-steps.

At the moment, two models are implemented in the package:

- [`Monteith`](@ref): the model found in Monteith and Unsworth (2013)
- `Missing`: if no computation of the energy balance is needed

# Arguments

- `object`: a [`Component`](@ref) struct ([`AbstractComponentModel`](@ref)), or a Dict/Array of.
- `meteo::Union{AbstractAtmosphere,Weather}`: meteorology structure, see [`Atmosphere`](@ref) or
[`Weather`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

Some models need input values for some variables. For example [`Monteith`](@ref) requires a
value for `Rₛ`, `d` and `skyFraction`. If you read the models from a file, you can
use [`init_status!`](@ref) (see examples).

# Examples

```julia
# ---Simple example---

meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

# Using the model of Monteith and Unsworth (2013) for energy, Farquhar et al. (1980) for
# photosynthesis, and Medlyn et al. (2011) for stomatal conductance:
leaf = LeafModels(energy = Monteith(),
            photosynthesis = Fvcb(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Rₛ = 13.747, skyFraction = 1.0, PPFD = 1500.0, d = 0.03)

energy_balance(leaf,meteo)

# ---Using several components---

leaf2 = copy(leaf)
leaf2.status.PPFD = 800.0

energy_balance([leaf,leaf2],meteo)

# You can use a Dict if you'd like to keep track of the leaf in the returned DataFrame:
energy_balance(Dict(:leaf1 => leaf, :leaf2 => leaf2), meteo)

# ---Using several meteo time-steps---

w = Weather([Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65),
             Atmosphere(T = 25.0, Wind = 1.5, P = 101.3, Rh = 0.55)],
             (site = "Test site",))

energy_balance(leaf, w)

# ---Using several meteo time-steps and several components---

energy_balance(Dict(:leaf1 => leaf, :leaf2 => leaf2), w)

# ---Using a model file---

model = read_model(joinpath(dirname(dirname(pathof(PlantBiophysics))),"test","inputs","models","plant_coffee.yml"))

# An example model file is available here:
# "https://raw.githubusercontent.com/VEZY/PlantBiophysics/main/test/inputs/models/plant_coffee.yml"

# Initialising the mandatory variables:
init_status!(model, Rₛ = 13.747, skyFraction = 1.0, PPFD = 1500.0, Tₗ = 25.0, d = 0.03)

# NB: To know which variables has to be initialised according to the models used, you can use
# `to_initialise(Component)`, *e.g.*:
to_initialise(model["Leaf"])

# Running a simulation for all component types in the same scene:
energy_balance!(model, meteo)

model["Leaf"].status.Rn
model["Leaf"].status.A
model["Leaf"].status.Cᵢ
```

# References

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5 (4):
919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Monteith, John L., et Mike H. Unsworth. 2013. « Chapter 13 - Steady-State Heat Balance: (i)
Water Surfaces, Soil, and Vegetation ». In Principles of Environmental Physics (Fourth Edition),
edited by John L. Monteith et Mike H. Unsworth, 217‑47. Boston: Academic Press.

Schymanski, Stanislaus J., et Dani Or. 2017. « Leaf-Scale Experiments Reveal an Important
Omission in the Penman–Monteith Equation ». Hydrology and Earth System Sciences 21 (2): 685‑706.
https://doi.org/10.5194/hess-21-685-2017.

Vezy, Rémi, Mathias Christina, Olivier Roupsard, Yann Nouvellon, Remko Duursma, Belinda Medlyn,
Maxime Soma, et al. 2018. « Measuring and modelling energy partitioning in canopies of varying
complexity using MAESPA model ». Agricultural and Forest Meteorology 253‑254 (printemps): 203‑17.
https://doi.org/10.1016/j.agrformet.2018.02.005.
"""
energy_balance!, energy_balance

function energy_balance!(object::AbstractComponentModel, meteo::AbstractAtmosphere, constants = Constants())
    is_init = is_initialised(object)
    !is_init && error("Some variables must be initialized before simulation")

    net_radiation!(object, meteo, constants)
    return nothing
end

# Same as above but non-mutating
function energy_balance(object::AbstractComponentModel, meteo::AbstractAtmosphere, constants = Constants())
    object_tmp = copy(object)
    energy_balance!(object_tmp, meteo, constants)

    return object_tmp
end

# energy_balance over several objects (e.g. all leaves of a plant) in an Array
function energy_balance!(object::O, meteo::AbstractAtmosphere, constants = Constants()) where O <: AbstractArray{<:AbstractComponentModel}

    for i in values(object)
        energy_balance!(i, meteo, constants)
    end

    return nothing
end

# energy_balance over several objects (e.g. all leaves of a plant) in a kind of Dict.
function energy_balance!(object::O, meteo::AbstractAtmosphere, constants = Constants()) where {O <: AbstractDict{N,<:AbstractComponentModel} where N}

    for (k, v) in object
        energy_balance!(v, meteo, constants)
    end

    return nothing
end

# same as the above but non-mutating
function energy_balance(
    object::O,
    meteo::AbstractAtmosphere,
    constants = Constants()
    ) where O <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N}

    # Copy the objects only once before the computation for performance reasons:
    object_tmp = copy(object)

    # Computation:
    energy_balance!(object_tmp, meteo, constants)

    return object_tmp
end

# energy_balance over several meteo time steps (called Weather) and possibly several components:
function energy_balance!(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N}

    # Check if the meteo data and the status have the same length (or length 1)
    check_status_wheather(object, meteo)

    # Computing for each time-steps:
    for (i, meteo_i) in enumerate(meteo.data)
        # Each object in a time-step:
        for obj in object
            energy_balance!(obj[i], meteo_i, constants)
        end
    end

end

# If we call weather with one component only, put it in an Array and call the function above
function energy_balance!(object::AbstractComponentModel, meteo::Weather, constants = Constants())
    energy_balance!([object], meteo, constants)
end

# energy_balance over several meteo time steps (same as above) but non-mutating
function energy_balance(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractComponentModel,AbstractDict{N,<:AbstractComponentModel} where N}

    object_tmp = copy(object)

    energy_balance!(object_tmp, meteo, constants)

    return object_tmp
end


text:
using BasePlus
using Base.Test

include("os.jl")
include("text.jl")


text:
using GigaSOM, DataFrames, XLSX, CSV, Test, Random, Distributed, SHA, JSON

owd = pwd()

checkDir()

@testset "GigaSOM test suite" begin
    #load and transform the .fcs data to be ready for computing
    include("io.jl")

    #apply the batch GigaSOM algorithm to the data, train it and test it
    include("batch.jl")

    #apply the parallel GigaSOM algorithm to the data, train it and test it
    include("parallel.jl")
end

cd(owd)


text:
"""
`init_subspec!(m::SmetsWoutersOrig)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::SmetsWoutersOrig)

    if subspec(m) == "ss0"
        return
    elseif subspec(m) == "ss1"
        # Fix γ at 0.4312
        ss1!(m)
    elseif subspec(m) =="ss2"
        # Diffuse prior
        ss2!(m)
    else
        error("This subspec is not defined.")
    end
end

function ss1!(m::SmetsWoutersOrig)
    m <= parameter(:γ, 0.4312, (-5.0, 5.0), (-5., 5.), Untransformed(), Normal(0.4, 0.1),
                   fixed=true, scaling = x -> 1 + x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="\\gamma")
end

# Diffuse prior
function ss2!(m::SmetsWoutersOrig)
 m <= parameter(:S′′, 6.3325, (-15., 15.), (-15., 15.), Untransformed(), Distributions.Normal(4., 4.5), fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.",
                   tex_label="S''")

    m <= parameter(:σ_c, 1.2312, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Distributions.Normal(1.5, 1.11), fixed=false,
                   tex_label="\\sigma_{c}")

    m <= parameter(:h, 0.7205, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="h: Consumption habit persistence.",
                   tex_label="h")

    m <= parameter(:λ_w, 1.5000, fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:ζ_w, 0.7937, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ϵ_w, 10.000, fixed=true,
                   description="ϵ_w: Curvature parameter in the Kimball aggregator for wages.",
                   tex_label="\\epsilon_{w}")

    m <= parameter(:ν_l, 2.8401, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Distributions.Normal(2.0, 2.25), fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:δ, 0.025, fixed=true,
                   description="δ: The capital depreciation rate.",
                   tex_label="\\delta" )

    m <= parameter(:ζ_p, 0.7813, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ϵ_p, 10.000, fixed=true,
                   description="ϵ_p: Curvature parameter in the Kimball aggregator for prices.",
                   tex_label="\\epsilon_{p}")

    m <= parameter(:ι_w, 0.4425, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ι_w: The weight attributed to last period's wage in wage indexation. (1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label="\\iota_w")
    m <= parameter(:ι_p, 0.3291, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:ppsi, 0.2648, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ppsi: Utilization costs.",
                   tex_label="\\psi")

    m <= parameter(:Φ, 1.4672, (1., 10.), (1.00, 10.00), ModelConstructors.Exponential(), Distributions.Normal(1.25, 0.36), fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:ψ1, 1.7985, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Distributions.Normal(1.5, 0.75), fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ρ, .8258, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")

    m <= parameter(:ψ2, 0.0893, (-0.5, 0.5), (-0.5, 0.5), Untransformed(), Distributions.Normal(0.12, 0.15), fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2239, (-0.5, 0.5), (-0.5, 0.5), Untransformed(), Distributions.Normal(0.12, 0.15), fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:π_star, 0.7000, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.62, 0.3), fixed=false, scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:β, 0.7420, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.3), fixed=false, scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="\\beta ")

    m <= parameter(:Lmean, -1., (-5., 5.), (-5., 5.), Untransformed(), Distributions.Normal(0.0, 6.0), fixed=false,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")

    m <= parameter(:γ, 0.4312, (-5.0, 5.0), (-5., 5.), Untransformed(), Distributions.Normal(0.4, 0.3), fixed=false, scaling = x -> 1 + x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="\\gamma")

    m <= parameter(:α, 0.24, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Normal(0.30, 0.15), fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's Cobb-Douglas production function.",
                   tex_label="\\alpha")

m <= parameter(:g_star, 0.1800, fixed=true,
                   description="g_star: 1 - (c_star + i_star)/y_star.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_z, 0.9676, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_b, 0.2703, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")
#added 4 9's to upper bound for the following parameter
    m <= parameter(:ρ_g, 0.9930, (1e-5, 0.9999999), (1e-5, 0.9999999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_μ, 0.5724, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_rm, 0.3000, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.", # CHECK THIS
                   tex_label="\\rho_{r^m}")

    m <= parameter(:ρ_λ_f, 0.8692, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w, 0.9546, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.", # CHECK THIS
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:η_λ_f, 0.7652, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w, 0.8936, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

    m <= parameter(:η_gz, 0.0500, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_z, 0.4618, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_b, 0.1818, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_g, 0.6090, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_μ, 0.4601, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_rm, 0.2397, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_r_m: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{r^m}")

    m <= parameter(:σ_λ_f, 0.1455, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w, 0.2089, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=false,
                   tex_label="\\sigma_{\\lambda_w}")

end


text:
 # api3.jl

include("../domain/domain.jl")

# test with different keyword values
include("./myfunctions_pids.jl")


text:
using Documenter, MonteCarloZoo

push!(LOAD_PATH,"../src/")
makedocs(sitename="MonteCarloZoo.jl Documentation",
         pages = [
            "Index" => "index.md",
         ],
         format = Documenter.HTML(prettyurls = false)
)

deploydocs(
    repo = "github.com/THargreaves/MonteCarloZoo.jl.git",
    devbranch = "master"
)


text:
module VlasovBenchmarks

using BenchmarkTools

const SUITE = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","VlasovSolver","LaxWendroff.jl"))
import .LaxWendroff
SUITE["LaxWendroff"] = BenchmarkGroup()
SUITE["LaxWendroff c"] = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","VlasovSolver","Upwind.jl"))
import .Upwind
SUITE["Upwind"] = BenchmarkGroup()
SUITE["Upwind с"] = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","VlasovSolver","Godunov.jl"))
import .Godunov
SUITE["Godunov constant"] = BenchmarkGroup()
SUITE["Godunov linear"] = BenchmarkGroup()
SUITE["Godunov linear VanLeer"] = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","VlasovSolver","SemiLagrangian.jl"))
import .SemiLagrangian
SUITE["SemiLagrangian linear"] = BenchmarkGroup()
SUITE["SemiLagrangian quadratic"] = BenchmarkGroup()
SUITE["SemiLagrangian cubic"] = BenchmarkGroup()

include(joinpath(dirname(@__FILE__),"..","src","VlasovSolver","PFC.jl"))
import .PFC
SUITE["PFC"] = BenchmarkGroup()

Δx = Dict()
Δt = Dict()
v = Dict()
f₀ = Dict()
f = Dict()
advect! = Dict()

for N in [100, 1000, 10000]
    Δx[N] = 0.01
    Δt[N] = 0.8*Δx[N]
    v[N] = 1
    f₀[N] = [1.0 + 0.01*sin(2π*i*Δx[N]) for i = 0:N]
    f[N] = similar(f₀[N])
    advect![N] = Dict()

    advect![N][:LaxWendroff] = LaxWendroff.generate_solver(f₀[N], f[N], v[N]*Δt[N]/Δx[N])
    SUITE["LaxWendroff"]["advect $N"] = @benchmarkable advect![$N][:LaxWendroff]()

    advect![N][:LaxWendroff_c] = LaxWendroff.generate_solver(f₀[N], f[N])
    SUITE["LaxWendroff c"]["advect $N"] = @benchmarkable advect![$N][:LaxWendroff_c]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:Upwind] = Upwind.generate_solver(f₀[N], f[N], v[N]*Δt[N]/Δx[N])
    SUITE["Upwind"]["advect $N"] = @benchmarkable advect![$N][:Upwind]()

    advect![N][:Upwind_с] = Upwind.generate_solver(f₀[N], f[N])
    SUITE["Upwind с"]["advect $N"] = @benchmarkable advect![$N][:Upwind_с]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:Godunov_constant] = Godunov.generate_solver(f₀[N], f[N], :Riemann_constant)
    SUITE["Godunov constant"]["advect $N"] = @benchmarkable advect![$N][:Godunov_constant]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:Godunov_linear] = Godunov.generate_solver(f₀[N], f[N], :Riemann_linear)
    SUITE["Godunov linear"]["advect $N"] = @benchmarkable advect![$N][:Godunov_linear]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:Godunov_linear_VanLeer] = Godunov.generate_solver(f₀[N], f[N], :Riemann_linear; flux_limiter = :VanLeer)
    SUITE["Godunov linear VanLeer"]["advect $N"] = @benchmarkable advect![$N][:Godunov_linear_VanLeer]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:SemiLagrangian_linear] = SemiLagrangian.generate_solver(f₀[N], f[N]; interpolation_order = :Linear)
    SUITE["SemiLagrangian linear"]["advect $N"] = @benchmarkable advect![$N][:SemiLagrangian_linear]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:SemiLagrangian_quadratic] = SemiLagrangian.generate_solver(f₀[N], f[N]; interpolation_order = :Quadratic)
    SUITE["SemiLagrangian quadratic"]["advect $N"] = @benchmarkable advect![$N][:SemiLagrangian_quadratic]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:SemiLagrangian_cubic] = SemiLagrangian.generate_solver(f₀[N], f[N]; interpolation_order = :Cubic)
    SUITE["SemiLagrangian cubic"]["advect $N"] = @benchmarkable advect![$N][:SemiLagrangian_cubic]($(v[N]*Δt[N]/Δx[N]))

    advect![N][:PFC] = PFC.generate_solver(f₀[N], f[N])
    SUITE["PFC"]["advect $N"] = @benchmarkable advect![$N][:PFC]($(v[N]*Δt[N]/Δx[N]))
end

end # module


text:
module Day14

using Advent: save_input

function solve_a(filename)
    memory = Dict{Int,UInt}()
    mask0 = nothing
    mask1 = nothing
    for line in eachline(filename)
        left, right = split(line, r" = ")
        if left == "mask"
            mask0 = reduce((x, y) -> (x << 1) | (y != '0'), right, init = zero(UInt64))
            mask1 = reduce((x, y) -> (x << 1) | (y == '1'), right, init = zero(UInt64))
        else
            address = parse(Int, left[5:end - 1])
            value = parse(UInt, right)
            memory[address] = value & mask0 | mask1
        end
    end
    memory |> values |> sum |> Int
end

function solve_b(filename)
    memory = Dict{Int,UInt}()
    mask1 = nothing
    maskf = nothing
    for line in eachline(filename)
        left, right = split(line, r" = ")
        if left == "mask"
            maskf = reduce((x, y) -> (x << 1) | (y == 'X'), right, init = zero(UInt64))
            mask1 = reduce((x, y) -> (x << 1) | (y == '1'), right, init = zero(UInt64))
        else
            address = parse(Int, left[5:end - 1]) & ~maskf | mask1
            value = parse(UInt, right)
            st = maskf
            while st != 0
                # cycle through bits in the mask
                memory[address | st] = value
                st = (st - 1) & maskf
            end
        end
    end
    memory |> values |> sum
end

solve(filename = save_input(2020, 14)) = solve_a(filename), solve_b(filename)
solve_big() = solve("day14_big.in")

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
abstract type AbstractPeriodicOrbitProblem end

# Periodic orbit computations by finite differences
abstract type AbstractPOFDProblem <: AbstractPeriodicOrbitProblem end

abstract type AbstractShootingProblem <: AbstractPeriodicOrbitProblem end

# get the number of time slices / sections
@inline getM(pb::AbstractPeriodicOrbitProblem) = pb.M

# update a problem with arguments
function update(prob::AbstractPeriodicOrbitProblem, F, dF, hopfpt, ζr, M, orbitguess_a, period) end

####################################################################################################
# if we use the same code as for newton (see below) in continuation, it is difficult to tell the eigensolver not to use the jacobian but instead the monodromy matrix. So we have to use a dedicated composite type for the jacobian to handle this case.

struct ShootingJacobian{Tpb <: AbstractShootingProblem, Torbitguess, Tp}
	pb::Tpb
	x::Torbitguess
	par::Tp
end

# evaluation of the jacobian
(shjac::ShootingJacobian)(dx) = shjac.pb(shjac.x, shjac.par, dx)

####################################################################################################
# newton wrapper
"""
$(SIGNATURES)

This is the Newton-Krylov Solver for computing a periodic orbit using (Standard / Poincaré) Shooting method.
Note that the linear solver has to be apropriately set up.

# Arguments

Similar as [`newton`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). These two problems have specific options to be tuned, we refer to their link for more information and to the tutorials.

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::AbstractShootingProblem, orbitguess, par, options::NewtonPar; kwargs...)
	return newton(prob,
			(x, p) -> (dx -> prob(x, p, dx)),
			orbitguess, par, options; kwargs...)
end

"""
$(SIGNATURES)

This is the deflated Newton-Krylov Solver for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar as [`newton`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref).

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::AbstractShootingProblem, orbitguess, options::NewtonPar, defOp::DeflationOperator; kwargs...)
	return newton(prob,
			(x, p) -> (dx -> prob(x, p, dx)),
			orbitguess, par,
			options, defOp; kwargs...)
end

####################################################################################################
# Continuation for shooting problems

"""
$(SIGNATURES)

This is the continuation method for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar as [`continuation`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref).

- `printPeriod` boolean to print the period of the solution. This is useful for `prob::PoincareShootingProblem` as this information is not easily available.
"""

function continuation(prob::AbstractShootingProblem, orbitguess, par, lens::Lens, _contParams::ContinuationPar, linearAlgo::AbstractBorderedLinearSolver; printPeriod = true, kwargs...)

	contParams = check(_contParams)

	options = contParams.newtonOptions

	if contParams.computeEigenValues
		contParams = @set contParams.newtonOptions.eigsolver = FloquetQaDShooting(contParams.newtonOptions.eigsolver)
	end

	if (prob isa PoincareShootingProblem)

		if printPeriod
			printSolutionPS = (x, p) -> getPeriod(prob, x, set(par, lens, p))
			return continuation(
				prob,
				(x, p) -> ShootingJacobian(prob, x, p),
				orbitguess, par, lens,
				contParams, linearAlgo;
				printSolution = printSolutionPS,
				kwargs...)
		end
	end

	return continuation(
		prob,
		(x, p) -> ShootingJacobian(prob, x, p),
		orbitguess, par, lens,
		contParams, linearAlgo;
		kwargs...)
end

"""
$(SIGNATURES)

This is the continuation routine for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar as [`continuation`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref).

- `printPeriod` boolean to print the period of the solution. This is useful for `prob::PoincareShootingProblem` as this information is not easily available.
"""
function continuation(prob::AbstractShootingProblem, orbitguess, par, lens::Lens, contParams::ContinuationPar; linearAlgo = BorderingBLS(), printPeriod = true, kwargs...)
	_linearAlgo = @set linearAlgo.solver = contParams.newtonOptions.linsolver
	return continuation(prob, orbitguess, par, lens, contParams, _linearAlgo; printPeriod = printPeriod, kwargs...)
end

####################################################################################################
"""
$(SIGNATURES)

Perform automatic branch switching from a Hopf bifurcation point labelled `ind_bif` in the list of the bifurcated points on a previously computed branch `br::ContResult`. It first compute a Hopf normal form.

# Arguments

- `F, dF, d2F, d3F`: function `(x,p) -> F(x,p)` and its differentials `(x,p,dx) -> d1F(x,p,dx)`, `(x,p,dx1,dx2) -> d2F(x,p,dx1,dx2)`... These are used to compute the Hopf normal form.
- `br` branch result from a call to `continuation`
- `ind_hopf` index of the bifurcation point in `br`
- `contParams` parameters for the call to `continuation`
- `prob` problem used to specify the way the periodic orbit is computed. It can be [`PeriodicOrbitTrapProblem`](@ref), [`ShootingProblem`](@ref) or [`PoincareShootingProblem`](@ref) .

# Optional arguments

- `linearPO` linear algorithm used for the computation of periodic orbits when `prob` is [`PeriodicOrbitTrapProblem`](@ref))
- `Jt` is the jacobian adjoint, used for computation of the eigen-elements of the jacobian adjoint, needed to compute the spectral projector for the Hopf normal form!!!!!! COMME NEWTON FOLD
- `δ = 1e-8` used for finite differences
- `δp` used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by `contParams.ds`. This allows to use a step larger than `contParams.dsmax`.
- `ampfactor = 1` factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
- `usedeflation = true` whether to use nonlinear deflation (see [Deflated problems](@ref)) to help finding the guess on the bifurcated branch

!!! note "Linear solver"
    You have to be carefull about the options `contParams.newtonOptions.linsolver`. In the case of Matrix-Free solver, you have to pass the right number of unknowns `N * M + 1`. Note that the options for the preconditioner are not accessible yet.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not `nothing`, is computed with Finite differences.
"""
function continuation(F, dF, d2F, d3F, br::ContResult, ind_bif::Int, _contParams::ContinuationPar, prob::AbstractPeriodicOrbitProblem ; Jt = nothing, δ = 1e-8, δp = nothing, ampfactor = 1, usedeflation = false, kwargs...)
	# compute the normal form of the branch point
	verbose = get(kwargs, :verbosity, 0) > 1 ? true : false
	cb = get(kwargs, :callbackN, (args...;kwargs...)->true)

	hopfpt = hopfNormalForm(F, dF, d2F, d3F, br, ind_bif ; Jt = Jt, δ = δ, nev = _contParams.nev, verbose = verbose)

	# compute predictor for point on new branch
	ds = isnothing(δp) ? _contParams.ds : δp
	Ty = typeof(ds)
	pred = predictor(hopfpt, ds; verbose = verbose, ampfactor = Ty(ampfactor))

	verbose && printstyled(color = :green, "#"^51*
			"\n--> Start Hopf branch switching.
			\n--> Bifurcation type = ", hopfpt.type,
			"\n----> newp = ", pred.p, ", δp = ", pred.p - br.bifpoint[ind_bif].param,
			"\n----> amplitude = ", pred.amp,
			"\n----> period = ", abs(2pi/pred.ω), "\n")

	# we compute a phase so that the constraint equation
	# < u(0) − u_hopf, ψ > is satisfied, i.e. equal to zero.
	ζr = real.(hopfpt.ζ); ζi = imag.(hopfpt.ζ)
	# this phase is for POTrap problem constraint to be satisfied
	ϕ = atan(dot(ζr, ζr), dot(ζi, ζr))
	verbose && printstyled(color = :green,"----> phase ϕ = ", ϕ/pi, "⋅π\n")

	M = getM(prob)
	orbitguess_a = [pred.orbit(t - ϕ) for t in LinRange(0, 2pi ,M+1)[1:M]]

	# build the variable to hold the functional for computing PO based on finite differences
	probPO, orbitguess = update(prob, F, dF, hopfpt, ζr, M, orbitguess_a, abs(2pi/pred.ω))

	if _contParams.newtonOptions.linsolver isa GMRESIterativeSolvers
		_contParams = @set _contParams.newtonOptions.linsolver.N = length(orbitguess)
	end

	# perform continuation
	branch, u, tau = continuation(probPO, orbitguess, set(br.params, br.param_lens, pred.p), br.param_lens, _contParams; kwargs...)

	return setproperties(branch; type = :PeriodicOrbit, functional = probPO), u, tau

end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ d307b722-47fa-4288-8291-a2d17fe276f9
#hide
using PlutoUI

# ╔═╡ c1c575f0-f8b3-11eb-3679-27a9b018100f
md"""

**Alpha Release**

## Nbdev
This project is an attempt to port fastai [Nbdev](https://github.com/fastai/nbdev) to Julia. Nbdev(both the oriignal version and the Julia version) are based on the idea of using Notebooks for literate programming.

On a very high level the base idea is this 👇

$(LocalResource("../docs/images/nbdev.png"))

While both projects try to achieve the same base goal i.e. generating source code and documentation from Notebooks, the Julia version aims to start from Pluto notebook while the python version starts from a Jupyter notebook.

On the documentation side, I am using Mkdocs as the document generator while the original version uses Jekyll as the static site generator of choice."""

# ╔═╡ e9a9f361-eaf1-4546-a83e-516eb189e3e5
md"## Getting-started"

# ╔═╡ c36b7a25-4b76-4c91-a77e-a677513158d0
md"Follow the [tutorial](https://sapal6.github.io/Nbdev.jl/tutorial/) to start experimenting with Nbdev."

# ╔═╡ bb78bb31-8a72-4e9e-a1c5-31457fd4326c
md"## First among the version

This is an alpha release of Nbdev so expect some breaks in plumbing here and there. Some features are still under development and are not included in this version.

* Feel free to log any issues, bugs (which I expect to come) [here](https://github.com/sapal6/Nbdev.jl/issues)

* Anything that you feel like is missing and want to discuss the same or havean idea then start a new discussion [here](https://github.com/sapal6/Nbdev.jl/discussions)"

# ╔═╡ 6ee57f9a-dfd7-4bea-802d-eba78f8fc99c
md"## What's more to come?

Well much more can be expected some of which are already thought of and some are not. However, the immediate requirements for future release are the following -->

* Backlinks to source code.
* Automatic document (markdown) generation using Github actions.
* Testing for including plots/graphs from notebook to markdown document."

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.11"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "0c3e067931708fa5641247affc1a1aceb53fff06"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.11"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═d307b722-47fa-4288-8291-a2d17fe276f9
# ╠═c1c575f0-f8b3-11eb-3679-27a9b018100f
# ╠═e9a9f361-eaf1-4546-a83e-516eb189e3e5
# ╠═c36b7a25-4b76-4c91-a77e-a677513158d0
# ╠═bb78bb31-8a72-4e9e-a1c5-31457fd4326c
# ╠═6ee57f9a-dfd7-4bea-802d-eba78f8fc99c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using LinearAlgebraicRepresentation
using Base.Test

@testset "LarGrid-0-1 Tests" begin
   @testset "grid_0" begin
      @test LinearAlgebraicRepresentation.grid_0(1) == [0  1]
      @test LinearAlgebraicRepresentation.grid_0(2) == [0  1  2]
      @test LinearAlgebraicRepresentation.grid_0(3) == [0  1  2  3]
   end

   @testset "Grid_1 Tests" begin
      @test repr(LinearAlgebraicRepresentation.grid_1(1)) == "[0; 1]"
      @test LinearAlgebraicRepresentation.grid_1(2) == [0 1; 1 2]
      @test LinearAlgebraicRepresentation.grid_1(3) == [0 1 2; 1 2 3]
   end

   @testset "LarGrid Tests" begin
      @test repr(LinearAlgebraicRepresentation.larGrid(1)(0)) == "[0 1]"
      @test repr(LinearAlgebraicRepresentation.larGrid(1)(1)) == "[0; 1]"
      @test LinearAlgebraicRepresentation.larGrid(2)(0) == [0  1  2]
      @test LinearAlgebraicRepresentation.larGrid(2)(1) == [0 1; 1 2]
      @test LinearAlgebraicRepresentation.larGrid(3)(0) == [0  1  2  3]
      @test LinearAlgebraicRepresentation.larGrid(3)(1) == [0  1  2; 1  2  3]
   end
end;

@testset "LarVertProd Tests" begin
	vertexDomain(n) = hcat([k for k in 0:n-1]...)
	@testset "LarVertProd 1D" begin
      shape = [3]
	vertLists = [vertexDomain(k+1) for k in shape]
      @test typeof(LinearAlgebraicRepresentation.larVertProd(vertLists))==Array{Int64,2}
      @test size(LinearAlgebraicRepresentation.larVertProd(vertLists))==(1, 4)
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,1]==[0]
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,4]==[3]
   end

   @testset "LarVertProd 2D" begin
      shape = [3,2]
	vertLists = [vertexDomain(k+1) for k in shape]
      @test typeof(LinearAlgebraicRepresentation.larVertProd(vertLists))==Array{Int64,2}
      @test size(LinearAlgebraicRepresentation.larVertProd(vertLists))==(2, 12)
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,1]==[0;0]
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,12]==[3;2]
   end

   @testset "LarVertProd 3D" begin
      shape = [3,2,1]
	vertLists = [vertexDomain(k+1) for k in shape]
      @test typeof(LinearAlgebraicRepresentation.larVertProd(vertLists))==Array{Int64,2}
      @test size(LinearAlgebraicRepresentation.larVertProd(vertLists))==(3, 24)
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,1]==[0;0;0]
      @test LinearAlgebraicRepresentation.larVertProd(vertLists)[:,24]==[3;2;1]
   end
end

@testset "Index2addr Tests" begin
   @testset "Shape 1D Tests" begin
      @test LinearAlgebraicRepresentation.index2addr([10])([0])==1
      @test LinearAlgebraicRepresentation.index2addr([10])([9])==10
      @test [LinearAlgebraicRepresentation.index2addr([10])([index]) for index in collect(0:9)]==collect(1:10)
   end

   @testset "Shape 2d Tests" begin
      aa = LinearAlgebraicRepresentation.cart([[0;1;2],[0;1]])
      bb = LinearAlgebraicRepresentation.cart([collect(0:9),collect(0:1)])
      cc = LinearAlgebraicRepresentation.cart([collect(0:2),collect(0:2)])
      dd = "Tuple{Int64,Int64}[(0,0),(0,1),(1,0),(1,1),(2,0),(2,1)]"
      ee = "Tuple{Int64,Int64}[(0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1)]"
      ff = repr( LinearAlgebraicRepresentation.cart([collect(0:2),collect(0:1)]) )
      @test [ LinearAlgebraicRepresentation.index2addr([3,2])(collect(index)) for index in aa ]==collect(1:6)
      @test [ LinearAlgebraicRepresentation.index2addr([10,2])(collect(index)) for index in bb ] == collect(1:20)
      @test [ LinearAlgebraicRepresentation.index2addr([3,3])(collect(index)) for index in cc] == collect(1:9)
   end

   @testset "Shape 3d Tests" begin
      aaa = LinearAlgebraicRepresentation.cart([collect(0:3),collect(0:2),collect(0:1)])
      bbb = LinearAlgebraicRepresentation.cart([[0;1;2],[0;1],[0;1]])
      ccc = LinearAlgebraicRepresentation.cart([[0;1],[0;1],[0;1;2]])
      @test LinearAlgebraicRepresentation.index2addr([3,2,1])([0,0,0]) == 1
      @test [LinearAlgebraicRepresentation.index2addr([4,3,2])(collect(index))
      	for index in aaa ] == collect(1:24)
      @test [LinearAlgebraicRepresentation.index2addr([3,2,2])(collect(index))
      	for index in bbb ] == collect(1:12)
      @test [LinearAlgebraicRepresentation.index2addr([2,2,3])(collect(index))
      	for index in ccc ] == collect(1:12)
   end
end

@testset "BinaryRange Tests" begin
	binaryRange(n) = bin.(range(0,2^n),n)
   @test binaryRange(1)==["0";"1"]
   @test binaryRange(2)==["00","01","10","11"]
   @test binaryRange(3)==["000","001","010","011","100","101","110","111"]
end

@testset "LarCellProd Tests" begin
   @testset "$shape" for shape in [[3,2,1],[3,2],[3],[1,1,1]]
      @testset "$d" for d in 0:length(shape)
         @test typeof(LinearAlgebraicRepresentation.larGridSkeleton(shape)(d)) == Array{Array{Int64,1},1}
      end
   end
end

@testset "FilterByOrder Tests" begin
	binaryRange(n) = bin.(range(0,2^n),n)
   term = "000"
   bit = '0'
   theTerm = convert(Array{Char,1},term)
   @test typeof(theTerm) == Array{Char,1}
   @test parse(Int8,bit) == 0
   @test [parse(Int8,bit) for bit in theTerm] == zeros(3)
   out = hcat([[parse(Int8,bit) for bit in term] for term in binaryRange(3)]...)
   @test typeof(out) == Array{Int8,2}
   @test size(out) == (3,8)
   @test repr(out) == "Int8[0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1]"

   @testset  "$n" for n in 1:4
      data = [LinearAlgebraicRepresentation.filterByOrder(n)[k] for (k,el) in enumerate(LinearAlgebraicRepresentation.filterByOrder(n))]
      @test sum(map(length,data)) == 2^n
   end

   @testset "$n,$k" for n in 1:4, k in 0:n
      @test length(LinearAlgebraicRepresentation.filterByOrder(n)[k+1]) == binomial(n,k)
   end
end

@testset "LarImageVerts Tests" "$shape" for shape in [[3,2,1],[3,2],[10,10,10]]
   @test size(LinearAlgebraicRepresentation.larImageVerts(shape)) == (length(shape),prod(shape + 1))
end
@testset "LarCuboids Tests" begin
    shape = [3,2,1]
    cubes = LinearAlgebraicRepresentation.larCuboids(shape,true)
    verts, cells = cubes
    @test typeof(verts) == Array{Float64,2}
    VV,EV,FV,CV = cells
    @testset "$basis" for basis in [VV,EV,FV,CV]
      @test typeof(basis) == Array{Array{Int64,1},1}
    end
end

@testset "LarImageVerts Tests" "$shape" for shape in [[3,2,1],[3,2],[10,10,10]]
   @test size(LinearAlgebraicRepresentation.larImageVerts(shape)) == (length(shape),prod(shape + 1))
end

@testset "LarCuboids Tests" begin
   @testset "LarCuboids Tests" "$shape" for shape in [[3,2,1],[1,1,1],[3,3,10]]
      cubes = LinearAlgebraicRepresentation.larCuboids(shape,true);
      verts, cells = cubes;
      VV,EV,FV,CV = cells;
      @test size(LinearAlgebraicRepresentation.larCuboids(shape)[1],2) == prod(collect(shape) + 1)
      @test length(LinearAlgebraicRepresentation.larCuboids(shape)[2]) == prod(shape)
      @test typeof(verts) == Array{Float64,2}
	  @test typeof(VV) == Array{Array{Int64,1},1}
	  @test typeof(EV) == Array{Array{Int64,1},1}
	  @test typeof(FV) == Array{Array{Int64,1},1}
	  @test typeof(FV) == Array{Array{Int64,1},1}
   end

   @testset "LarCuboids Tests" "$shape" for shape in [[3,2],[1,1],[3,10]]
      @test size(LinearAlgebraicRepresentation.larCuboids(shape)[1],2) == prod(collect(shape) + 1)
      @test length(LinearAlgebraicRepresentation.larCuboids(shape)[2]) == prod(shape)
      cubes = LinearAlgebraicRepresentation.larCuboids(shape,true)
      verts, cells = cubes
      @test typeof(verts) == Array{Float64,2}
      VV,EV,FV = cells
	  @test typeof(VV) == Array{Array{Int64,1},1}
	  @test typeof(EV) == Array{Array{Int64,1},1}
	  @test typeof(FV) == Array{Array{Int64,1},1}
   end

   @testset "LarCuboids Tests" "$shape" for shape in [[3],[1],[10]]
      @test size(LinearAlgebraicRepresentation.larCuboids(shape)[1],2) == prod(collect(shape) + 1)
      @test length(LinearAlgebraicRepresentation.larCuboids(shape)[2]) == prod(shape)
      cubes = LinearAlgebraicRepresentation.larCuboids(shape,true)
      verts, cells = cubes
      @test typeof(verts) == Array{Float64,2}
      VV,EV = cells
	  @test typeof(VV) == Array{Array{Int64,1},1}
	  @test typeof(EV) == Array{Array{Int64,1},1}
   end
end

@testset "LinearAlgebraicRepresentation.LarModelProduct Tests" begin
   geom_0 = hcat([[x] for x=0.:10]...)
   topol_0 = [[i,i+1] for i=1:9]
   geom_1 = hcat([[0.],[1.],[2.]]...)
   topol_1 = [[1,2],[2,3]]
   model_0 = (geom_0, topol_0)
   model_1 = (geom_1, topol_1)
   model_2 = LinearAlgebraicRepresentation.larModelProduct(model_0, model_1)
   model_3 = LinearAlgebraicRepresentation.larModelProduct(model_2, model_1)
   @test size(model_2[1],2) == size(model_0[1],2) * size(model_1[1],2)
   @test length(model_2[2]) == length(model_0[2]) * length(model_1[2])
   @test size(model_3[1],2) == size(model_2[1],2) * size(model_1[1],2)
   @test length(model_3[2]) == length(model_2[2]) * length(model_1[2])
   V,EV = model_0
   @test typeof(V) == Array{Float64,2}
   @test typeof(EV) == Array{Array{Int64,1},1}
   V,EV = model_1
   @test typeof(V) == Array{Float64,2}
   @test typeof(EV) == Array{Array{Int64,1},1}
   V,FV = model_2
   @test typeof(V) ==  Array{Float64,2}
   @test typeof(FV) == Array{Array{Int64,1},1}
   V,CV = model_3
   @test typeof(V) ==  Array{Float64,2}
   @test typeof(CV) == Array{Array{Int64,1},1}

   modelOne = ([0 0 1 1; 0 1 0 1], Array{Int64,1}[[1, 2, 3, 4]])
   modelTwo = ([0 1], Array{Int64,1}[[1, 2]])
   @test LinearAlgebraicRepresentation.larModelProduct(modelOne, modelTwo) == LinearAlgebraicRepresentation.larModelProduct(modelTwo, modelOne)
   @test LinearAlgebraicRepresentation.larModelProduct(modelTwo, modelTwo)[1] == [0  0  1  1; 0  1  0  1]
   @test LinearAlgebraicRepresentation.larModelProduct(modelOne,modelOne)[2][1] == Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
end


text:
module TestMMTF

using Test
using MMTF

@testset "Codecs - Decoding Test" begin
    #runlengthdecode
    test_data = Int32[15,3,100,2,111,4,10000,6]
    output_data = Int32[15,15,15,100,100,111,111,111,111,10000,10000,10000,10000,10000,10000]
    eval_data = MMTF.runlengthdecode(test_data)
    @test eval_data==output_data

    #deltadecode
    test_data = Int32[15,3,100,-1,11,4]
    output_data = Int32[15,18,118,117,128,132]
    eval_data = MMTF.deltadecode(test_data)
    @test eval_data==output_data

    #deltarecursivefloatdecode
    test_data =  UInt8[0x7f,0xff,0x44,0xab,0x01,0x8f,0xff,0xca]
    output_data = Float32[50.346, 50.745, 50.691]
    eval_data = MMTF.deltarecursivefloatdecode(test_data,1000)
    @test eval_data==output_data

    #runlengthfloatdecode
    test_data = UInt8[ 0x00,0x00,0x00,0x64,0x00,0x00,0x00,0x03]
    output_data = Float32[1.00,1.00,1.00]
    eval_data = MMTF.runlengthfloatdecode(test_data,100)
    @test eval_data==output_data

    #runlengthdeltaintdecode
    test_data = UInt8[ 0x00,0x00,0x00,0x01,0x00,0x00,0x00,0x07]
    output_data = Int32[1,2,3,4,5,6,7]
    eval_data = MMTF.runlengthdeltaintdecode(test_data)
    @test eval_data==output_data

    #runlengthchardecode
    test_data = UInt8[ 0x00,0x00,0x00,0x41,0x00,0x00,0x00,0x04]
    output_data = Char['A','A','A','A']
    eval_data = MMTF.runlengthchardecode(test_data)
    @test eval_data==output_data

    #stringdecode
    test_data = UInt8[0x42, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x43, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00]
    output_data = String["B","A","C","A","A","A"]
    eval_data = MMTF.stringdecode(test_data)
    @test eval_data==output_data

    #intdecode
    test_data = UInt8[0x07, 0x06, 0x06, 0x07, 0x07]
    output_data = Int8[7,6,6,7,7]
    eval_data = MMTF.intdecode(test_data)
    @test eval_data==output_data

    #fourbyteintdecode
    test_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    output_data = Int32[1, 131073, 0, 2]
    eval_data = MMTF.fourbyteintdecode(test_data)
    @test eval_data==output_data
end

@testset "Codecs - Encoding Test" begin
    #runlengthencode
    output_data = Int32[15,3,100,2,111,4,10000,6]
    test_data = Int32[15,15,15,100,100,111,111,111,111,10000,10000,10000,10000,10000,10000]
    eval_data = MMTF.runlengthencode(test_data)
    @test eval_data==output_data

    #deltaencode
    output_data = Int32[15,3,100,-1,11,4]
    test_data = Int32[15,18,118,117,128,132]
    eval_data = MMTF.deltaencode(test_data)
    @test eval_data==output_data

    #deltarecursivefloatencode
    output_data = UInt8[0x7f, 0xff, 0x44, 0xab, 0x01, 0x8f, 0xff, 0xca]
    test_data = Float32[50.346, 50.745, 50.691]
    eval_data = MMTF.deltarecursivefloatencode(test_data,1000)
    @test eval_data==output_data

    #runlengthfloatencode
    output_data = UInt8[0x00, 0x00, 0x00, 0x64, 0x00, 0x00, 0x00, 0x03]
    test_data = Float32[1.00,1.00,1.00]
    eval_data = MMTF.runlengthfloatencode(test_data,100)
    @test eval_data==output_data

    #runlengthdeltaintencode
    output_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x07]
    test_data = Int32[1,2,3,4,5,6,7]
    eval_data = MMTF.runlengthdeltaintencode(test_data)
    @test eval_data==output_data

    #runlengthcharencode
    output_data = UInt8[0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x04]
    test_data = Char['A','A','A','A']
    eval_data = MMTF.runlengthcharencode(test_data)
    @test eval_data==output_data

    #stringencode
    output_data = UInt8[0x42, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x43, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00]
    test_data = String["B","A","C","A","A","A"]
    eval_data = MMTF.stringencode(test_data)
    @test eval_data==output_data

    #intencode
    output_data = UInt8[0x07, 0x06, 0x06, 0x07, 0x07]
    test_data = Int8[7,6,6,7,7]
    eval_data = MMTF.intencode(test_data)
    @test eval_data==output_data

    #fourbyteintencode
    output_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    test_data = Int32[1, 131073, 0, 2]
    eval_data = MMTF.fourbyteintencode(test_data)
    @test eval_data==output_data
end

@testset "Converters Test" begin
    #Bytes to 1 byte ints
    test_data = UInt8[0x07, 0x06, 0x06, 0x07, 0x07]
    output_data = Int8[7,6,6,7,7]
    eval_data = MMTF.ints(test_data,1)
    @test eval_data==output_data

    #Bytes to 2 byte ints
    test_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    output_data = Int16[0,1,2,1,0,0,0,2]
    eval_data = MMTF.ints(test_data,2)
    @test eval_data==output_data

    #Bytes to 4 byte ints
    test_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    output_data = Int32[1, 131073, 0, 2]
    eval_data = MMTF.ints(test_data,4)
    @test eval_data==output_data

    #Floats to ints
    test_data = Float32[10.001,100.203,124.542]
    output_data = Int32[10001,100203,124542]
    eval_data = MMTF.ints(test_data,1000)
    @test eval_data==output_data

    #Ints to floats
    output_data= Float32[10.001,100.203,124.542]
    test_data = Int32[10001,100203,124542]
    eval_data = MMTF.floats(test_data,1000)
    @test eval_data==output_data

    #Ints to bytes
    output_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    test_data = Int32[1, 131073, 0, 2]
    eval_data = MMTF.bytes(test_data)
    @test eval_data==output_data

    #Ints to Chars
    test_data = Int32[66,63,67]
    output_data = Char['B', '?','C']
    eval_data = MMTF.chars(test_data)
    @test eval_data==output_data

    #decodechainlist
    test_data = UInt8[0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00]
    output_data = String["A", "A","A","A","A","A"]
    eval_data = MMTF.decodechainlist(test_data)
    @test eval_data==output_data

    #encodechainlist
    output_data = UInt8[0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00, 0x41, 0x00, 0x00, 0x00]
    test_data = String["A", "A","A","A","A","A"]
    eval_data = MMTF.encodechainlist(test_data)
    @test eval_data==output_data

    #recursiveindexdecode
    test_data = Int16[1,420,32767,0,120,-32768,0,32767,2000]
    output_data = Int32[1,420,32767,120,-32768,34767]
    eval_data = MMTF.recursiveindexdecode(test_data)
    @test eval_data==output_data

    #recursiveindexencode
    output_data = Int16[1,420,32767,0,120,-32768,0,32767,2000]
    test_data = Int32[1,420,32767,120,-32768,34767]
    eval_data = MMTF.recursiveindexencode(test_data)
    @test eval_data==output_data
end

@testset "Utils Test" begin
    #parseheader
    test_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    codec,length,param,bytearray = MMTF.parseheader(test_data)
    @test codec==1
    @test length==131073
    @test param==0
    @test size(bytearray,1)==4

    #addheader
    test_data = UInt8[0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02]
    codec=1
    length=131073
    param=0
    bytearray=UInt8[0x00, 0x00, 0x00, 0x02]
    eval_data = MMTF.addheader(bytearray,codec,length,param)
    @test eval_data==test_data
end

@testset "Round Trip Test" begin
id_list = [
    #
    "1a1q",
    # // Just added to check
    "9pti",
    # // An entity that has no chain
    "2ja5",
    # // A couple of examples of multiple disulpgide bonds being formed.
    "3zxw",
    "1nty",
    # // A weird residue case
    "2eax",
    # // A Deuterated Structure
    "4pdj",
    # // Weird bioassembly
    "4a1i",
    # // Multi model structure
    "1cdr",
    # // Another weird structure
    "3zyb",
    # //Standard structure
    "4cup",
    # // Weird NMR structure
    "1o2f",
    # // B-DNA structure
    "1bna",
    # // DNA structure
    "4y60",
    # // Sugar structure
    "1skm",
    # // Calpha atom is missing (not marked as calpha)
    "1lpv",
    # // NMR structure with multiple models - one of which has chain missing
    "1msh",
    # // No ATOM records just HETATM records (in PDB). Opposite true for MMCif. It's a D-Peptide.
    "1r9v",
    # // Biosynthetic protein
    "5emg",
    # // Micro heterogenity
    "4ck4",
    # // Ribosome
    "4v5a",
    # // Negative residue numbers
    "5esw",
    # // A tiny example case
    "3njw",
    # // A GFP example with weird seqres records
    "1ema"]

    for pdbid in id_list
        println("PDB ID - ",pdbid)
        data_in = MMTF.fetchmmtf(pdbid)
        MMTF.writemmtf(data_in,"$(pdbid).mmtf")
        data_rt = MMTF.parsemmtf("$(pdbid).mmtf")
        @test data_in==data_rt
    end
end
end


text:
function bin2dec(str)
  x = 0
  for c in str
    x *= 2
    x += Int(c)-Int('0')
  end
  return x
end

function stringbit2num(str)
  bias = 2^10 -1
  m = 1 + bin2dec(str[13:end]) / 2^52
  ex = bin2dec(str[2:12]) - bias
  s = bin2dec(str[1])
  if s==1
    s=-1
  else
    s=1
  end
  val = s * m * 2^ex
  return val
end

println(bitstring(4.0))
println(stringbit2num(bitstring(4.0)))

max_stars_repo_path: II rok/I semestr/ANM/Cwiczenia 01/zad7_2.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
# TODO:
# - use keyword argument to enable/disable showing indices?
# - parse the "status" flag for every node and edge? needed for `des_pipes` (design
#   problem); not sure of other use cases
# * for "epanet-data" objects, identify shutoff valves; might want to do that in
#   WaterModels?

"""
Build pygraphviz graph object from a WaterModels network dictionary parsed from an EPANET
file. If a `solution` dict is provided, it should be for the same time as that of the `data`
dict. The `layout` option equates to the layout functions of graphviz (dot, neato, etc.) and
any provided `args` will also be passed to graphviz. Provide `layout=neato` and `args=-n` to
use the provided graph coordinates for the layout.
"""
function build_graph(data::Dict{String, <:Any},
                     solution::Union{Nothing, Dict{String, <:Any}} = nothing;
                     layout::String="dot", args::String="")
    # presumes input is data dict from WaterModels.parse_file() OR WaterModels.parse_epanet()

    # TODO:
    # - check that input dictionary has the expected fields

    G = pgv.AGraph(strict = false, directed = true)

    # Add nodal data to the directed graph object.
    add_nodes!(G, data["node"])
    node_labels!(G, data["demand"])
    node_labels!(G, data["reservoir"])
    node_labels!(G, data["tank"])

    # Add link data to the directed graph object.
    max_pipe_diameter = maximum([x.second["diameter"] for x in data["pipe"]])
    _add_pipe_to_graph!.(Ref(G), values(data["pipe"]), max_pipe_diameter)
    _add_des_pipe_to_graph!.(Ref(G), values(data["des_pipe"]))
    _add_short_pipe_to_graph!.(Ref(G), values(data["short_pipe"]))
    _add_pump_to_graph!.(Ref(G), values(data["pump"]))
    _add_regulator_to_graph!.(Ref(G), values(data["regulator"]))
    _add_valve_to_graph!.(Ref(G), values(data["valve"]))

    # if solution dict provided
    if solution !== nothing
        add_solution!(G, data, solution)
    end

    # set the layout here--can be reused for displaying different solutions
    try
        G.layout(prog=layout, args=(args))
    catch
        G.layout(prog="dot")
        # by putting this warning after, it should display only if `G.draw()` was actually
        # successful using "dot"; if it fails for another reason, that error should be shown
        @warn "$layout $args is not supported; dot without args was used instead"
    end
    # save the layout information for the user to inspect later
    PyCall.set!(G."graph_attr", "layout_args", "$layout $args")
    # get(G.graph_attr, "layout_args") # this is how to inspect

    return G
end # function build_graph

"""
Update a graph for a different time. Updates reservoir elevation (aka head), demands, and
solution information (if provided).
"""
function update_graph!(G::PyCall.PyObject, data::Dict{String,Any},
                       solution::Union{Nothing, Dict{String, <:Any}} = nothing)
    # TODO:
    # - check that reservoir elevation actually changes for a network that has a reservoir
    #   pattern
    # - check that there aren't any other time-dependent information to update

    # update reservoir elevation
    for (key, res) in data["reservoir"]
        nodekey = string(res["node"])
        node = data["node"][nodekey]
        nodeobj = @pycall G."get_node"(nodekey)::PyObject
        # Color according to the node's elevation.
        elmin = parse(Float64, get(G.graph_attr, "elmin", 0.0))
        elmax = parse(Float64, get(G.graph_attr, "elmax", 1.0))
        clr = HSV(get(viridis, node["elevation"], (elmin, elmax)))
        clrstr = string(clr.h / 360) * " " * string(clr.s) * " " * string(clr.v)
        fntclr = clr.v < 0.6 ? "white" : "black"
        PyCall.set!(nodeobj.attr, "fillcolor", clrstr)
        PyCall.set!(nodeobj.attr, "fontcolor", fntclr)
    end

    # update demands
    node_labels!(G, data["demand"])

    # update solution
    if solution !== nothing
        add_solution!(G, data, solution)
    end
end

"""
Add nodes to the pygraphviz graph object, including node attributes for name label,
elevation, and coordinates
"""
function add_nodes!(G::PyCall.PyObject, nodes::Dict{String,Any})
    # Determine minimum and maximum elevations.
    elevation_min = minimum([x["elevation"] for x in values(nodes)])
    elevation_max = maximum([x["elevation"] for x in values(nodes)])

    # Determine minimum and maximum coordinate values.
    if all(map(x -> haskey(x, "coordinates"), values(nodes)))
        x_min = minimum([x["coordinates"][1] for x in values(nodes)])
        x_max = maximum([x["coordinates"][1] for x in values(nodes)])
        y_min = minimum([x["coordinates"][2] for x in values(nodes)])
        y_max = maximum([x["coordinates"][2] for x in values(nodes)])
    else
        x_min, x_max, y_min, y_max = Inf, -Inf, Inf, -Inf
    end

    # scale to use to achieve "good" relative positioning with graphviz output --
    # calculate this systematically from the number of nodes?
    #scale = 20.0 # works with exclamation points as part of `pos` argument
    scale = 2e3 # no `!` but `-n` flag

    # add elmin and elmax as graph attributes for use in creating a colorbar in
    # colorbar (or via write_visualization)
    PyCall.set!(G."graph_attr", "elmin", elevation_min)
    PyCall.set!(G."graph_attr", "elmax", elevation_max)

    # create and populate pygraphviz nodes
    for (key, node) in nodes
        # Get the label to be used for the node.
        label = _get_comp_label(node)

        # Color according to the node's elevation.
        clr = HSV(get(viridis, node["elevation"], (elevation_min, elevation_max)))

        # Convert clr to a string that is in a form suitable for graphviz.
        clrstr = string(clr.h / 360) * " " * string(clr.s) * " " * string(clr.v)

        # Change font color depending on the background color.
        fntclr = clr.v < 0.6 ? "white" : "black"

        if all(map(x -> haskey(x, "coordinates"), values(nodes)))
            # Get scaled version of the coordinates.
            x = scale * (node["coordinates"][1] - x_min) / (x_max - x_min)
            y = scale * (node["coordinates"][2] - y_min) / (y_max - y_min)

            # Exclamation point forces exact positioning using Neato. But probably better to
            # not have it and use "-n" argument instead with `draw` command
            #pos = string(x) * "," * string(y) * "!"
            pos = string(x) * "," * string(y)

            # Add the node to the graph object with specified attributes.
            G.add_node(node["index"], label = label, elevation = node["elevation"],
                pos = pos, style = "filled", fillcolor = clrstr, fontcolor = fntclr)
        else
            # Add the node to the graph object with specified attributes.
            G.add_node(node["index"], label = label, elevation = node["elevation"],
                style = "filled", fillcolor = clrstr, fontcolor = fntclr)
        end
    end
end

"""Add labels for junctions, tanks, and reservoirs."""
function node_labels!(G::PyCall.PyObject, nodes::Dict{String,Any})
    for (key, node) in nodes
        nodekey = node["node"]
        nodeobj = @pycall G."get_node"(nodekey)::PyObject
        name = node["name"]

        # add node index to the label if different from the name
        if name==string(nodekey)
            label = name
        else
            label = name*" ("*string(nodekey)*")"
        end

        # will the node type ever change from the original source_id description? some link
        # types do (pipes -> valves)
        node_type = node["source_id"][1]

        if node_type == "reservoir"
            PyCall.set!(nodeobj.attr, "label", "Rsvr\\n"*label)
            PyCall.set!(nodeobj.attr, "shape", "diamond")
        elseif node_type == "tank"
            PyCall.set!(nodeobj.attr, "label", "Tank\\n"*label)
            PyCall.set!(nodeobj.attr, "shape", "rectangle")
        else
            dem = node["flow_nominal"]

            if dem == 0
                PyCall.set!(nodeobj.attr, "label", label)
            else
                dem = @sprintf("%2.2g", dem)
                PyCall.set!(nodeobj.attr, "label", label*"\\nd: "*dem)
            end
        end
    end
end

function _get_comp_label(comp::Dict{String, <:Any})
    source_name = comp["source_id"][2]

    if string(comp["index"]) == source_name
        return source_name
    else
        return "$(source_name) ($(string(comp["index"])))"
    end
end

function _get_link_arrowhead(link::Dict{String, <:Any})
    return "normal"
end

function _get_link_dir(link::Dict{String, <:Any})
    if link["flow_direction"] == _WM.UNKNOWN
        return "none"
    elseif link["flow_direction"] == _WM.POSITIVE
        return "forward"
    elseif link["flow_direction"] == _WM.NEGATIVE
        return "back"
    end
end

function _add_pipe_to_graph!(graph::PyCall.PyObject, pipe::Dict{String, <:Any}, max_diameter::Float64)
    index, label = string(pipe["index"]), _get_comp_label(pipe)
    penwidth = pipe["diameter"] * inv(max_diameter) * 10.0
    if penwidth > 8.0
        pad = "  "
    elseif penwidth > 4.0
        pad = " "
    else
        pad = ""
    end
    dir, arrowhead = _get_link_dir(pipe), _get_link_arrowhead(pipe)
    length = @sprintf("%.5g m", pipe["length"])
    # save `pad` for use with adding solution labels
    graph.add_edge(pipe["node_fr"], pipe["node_to"], index, dir = dir, headclip = "true",
                   arrowhead = arrowhead, penwidth = penwidth, pad = pad,
                   label = "$(pad)$(label)\\n$(pad)$(length)")
end

function _add_des_pipe_to_graph!(graph::PyCall.PyObject, des_pipe::Dict{String, <:Any})
    _add_pipe_to_graph!(graph, des_pipe, 1.0)
end

function _add_short_pipe_to_graph!(graph::PyCall.PyObject, short_pipe::Dict{String, <:Any})
    index, label = string(short_pipe["index"]), _get_comp_label(short_pipe)
    dir, arrowhead = _get_link_dir(short_pipe), _get_link_arrowhead(short_pipe)
    graph.add_edge(short_pipe["node_fr"], short_pipe["node_to"], index,
        dir = dir, arrowhead = arrowhead, label = "$(label)")
end

function _add_pump_to_graph!(graph::PyCall.PyObject, pump::Dict{String, <:Any})
    index, label = string(pump["index"]), _get_comp_label(pump)
    dir, arrowhead = _get_link_dir(pump), _get_link_arrowhead(pump)
    graph.add_edge(pump["node_fr"], pump["node_to"], index, dir = dir,
        arrowhead = arrowhead, label = "Pmp\\n$(label)", color = "red", style = "bold")
end

function _add_regulator_to_graph!(graph::PyCall.PyObject, regulator::Dict{String, <:Any})
    index, label = string(regulator["index"]), _get_comp_label(regulator)
    dir, arrowhead = _get_link_dir(regulator), _get_link_arrowhead(regulator)
    graph.add_edge(regulator["node_fr"], regulator["node_to"], index, dir = dir,
        arrowhead = arrowhead, label = "Reg\\n$(label)", color = "purple", style = "bold")
end

function _add_valve_to_graph!(graph::PyCall.PyObject, valve::Dict{String, <:Any})
    index, label = string(valve["index"]), _get_comp_label(valve)
    dir, arrowhead = _get_link_dir(valve), _get_link_arrowhead(valve)
    if dir != "none"
        label = "CV\\n$(label)"
    else
        label = "SV\\n$(label)"
    end
    graph.add_edge(valve["node_fr"], valve["node_to"], index, dir = dir,
        arrowhead = arrowhead, label = label, color = "blue", style = "bold")
end

"""
add solution values to the node and link labels
"""
function add_solution!(G::PyCall.PyObject, data::Dict{String,Any},
                       solution::Dict{String,Any})
    # add head to the node labels (could alternatively show pressure)
    for (key,node) in solution["node"]
        head = @sprintf("%2.2g", node["h"])
        nodeobj = @pycall G."get_node"(key)::PyObject
        label = get(nodeobj.attr, "label")
        # remove existing solution string (if it exists)
        idx = findlast("\\n", label)
        if idx != nothing
            idx = idx[1]
            #if contains(label[idx:end], "h:") # requires Julia >=1.5
            if occursin("h:", label[idx:end])
                label = label[1:idx-1]
            end
        end
        PyCall.set!(nodeobj.attr, "label", label*"\\nh: "*head)
    end
    # add flow to the labels for pipes and valves
    # Byron added this set, but not all of these fields exist in the solution object
    #pipesplus = ["pipe", "des_pipe", "pump", "regulator", "short_pipe", "valve"]
    links = ["pump", "pipe", "short_pipe", "valve"]
    for linktype in links
        for (key,pipesol) in solution[linktype]
            flow = _val_string_cut(pipesol["q"], 1e-10)
            link = data[linktype][key]
            # may also need to use `key` if multiple pipes between nodes
            edgeobj = @pycall G."get_edge"(link["node_fr"], link["node_to"])::PyObject
            label = get(edgeobj.attr, "label")
            # remove existing solution string (if it exists)
            idx = findlast("\\n", label)
            if idx != nothing
                idx = idx[1]
                #if contains(label[idx:end], "q:") # requires Julia >=1.5
                if occursin("q:", label[idx:end])
                    label = label[1:idx-1]
                end
            end
            pad = get(edgeobj.attr, "pad")
            PyCall.set!(edgeobj.attr, "label", label*"\\n$(pad)q: "*flow)
        end
    end
end

function _val_string_cut(val::Real, cut::Real)
    if val < cut
        return "0"
    else
        return @sprintf("%2.2g", val)
    end
end