text:
export produce_or_load, tagsave, @tagsave, safesave

"""
    produce_or_load([path="",] c, f; kwargs...) -> file, s
Let `s = joinpath(path, savename(prefix, c, suffix))`.
If a file named `s` exists then load it and return it, along
with the global path that it is saved at (`s`).

If the file does not exist then call `file = f(c)`, with `f` your function
that produces your data. Then save `file` as `s` and then return `file, s`.
The function `f` must return a dictionary,
the macros [`@dict`](@ref) and [`@strdict`](@ref) can help with that.

## Keywords
* `tag = true` : Save the file using [`tagsave`](@ref).
* `gitpath = projectdir()` : Path to search for a Git repo.
* `suffix = "bson", prefix = default_prefix(c)` : Used in `savename`.
* `force = false` : If `true` then don't check if file `s` exists and produce
  it and save it anyway.
* `loadfile = true` : If `false`, this function does not actually load the
  file, but only checks if it exists. The return value in this case is always
  `nothing, s`, regardless of whether the file exists or not. If it doesn't
  exist it is still produced and saved.
* `verbose = true` : print info about the process, if the file doesn't exist.
* `kwargs...` : All other keywords are propagated to `savename`.

See also [`savename`](@ref).
"""
produce_or_load(c, f; kwargs...) = produce_or_load("", c, f; kwargs...)
function produce_or_load(path::String, c, f;
    tag::Bool = true, gitpath = projectdir(), loadfile = true,
    suffix = "bson", prefix = default_prefix(c),
    force = false, verbose = true, kwargs...)

    s = joinpath(path, savename(prefix, c, suffix; kwargs...))

    if !force && isfile(s)
        if loadfile
            file = wload(s)
            return file, s
        else
            return nothing, s
        end
    else
        if force
            verbose && @info "Producing file $s now..."
        else
            verbose && @info "File $s does not exist. Producing it now..."
        end
        file = f(c)
        try
            if tag
                tagsave(s, file; safe = false, gitpath = gitpath)
            else
                wsave(s, copy(file))
            end
            verbose && @info "File $s saved."
        catch er
            @warn "Could not save file. Error stacktrace:"*
            Base.showerror(stderr, er, stacktrace(catch_backtrace()))
        end
        if loadfile
            return file, s
        else
            return nothing, s
        end
    end
end

################################################################################
#                             tag saving                                       #
################################################################################
"""
    tagsave(file::String, d::Dict; safe = false, gitpath = projectdir(), storepatch = true, force = false)
First [`tag!`](@ref) dictionary `d` and then save `d` in `file`.
If `safe = true` save the file using [`safesave`](@ref).

"Tagging" means that when saving the dictionary, an extra field
`:gitcommit` is added to establish reproducibility of results using
Git. If the Git repository is dirty, one more field `:gitpatch` is
added that stores the difference string.  If a dictionary already
contains a key `:gitcommit`, it is not overwritten, unless,
`force=true`. For more details, see [`tag!`](@ref).
"""
function tagsave(file, d; safe::Bool = false, gitpath = projectdir(), storepatch = true, force = false, source = nothing)
    d2 = tag!(d, gitpath=gitpath, storepatch=storepatch, force=force, source=source)
    if safe
        safesave(file, copy(d2))
    else
        wsave(file, copy(d2))
    end
    return d2
end

"""
    @tagsave(file::String, d::Dict; kwargs...)
Same as [`tagsave`](@ref) but one more field `:script` is added that records
the local path of the script and line number that called `@tagsave`, see [`@tag!`](@ref).
"""
macro tagsave(file,d,args...)
    s = QuoteNode(__source__)
    N = length(args)
    (N == 0 || all(iskwdefinition.(args))) &&
        return :(tagsave($(esc(file)), $(esc(d)),$(esc.(convert_to_kw.(args))...),source=$s))
    # First optional arg. is not needed as if it's not provided dispatch
    # is done through the kw-version of the function (ie. this line is
    # never reached)
    default = [
        :(projectdir()), #gitpath
        :(true),         #storepatch
    ]
    :(tagsave($(esc(file)),$(esc(d)), $(esc.(args)...), $(esc.(default[N:end])...), $s))
    # Use this code after the deprecation warning for the non-kw version
    # is removed.
    # throw(MethodError(@tagsave,args...))
end

################################################################################
#                          Backup files before saving                          #
################################################################################

# Implementation inspired by behavior of GROMACS
"""
    safesave(filename, data)

Safely save `data` in `filename` by ensuring that no existing files
are overwritten. Do this by renaming already existing data with a backup-number
ending like `#1, #2, ...`. For example if `filename = test.bson`, the first
time you `safesave` it, the file is saved normally. The second time
the existing save is renamed to `test_#1.bson` and a new file `test.bson`
is then saved.

If a backup file already exists then its backup-number is incremented
(e.g. going from `#2` to `#3`). For example safesaving `test.bson` a third time
will rename the old `test_#1.bson` to `test_#2.bson`, rename the old
`test.bson` to `test_#1.bson` and then save a new `test.bson` with the latest
`data`.

See also [`tagsave`](@ref).
"""
function safesave(f, data)
    recursively_clear_path(f)
    wsave(f, data)
end

#take a path of a results file and increment its prefix backup number
function increment_backup_num(filepath)
    path, filename = splitdir(filepath)
    fname, suffix = splitext(filename)
    m = match(r"^(.*)_#([0-9]+)$", fname)
    if m == nothing
        return joinpath(path, "$(fname)_#1$(suffix)")
    end
    newnum = string(parse(Int, m.captures[2]) +1)
    return joinpath(path, "$(m.captures[1])_#$newnum$(suffix)")
end

#recursively move files to increased backup number
function recursively_clear_path(cur_path)
    isfile(cur_path) || return
    new_path=increment_backup_num(cur_path)
    if isfile(new_path)
        recursively_clear_path(new_path)
    end
    mv(cur_path, new_path)
end

################################################################################
#                    Compliment to dict_list: tmpsave                          #
################################################################################
export tmpsave
using Random
"""
    tmpsave(dicts::Vector{Dict} [, tmp]; kwargs...) -> r
Save each entry in `dicts` into a unique temporary file in the directory `tmp`.
Then return the list of file names (relative to `tmp`) that were used
for saving each dictionary.

`tmp` defaults to `projectdir("_research", "tmp")`.

See also [`dict_list`](@ref).

## Keywords
* `l = 8` : number of characters in the random string.
* `prefix = ""` : prefix each temporary name will have.
* `suffix = "bson"` : ending of the temporary names (no need for the dot).
"""
function tmpsave(dicts, tmp = projectdir("_research", "tmp");
    l = 8, suffix = "bson", prefix = "")

    mkpath(tmp)
    n = length(dicts)
    existing = readdir(tmp)
    r = String[]
    i = 0
    while i < n
        x = prefix*randstring(l)*"."*suffix
        while x ∈ r || x ∈ existing
            x = prefix*randstring(l)*"."*suffix
        end
        i += 1
        push!(r, x)
        wsave(joinpath(tmp, x), copy(dicts[i]))
    end
    r
end

@deprecate tagsave(file, d, p::String) tagsave(file, d, gitpath=p)
@deprecate tagsave(file, d, safe::Bool, gitpath = projectdir(), storepatch = true, source = nothing) tagsave(file,d,safe=safe,gitpath=gitpath,storepatch=storepatch,source=source)
# TODO: When removing the deprecation warning, the tests must be adapted
# to only use the kw-version of this function. Also the code from the
# macro version for parsing the non-kw arguments can be replaced.


text:
module GeoJSONEx
using DataFrames

export to_features, to_featurecollection

function to_features(data::DataFrame, geofield::Symbol)
    propertykeys = [k for k in names(data) if k != geofield]
    [Dict("type" => "Feature",
          "geometry" => data[geofield][i],
          "properties" => Dict(String(k) => data[k][i] for k in propertykeys),
    ) for i in 1:size(data, 1)]
end

function to_featurecollection(data::DataFrame, geofield::Symbol)
    Dict("type" => "FeatureCollection",
        "features" => to_features(data, geofield))
end

end

ext: jl

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
"""
    AbstractPathState

An abstract type that provides information from shortest paths calculations.
"""
abstract type AbstractPathState end

"""
    is_ordered(e)
Return true if the source vertex of edge `e` is less than or equal to
the destination vertex.
"""
is_ordered(e::AbstractEdge) = src(e) <= dst(e)

"""
    add_vertices!(g, n)
Add `n` new vertices to the graph `g`.
Return the number of vertices that were added successfully.
"""
add_vertices!(g::AbstractGraph, n::Integer) = sum([add_vertex!(g) for i = 1:n])

"""
    indegree(g[, v])

Return a vector corresponding to the number of edges which end at each vertex in
graph `g`. If `v` is specified, only return degrees for vertices in `v`.
"""
indegree(g::AbstractGraph, v::Integer) = length(inneighbors(g, v))
indegree(g::AbstractGraph, v::AbstractVector = vertices(g)) = [indegree(g, x) for x in v]

"""
    outdegree(g[, v])

Return a vector corresponding to the number of edges which start at each vertex in
graph `g`. If `v` is specified, only return degrees for vertices in `v`.
"""
outdegree(g::AbstractGraph, v::Integer) = length(outneighbors(g, v))
outdegree(g::AbstractGraph, v::AbstractVector = vertices(g)) = [outdegree(g, x) for x in v]

"""
    degree(g[, v])
Return a vector corresponding to the number of edges which start or end at each
vertex in graph `g`. If `v` is specified, only return degrees for vertices in `v`.
For directed graphs, this value equals the incoming plus outgoing edges.
For undirected graphs, it equals the connected edges.
"""
function degree end
@traitfn degree(g::::IsDirected, v::Integer) = indegree(g, v) + outdegree(g, v)
@traitfn degree(g::::(!IsDirected), v::Integer) = indegree(g, v)

degree(g::AbstractGraph, v::AbstractVector = vertices(g)) = [degree(g, x) for x in v]

"""
    Δout(g)

Return the maximum [`outdegree`](@ref) of vertices in `g`.
"""
Δout(g) = noallocextreme(outdegree, (>), typemin(Int), g)
"""
    δout(g)

Return the minimum [`outdegree`](@ref) of vertices in `g`.
"""
δout(g) = noallocextreme(outdegree, (<), typemax(Int), g)

"""
    Δin(g)

Return the maximum [`indegree`](@ref) of vertices in `g`.
"""
Δin(g) = noallocextreme(indegree, (>), typemin(Int), g)

"""
    δin(g)

Return the minimum [`indegree`](ref) of vertices in `g`.
"""
δin(g) = noallocextreme(indegree, (<), typemax(Int), g)

"""
    Δ(g)

Return the maximum [`degree`](@ref) of vertices in `g`.
"""
Δ(g) = noallocextreme(degree, (>), typemin(Int), g)

"""
    δ(g)
Return the minimum [`degree`](@ref) of vertices in `g`.
"""
δ(g) = noallocextreme(degree, (<), typemax(Int), g)

"""
    noallocextreme(f, comparison, initial, g)
Compute the extreme value of `[f(g,i) for i=i:nv(g)]` without gathering them all
"""
function noallocextreme(f, comparison, initial, g)
    value = initial
    for i in vertices(g)
        funci = f(g, i)
        if comparison(funci, value)
            value = funci
        end
    end
    return value
end

"""
    degree_histogram(g, degfn=degree)

Return a Dict with values representing the number of vertices that have degree
represented by the key.

Degree function (for example, `indegree` or `outdegree`) may be specified by
overriding `degfn`.
"""
function degree_histogram(g::AbstractGraph{T}, degfn=degree) where T
    hist = Dict{T,Int}()
    for v in vertices(g)        # minimize allocations by
        for d in degfn(g, v)    # iterating over vertices
            hist[d] = get(hist, d, 0) + 1
        end
    end
    return hist
end

"""
    neighbors(g, v)

Return a list of all neighbors reachable from vertex `v` in `g`.
For directed graphs, the default is equivalent to [`outneighbors`](@ref);
use [`all_neighbors`](@ref) to list inbound and outbound neighbors.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
neighbors(g::AbstractGraph, v::Integer) = outneighbors(g, v)

"""
    all_neighbors(g, v)
Return a list of all inbound and outbound neighbors of `v` in `g`.
For undirected graphs, this is equivalent to both [`outneighbors`](@ref)
and [`inneighbors`](@ref).

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
function all_neighbors end
@traitfn all_neighbors(g::::IsDirected, v::Integer) =
    union(outneighbors(g, v), inneighbors(g, v))
@traitfn all_neighbors(g::::(!IsDirected), v::Integer) =
    neighbors(g, v)

"""
    common_neighbors(g, u, v)

Return the neighbors common to vertices `u` and `v` in `g`.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
common_neighbors(g::AbstractGraph, u::Integer, v::Integer) =
    intersect(neighbors(g, u), neighbors(g, v))

"""
    has_self_loops(g)

Return true if `g` has any self loops.
"""
has_self_loops(g::AbstractGraph) = nv(g) == 0 ? false : any(v -> has_edge(g, v, v), vertices(g))

"""
    num_self_loops(g)

Return the number of self loops in `g`.
"""
num_self_loops(g::AbstractGraph) = nv(g) == 0 ? 0 : sum(v -> has_edge(g, v, v), vertices(g))

"""
    density(g)
Return the density of `g`.
Density is defined as the ratio of the number of actual edges to the
number of possible edges (``|V|×(|V|-1)`` for directed graphs and
``\\frac{|V|×(|V|-1)}{2}`` for undirected graphs).
"""
function density end
@traitfn density(g::::IsDirected) =
ne(g) / (nv(g) * (nv(g) - 1))
@traitfn density(g::::(!IsDirected)) =
(2 * ne(g)) / (nv(g) * (nv(g) - 1))

"""
    squash(g)

Return a copy of a graph with the smallest practical type that
can accommodate all vertices.
"""
function squash(g::AbstractGraph)
    gtype = is_directed(g) ? DiGraph : Graph
    validtypes = [UInt8, UInt16, UInt32, UInt64, Int]
    nvg = nv(g)
    for T in validtypes
        nvg < typemax(T) && return gtype{T}(g)
    end
end

"""
    weights(g)

Return the weights of the edges of a graph `g` as a matrix. Defaults
to [`LightGraphs.DefaultDistance`](@ref).
"""
weights(g::AbstractGraph) = DefaultDistance(nv(g))


text:
## Deprecated in v0.2

function compute_message end
function update_node end
function update_edge end

compute_message(l, xi, xj, e) = compute_message(l, xi, xj)
update_node(l, x, m̄) = m̄
update_edge(l, e, m) = e

function propagate(l::GNNLayer, g::GNNGraph, aggr, x, e=nothing)
    @warn """
          Passing a GNNLayer to propagate is deprecated,
          you should pass the message function directly.
          The new signature is `propagate(f, g, aggr; [xi, xj, e])`.

          The functions `compute_message`, `update_node`,
          and `update_edge` have been deprecated as well. Please
          refer to the documentation.
          """
    m = apply_edges((a...) -> compute_message(l, a...), g, x, x, e)
    m̄ = aggregate_neighbors(g, aggr, m)
    x = update_node(l, x, m̄)
    e = update_edge(l, e, m)
    return x, e
end

## Deprecated in v0.3

@deprecate copyxj(xi, xj, e)  copy_xj(xi, xj, e)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# Mapping NDArray functions to Base-like API

const _ndsig = Dict{Symbol,Expr}()
const _nddoc = Dict{Symbol,Any}()

_isinplace(name::Symbol) = endswith(string(name), "!")

_writable(name::Symbol, x) =
  _isinplace(name) ? :(@assert $x.writable "this NDArray isn't writable") : :()

function _outexpr(name::Symbol, x #= the first arg of `sig` =#)
  if _isinplace(name)  # `func!`
    Ptr, 1, :([[MX_handle(x.handle)]]), :($x)
  else
    retexpr = :(NDArray(MX_NDArrayHandle(unsafe_load(hdls_ref[], 1))))
    Ref, 0, :(Ref{Ptr{MX_handle}}(C_NULL)), retexpr
  end
end

_broadcast_target(sig::Expr) = sig.args[2].args[].args[end]

"""
Generate docstring from function signature
"""
function _docsig(fname::Symbol, sig::Expr, opname::String)
  if fname !== :broadcasted
    get(_nddoc, fname, "    $sig") * "\n" * _getdocdefine(opname)
  else
    name = _broadcast_target(sig)
    str = get(_nddoc, name, "")
    _nddoc[name] = false  # change to false, denote docstring has been set up
    if isempty(str)
      sig_ = Expr(:call, Symbol(name, "."), sig.args[3:end]...)
      str = "    $sig_"
    end
    if str ≠ false
      # append "Defined in ..."
      def = _getdocdefine(opname)
      str = if str isa Markdown.MD
        str = Markdown.MD(copy(str.content), copy(str.meta))
        push!(str, Markdown.Paragraph(def))
        str
      else
        str * def
      end

      @eval @doc $str $name
    end
    ""
  end
end

"""
    @_remap(sig::Expr, imp::Expr)

Creating a function in signature `sig` with the function implementation `imp`.

## Arguments
- `sig` is the function signature.
  If the function name ends with `!`, it will invoke the corresponding inplace
  call.
- `imp` is the underlying libmxnet API call

"""
macro _remap(sig::Expr, imp::Expr)
  d = splitdef(:($sig = $imp))
  @capture d[:name] (M_.fname_|fname_)

  opname = string(imp.args[1])

  if isa(imp.args[2], Expr) && imp.args[2].head == :parameters
    ndin = imp.args[3:end]
    mxargs = imp.args[2].args
  else  # no keyword arguments
    ndin = imp.args[2:end]
    mxargs = []
  end

  mxkeys = map(x -> string(x.args[1]), mxargs)
  mxvals = Expr(:vect, map(x -> :(dump_mx_param($(x.args[2]))), mxargs)...)
  ndhlds = Expr(:vect, map(x -> :($(x).handle), ndin)...)

  # handler for `func!` which has side effect on first argument.
  T, n_output, hdls_ref, retexpr = _outexpr(fname, _firstarg(sig))

  assert_expr = _writable(fname, _firstarg(sig))

  func_body = quote
    $assert_expr
    op_handle = _get_cached_libmx_op_handle($opname)
    n_output = Ref(Cint($n_output))
    hdls_ref = $hdls_ref
    @mxcall(:MXImperativeInvoke,
            (MX_handle,
             Cint,
             Ptr{MX_handle},
             Ref{Cint},
             $T{Ptr{MX_handle}},
             Cint,
             char_pp,
             char_pp),
            op_handle,
            $(length(ndin)),
            $(ndhlds),
            n_output,
            hdls_ref,
            $(length(mxargs)),
            $mxkeys,
            $mxvals)
    $retexpr
  end

  docstr = _docsig(fname, sig, opname)
  func_def = Expr(:function, sig, func_body)

  esc(quote
    @doc $docstr
    $func_def
  end)
end

macro _remap(sig::Expr, imp::Symbol)
  imp = _ndsig[imp]

  esc(quote
    @_remap($sig, $imp)
  end)
end


text:
DIAGRAMS = Dict()
include("lotvol.jl")
include("orego.jl")
include("hires.jl")
include("pollution.jl")


text:
module MicroscopyTools

using Tullio
using FFTW
using NDTools
using Statistics
using IndexFunArrays

include("binning.jl")
include("calculation_tools.jl")
include("utils.jl")
include("grid_generation.jl")

end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
using Logging
using JuMP, Gurobi
using DecisionProgramming

const O = 1  # Chance node: lemon or peach
const T = 2  # Decision node: pay stranger for advice
const R = 3  # Chance node: observation of state of the car
const A = 4  # Decision node: purchase alternative
const O_states = ["lemon", "peach"]
const T_states = ["no test", "test"]
const R_states = ["no test", "lemon", "peach"]
const A_states = ["buy without guarantee", "buy with guarantee", "don't buy"]

@info("Creating the influence diagram.")
S = States([
    (length(O_states), [O]),
    (length(T_states), [T]),
    (length(R_states), [R]),
    (length(A_states), [A]),
])
C = Vector{ChanceNode}()
D = Vector{DecisionNode}()
V = Vector{ValueNode}()
X = Vector{Probabilities}()
Y = Vector{Consequences}()

I_O = Vector{Node}()
X_O = [0.2, 0.8]
push!(C, ChanceNode(O, I_O))
push!(X, Probabilities(O, X_O))

I_T = Vector{Node}()
push!(D, DecisionNode(T, I_T))

I_R = [O, T]
X_R = zeros(S[O], S[T], S[R])
X_R[1, 1, :] = [1,0,0]
X_R[1, 2, :] = [0,1,0]
X_R[2, 1, :] = [1,0,0]
X_R[2, 2, :] = [0,0,1]
push!(C, ChanceNode(R, I_R))
push!(X, Probabilities(R, X_R))

I_A = [R]
push!(D, DecisionNode(A, I_A))

I_V1 = [T]
Y_V1 = [0.0, -25.0]
push!(V, ValueNode(5, I_V1))
push!(Y, Consequences(5, Y_V1))

I_V2 = [A]
Y_V2 = [100.0, 40.0, 0.0]
push!(V, ValueNode(6, I_V2))
push!(Y, Consequences(6, Y_V2))

I_V3 = [O, A]
Y_V3 = [-200.0 0.0 0.0;
        -40.0 -20.0 0.0]
push!(V, ValueNode(7, I_V3))
push!(Y, Consequences(7, Y_V3))

validate_influence_diagram(S, C, D, V)
sort!.((C, D, V, X, Y), by = x -> x.j)

P = DefaultPathProbability(C, X)
U = DefaultPathUtility(V, Y)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, S, D)
x_s = PathCompatibilityVariables(model, z, S, P)
EV = expected_value(model, x_s, U, P)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)

@info("Printing decision strategy:")
print_decision_strategy(S, Z)

@info("Computing utility distribution.")
udist = UtilityDistribution(S, P, U, Z)

@info("Printing utility distribution.")
print_utility_distribution(udist)

@info("Printing expected utility.")
print_statistics(udist)


text:
using Allnix
function run_condition(
    n::Int64,
    x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64})
    @inbounds for i = 1:n
         if z[i] < 1500.
             y[i] = x[i]/13. + z[i] - 10.
         else
             y[i] = x[i]/13. + z[i] + 10.
         end
    end
end

function run_xyz(
    n::Int64,
    x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}
    )

    @inbounds for i = 1:n
         y[i] = x[i]/13. + z[i] - 10.
    end
end

@testset "test_xyz" begin
n = 1000000
count = 1000
a::Float64 = 0.5
x = Array{Float64, 1}(n)
y = Array{Float64, 1}(n)
z = Array{Float64, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 0.0
    z[i] = 1.5 * (i+1)
end
#val, t, bytes, gctime, memallocs = @timed run_xyz(1, n, x, y, z)
#println("test_xyz: $t, y = $(y[n])")
#val, t, bytes, gctime, memallocs = @timed run_xyz(count, n, x, y, z)

# warm up
val, t, bytes, gctime, memallocs = @timed run_xyz(n, x, y, z)
println("test_xyz: $t, y = $(y[n])")

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    run_xyz(n, x, y, z)
end
println("test_xyz: $t, y = $(y[n])")
end

@testset "test_condition" begin
n = 1000000
count = 1000
a::Float64 = 0.5
x = Array{Float64, 1}(n)
y = Array{Float64, 1}(n)
z = Array{Float64, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 0.0
    z[i] = 1.5 * (i+1)
end
#val, t, bytes, gctime, memallocs = @timed run_condition(1, n, x, y, z)
#println("test_condition: $t, y = $(y[n])")
#val, t, bytes, gctime, memallocs = @timed run_condition(count, n, x, y, z)

# warm up
val, t, bytes, gctime, memallocs = @timed run_condition(n, x, y, z)
println("test_condition: $t, y = $(y[n])")
val, t, bytes, gctime, memallocs = @timed for i in 1:count
    run_condition(n, x, y, z)
end
println("test_condition: $t, y = $(y[n])")
end

# Test daxpy
@testset "daxpy" begin
n = 1000000
count = 1000
a::Float64 = 0.5
x = Array{Float64, 1}(n)
y = Array{Float64, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 1.0
end

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    Allnix.daxpy(a,x,y)
end

println("daxpy time: $t, y = $(y[n])")
# my Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz takes less than one second
@test t <= 1.0
end

@testset "axpy" begin
# - Float64 axpy - #
@testset "Float64" begin
n = 1000000
count = 1000
a::Float64 = 0.5
x = Array{Float64, 1}(n)
y = Array{Float64, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 1.0
end

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    Allnix.axpy(a,x,y)
end

println("axpy{Float64} time: $t")
@test t <= 1.0
end
# - Float32 axpy - #
@testset "Float32" begin
n = 1000000
count = 1000
a::Float32 = 0.5
x = Array{Float32, 1}(n)
y = Array{Float32, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 1.0
end

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    Allnix.axpy(a,x,y)
end

println("axpy{Float32} time: $t")
# I guess this is slower than Float64 because Julia uses
# Float64 internally
@test t <= 3.0
end

# - Int64 axpy - #
@testset "Int64" begin
n = 1000000
count = 1000
a::Int64 = 2
x = Array{Int64, 1}(n)
y = Array{Int64, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 1
end

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    Allnix.axpy(a,x,y)
end

println("axpy{Int64} time: $t")
@test t <= 1.2
end
# - Int32 axpy - #
@testset "Int32" begin
n = 1000000
count = 1000
a::Int32 = 2
x = Array{Int32, 1}(n)
y = Array{Int32, 1}(n)
for i in 1:n
    x[i] = i - 1
    y[i] = 1
end

val, t, bytes, gctime, memallocs = @timed for i in 1:count
    Allnix.axpy(a,x,y)
end

println("axpy{Int32} time: $t")
@test t <= 1.0
end
end


text:
# test of calculateExpCF function
# Claudia December 2014

# Case G -----------------

#println("------ Case G ----------")
include("../examples/case_g_example.jl")
# include(joinpath(dirname(pathof(PhyloNetworks)),  "..","examples","case_g_example.jl"))
error1 = false
ind = 0

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q1);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 3 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 2 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    qnet.hybrid[1].prev.number != -7 ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != 0.2-log(1-0.1*(1-exp(-1.1))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 1")
    global error1 |= true
    global ind = 1
end

q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q2);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should have 1 hybrid nodes") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [0.9*(1-2/3*exp(-0.1))+0.1*1/3*exp(-1.0), 0.9*(1/3*exp(-0.1))+0.1*(1-2/3*exp(-1.0)), 0.9*(1/3*exp(-0.1))+0.1*(1/3*exp(-1.0))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 2")
    global error1 |= true
    global ind = 2
end

q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q3);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should have 1 hybrid nodes") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [0.9*(1/3*exp(-0.1))+0.1*1/3*exp(-1.0), 0.9*(1/3*exp(-0.1))+0.1*(1-2/3*exp(-1.0)), 0.9*(1-2/3*exp(-0.1))+0.1*(1/3*exp(-1.0))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 3")
    global error1 |= true
    global ind = 3
end

q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q4);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 > 0.30001 || qnet.t1 < 0.2999999999 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 4")
    global error1 |= true
    global ind = 4
end

q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q5);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 3 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 2 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    qnet.hybrid[1].prev.number != -3 ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != 0.2-log(1-0.1*(1-exp(-0.1))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 5")
    global error1 |= true
    global ind = 5
end

if error1
    throw("error Case G in quartet $(ind)")
end

# Case F Bad Diamond I -----------------

#println("------ Case F Bad diamond I ----------")
include("../examples/case_f_example.jl");
# include(joinpath(dirname(pathof(PhyloNetworks)),  "..","examples","case_f_example.jl"))
error1 = false
ind = 0
parameters!(net)

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q1);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [(1-(0.7*(1-exp(-0.2)))-(0.3*(1-exp(-0.1))))/3, (1+2*(0.7*(1-exp(-0.2)))-(0.3*(1-exp(-0.1))))/3,(1-(0.7*(1-exp(-0.2)))+2*(0.3*(1-exp(-0.1))))/3] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 1")
    global error1 |= true
    global ind = 1
end

q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q2);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    qnet.t1 != 0.1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [1,2,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1-2/3*(exp(-0.1)),1/3*(exp(-0.1)),1/3*(exp(-0.1))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 2")
    global error1 |= true
    global ind = 2
end

q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q3);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    qnet.t1 != 0.1-log(1-(0.3*(1-exp(-0.1)))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,2,1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 3")
    global error1 |= true
    global ind = 3
end

q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q4);

try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    qnet.t1 != 0.1-log(1-(0.7*(1-exp(-0.2)))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 4")
    global error1 |= true
    global ind = 4
end

q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q5);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [(1-0.7*(1-exp(-0.2))-0.3*(1-exp(-0.1)))/3,(1+2*0.7*(1-exp(-0.2))-0.3*(1-exp(-0.1)))/3,(1-0.7*(1-exp(-0.2))+2*0.3*(1-exp(-0.1)))/3] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 5")
    global error1 |= true
    global ind = 5
end

if error1
    throw("error Case F in quartet $(ind)")
end

# Case I Bad Diamond II -----------------

#println("------ Case I Bad diamond II ----------")
include("../examples/case_i_example.jl");
# include(joinpath(dirname(pathof(PhyloNetworks)),  "..","examples","case_i_example.jl"))
error1 = false
ind = 0

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q1);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [0.9*(1/3*exp(-1))+0.1*(1-2/3*exp(-1)),0.9*(1-2/3*exp(-1))+0.1*1/3*exp(-1),0.9*(1/3*exp(-1))+0.1*(1/3*exp(-1))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 1")
    global error1 |= true
    global ind = 1
end

q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q2);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 3 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 4 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    (qnet.hybrid[1].prev.number != -2 && qnet.hybrid[1].prev.number != -3) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    !approxEq(qnet.t1,-log(1+0.1*(1-exp(-1))-0.1*0.1*(1-exp(-1-1))-0.1*0.1*(1-exp(-1))-0.9*0.9*(1-exp(-2.)))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [1,2,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1-2/3*(exp(-qnet.t1)),1/3*(exp(-qnet.t1)),1/3*(exp(-qnet.t1))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 2")
    global error1 |= true
    global ind = 2
end

q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q3);
try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 3 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 4 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    (qnet.hybrid[1].prev.number != -6 && qnet.hybrid[1].prev.number != -3) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    !approxEq(qnet.t1,-log(1+0.1*(1-exp(-2))-0.1*0.1*(1-exp(-1))-0.1*0.1*(1-exp(-2))-0.9*0.9*(1-exp(-2.)))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,2,1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 3")
    global error1 |= true
    global ind = 3
end

q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q4);

try
    identifyQuartet!(qnet)
    qnet.which != 1 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 3 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 4 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    (qnet.hybrid[1].prev.number != -6 && qnet.hybrid[1].prev.number != -3) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 0 || qnet.numHybrids != 0 ? error("qnet should have 0 hybrid nodes") : nothing
    !approxEq(qnet.t1,-log(1+0.1*(1-exp(-1))-0.1*0.1*(1-exp(-1))-0.1*0.1*(1-exp(-1))-0.9*0.9*(1-exp(-3)))) ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [1,1,2,2] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [2,1,2] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [1/3*exp(-qnet.t1),1-2/3*exp(-qnet.t1),1/3*exp(-qnet.t1)] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 4")
    global error1 |= true
    global ind = 4
end

q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);
qnet = extractQuartet!(net,q5);
try
    identifyQuartet!(qnet)
    qnet.which != 2 ? error("qnet which not correctly assigned") : nothing
    qnet.hybrid[1].k != 4 ? error("qnet.hybrid[1].k not correctly assigned") : nothing
    qnet.hybrid[1].typeHyb != 5 ? error("qnet.hybrid[1].typeHyb not correctly assigned") : nothing
    !isa(qnet.hybrid[1].prev,Nothing) ? error("qnet.hybrid[1].prev not correctly assigned") : nothing

    eliminateHybridization!(qnet)
    size(qnet.hybrid,1) != 1 || qnet.numHybrids != 1 ? error("qnet should not have hybrid nodes anymore") : nothing
    qnet.t1 != -1 ? error("internal edge length not correctly updated") : nothing

    updateSplit!(qnet)
    qnet.split != [-1,-1,-1,-1] ? error("qnet.split not correctly assigned") : nothing

    updateFormula!(qnet)
    qnet.formula != [-1,-1,-1] ? error("qnet.formula not correctly assigned") : nothing

    calculateExpCF!(qnet)
    qnet.expCF != [0.9*(1/3*exp(-1))+0.1*(1-2/3*exp(-1)),0.9*(1-2/3*exp(-1))+0.1*1/3*exp(-1),0.9*(1/3*exp(-1))+0.1*(1/3*exp(-1))] ? error("qnet.expCF wrongly calculated") : nothing
catch
    println("errors in quartet 5")
    global error1 |= true
    global ind = 5
end

if error1
    throw("error Case I in quartet $(ind)")
end


text:
# Autogenerated wrapper script for QD_jll for x86_64-unknown-freebsd11.1-cxx11
export libqd_f_main, libqd, libqdmod

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libqd_f_main`
const libqd_f_main_splitpath = ["lib", "libqd_f_main.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqd_f_main_path = ""

# libqd_f_main-specific global declaration
# This will be filled out by __init__()
libqd_f_main_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqd_f_main = "libqd_f_main.so.0"

# Relative path to `libqd`
const libqd_splitpath = ["lib", "libqd.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqd_path = ""

# libqd-specific global declaration
# This will be filled out by __init__()
libqd_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqd = "libqd.so.0"

# Relative path to `libqdmod`
const libqdmod_splitpath = ["lib", "libqdmod.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqdmod_path = ""

# libqdmod-specific global declaration
# This will be filled out by __init__()
libqdmod_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqdmod = "libqdmod.so.0"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"QD")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global libqd_f_main_path = normpath(joinpath(artifact_dir, libqd_f_main_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqd_f_main_handle = dlopen(libqd_f_main_path)
    push!(LIBPATH_list, dirname(libqd_f_main_path))

    global libqd_path = normpath(joinpath(artifact_dir, libqd_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqd_handle = dlopen(libqd_path)
    push!(LIBPATH_list, dirname(libqd_path))

    global libqdmod_path = normpath(joinpath(artifact_dir, libqdmod_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqdmod_handle = dlopen(libqdmod_path)
    push!(LIBPATH_list, dirname(libqdmod_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
mutable struct DiffMethods
  closurell::Union{Function, Void}
  closurelp::Union{Function, Void}
  closurelt::Union{Function, Void}
  tapegll::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}
  tapeglp::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}
  tapeglt::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}
  tapetll::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}
  tapetlp::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}
  tapetlt::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}
end

DiffMethods(;
  closurell::Union{Function, Void}=nothing,
  closurelp::Union{Function, Void}=nothing,
  closurelt::Union{Function, Void}=nothing,
  tapegll::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}=nothing,
  tapeglp::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}=nothing,
  tapeglt::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape, Void}=nothing,
  tapetll::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}=nothing,
  tapetlp::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}=nothing,
  tapetlt::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape, Void}=nothing
) =
 DiffMethods(closurell, closurelp, closurelt, tapegll, tapeglp, tapeglt, tapetll, tapetlp, tapetlt)

mutable struct DiffState
  resultll::Union{DiffResults.DiffResult, Void}
  resultlp::Union{DiffResults.DiffResult, Void}
  resultlt::Union{DiffResults.DiffResult, Void}
  cfggll::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}
  cfgglp::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}
  cfgglt::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}
  cfgtll::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}
  cfgtlp::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}
  cfgtlt::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}
end

DiffState(;
  resultll::Union{DiffResults.DiffResult, Void}=nothing,
  resultlp::Union{DiffResults.DiffResult, Void}=nothing,
  resultlt::Union{DiffResults.DiffResult, Void}=nothing,
  cfggll::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}=nothing,
  cfgglp::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}=nothing,
  cfgglt::Union{ReverseDiff.GradientConfig, ForwardDiff.GradientConfig, Void}=nothing,
  cfgtll::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}=nothing,
  cfgtlp::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}=nothing,
  cfgtlt::Union{ReverseDiff.HessianConfig, ForwardDiff.HessianConfig, Void}=nothing
) =
  DiffState(resultll, resultlp, resultlt, cfggll, cfgglp, cfgglt, cfgtll, cfgtlp, cfgtlt)

==(z::DiffState, w::DiffState) = reduce(&, [getfield(z, n) == getfield(w, n) for n in fieldnames(DiffState)])

isequal(z::DiffState, w::DiffState) = reduce(&, [isequal(getfield(z, n), getfield(w, n)) for n in fieldnames(DiffState)])

mutable struct DiffOptions
  mode::Symbol
  order::Integer
  targets::Vector{Bool}
  chunksize::Integer
  compiled::Bool

  function DiffOptions(mode::Symbol, order::Integer, targets::Vector{Bool}, chunksize::Integer, compiled::Bool)
    @assert (mode == :reverse || mode == :forward) "Mode of automatic differentation must be :reverse or :forward, got $mode"
    @assert (order == 1 || order == 2) "Order of differentiation must be 1 or 2, got order=$order"
    @assert (length(targets) == 3) "Length of targets must be 3, got $length(targets)-length vector"
    @assert chunksize >= 0 "chunksize can not be negative, got chunksize=$chunksize"
    new(mode, order, targets, chunksize, compiled)
  end
end

DiffOptions(mode::Symbol, order::Integer) = DiffOptions(mode, order, fill(false, 3), 0, mode == :reverse ? true : false)

DiffOptions(;
  mode::Symbol=:reverse, order::Integer=1, targets::Vector{Bool}=fill(false, 3), chunksize::Integer=0, compiled::Bool=true
) =
  DiffOptions(mode, order, targets, chunksize, compiled)


text:
using RecyclePermanent
using JLD2

include("../models/bpModel.jl")

@load "data/bpData.jld" ys θ
model = makeBPModel(θ, ys)

N = 40 * model.n

println("Single-threaded")
rpio = RPIO{model.Tx}(N, model.n, 1)
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimateDiscrete(model, rpio)
@time biasedEstimateDiscrete(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimateDiscrete(model, rpio)
@time rpEstimateDiscrete(model, rpio)

println("Multi-threaded")
rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimateDiscrete(model, rpio)
@time biasedEstimateDiscrete(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimateDiscrete(model, rpio)
@time rpEstimateDiscrete(model, rpio)


text:
nums = map(str -> [parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)], readlines("inputs/6-10/07.txt"))[1]
using Statistics

function part1(nums)
  med = Statistics.median(nums)
  t = 0
  for n in nums
    t += abs(n-med)
  end
  t
end

sum_to(n) = n*(n+1) / 2

function part2(nums)
  best = typemax(Int)
  for n in minimum(nums):maximum(nums)
    t = 0
    for crab in nums
      t += sum_to(abs(n - crab))
    end
    best = t < best ? t : best
  end
  round(Int, best)
end

function part2faster(nums)
  med = floor(mean(nums) - 0.5)
  t = 0
  for n in nums
    t += sum_to(abs(n-med))
  end
  t
end

day7part2(nums) = map(n->(n->n*(n+1)/2)(n-floor(mean(nums))),nums)|>sum

sum_to(n) = n*(n+1) / 2


text:
input = open("../inputs/day02.txt") do f
    read(f, String)
end

struct Record
    from::Int64
    to::Int64
    chr::Char
    pwd::String
end

format_re = r"^(\d+)-(\d+) ([a-z]): ([a-z]+)$"

function parseLine(ln)
    # 1-3 a: abcde
    capts = match(format_re, ln).captures
    return Record(
        parse(Int64, capts[1]),
        parse(Int64, capts[2]),
        capts[3][1],
        capts[4],
    )

end

function count_valid(is_valid::Function, records::Array{Record})::Int64
    return sum([
        1 for record in records if is_valid(record)
    ])
end

records = [ parseLine(ln) for ln in split(input, "\n") ]

valid1 = count_valid(records) do record
    matches = count(c -> c == record.chr, record.pwd)
    record.from <= matches <= record.to
end

valid2 = count_valid(records) do record
    ok1 = record.pwd[record.from] == record.chr
    ok2 = record.pwd[record.to] == record.chr
    ok1 != ok2
end

println("Part 1: $(valid1)")
println("Part 2: $(valid2)")


text:
######################################################################
# gif.jl
#
# This file makes gifs of simulations
######################################################################

"""
Creates a gif.

`gif(m, x, f, p, num_steps=10, filename="out.gif")`

where
* `m` is a `SearchDomain`
* `x` is a `Vehicle`
* `f` is a subtype of `AbstractFilter`
* `p` is a subtype of `Policy`
"""
function gif(m::SearchDomain,
             x::Vehicle,
             f::AbstractFilter,
             p::Policy,
             tc::TerminationCondition=StepThreshold(10),
             filename="out.gif";
             seconds_per_step=0.5,
             show_mean=false,
             show_cov=false,
             show_path=false
            )

    frames = Frames(MIME("image/png"), fps=20)
    path_x = Float64[]
    path_y = Float64[]

    # Plot the original scene
    plot(m, f, x)
    title("\$t\$ = 0 s")
    push!(frames, gcf())
    push!(path_x, x.x)
    push!(path_y, x.y)
    close()

    # 20 is the fps
    frames_per_step = round(Int, seconds_per_step * 20)

    step_count = 0
    while !is_complete(f, tc, step_count)
        old_pose = (x.x, x.y, x.heading)
        o = observe(m,x)
        update!(f, x, o)
        a = action(m, x, o, f, p)
        act!(m,x,a)

        move_target!(m)

        new_pose = (x.x, x.y, x.heading)
        # Plot everything in between old pose and new pose
        dx = (new_pose[1] - old_pose[1]) / frames_per_step
        dy = (new_pose[2] - old_pose[2]) / frames_per_step
        dh = a[3] / frames_per_step
        for j = 1:frames_per_step
            #figure()
            # determine intermediate pose
            new_h = mod(old_pose[3] + dh, 360.0)
            old_pose = (old_pose[1] + dx, old_pose[2] + dy, new_h)

            push!(path_x, old_pose[1])
            push!(path_y, old_pose[2])
            plot(m, f, old_pose;show_mean=show_mean, show_cov=show_cov)
            if show_path
                plot(path_x, path_y)
            end
            title("\$t\$ = $step_count s")
            push!(frames, gcf())
            close()
        end
        step_count += 1
    end
    #write(filename, frames)
    write("temp.mp4", frames)
    write("temp.gif", frames)
end

function gif(m::SearchDomain,
             xarr::Vector{Vehicle},
             farr,
             parr,
             num_steps::Int=10,
             filename="out.gif";
             seconds_per_step=0.5,
             show_mean=false,
             show_cov=false,
             show_path=false,
             alpha=1.0
            )

	num_vehicles = length(xarr)
	obsarr = zeros(num_vehicles)
	new_poses = Array(Pose, num_vehicles)
	old_poses = Array(Pose, num_vehicles)
	diffs = Array(Pose, num_vehicles)

	frames = Frames(MIME("image/png"), fps=20)

	# Set up the paths
	x_paths = Array(Vector{Float64}, num_vehicles)
	y_paths = Array(Vector{Float64}, num_vehicles)
	for (xi,x) in enumerate(xarr)
		x_paths[xi] = [x.x]
		y_paths[xi] = [x.y]
	end

	# Plot the original scene
	plot(m, farr, xarr)
	push!(frames, gcf())
	close()

	# 20 is the fps
	frames_per_step = round(Int, seconds_per_step * 20)

	colors = ["b", "r"]
	if num_vehicles == 1
		colors = ["b"]
	elseif num_vehicles == 2
		colors = ["b", "r"]
	elseif num_vehicles == 3
		colors = ["b", "r", "g"]
	elseif num_vehicles == 4
		colors = ["b", "r", "g", "k"]
	end

	for i = 1:num_steps
		for (xi,x) in enumerate(xarr)
			f = farr[xi]
			p = parr[xi]
			old_pose = (x.x, x.y, x.heading)
			old_poses[xi] = old_pose
			o = observe(m,x)
			update!(f, x, o)
			a = action(m, x, o, f, p)
			act!(m,x,a)
			new_pose = (x.x, x.y, x.heading)
			new_poses[xi] = new_pose
			# Plot everything in between old pose and new pose
			dx = (new_pose[1] - old_pose[1]) / frames_per_step
			dy = (new_pose[2] - old_pose[2]) / frames_per_step
			dh = a[3] / frames_per_step
			diffs[xi] = (dx,dy,dh)
		end
		for j = 1:frames_per_step
			#figure()
			# determine intermediate pose
			for (xi, x) in enumerate(xarr)
				dx = diffs[xi][1]
				dy = diffs[xi][2]
				dh = diffs[xi][3]
				old_pose = old_poses[xi]
				new_h = mod(old_pose[3] + dh, 360.0)
				old_pose = (old_pose[1] + dx, old_pose[2] + dy, new_h)
				old_poses[xi] = (old_pose)

				push!(x_paths[xi], old_pose[1])
				push!(y_paths[xi], old_pose[2])
			end
			plot(m, farr, old_poses; show_mean=show_mean, show_cov=show_cov, alpha=alpha)
			if show_path
				for (xi, x) in enumerate(xarr)
					plot(x_paths[xi], y_paths[xi], colors[xi])
				end
			end
			push!(frames, gcf())
			close()
		end
	end
	write(filename, frames)
end

# TODO: find a place for this stuff.... probably not here

# really, this is for white sands results
function simgif(j::LocTuple, actions::Vector{Pose}, observations::Vector{Float64}, b::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false, show_path=false)
	# always start the vehicle in the center
	# NO, don't do that
	#x.x = m.length / 2.0
	#x.y = m.length / 2.0
	theta!(m, j[1], j[2])

	num_vehicles = 1
	new_poses = Array(Pose, num_vehicles)
	old_poses = Array(Pose, num_vehicles)
	diffs = Array(Pose, num_vehicles)

	frames = Frames(MIME("image/png"), fps=20)

	# Set up the paths
	x_path = [x.x]
	y_path = [x.y]

	# Plot the original scene
	#plot(m, b[1], x, show_mean=show_mean, show_cov=show_cov, obs=observations[1])
	plot(m, b[1], x, show_mean=show_mean, show_cov=show_cov)
	push!(frames, gcf())
	close()

	# 20 is the fps
	seconds_per_step = 2.0
	frames_per_step = round(Int, seconds_per_step * 20)

	# warn user if jammer is not where it should be
	#if abs(dx) > x.x || abs(dy) > x.y
	#	println("WARNING: jammer outside search domain.")
	#end
	colors = ["b", "r"]
	if num_vehicles == 1
		colors = ["b"]
	elseif num_vehicles == 2
		colors = ["b", "r"]
	elseif num_vehicles == 3
		colors = ["b", "r", "g"]
	elseif num_vehicles == 4
		colors = ["b", "r", "g", "k"]
	end

	# loop through all observations...
	num_obs = length(observations)
	for (oi,o) in enumerate(observations)
		println("oi = ", oi)
		xi = 1
		old_pose = (x.x, x.y, x.heading)
		old_poses[xi] = old_pose

		a = actions[oi]
		act!(m, x, actions[oi])
		new_pose = (x.x, x.y, x.heading)
		new_poses[xi] = new_pose
		# Plot everything in between old pose and new pose
		dx = (new_pose[1] - old_pose[1]) / frames_per_step
		dy = (new_pose[2] - old_pose[2]) / frames_per_step
		dh = a[3] / frames_per_step
		diffs = (dx, dy, dh)

		#plot(m, b[oi], x, show_mean=true, show_cov=true)
		#println("oi = ", oi)
		FEBOL.plot(m, b[oi], x, show_mean=show_mean, show_cov=show_cov)
		for j = 1:frames_per_step
			xi = 1
			figure()
			# determine intermediate pose
			dx, dy, dh = diffs
			old_pose = old_poses[xi]
			new_h = mod(old_pose[3] + dh, 360.0)
			old_pose = (old_pose[1] + dx, old_pose[2] + dy, new_h)
			old_poses[xi] = (old_pose)
			push!(x_path, old_pose[1])
			push!(y_path, old_pose[2])
			#plot(m, farr, old_poses; show_mean=show_mean, show_cov=show_cov, alpha=alpha)
			#FEBOL.plot(m, b[oi], old_pose, show_mean=show_mean, show_cov=show_cov, obs=o)
			FEBOL.plot(m, b[oi], old_pose, show_mean=show_mean, show_cov=show_cov)
			if show_path
				plot(x_path, y_path, colors[1])
			end

			push!(frames, gcf())
			close()
		end
		#savefig("temp_$(oi).png", format="png")
		#hold(false)
		#act!(m, x, actions[oi])
	end
	filename="temp.gif"
	write(filename, frames)
end

# really, this is for white sands results
# This shouldn't be so hard...
export simgif2
function simgif2(states::Vector{Pose}, beliefs::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false, show_path=false, j::LocTuple=m.theta)

	theta!(m, j[1], j[2])
	frames_per_second = 4
	frames = Frames(MIME("image/png"), fps=frames_per_second)

	# Set up the paths
	x_path = [x.x]
	y_path = [x.y]

	# Plot the original scene
	plot(m, beliefs[1], x, show_mean=show_mean, show_cov=show_cov,obs=0)
	push!(frames, gcf())
	close()

	# 20 is the fps
	seconds_per_step = 1.0
	frames_per_step = round(Int, seconds_per_step * frames_per_second)

	colors = ["b"]

	# loop through all beliefs...
	num_b = length(beliefs)
	s = states[1]
	for bi = 2:(num_b+1)
		println("bi = ", bi)
		sp = states[num_b]
		if bi < (num_b + 1)
			sp = states[bi]
		end

		# Plot everything in between old pose and new pose
		dx = (sp[1] - s[1]) / frames_per_step
		dy = (sp[2] - s[2]) / frames_per_step
		dh = fit_180(sp[3] - s[3]) / frames_per_step
		diffs = (dx, dy, dh)

		#plot(m, b[oi], x, show_mean=true, show_cov=true)
		#println("oi = ", oi)
		#FEBOL.plot(m, beliefs[bi-1], x,show_mean=show_mean, show_cov=show_cov)
		for j = 1:frames_per_step
			figure()
			# determine intermediate pose
			dx, dy, dh = diffs
			new_h = mod(s[3] + dh, 360.0)
			s = (s[1] + dx, s[2] + dy, new_h)

			push!(x_path, s[1])
			push!(y_path, s[2])

			FEBOL.plot(m, beliefs[bi-1], s, show_mean=show_mean, show_cov=show_cov, obs=bi-1)
			if show_path
				# changed by LD for rss submission
				#plot(x_path, y_path, colors[1])
				plot(x_path, y_path, "k")
			end

			push!(frames, gcf())
			close()
		end
	end
	#filename="temp.mp4"
	filename="temp.gif"
	write(filename, frames)
end


text:
function testplot_line1()
  plot(rand(100,3))
end

function testplot_fn1()
  plot(0:0.01:4π, [sin,cos])
end

function testplot_guides1()
  plot(rand(10); title="TITLE", xlabel="XLABEL", ylabel="YLABEL", background_color=:red)
end

function testplot_points1()
  plot(Vector[rand(10), rand(20)]; marker=:ellipse, markersize=8)
end

function testplot_points2()
  plot(Vector[rand(10), rand(20)]; marker=:ellipse, markersize=8, markercolors=[:red,:blue])
end


text:
using Test
using SparseIR

include("_util.jl")

@testset "bessels.jl" begin
    @test typestable(SparseIR.sphericalbesselj, [Int, Float64])

    @testset "small_x" begin
        n = 11
        x = 1e0
        # Mathematica
        ref = 3.099551854790080034457495747083911933213405593516888829346e-12
        @test SparseIR.sphericalbesselj(n, x) ≈ ref
    end

    @testset "large_x" begin
        n = 11
        x = 2e9
        # Mathematica
        ref = 2.020549136012222873136295240464724116017069579717738807316e-10
        @test SparseIR.sphericalbesselj(n, x) ≈ ref
    end
end


text:
module StringRepFilter

## ------------------------------------------------------------------
export has_match
export foreach_match, findall_match, findfirst_match, filter_match
export strrep
export show_strrep, show_elstrrep

## ------------------------------------------------------------------
_match(rstr::String, str::String) = match(Regex(rstr), str)
_match(reg::Regex, str::String) = match(reg, str)
_match(r, str) = _match(string(r), string(str))
# _eltype(obj) = try; eltype(obj) catch; Any end

## ------------------------------------------------------------------
"""
    strrep(obj)::String

Return an string representation of the object.
Overwrite it for custom behavior.
By default it uses `Base.string`

"""
strrep(obj) = string(obj)

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

"""
Check if the query match the object
"""
has_match

"""
    has_match(dict::Union{AbstractDict, NamedTuple}, qp::Pair)::Bool

The query pair will try to match the key and the value in the dictionary.

# Examples
```julia-repl
julia> d = Dict(:A => "AAB")
Dict{Symbol, String} with 1 entry:
  :A => "AAB"

julia> has_match(d, :A => r"^A")
true

julia> has_match(d, :B => r"^A")
false
```
"""
function has_match(dict::AbstractDict, qp::Pair)
    rkey, q = qp
    !haskey(dict, rkey) && return false
    has_match(dict[rkey], q)
end

"""
    has_match(pair::Pair, qp::Pair)::Bool

The query pair will try to match the first and secund element in the target pair.

# Examples
```julia-repl
julia> p = Pair(:A ,"AAB")
:A => "AAB"

julia> has_match(p, :A => r"^A")
true

julia> has_match(p, :B => r"^A")
false

julia> has_match(p, "B" => r"^A")
false
```
"""
function has_match(pair::Pair, qp::Pair)
    rkey, q = qp
    pkey, val = pair
    !isequal(rkey, pkey) && return false
    has_match(val, q)
end

"""
    has_match(obj, qp::Pair)::Bool

The query pair will try to match the property and the value in the object.

# Examples
```julia-repl
julia> struct Obj; A end

julia> o = Obj("AAB")
Obj("AAB")

julia> has_match(o, :A => r"^A")
true

julia> has_match(o, :B => r"^A")
false
```
"""
function has_match(obj, qp::Pair)
    rkey, q = qp
    !hasproperty(obj, rkey) && return false
    has_match(getproperty(obj, rkey), q)
end

"""
    has_match(obj, q::Symbol)::Bool

Equivalent hasproperty(obj, q)

"""
function has_match(obj, q::Symbol)
    hasproperty(obj, q)
end

"""
    has_match(obj, q::Type)::Bool

Equivalent obj isa T

"""
function has_match(obj, T::Type)
    obj isa T
end

"""
    has_match(obj, q::Type)::Bool

Equivalent obj isa T

"""
function has_match(obj, f::Function)
    f(obj)::Bool
end

"""
    has_match(obj, qp::Union{Regex, String})

Check if the query match the string representation of the object.
"""
has_match(obj, p::Union{Regex, String}) = !isnothing(_match(p, strrep(obj)))
has_match(obj, p) = !isnothing(_match(p, strrep(obj)))

"""
    has_match(obj, qps::Vector)

Check if all queries match the string representation of the object.
"""
function has_match(obj, qs::Vector)
    for qi in qs
        ismatch = has_match(obj, qi)
        !ismatch && return ismatch
    end
    return true
end

"""
    has_match(obj, qs::Tuple)::Vector

It negates the query match. Equivalent to !has_match(obj, collect(qs))
"""
has_match(obj, qs::Tuple) = !has_match(obj, collect(qs))

"""
    has_match(obj, q, qs...)::Bool

Check if any queries match the string representation of the object.
"""
function has_match(obj, q, qs...)
    ismatch = has_match(obj, q)
    ismatch && return true
    for qi in qs
        ismatch = has_match(obj, qi)
        ismatch && return true
    end
    return false
end

## ------------------------------------------------------------------
"""
    foreach_match(f::Function, col, qp, qps...)::Nothing

Apply `f(i::Int, elm)` to each element (`elm`) in `col` which match the queries (`i` is the element index).
The index is based on the `iterator` interface (e.g. index 1 means the element yield by the first iteration).
If the return value of `f` is `=== true` the iteration will stop.

# Examples
```julia-repl
julia> d = Dict("A" => "AAB", :B => "BBC", "C" => "CCB")
Dict{Any, String} with 3 entries:
  "A" => "AAB"
  "C" => "CCB"
  :B  => "BBC"

julia> foreach_match(d, r":B") do i, elm
            println(i, ": ", elm)
        end
3: Pair{Any, String}(:B, "BBC")

julia> foreach_match(d, r"") do i, elm
            println(i, ": ", elm)
            i == 1
        end
1: Pair{Any, String}("A", "AAB")
```
"""
function foreach_match(f::Function, col, qp, qps...)
    for (i, elm) in enumerate(col)
        match_flag = has_match(elm, qp)
        for qpi in qps
            match_flag |= has_match(elm, qpi)
        end
        if match_flag
            flag = f(i, elm)
            (flag === true) && return nothing
        end
    end
    return nothing
end

## ------------------------------------------------------------------
"""
Print the string representation (yield by `strrep`) of the object.

"""
show_strrep(obj) = println(strrep(obj))

"""
    show_elstrrep(col, n = 10)

Print the string representation (yield by `strrep`) of the first `n` items in the collection.
This method is intended mainly for debugging proposes.

# Examples
```julia-repl
julia> d = Dict("A" => "AAB", :B => "BBC", "C" => "CCB")
Dict{Any, String} with 3 entries:
    "A" => "AAB"
    "C" => "CCB"
    :B  => "BBC"

julia> show_elstrrep(d)
Pair{Any, String}("A", "AAB")
Pair{Any, String}("C", "CCB")
Pair{Any, String}(:B, "BBC")
```
"""
function show_elstrrep(col, n = 10)
    foreach_match(col, r"") do i, elm
        show_strrep(elm)
        if (i == n)
            println("[...]")
            return true
        end
        return false
    end
    return nothing
end

## ------------------------------------------------------------------
"""
    findall_match(col, qp, qps...)::Vector{Int}

Returns all the indexes of the elements in the collection which match the queries.
The index is based on the `iterator` interface (e.g. index 1 means the element yield by the first iteration).
For capturing the elements instead use `find_match`.

# Examples
```julia-repl
julia> d = Dict("A" => "AAB", :B => "BBA")
Dict{Any, String} with 2 entries:
  "A" => "AAB"
  :B  => "BBA"

julia> show_elstrrep(d)
Pair{Any, String}("A", "AAB")
Pair{Any, String}(:B, "BBA")

julia> findall_match(d, "A" => r"^A")
1-element Vector{Int64}:
  1

julia> findall_match(d, :A => r"^A")
Int64[]
```
"""
function findall_match(col, qp, qps...)
    founds = Int[]
    foreach_match(col, qp, qps...) do i, elm
        push!(founds, i)
        return nothing
    end
    return founds
end

"""
    findfirst_match(col, qp, qps...)::Union{Nothing, Tuple}

Returns the tuple (index, item) for the first item in the collection which match the queries.
If no match is found returns nothing.
The index is based on the `iterator` interface (e.g. index 1 means the element yield by the first iteration).

# Examples
```julia-repl
julia> d = Dict("A" => "AAB", :B => "BBA")
Dict{Any, String} with 2 entries:
  "A" => "AAB"
  :B  => "BBA"

julia> show_elstrrep(d)
Pair{Any, String}("A", "AAB")
Pair{Any, String}(:B, "BBA")

julia> findfirst_match(d, "A" => r"^A")
(1, Pair{Any, String}("A", "AAB"))

julia> findfirst_match(d, r"\".+B\"")
(1, Pair{Any, String}("A", "AAB"))

julia> findfirst_match(d, "ZZ") |> typeof
Nothing

julia> all(filter_match(d, r"[A|B]") .== filter_match(d, "A", "B"))
true
```
"""
function findfirst_match(col, qp, qps...)
    found = nothing
    foreach_match(col, qp, qps...) do i, elm
        found = (i, elm)
        return true
    end
    return found
end

"""
    filter_match(col, qp, qps...)::Vector

Returns the elements which match the queries.
If no match is found returns an empty vector.

# Examples
```julia-repl
julia> d = Dict("A" => "AAB", "B" => "BBC", "C" => "CCD")
Dict{String, String} with 3 entries:
  "B" => "BBC"
  "A" => "AAB"
  "C" => "CCD"

julia> show_elstrrep(values(d))
BBC
AAB
CCD

julia> filter_match(values(d), r"")
3-element Vector{String}:
 "BBC"
 "AAB"
 "CCD"

julia> filter_match(values(d), r"^[^B]*\$")
1-element Vector{String}:
 "CCD"

julia> filter_match(values(d), r"Z")
String[]
```
"""
function filter_match(col, qp, qps...)
    # T = _eltype(col)
    T = eltype(col)
    founds = T[]
    foreach_match(col, qp, qps...) do i, elm
        push!(founds, elm)
    end
    return founds
end

##
end # module


text:
using Symata
using Test

@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1.0
          den = 1.0
          For( j=1, j <= N, Increment(j),
            begin
              num = If(j != i, num = num * (x-xi[j]), num)
              den = If(j != i, den = den * (xi[i]-xi[j]), den)
            end
          )
          sum = sum + yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end

@sym begin
  CollocationWeightedResidualMethod(xi_, yi_, doprint_) := Module([Y, R],
    begin
      Y(x_) := LagrangePolynomial(xi, yi)
      If(doprint, Println("\nY(x) = ", Y(x), "\n"))
      #
      # Can be formulated as ytile(x) = F(x) + C1(a) * Ψ(x)
      #
      # F(x_) := 2*x^2 - x
      # C(a_) := -4a
      # Ψ(x_) := x^2 - x
      #
      #Y(x_) := F(x) + C(a)*Ψ(x)
      #
      R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
      R(x_) = Simplify(R(x) ./ (a => -C1/4))
      If(doprint, Println("R(x) = ", R(x), "\n"))
      Return(R(x))
    end
  )
end

function wrm(x1,y1, doprint=true)
    setsymata(:xin, List(x1...))
    setsymata(:yin, List(y1...))
    setsymata(:doprint, doprint ? true : false)
    #symtranseval(Meta.parse("Compile(Evaluate(CollocationWeightedResidualMethod(xin,yin,doprint)))"))
    symparseeval("Compile(Evaluate(CollocationWeightedResidualMethod(xin,yin,doprint)))")
end

println("\nExample 7.13: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("Residual for Weighted Residual Method using 1 point Lagragian Polynomial\n")

#
# A consequence of the Julia function based approach is that Y(x) and R(x) are local
# to the Symata function and are not available for examination, e.g. @sym R(x)
#

lp_eq = @sym ToString(LagrangePolynomial([0.0,0.5,1.0], [0.0,a,1.0]))
println(lp_eq)
@test lp_eq == "x*(-1.0 + 2.0*x + (-4.0)*a*(-1.0 + x))"
wrm_eq = @sym ToString(CollocationWeightedResidualMethod([0.0,0.5,1.0], [0.0,a,1.0], false))
println(wrm_eq)
@test wrm_eq == "4.0 + 2.0*C1 - 3x - 4x*(-1.0 + 2.0*x + C1*(-1.0 + x))"

#
# Another consideration is that below yt is yt(C1::Any, x::Any)
#

yt = wrm([0.0,0.5,1.0], [0.0,:a,1.0], true)

println("( Example 7.13 gives: R = -4x^2*(2 + C1) + x*(1 + 4C1) + 2*(2 + C1) )", "\n")

"
```
The function ytilde(x, C1) contains the 1-point Weighted Residual Method
polynomial with parameter C1.
```
Arguments to ytilde(x, C1):

x::Float64          : x value\n
C1::Float64         : Weighted Residual Method constant.

C1 can be estimated by e.g. collocation, subdomain, least squares or Galerkin.
See examples 7.15 through 7.18 and Fig.7.14.jl
"
ytilde(x, C1) = yt(C1, x)
println("ytilde(1.0, 0.0) = $(ytilde(1.0, 0.0))")
@test ytilde(1.0, 0.0) == -3.0

println("ytilde(0.0, 1.0) = $(ytilde(0.0, 1.0))")
@test ytilde(0.0, 1.0) == 6.0

#
# Tests with @symExpr
#

setsymata(:x1, [0.0,0.5,1.0])
wrm_1pt_col(x, C1) = @symExpr CollocationWeightedResidualMethod(x1, [0.0,a,1.0], false)

#
# The order of the arguments can be chosen. This is nice.
#

println("wrm_1pt_col(1.0, 0.0) = $(wrm_1pt_col(1.0, 0.0))")
@test wrm_1pt_col(1.0, 0.0) == ytilde(1.0, 0.0)

#
# Like @sym, I think @symExpr is not intended to return an expression from a Julia function.
#

function wrm_1(xin, xinval, C1inval)
  setsymata(:x1, xin)
  f(x, C1) = @symExpr CollocationWeightedResidualMethod(x1, [0.0,a,1.0], false)
  f(xinval, C1inval)
end

@test wrm_1([0.0,0.5,1.0], 1.0, 0.0) == ytilde(1.0, 0.0)


text:
module MCVI

using JSON
using POMDPs
using ParticleFilters
using Random
using Printf
using Distributed
using LinearAlgebra
import Statistics
import POMDPs: solve, action, rand, simulate, updater, initialize_belief, update

using POMDPLinter: @POMDP_require, @implemented

# Currently the only method of ParticleFilters.obs_weight is obs_weight(m, s, a, sp, o)
# Since MCVI uses obs_weight(m, a, sp, o) and obs_weight(m, sp, o), we define them here
# MCVI should be fixed to use the s, a, sp version.
ParticleFilters.obs_weight(m::POMDP, a, sp, o) = pdf(observation(m, a, sp), o)
ParticleFilters.obs_weight(m::POMDP, sp, o) = pdf(observation(m, sp), o)

const Reward = Float64

# implementation warning handled by POMDPrequire
function init_lower_action end

include("subspace.jl")
include("belief.jl")
include("policy.jl")
include("simulate.jl")
include("alphaedge.jl")
include("updater.jl")
include("solver.jl")

export MCVISolver,
       MCVIPolicy,

       solve,
       action,
       MCVIUpdater,
       MCVISimulator,
       updater,
       dump_json,
       upper_bound,
       lower_bound

end # module


text:
using YaoBlocks
export LabelBlock

"""
    LabelBlock{BT,N} <: TagBlock{N}

A marker to mark a circuit applying on a continous block for better plotting.
"""
struct LabelBlock{BT<:AbstractBlock,N} <: TagBlock{BT,N}
    content::BT
    name::String
end

YaoBlocks.content(cb::LabelBlock) = cb.content
function LabelBlock(x::BT, name::String) where {N,BT<:AbstractBlock{N}}
    LabelBlock{BT,N}(x, name)
end

function is_continuous_chunk(x)
    length(x) == 0 && return true
    return length(x) == maximum(x)-minimum(x)+1
end

YaoBlocks.PropertyTrait(::LabelBlock) = YaoBlocks.PreserveAll()
YaoBlocks.mat(::Type{T}, blk::LabelBlock) where {T} = mat(T, content(blk))
YaoBlocks.apply!(reg::YaoBlocks.AbstractRegister, blk::LabelBlock) = apply!(reg, content(blk))
YaoBlocks.chsubblocks(blk::LabelBlock, target::AbstractBlock) = LabelBlock(target, blk.name)

Base.adjoint(x::LabelBlock) = LabelBlock(adjoint(content(x)), endswith(x.name, "†") ? x.name[1:end-1] : x.name*"†")
Base.copy(x::LabelBlock) = LabelBlock(copy(content(x)), x.name)
YaoBlocks.Optimise.to_basictypes(block::LabelBlock) = block

export label
label(b::AbstractBlock, str::String) = LabelBlock(b, str)


text:
using Contextual
using TinyCassette

using Test: @testset, @test, @test_throws

@testset "Contextual" begin
    include("test_inline_contextual.jl")
    include("test_outline_contextual.jl")
end