text:
"""
    Clock{T<:Number}

    Clock{T}(time, iteration)

Keeps track of the current `time` and `iteration` number.
"""
mutable struct Clock{T<:Number}
       time :: T
  iteration :: Int
end


text:
"""
Tree evolution
Author: Alex Mallen (atmallen@uw.edu)

Can trees grow down?

Trees are tall because they compete for light.
    In a forest where all trees are at least 10 meters tall, with 8 meters of trunk, each tree
    could be 8 meters shorter and recieve the same amount of light, but
    spend much less energy on growth.

    However, the dominant strategy, paradoxically, is to grow tall because regardless
    of how tall the neighboring trees are, it is more advantageous to grow tall. This
    is because for every millimeter of height gained by the tree, y units of energy
    are used to grow, while z units of energy are gained due to increased light,
    where y < z.

    If trees were rational human beings, they could all be better off by
    coming to an agreement that changes their dominant strategy: A tree gives its
    neighbor (in the direction of the sun) x units of energy for every millimeter shorter it becomes, where
    z - y < x < y < z. The reason x < y is because the x delivered to one tree
    must come from the y gained by its neighbor which is now able to be shorter

    z - y < y ⟹ y > z/2 is one important constraint on the system.

    More precisely, "advantage" for a tree (or its genes) is defined as likelihood
    of persisting long-term, and we assume that energy is a good proxy for this.

    (Is this green-beard altruism?) I don't think so because it doesn't distinguish like from unlike
    The tree actually needs to shrink when it donates
    The difference between economics and evolution is that capitalists are profit maximizing while genes are *relative* profit maximizing
"""

cd(@__DIR__)
cd("..")
using Pkg
Pkg.activate(".")

using Statistics
using Plots
using Random

seed = rand(1:1000)
Random.seed!(seed)
println(string("seed", seed))

abstract type Tree end

mutable struct AltruisticTree <: Tree
    id::Int32
    height::Int32
    donation_rate::Float32  # x
    growth_cost::Float32  # y
    growth_gain::Float32  # z
    energy::Float32
    age::Int64
end

modp1(x::Integer, m::Integer) = x < 1 ? m + x % m : 1 + (x - 1) % m

function both_animate!(trees::Vector{Tree}, grid::Matrix{Int32}, trees_cache::Vector{Tree}, generations::Int64, name::String)
    ydim, xdim = size(grid)
    avg_energies = Vector{Float32}()
    avg_heights = Vector{Float32}()
    anim = @animate for gen in 1:generations
        # first update each tree's energy
        avg_energy = 0
        for i in 1:length(trees)
            neigh1, neigh2 = trees[modp1(i - 1, length(trees))], trees[modp1(i + 1, length(trees))]  # neighbor in the direction of the sun
            tree = trees[i]
            tree.energy = tree.growth_gain * (tree.height - neigh1.height + 45) - tree.growth_cost * tree.height
            if typeof(neigh2) === AltruisticTree && neigh2.height > tree.height
                tree.energy += neigh2.donation_rate * (neigh2.height - tree.height)  # TODO this 45 ÷ 2 might need to be changed to avg_height or neigh2.height
            end
            if typeof(tree) == AltruisticTree && tree.height > neigh1.height
                tree.energy -= tree.donation_rate * (tree.height - neigh1.height)  # TODO this 45 ÷ 2 might need to be changed to avg_height or neigh2.height
            end
            avg_energy += tree.energy
            if tree.energy <= 0
                println(tree.energy)
            end
        end
        avg_energy = round(avg_energy / length(trees), digits=1)
        push!(avg_energies, avg_energy)

        avg_height = 0
        for i in 1:length(trees)
            # draw trees[i] onto current grid
            # grid[2:trees[i].height + 1, 4 * (i - 1) + 2(gen % 2 + 1)] = trees[i].id
            grid[2:ydim - 4, 4 * (i - 1) + 2] .= 0
            grid[2:trees[i].height + 1, 4 * (i - 1) + 2] .= trees[i].id
            avg_height += trees[i].height

            # now update trees
            j1, j2 = modp1(i - 1, length(trees)), modp1(i + 1, length(trees))
            trees_cache[i] = rand() < trees[j1].energy / (trees[j1].energy + trees[j2].energy) ? trees[j1] : trees[j2]
            trees_cache[i].height = min(max(trees_cache[i].height + rand(-1:1), 1), 45)  # random mutation: height must be between 1 and 45
        end
        avg_height = round(avg_height / length(trees), digits=1)
        push!(avg_heights, avg_height)
        trees, trees_cache = trees_cache, trees
        # plot
        heatmap(1:xdim, 1:ydim, grid,
                colorbar=false, grid=false, ticks=false,
                showaxis=false, title=string(gen, " ", avg_height, " both; ", name), dpi=150)
    end
    anim, avg_energies, avg_heights
end

function continuous_time_evolve!(trees::Vector{AltruisticTree}, occ_idxs::Set{Int64}, free_idxs::Set{Int64}, generations::Int64, age_cap::Int64=100, spawnrate::Float64=0.01)
    avg_energies = Vector{Float32}()
    avg_heights = Vector{Float32}()
    avg_xs = Vector{Float32}()
    counts = Vector{Int64}()

        anim = @animate for gen in 1:generations
        heights = Vector{Float32}()
        energies = Vector{Float32}()
        xs = Vector{Float32}()
        avg_height = 0
        avg_energy = 0
        count = 0
        avg_x = 0

        # update energies and ages
        for idx in occ_idxs
            tree = trees[idx]
            idx_left = get_left_neigh_idx(trees, occ_idxs, idx)
            neigh_left = trees[idx_left]
            idx_right = get_right_neigh_idx(trees, occ_idxs, idx)
            neigh_right = trees[idx_right]
            tree.energy = tree.growth_gain * (tree.height - neigh_left.height / modp1(idx - idx_left, length(trees)) + 10) - tree.growth_cost * tree.height
            offset = 0  # diff between tree heights at which altruistic trees will start giving
            if neigh_right.height > tree.height + offset
                tree.energy += neigh_right.donation_rate * (neigh_right.height - tree.height + offset)  # TODO this 45 ÷ 2 might need to be changed to avg_height or neigh2.height
            end
            if tree.height > neigh_left.height + offset
                tree.energy -= tree.donation_rate * (tree.height - neigh_left.height + offset)  # TODO this 45 ÷ 2 might need to be changed to avg_height or neigh2.height
            end
            count += 1
            tree.age += 1
            push!(heights, tree.height)
            push!(xs, tree.donation_rate)
            push!(energies, tree.energy)
            avg_height += tree.height
            avg_energy += tree.energy
            avg_x += tree.donation_rate
        end

        avg_height = round(avg_height / length(occ_idxs), digits=1)
        push!(avg_heights, avg_height)
        avg_energy = round(avg_energy / length(occ_idxs), digits=1)
        push!(avg_energies, avg_energy)
        avg_x = round(avg_x / length(occ_idxs), digits=2)
        push!(avg_xs, avg_x)
        push!(counts, count)

        print("correlation between height and energy: ")
        println(round(cor(heights, energies), digits=3))
        print("between x and energy: ")
        println(round(cor(xs, energies), digits=3))
        p1 = histogram(heights, bins=0:2:46, title=gen, xaxis="height", ylims=(-0.1, length(trees) / 4), leg=false)
        p2 = histogram(xs, bins=0:0.2:5, title=gen, xaxis="donation rate (x)", ylims=(-0.1, length(trees) / 4), leg=false)
        p3 = histogram(energies, bins=-50:4:60, title=gen, xaxis="energy", ylims=(-0.1, length(trees) / 4), leg=false)
        plot(p1, p2, p3)

        # birth and death
        for idx in occ_idxs  # concurrent modification seems to be ok in this case
            tree = trees[idx]
            if tree.energy <= 0 || tree.age > age_cap  # death
                pop!(occ_idxs, idx)
                push!(free_idxs, idx)
            elseif length(free_idxs) > 0 && rand() < spawnrate * tree.energy  # birth
                new_height = min(45, max(1, tree.height + rand(-1:1)))
                new_x = max(0, tree.donation_rate + randn() / 10)
                new_tree = AltruisticTree(tree.id, new_height, new_x, tree.growth_cost, tree.growth_gain, 0, 0)
                new_idx = rand(free_idxs)
                pop!(free_idxs, new_idx)
                push!(occ_idxs, new_idx)
                trees[new_idx] = new_tree
            end
        end
    end
    anim, avg_energies, avg_heights, counts, avg_xs
end

function continuous_time_evolve_no_animate!(trees::Vector{AltruisticTree}, occ_idxs::Set{Int64}, free_idxs::Set{Int64}, generations::Int64, age_cap::Int64=100, spawnrate::Float64=0.01)
    for gen in 1:generations

        # update energies and ages
        for idx in occ_idxs
            tree = trees[idx]
            idx_left = get_left_neigh_idx(trees, occ_idxs, idx)
            neigh_left = trees[idx_left]
            idx_right = get_right_neigh_idx(trees, occ_idxs, idx)
            neigh_right = trees[idx_right]
            tree.energy = tree.growth_gain * (tree.height - neigh_left.height / modp1(idx - idx_left, length(trees)) + 10) - tree.growth_cost * tree.height
            offset = 0  # diff between tree heights at which altruistic trees will start giving #TODO make this evolve
            if neigh_right.height > tree.height + offset
                tree.energy += neigh_right.donation_rate * (neigh_right.height - tree.height + offset)
                # TODO decrease this tree so that its trunk height is its neighbor's height
            end
            if tree.height > neigh_left.height + offset
                tree.energy -= tree.donation_rate * (tree.height - neigh_left.height + offset)
                # TODO decrease height of tree --I think this needs a trunk_height parameter
            end
            tree.age += 1
        end

        # birth and death
        for idx in occ_idxs  # concurrent modification seems to be ok in this case
            tree = trees[idx]
            if tree.energy <= 0 || tree.age > age_cap  # death
                pop!(occ_idxs, idx)
                push!(free_idxs, idx)
            elseif length(free_idxs) > 0 && rand() < spawnrate * tree.energy  # birth
                new_height = min(45, max(1, tree.height + rand(-1:1)))
                new_x = max(0, tree.donation_rate + randn() / 10)
                new_tree = AltruisticTree(tree.id, new_height, new_x, tree.growth_cost, tree.growth_gain, 0, 0)
                new_idx = rand(free_idxs)
                pop!(free_idxs, new_idx)
                push!(occ_idxs, new_idx)
                trees[new_idx] = new_tree
            end
        end
    end

    heights = Vector{Float32}()
    energies = Vector{Float32}()
    xs = Vector{Float32}()
    count = 0
    for idx in occ_idxs
        tree = trees[idx]
        push!(heights, tree.height)
        push!(energies, tree.energy)
        push!(xs, tree.donation_rate)
        count += 1
    end
    heights, energies, xs, count
end

function get_left_neigh_idx(trees::Vector{AltruisticTree}, occ_idxs::Set{Int64}, idx::Int64)
    curr = modp1(idx - 1, length(trees))
    while !(curr in occ_idxs)
        curr =  modp1(curr - 1, length(trees))
    end
    curr
end

function get_right_neigh_idx(trees::Vector{AltruisticTree}, occ_idxs::Set{Int64}, idx::Int64)
    curr = modp1(idx + 1, length(trees))
    while !(curr in occ_idxs)
        curr =  modp1(curr + 1, length(trees))
    end
    curr
end

# 1 normal trees
# 2 altruistic trees
# 3 normal vs altruistic trees
x, y, z = 1, 3, 4
num_trees = 200
xdim = 4 * num_trees
ydim = 50  # ⟹ 45 is max height
grid = zeros(Int32, ydim, xdim)
grid[1, :] .= -2 * num_trees  # ground
grid[ydim - 3:ydim - 1, 2:2 + num_trees ÷ 7] .= 2 * num_trees  # sunshine!
generations = 10_000
age_cap = 100
spawnrate = 0.1
name = string("gens=", generations, " trees=", num_trees, " x,y,z=", x, ",", y, ",", z,
              " sr=", spawnrate, " ac=", age_cap, " seed=", seed)

function both_trees_init()
    trees = Vector{Tree}()  # from l to r
    for i in 1:num_trees
        height = rand(10:20)
        if rand() < ratio
            push!(trees, NormalTree(i + 10, height, y, z, 0, 0))
        else
            push!(trees, AltruisticTree(-i - 10, height, x, y, z, 0, 0))
        end
    end
    results = both_animate!(trees, grid, similar(trees), generations, name)
    results
end

function continuous_time_init()
    trees = Vector{AltruisticTree}(undef, 2 * num_trees)
    perm = randperm(length(trees))
    idxs = Set(perm[1:num_trees])
    free_idxs = Set(perm[num_trees + 1:end])
    for i in idxs
        height = rand(10:20)
        age = rand(0:age_cap - 1)
        trees[i] = AltruisticTree(-i - 10, height, max(x + randn(), 0), y, z, 0, age)
    end
    results = continuous_time_evolve!(trees, idxs, free_idxs, generations, age_cap, spawnrate)
    results
end

function init_test_x()
    n_trials = 10
    n_quantiles = 5
    std = 0.1
    d = 1 / (n_quantiles + 1)
    quantiles = d:d:1 - d
    x_inits = 0:0.25:3
    s = (length(x_inits), n_quantiles)
    q_heights = Matrix{Float32}(undef, s)
    q_energies = Matrix{Float32}(undef, s)
    q_end_xs = Matrix{Float32}(undef, s)
    q_counts = Matrix{Float32}(undef, s)

    for pair in enumerate(x_inits)
        i, x = pair
        println(x)
        heights = Vector{Float32}()
        energies = Vector{Float32}()
        end_xs = Vector{Float32}()
        counts = Vector{Float32}()
        for trial in 1:n_trials
            trees = Vector{AltruisticTree}(undef, 2 * num_trees)
            perm = randperm(length(trees))
            idxs = Set(perm[1:num_trees])
            free_idxs = Set(perm[num_trees + 1:end])
            for i in idxs
                height = rand(10:20)
                age = rand(0:age_cap - 1)
                trees[i] = AltruisticTree(-i - 10, height, max(x + randn() * std, 0), y, z, 0, age)
            end
            hs, es, xs, count = continuous_time_evolve_no_animate!(trees, idxs, free_idxs, generations, age_cap, spawnrate)
            push!(heights, mean(hs))
            push!(energies, mean(es))
            push!(end_xs, mean(xs))
            push!(counts, count)
        end
        sort!(heights)
        sort!(energies)
        sort!(end_xs)
        sort!(counts)
        for pair in enumerate(quantiles)
            j, q = pair
            q_heights[i, j] = get_quantile(heights, q)
            q_energies[i, j] = get_quantile(energies, q)
            q_end_xs[i, j] = get_quantile(end_xs, q)
            q_counts[i, j] = get_quantile(counts, q)

        end
    end
    labels = reshape([string("x--", round(q, digits=2)) for q in quantiles], 1, n_quantiles)
    p = plot(x_inits, q_end_xs, label=labels, xaxis="initial x", legend=:top)
    savefig(p, string("gifs\\", name, ".png"))
    p, q_heights, q_energies, q_end_xs, q_counts
end

get_quantile(list::Vector{Float32}, q::Float64) = list[Int(floor(length(list) * q + 1))]

# anim, avg_energies, avg_heights, counts, avg_xs = continuous_time_init()
# g = gif(anim, string("gifs\\", name, ".gif"), fps=15)
# p = plot(avg_xs, lab="x")
# plot!(p, counts ./ maximum(counts), lab="# trees")
# plot!(p, avg_energies ./ maximum(avg_energies), lab="energy/tree")
# plot!(p, avg_heights ./ maximum(avg_heights), lab="height", xaxis="time", legend=:top, ylims=(-0.1, maximum(avg_xs) + 0.1))
# savefig(p, string("gifs\\", name, ".png"))
p, q_heights, q_energies, q_end_xs, q_counts = init_test_x()
p

ext: jl

max_stars_repo_name: AlexTMallen/the-multiverse

max_issues_repo_name: AlexTMallen/the-multiverse

max_forks_repo_name: AlexTMallen/the-multiverse


text:
@doc raw"""
```
t = find_threshold(MinimumError(), histogram, edges)
```

Under the assumption that the histogram is a mixture of two Gaussian
distributions the threshold is chosen such that the expected misclassification
error rate is minimised.

# Output

Returns a real number `t` in `edges`. The `edges` parameter represents an
`AbstractRange` which specifies the intervals associated with the histogram bins.

# Details

Let ``f_i`` ``(i=1 \ldots I)`` denote the number of observations in the
``i``th bin of the histogram. Then the probability that an observation
belongs to the ``i``th bin is given by  ``p_i = \frac{f_i}{N}`` (``i = 1,
\ldots, I``), where ``N = \sum_{i=1}^{I}f_i``.

The minimum error thresholding method assumes that one can find a threshold
``T`` which partitions the data into two categories,  ``C_0`` and ``C_1``, such that
the data can be modelled by a mixture of two Gaussian distribution. Let
```math
P_0(T) = \sum_{i = 1}^T p_i \quad \text{and} \quad P_1(T) = \sum_{i = T+1}^I p_i
```
denote the cumulative probabilities,
```math
\mu_0(T) = \sum_{i = 1}^T i \frac{p_i}{P_0(T)} \quad \text{and} \quad \mu_1(T) = \sum_{i = T+1}^I i \frac{p_i}{P_1(T)}
```
denote the means, and
```math
\sigma_0^2(T) = \sum_{i = 1}^T (i-\mu_0(T))^2 \frac{p_i}{P_0(T)} \quad \text{and} \quad \sigma_1^2(T) = \sum_{i = T+1}^I (i-\mu_1(T))^2 \frac{p_i}{P_1(T)}
```
denote the variances of categories ``C_0`` and ``C_1``, respectively.

Kittler and Illingworth proposed to use the minimum error criterion function
```math
J(T) = 1 + 2 \left[ P_0(T) \ln \sigma_0(T) + P_1(T) \ln \sigma_1(T) \right] - 2 \left[P_0(T) \ln P_0(T) + P_1(T) \ln P_1(T) \right]
```
to assess the discreprancy between the mixture of Gaussians implied by a particular threshold ``T``,
and the piecewise-constant probability density function represented by the histogram.
The discrete value ``T`` which minimizes the function ``J(T)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(MinimumError(), counts[1:end], edges)
```

# References

1. J. Kittler and J. Illingworth, “Minimum error thresholding,” Pattern Recognition, vol. 19, no. 1, pp. 41–47, Jan. 1986. [doi:10.1016/0031-3203(86)90030-0](https://doi.org/10.1016/0031-3203%2886%2990030-0)
2. Q.-Z. Ye and P.-E. Danielsson, “On minimum error thresholding and its implementations,” Pattern Recognition Letters, vol. 7, no. 4, pp. 201–206, Apr. 1988. [doi:10.1016/0167-8655(88)90103-1](https://doi.org/10.1016/0167-8655%2888%2990103-1)
"""
function find_threshold(algorithm::MinimumError, histogram::AbstractArray, edges::AbstractRange)
    frequencies = cumsum(histogram)
    X = cumsum(edges .* histogram)
    X² = cumsum(edges.^2 .* histogram)
    𝔼X₁, 𝔼X₁², 𝔼X₂, 𝔼X₂² = compute_expectations(frequencies, X, X²)
    minJ = typemax(Float64)
    index = firstindex(frequencies)
    for t =  firstindex(frequencies):lastindex(frequencies)-1
        P₁ = max(frequencies[t] / frequencies[end], eps())
        P₂ = 1 - P₁
        # Var(X) = E[X²] - E[X]² = E[X²] - μ²
        σ₁ = sqrt(max(𝔼X₁²[t] - 𝔼X₁[t]^2, eps()))
        σ₂ = sqrt(max(𝔼X₂²[t] - 𝔼X₂[t]^2, eps()))
        J = 1 + 2*(P₁*log(σ₁) + P₂*log(σ₂)) - 2*(P₁*log(P₁)+ P₂*log(P₂))
        if (J < minJ && isfinite(J))
            minJ = J
            index = t
        end
    end
    edges[index]
end

#=
  If the histogram consists of integers then we risk integer overflow when we
  take the cumulative sums. Hence, we explicitlty convert the histogram to
  Float.
=#
function find_threshold(algorithm::MinimumError, histogram::AbstractArray{T}, edges::AbstractRange) where T <: Int
    find_threshold(algorithm, convert.(Float64, histogram), edges)
end

function compute_expectations(frequencies, X, X²)
    dimensions = axes(frequencies)
    𝔼X₁ = zeros(dimensions)
    𝔼X₁² = zeros(dimensions)
    𝔼X₂ = zeros(dimensions)
    𝔼X₂² = zeros(dimensions)
    for i = firstindex(X):lastindex(X)-1
        𝔼X₁[i] = X[i] / frequencies[i]
        𝔼X₁²[i] = X²[i] / frequencies[i]
        𝔼X₂[i] = (X[end] - X[i]) /  (frequencies[end] - frequencies[i])
        𝔼X₂²[i] = (X²[end] - X²[i]) /  (frequencies[end] - frequencies[i])
    end
    𝔼X₁, 𝔼X₁², 𝔼X₂, 𝔼X₂²
end


text:
"""
A `ClosedInterval(left, right)` is an interval set that includes both its upper and lower bounds. In
mathematical notation, the constructed range is `[left, right]`.
"""
immutable ClosedInterval{T} <: AbstractInterval{T}
    left::T
    right::T

    ClosedInterval(l::T, r::T) = new(l, r)
end

ClosedInterval{T}(left::T, right::T) = ClosedInterval{T}(left, right)
(::Type{ClosedInterval{T}}){T}(left, right) =
    ClosedInterval{T}(checked_conversion(T, left, right)...)

function ClosedInterval(left, right)
    # Defining this as ClosedInterval(promote(left, right)...) has one problem:
    # if left and right do not promote to a common type, it triggers a StackOverflow.
    T = promote_type(typeof(left), typeof(right))
    ClosedInterval{T}(checked_conversion(T, left, right)...)
end

ClosedInterval(i::AbstractInterval) = convert(ClosedInterval{eltype(i)}, i)
(::Type{ClosedInterval{T}}){T}(i::AbstractInterval) = convert(ClosedInterval{T}, i)

..(x, y) = ClosedInterval(x, y)

±(x, y) = ClosedInterval(x - y, x + y)
±(x::CartesianIndex, y) = map(ClosedInterval, (x - y).I, (x + y).I)

show(io::IO, I::ClosedInterval) = print(io, I.left, "..", I.right)

in(v, I::ClosedInterval) = I.left <= v <= I.right
in(a::ClosedInterval, b::ClosedInterval) = (b.left <= a.left) & (a.right <= b.right)

isempty(A::ClosedInterval) = A.left > A.right

isequal(A::ClosedInterval, B::ClosedInterval) = (isequal(A.left, B.left) & isequal(A.right, B.right)) | (isempty(A) & isempty(B))

==(A::ClosedInterval, B::ClosedInterval) = (A.left == B.left && A.right == B.right) || (isempty(A) && isempty(B))

const _closed_interval_hash = UInt == UInt64 ? 0x1588c274e0a33ad4 : 0x1e3f7252

hash(I::ClosedInterval, h::UInt) = hash(I.left, hash(I.right, hash(_closed_interval_hash, h)))

minimum(I::ClosedInterval) = I.left
maximum(I::ClosedInterval) = I.right

function intersect(A::ClosedInterval, B::ClosedInterval)
    left = max(A.left, B.left)
    right = min(A.right, B.right)
    ClosedInterval(left, right)
end

function union{T<:AbstractFloat}(A::ClosedInterval{T}, B::ClosedInterval{T})
    max(A.left, B.left) <= nextfloat(min(A.right, B.right)) || throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(A, B)
end

function union(A::ClosedInterval, B::ClosedInterval)
    max(A.left, B.left) <= min(A.right, B.right) || throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(A, B)
end

function _union(A::ClosedInterval, B::ClosedInterval)
    left = min(A.left, B.left)
    right = max(A.right, B.right)
    ClosedInterval(left, right)
end

issubset(A::ClosedInterval, B::ClosedInterval) = ((A.left in B) && (A.right in B)) || isempty(A)

⊇(A::ClosedInterval, B::ClosedInterval) = issubset(B, A)


text:
# Same as example/from_csv.jl, but without the time-consuming plot.

using SellmeierFit
using SimpleConstants  # for micro = 1e-6

## Main section
# Load a set of measured refractive index data from Malitson.csv
dir = @__DIR__
path = joinpath(dir, "data", "Malitson.csv")
(; λ, ε) = SellmeierFit.read(path)

# Fit the loaded data to the Sellmeier equation.
(; mdl, err) = fit_sellmeier(λ, ε)

## Analysis section
# The calculated fit parameters are the same as those shown in RefractiveIndex.info's
# Dispersion formula section.
println("B = $(mdl.str)")
println("√C (µm) = $(mdl.λres ./ micro)")
println("Error between data and fit = $err")


text:
using HDF5
include("./karger.jl")
using .Karger

# USAGE: julia src/julia/calculate_watershed.jl filename
# where filename is just the name (without path or .h5).
# The file has to be an HDF5 file with fields n, edges, weights and seeds. The output
# will be written in a file with the same name in the "watershed" directory

graph_path = "results/graphs/" * ARGS[1] * ".h5"
result_path = "results/watershed/" * ARGS[1] * ".h5"

n = h5read(graph_path, "n")
edges = h5read(graph_path, "edges")
weights = h5read(graph_path, "weights")
seeds = h5read(graph_path, "seeds")
g = Graph(n, transpose(edges) .+ 1, weights)
mkpath(dirname(result_path))
cp(graph_path, result_path, force=true)
h5write(result_path, "potential", watershed(g, seeds))


text:
function add_gigasecond(date::DateTime)

end


text:
module Crypto

using Pkg

# TODO
#  Warning: `Pkg.dir(pkgname, paths...)` is deprecated; instead, do `import ColonyRuntime; joinpath(dirname(pathof(ColonyRuntime)), "..", paths.
#  But, pathof just say LoadError: UndefVarError: ColonyRuntime not defined
crypto_lib_path = Pkg.dir("ColonyRuntime") * "lib/cryptolib.so"

function prvkey()
  prvkey_cstr = ccall((:prvkey, crypto_lib_path), Cstring, ())
  return unsafe_string(prvkey_cstr)
end

function id(prvkey::AbstractString)
  id_cstr = ccall((:id, crypto_lib_path), Cstring, (Cstring,), prvkey)
  return unsafe_string(id_cstr)
end

function sign(msg::AbstractString, prvkey::AbstractString)
  sig_cstr = ccall((:sign, crypto_lib_path), Cstring, (Cstring,Cstring), msg, prvkey)
  return unsafe_string(sig_cstr)
end

function hash(msg::AbstractString)
  msg_cstr = ccall((:hash, crypto_lib_path), Cstring, (Cstring,), msg)
  return unsafe_string(msg_cstr)
end

function recoverid(msg::AbstractString, sig::AbstractString)
  id_cstr = ccall((:recoverid, crypto_lib_path), Cstring, (Cstring,Cstring), msg, sig)
  return unsafe_string(id_cstr)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 98d66a30-c818-4e3d-ab62-69ebe6628815
using CairoMakie

# ╔═╡ 56ac4db9-96b4-4dfa-bb2b-0c1aa3cd551d
using Statistics

# ╔═╡ dfc116e7-606b-41cc-b6e9-c3f9ea336496
import PWS

# ╔═╡ 2737ac8b-c736-4a8f-a965-215b56788c05
system = PWS.gene_expression_system(dtimes=0:0.1:5)

# ╔═╡ 06ecfe25-740c-4267-9a2e-8ab512ba0562
md"""
# Direct Monte Carlo Estimate
"""

# ╔═╡ 86b9980f-3da8-4c49-8791-4045fff3b283
num_responses = 10^3

# ╔═╡ e6659282-7226-4b8b-a747-6478c1e2eac0
directmc = PWS.DirectMCEstimate(100)

# ╔═╡ 02932bf1-28dd-45e1-89d3-6e5b63d8f34a
directmc_result = PWS.mutual_information(system, directmc; num_responses)

# ╔═╡ 77d74d1d-59f6-4dc4-a79e-8a49954de27d
md"""
# Particle Filter Estimate
"""

# ╔═╡ 9bbc49c1-6b3a-4a3a-9f7b-b149bb00b44f
smc = PWS.SMCEstimate(100)

# ╔═╡ 2de7c4c7-906c-4035-a0b5-d459e11229a4
smc_result = PWS.mutual_information(system, smc; num_responses)

# ╔═╡ 3a8f3d25-01b2-40fb-a367-0245a0e120d4
md"# Thermodynamic Integration Estimate"

# ╔═╡ 52f1d352-7087-459d-b5c2-ec2fffb3f4fb
ti = PWS.TIEstimate(0, 20, 20)

# ╔═╡ f4262c1a-4df1-475c-8c3e-0a6a6f2c01b6
ti_result = PWS.mutual_information(system, ti; num_responses)

# ╔═╡ 5246f4b7-54c5-447b-b600-226357019c1e
md"# Comparison"

# ╔═╡ 5759fff6-d196-4e25-af7f-1d05b88ea87f
function plot_comparison(results...)
	fig = Figure()
	ax1 = Axis(fig[1, 1])
	ax2 = Axis(fig[2, 1])
	for r in results
		mi = mean(r.MutualInformation)
		error = 3*(std(r.MutualInformation) / sqrt(num_responses))
		# band!(system.dtimes, mi - error/2, mi + error/2)
		lines!(ax1, system.dtimes, mi)
		lines!(ax2, system.dtimes, mi ./ system.dtimes)
	end
	fig
end

# ╔═╡ a403bac4-2dec-431a-806f-5444a577d931
plot_comparison(directmc_result, smc_result)

# ╔═╡ fb71e194-8f76-46a9-8717-02de81272b3f
md"# Helper Functions"

# ╔═╡ 6ad0591a-15b5-436b-a812-1834ae68f503
function plot_result(result)
	fig = Figure()
	Axis(fig[1, 1])
	for mi in result.MutualInformation
		lines!(system.dtimes, mi, color=:grey, linewidth=0.3)
	end
	lines!(system.dtimes, mean(result.MutualInformation), linewidth=5)
	fig
end

# ╔═╡ 09335008-07b8-4905-a028-218d5d5020bb
plot_result(directmc_result)

# ╔═╡ d41eb941-f1b2-40bc-a0c6-0a3f82651796
plot_result(smc_result)

# ╔═╡ ad6f634d-a4f6-4b34-97e1-bad2c688ef98
plot_result(ti_result)

# ╔═╡ Cell order:
# ╠═dfc116e7-606b-41cc-b6e9-c3f9ea336496
# ╠═2737ac8b-c736-4a8f-a965-215b56788c05
# ╟─06ecfe25-740c-4267-9a2e-8ab512ba0562
# ╠═86b9980f-3da8-4c49-8791-4045fff3b283
# ╠═e6659282-7226-4b8b-a747-6478c1e2eac0
# ╠═02932bf1-28dd-45e1-89d3-6e5b63d8f34a
# ╠═09335008-07b8-4905-a028-218d5d5020bb
# ╟─77d74d1d-59f6-4dc4-a79e-8a49954de27d
# ╠═9bbc49c1-6b3a-4a3a-9f7b-b149bb00b44f
# ╠═2de7c4c7-906c-4035-a0b5-d459e11229a4
# ╠═d41eb941-f1b2-40bc-a0c6-0a3f82651796
# ╟─3a8f3d25-01b2-40fb-a367-0245a0e120d4
# ╠═52f1d352-7087-459d-b5c2-ec2fffb3f4fb
# ╠═f4262c1a-4df1-475c-8c3e-0a6a6f2c01b6
# ╠═ad6f634d-a4f6-4b34-97e1-bad2c688ef98
# ╟─5246f4b7-54c5-447b-b600-226357019c1e
# ╠═5759fff6-d196-4e25-af7f-1d05b88ea87f
# ╠═a403bac4-2dec-431a-806f-5444a577d931
# ╟─fb71e194-8f76-46a9-8717-02de81272b3f
# ╠═98d66a30-c818-4e3d-ab62-69ebe6628815
# ╠═56ac4db9-96b4-4dfa-bb2b-0c1aa3cd551d
# ╠═6ad0591a-15b5-436b-a812-1834ae68f503


text:
using Photometry
using Documenter

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

makedocs(;
    modules = [Photometry],
    authors = "Miles Lucas <mdlucas@hawaii.edu>",
    repo = "https://github.com/JuliaAstro/Photometry.jl/blob/{commit}{path}#L{line}",
    sitename = "Photometry.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliaastro.github.io/Photometry.jl",
        assets = String[],),
    pages = [
        "Home" => "index.md",
        "Aperture Photometry" => [
            "Getting Started" => "apertures/index.md",
            "Apertures" => "apertures/apertures.md",
            "Examples" => "apertures/examples.md"
        ]
    ],)

deploydocs(;
    repo = "github.com/JuliaAstro/Photometry.jl",)


text:
module CUSPARSE

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcusparse, unsafe_free!, @retry_reclaim

using CEnum

using Memoize

using LinearAlgebra
using LinearAlgebra: HermOrSym

using Adapt

using DataStructures

using SparseArrays

const SparseChar = Char

# core library
include("libcusparse_common.jl")
include("error.jl")
include("libcusparse.jl")
include("libcusparse_deprecated.jl")

include("array.jl")
include("util.jl")
include("types.jl")

# low-level wrappers
include("helpers.jl")
include("management.jl")
include("level1.jl")
include("level2.jl")
include("level3.jl")
include("preconditioners.jl")
include("conversions.jl")
include("generic.jl")

# high-level integrations
include("interfaces.jl")

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,cusparseHandle_t}}()

# cache for created, but unused handles
const old_handles = DefaultDict{CuContext,Vector{cusparseHandle_t}}(()->cusparseHandle_t[])

function handle()
    CUDA.detect_state_changes()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUSPARSE, ctx)) do
            handle = if isempty(old_handles[ctx])
                cusparseCreate()
            else
                pop!(old_handles[ctx])
            end

            finalizer(current_task()) do task
                push!(old_handles[ctx], handle)
            end
            # TODO: cusparseDestroy to preserve memory, or at exit?

            cusparseSetStream(handle, stream())

            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

@inline function set_stream(stream::CuStream)
    ctx = context()
    tls = task_local_storage()
    handle = get(tls, (:CUSPARSE, ctx), nothing)
    if handle !== nothing
        cusparseSetStream(handle, stream)
    end
    return
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end
end

end


text:
module PkgA
using Distributions

randchisq1() = randn() ^ 2

end # module


text:
"""
    convert2image(array) -> Array{RGB}

Convert the given SVHN tensor in WHCN format (or feature vector/matrix) to a
`RGB` array in HWN format.

```julia
julia> SVHN2.convert2image(SVHN2.traindata()[1]) # full training dataset
32×32×50000 Array{RGB{N0f8},3}:
[...]

julia> SVHN2.convert2image(SVHN2.traindata(1)[1]) # first training image
32×32 Array{RGB{N0f8},2}:
[...]
```
"""
function convert2image(array::AbstractArray{T}) where {T<:Number}
    nlast = size(array)[end]
    array = reshape(array, 32, 32, 3, :)
    array = permutedims(array, (3, 2, 1, 4))
    if size(array)[end] == 1 && nlast != 1
        array = dropdims(array, dims=4)
    end
    img = _colorview(RGB, _norm_array(array))

    img
end
_norm_array(array::AbstractArray) = array
_norm_array(array::AbstractArray{<:Integer}) = reinterpret(N0f8, convert(Array{UInt8}, array))


text:
f = open("test_file.txt")
while !eof(f)
    str = readline(f)
    display(str)
end
close(f)

for line in eachline("test_file.txt")
    display(line)
end

str = read("test_file.txt", String)

lines = readlines("test_file.txt")

println("Line #2 says: ", lines[2])
println()
println("Here are all the lines which have between 1 and 18 characters:\n")
for line in lines
    if 1 ≤ length(line) ≤ 18
        println(line)
    end
end

f = open("created_data.txt", "w")
for i = 1:5
     # Create random strings of letters
    str = String(rand('a':'z', 50))
    write(f, str * "\n")  # Write string to stream f
end
println(f) # println can be used with streams too

# Print Fibonacci numbers to file
x = y = 1
print(f, "$x $y")
for i = 1:50
    z = x + y
    x = y
    y = z
    print(f, " $z")
end
println(f)
close(f)

# Read file and print each line
for line in eachline("created_data.txt")
    println(line)
end

using DelimitedFiles

# Write file
A = rand(-100:100, 8,3)    # Sample data
writedlm("created_data.txt", A, ',')

# Print file
for line in eachline("created_data.txt")
    println(line)
end

# Read into array
B = readdlm("created_data.txt", ',')

isequal(A,B) # Check identical

function word_value(word)
    return sum(collect(word) .- 'A' .+ 1)
end
word_value("SKY")

trinums = [n*(n+1)÷2 for n = 1:50]
words = readdlm("p042_words.txt", ',', String)
nbrtriwords = count([word_value(word) ∈ trinums for word in words])
println("There are $nbrtriwords triangle words in the list")

function interpret_3digit_file(filename)
    array = zeros(Int64, 3, countlines(filename))
    for (iteration,line) in enumerate(eachline(filename))
        # Capture all 3 digits in this line
        pattern = Regex("(-?\\d+),(-?\\d+),(-?\\d+)") # easily written as r"(-?\d+),(-?\d+),(-?\d+)"
        m = match(pattern, line)
        # Convert all 3 captured strings into integers
        array[:,iteration] = parse.(Int, m.captures)
    end
    return array
end

interpret_3digit_file("created_data.txt")

function interpret_words_file(filename)
    sentences = Array[]
    for line in eachline(filename)
        # Capture all words in this line
        pattern = Regex("(\\w+)") # easily written as r"(\w+)"
        ms = eachmatch(pattern, line)
        # Collect all words into an array
        sentence = [m.match for m in ms]
        # Push current array into total array
        if !isempty(sentence)
            push!(sentences, sentence)
        end
    end
    return sentences
end

interpret_words_file("test_file.txt")


text:
module Rectangles

using ..Ahorn, Maple

displayName = "Rectangles"
group = "Brushes"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing

material = '0'
filled = true

selection = nothing
rectangleBrush = nothing

function drawRectangle(rect::Ahorn.Rectangle, layer::Ahorn.Layer, filled=filled)
    pixels = fill(true, rect.h, rect.w)
    if !filled
        pixels[2:end - 1, 2:end - 1] .= false
    end

    global rectangleBrush = Ahorn.Brush("Rectangle", pixels)

    Ahorn.drawBrush(rectangleBrush, layer, rect.x, rect.y)
end

function drawRectangles(layer::Ahorn.Layer, room::Ahorn.DrawableRoom, camera::Ahorn.Camera)
    if selection !== nothing
        drawRectangle(selection, toolsLayer)
    end
end

function applyRectangle!(rect::Ahorn.Rectangle, material::Char, filled)
    Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Rectangle $(material)", Ahorn.loadedState.room))

    Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    Ahorn.applyBrush!(rectangleBrush, tiles, material, rect.x, rect.y)
end

function cleanup()
    global selection = nothing

    Ahorn.redrawLayer!(toolsLayer)
end

function setMaterials!(layer::Ahorn.Layer)
    tileNames = Ahorn.tileNames(layer)
    displayNames = Ahorn.getBrushMaterialsNames(layer, get(Ahorn.config, "sort_brush_tiles_by_name", true))

    Ahorn.setMaterialList!(displayNames, row -> row[1] == get(tileNames, material, nothing))
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    layerName = Ahorn.layerName(targetLayer)
    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(layerName)", tileNames["Air"])[1]

    wantedMode = get(Ahorn.persistence, "brushes_rectangle_mode", "Filled")
    Ahorn.updateTreeView!(subTools, ["Filled", "Hollow"], row -> row[1] == wantedMode)

    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")
    Ahorn.updateLayerList!(["fgTiles", "bgTiles"], row -> row[1] == wantedLayer)

    Ahorn.redrawingFuncs["tools"] = drawRectangles
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)
    Ahorn.persistence["brushes_layer"] = selected

    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(selected)", tileNames["Air"])[1]
    setMaterials!(targetLayer)
end

function subToolSelected(list::Ahorn.ListContainer, selected::String)
    Ahorn.persistence["brushes_rectangle_mode"] = selected
    global filled = selected == "Filled"
end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    tileNames = Ahorn.tileNames(targetLayer)
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = tileNames[selected]
    global material = tileNames[selected]
end

function materialFiltered(list::Ahorn.ListContainer)
    tileNames = Ahorn.tileNames(targetLayer)
    Ahorn.selectRow!(list, row -> row[1] == get(tileNames, material, '0'))
end

function getFavorites()
    if targetLayer !== nothing
        key = "favorites_brushes_$(targetLayer.name)"

        return Ahorn.getFavorites(Ahorn.persistence, key)
    end

    return []
end

function materialDoubleClicked(material::String)
    if targetLayer !== nothing
        key = "favorites_brushes_$(targetLayer.name)"
        tileNames = Ahorn.tileNames(targetLayer)

        Ahorn.toggleFavorite(Ahorn.persistence, key, material)
        setMaterials!(targetLayer)
        Ahorn.updateMaterialFilter!(targetLayer.name)
    end
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "fgTiles")
end

function leftClick(x::Number, y::Number)
    if selection !== nothing
        applyRectangle!(selection, material, filled)

        global selection = nothing

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

function middleClick(x::Number, y::Number)
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    tileNames = Ahorn.tileNames(targetLayer)
    target = get(tiles.data, (y, x), '0')

    global material = target
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = material
    Ahorn.selectMaterialList!(tileNames[target])
end

function mouseMotion(x::Number, y::Number)
    # Only make a preview square if we aren't dragging
    if !Ahorn.mouseButtonHeld(0x1)
        newRect = Ahorn.Rectangle(x, y, 1, 1)

        if newRect != selection
            global selection = newRect

            Ahorn.redrawLayer!(toolsLayer)
        end
    end
end

function selectionMotion(rect::Ahorn.Rectangle)
    if rect != selection
        global selection = rect

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function selectionFinish(rect::Ahorn.Rectangle)
    if selection !== nothing
        applyRectangle!(selection, material, filled)

        global selection = nothing

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test, Statistics

### temperature patch tests, sl tet4, dl tet4, sl tet10, dl tet 10

tet4_meshfile = joinpath("test_problems_mortar_3d", "tet4.inp")
tet10_meshfile = joinpath("test_problems_mortar_3d", "tet10.inp")

# patch test temperature + abaqus inp + tet4
mesh = abaqus_read_mesh(tet4_meshfile)

upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)

lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)

interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

JuliaFEM.diagnose_interface(interface, 0.0)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("sl_lin_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
@info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)

#=
initialize!(solver)
assemble!(solver)
M, K, Kg, f, fg = get_field_assembly(solver)
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
K = K + Kg + Kb
f = f + fg + fb
K = 1/2*(K + K')
M = 1/2*(M + M')
=#

# patch test temperature + abaqus inp + tet4 + dual basis + adjust
mesh = abaqus_read_mesh(tet4_meshfile)

upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)

lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)

interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_master_elements, interface_slave_elements)
interface.properties.dual_basis = true
#interface.properties.adjust = true

JuliaFEM.diagnose_interface(interface, 0.0)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("dl_lin_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
@debug("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)

# patch test temperature + abaqus inp + tet10, quadratic surface elements
mesh = abaqus_read_mesh(tet10_meshfile)

upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)

lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)

interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_master_elements, interface_slave_elements)

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.alpha = 0.0
#   JuliaFEM.diagnose_interface(interface, 0.0)

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("sl_quad_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
@info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-6
@test isapprox(stdT, 0.0; atol=1.0e-6)

# patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2
mesh = abaqus_read_mesh(tet10_meshfile)

upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)

lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)

interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.dual_basis = true
interface.properties.alpha = 0.2

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("dl_quad_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
@info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-10
@test isapprox(stdT, 0.0; atol=1.0e-10)

### displacement patch tests, sl tet4, dl tet4, sl tet10, dl tet 10

# patch test displacement + abaqus inp + tet4 + adjust

mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
    mesh.nodes[nid][3] += 0.2
end

upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)

lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)

bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)

bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)

interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("sl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)

# patch test displacement + abaqus inp + tet4 + adjust + dual basis

mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
    mesh.nodes[nid][3] += 0.2
end

upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)

lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)

bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)

bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)

interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("dl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)

# patch test displacement + abaqus inp + tet10 + adjust

mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
    mesh.nodes[nid][3] += 0.2
end

upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)

lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)

bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)

bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)

interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("sl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)

# patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2

mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
    mesh.nodes[nid][3] += 0.2
end

upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)

lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)

bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)

bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)

bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)

bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)

interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)

removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)

interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
interface.properties.alpha = 0.2

analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("dl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)

displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: savingsplans
using AWS.Compat
using AWS.UUIDs

"""
    create_savings_plan(commitment, savings_plan_offering_id)
    create_savings_plan(commitment, savings_plan_offering_id, params::Dict{String,<:Any})

Creates a Savings Plan.

# Arguments
- `commitment`: The hourly commitment, in USD. This is a value between 0.001 and 1 million.
  You cannot specify more than three digits after the decimal point.
- `savings_plan_offering_id`: The ID of the offering.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Unique, case-sensitive identifier that you provide to ensure the
  idempotency of the request.
- `"purchaseTime"`: The time at which to purchase the Savings Plan, in UTC format
  (YYYY-MM-DDTHH:MM:SSZ).
- `"tags"`: One or more tags.
- `"upfrontPaymentAmount"`: The up-front payment amount. This is a whole number between 50
  and 99 percent of the total value of the Savings Plan. This parameter is supported only if
  the payment option is Partial Upfront.
"""
create_savings_plan(commitment, savingsPlanOfferingId; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/CreateSavingsPlan", Dict{String, Any}("commitment"=>commitment, "savingsPlanOfferingId"=>savingsPlanOfferingId, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_savings_plan(commitment, savingsPlanOfferingId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/CreateSavingsPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("commitment"=>commitment, "savingsPlanOfferingId"=>savingsPlanOfferingId, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_queued_savings_plan(savings_plan_id)
    delete_queued_savings_plan(savings_plan_id, params::Dict{String,<:Any})

Deletes the queued purchase for the specified Savings Plan.

# Arguments
- `savings_plan_id`: The ID of the Savings Plan.

"""
delete_queued_savings_plan(savingsPlanId; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DeleteQueuedSavingsPlan", Dict{String, Any}("savingsPlanId"=>savingsPlanId); aws_config=aws_config)
delete_queued_savings_plan(savingsPlanId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DeleteQueuedSavingsPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("savingsPlanId"=>savingsPlanId), params)); aws_config=aws_config)

"""
    describe_savings_plan_rates(savings_plan_id)
    describe_savings_plan_rates(savings_plan_id, params::Dict{String,<:Any})

Describes the specified Savings Plans rates.

# Arguments
- `savings_plan_id`: The ID of the Savings Plan.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filters"`: The filters.
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  additional results, make another call with the returned token value.
- `"nextToken"`: The token for the next page of results.
"""
describe_savings_plan_rates(savingsPlanId; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlanRates", Dict{String, Any}("savingsPlanId"=>savingsPlanId); aws_config=aws_config)
describe_savings_plan_rates(savingsPlanId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlanRates", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("savingsPlanId"=>savingsPlanId), params)); aws_config=aws_config)

"""
    describe_savings_plans()
    describe_savings_plans(params::Dict{String,<:Any})

Describes the specified Savings Plans.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filters"`: The filters.
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  additional results, make another call with the returned token value.
- `"nextToken"`: The token for the next page of results.
- `"savingsPlanArns"`: The Amazon Resource Names (ARN) of the Savings Plans.
- `"savingsPlanIds"`: The IDs of the Savings Plans.
- `"states"`: The states.
"""
describe_savings_plans(; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlans"; aws_config=aws_config)
describe_savings_plans(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlans", params; aws_config=aws_config)

"""
    describe_savings_plans_offering_rates()
    describe_savings_plans_offering_rates(params::Dict{String,<:Any})

Describes the specified Savings Plans offering rates.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filters"`: The filters.
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  additional results, make another call with the returned token value.
- `"nextToken"`: The token for the next page of results.
- `"operations"`: The specific AWS operation for the line item in the billing report.
- `"products"`: The AWS products.
- `"savingsPlanOfferingIds"`: The IDs of the offerings.
- `"savingsPlanPaymentOptions"`: The payment options.
- `"savingsPlanTypes"`: The plan types.
- `"serviceCodes"`: The services.
- `"usageTypes"`: The usage details of the line item in the billing report.
"""
describe_savings_plans_offering_rates(; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlansOfferingRates"; aws_config=aws_config)
describe_savings_plans_offering_rates(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlansOfferingRates", params; aws_config=aws_config)

"""
    describe_savings_plans_offerings()
    describe_savings_plans_offerings(params::Dict{String,<:Any})

Describes the specified Savings Plans offerings.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"currencies"`: The currencies.
- `"descriptions"`: The descriptions.
- `"durations"`: The durations, in seconds.
- `"filters"`: The filters.
- `"maxResults"`: The maximum number of results to return with a single call. To retrieve
  additional results, make another call with the returned token value.
- `"nextToken"`: The token for the next page of results.
- `"offeringIds"`: The IDs of the offerings.
- `"operations"`: The specific AWS operation for the line item in the billing report.
- `"paymentOptions"`: The payment options.
- `"planTypes"`: The plan type.
- `"productType"`: The product type.
- `"serviceCodes"`: The services.
- `"usageTypes"`: The usage details of the line item in the billing report.
"""
describe_savings_plans_offerings(; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlansOfferings"; aws_config=aws_config)
describe_savings_plans_offerings(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/DescribeSavingsPlansOfferings", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/ListTagsForResource", Dict{String, Any}("resourceArn"=>resourceArn); aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds the specified tags to the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tags`: One or more tags. For example, { \"tags\": {\"key1\":\"value1\",
  \"key2\":\"value2\"} }.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/TagResource", Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the specified tags from the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tag_keys`: The tag keys.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/UntagResource", Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = savingsplans("POST", "/UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys), params)); aws_config=aws_config)


text:
function read_map(map::String)
    map_dir = joinpath(dirname(@__DIR__), "maps")
    datac = RGB24.(PNGFiles.load(joinpath(map_dir,map)*".png"))
    datah = Float32.(reinterpret(UInt8,PNGFiles.load(joinpath(map_dir,MAP_LIST[map])*".png")))
    if size(datac,1) == 2*size(datah,1) # some height maps need to be upsampled
        datah = upsample_map_2x(datah)
    end
    return datac, datah
end

function upsample_map_2x(data::Matrix)
    # simple linear interpolation
    w,h = size(data)
    data_2x = similar(data,2w,2h)
    @inbounds for j = 1:h, i = 1:w
        ii,jj = 2*(i-1)+1, 2*(j-1)+1
        data_2x[ii,  jj]   = data[i,j]
        data_2x[ii+1,jj]   = (data[i,j] + data[mod1(i+1,w),j])/2
        data_2x[ii,  jj+1] = (data[i,j] + data[i,mod1(j+1,h)])/2
        data_2x[ii+1,jj+1] = (data[i,j] + data[mod1(i+1,w),mod1(j+1,h)])/2
    end
    return data_2x
end

const MAP_LIST = Dict{String,String}(
    "C1W" => "D1",
    "C2W" => "D2",
    "C3" => "D3",
    "C4" => "D4",
    "C5W" => "D5",
    "C6W" => "D6",
    "C7W" => "D7",
    "C8" => "D6",
    "C9W" => "D9",
    "C10W" => "D10",
    "C11W" => "D11",
    "C12W" => "D11",
    "C13" => "D13",
    "C14" => "D14",
    "C14W" => "D14",
    "C15" => "D15",
    "C16W" => "D16",
    "C17W" => "D17",
    "C18W" => "D18",
    "C19W" => "D19",
    "C20W" => "D20",
    "C21" => "D21",
    "C22W" => "D22",
    "C23W" => "D21",
    "C24W" => "D24",
    "C25W" => "D25",
    "C26W" => "D18",
    "C27W" => "D15",
    "C28W" => "D25",
    "C29W" => "D16",
    )

max_stars_repo_name: jebej/VoxelSpace

max_issues_repo_name: jebej/VoxelSpace

max_forks_repo_name: jebej/VoxelSpace


text:
# the number of masked elements in each set of mosaic
function nmasked_perset_mask!(nmasked::AbstractVector{Int}, mosaic::SetMosaic, elmask::AbstractVector{Bool},
                              set2index::Union{Nothing,Dict{Int,Int}} = nothing)
    @inbounds for elix in eachindex(elmask)
        elmask[elix] || continue
        for setix in view(mosaic.setXelm, :, elix)
            if set2index === nothing
                nmasked[setix] += 1
            else
                ix = get(set2index, setix, 0)
                if ix > 0
                    nmasked[ix] += 1
                end
            end
        end
    end
    return nmasked
end

# the number of masked elements in each set of mosaic
function nmasked_perset(mosaic::SetMosaic, elmasks::AbstractVector,
                        set2index::Union{Nothing,Dict{Int,Int}} = nothing)
    res = zeros(Int, set2index === nothing ? nsets(mosaic) : length(set2index), length(elmasks))
    @inbounds for (expix, elmask) in enumerate(elmasks)
        @assert length(elmask) == nelements(mosaic)
        nmasked_perset_mask!(view(res, :, expix), mosaic, elmask, set2index)
    end
    return res
end

# the number of masked elements in each set of mosaic
function nmasked_perset(mosaic::SetMosaic, elmasks::AbstractMatrix{Bool},
                        set2index::Union{Nothing,Dict{Int,Int}} = nothing)
    @assert size(elmasks, 1) == nelements(mosaic)
    res = zeros(Int, set2index === nothing ? nsets(mosaic) : length(set2index), size(elmasks, 2))
    @inbounds for expix in 1:size(elmasks, 2)
        nmasked_perset_mask!(view(res, :, expix), mosaic, view(elmasks, :, expix), set2index)
    end
    return res
end

# info for the set × mask overlap for the MaskedSetMosaic
struct MaskOverlap
    nmasked::Int    # number of masked elements
    nunmasked::Int  # number of unmasked elements
    used::Bool      # true if used for the cover problem
end

"""
[`SetMosaic`](@ref) with the elements masks (selections) on top.
Sets that are not overlapping with the masks are excluded(skipped) from `MaskedSetMosaic`.
Optionally, the filtering can include testing for the minimal overlap significance P-value.

The tiles of non-overlapped sets are removed, the tiles that have identical membership
for all the masked sets are squashed into a single tile.

## Type parameters
* `T`: the type of elements
* `S`: the type of set ids
* `E`: the type of experiment ids
"""
mutable struct MaskedSetMosaic{T, S, E} <: AbstractWeightedSetMosaic{T, S, E, Float64}
    original::SetMosaic{T, S}       # original mosaic
    elmasks::BitMatrix              # elements×masks
    elunionmask::BitVector          # all masked elements
    total_masked::Vector{Int}       # total number of masked elements for each mask

    # info for overlaps with all experiment masks
    # experiments are ordered by experiment index
    # if there's no overlap, the value is "missing"
    setXexp_olaps::Matrix{MaskOverlap}

    loc2glob_setix::Vector{Int}         # local set index to global set index in original mosaic
    glob2loc_setix::SparseVector{Int}   # global set index (in orig. mosaic) to local set index

    ix2experiment::Vector{E}        # experiment index to the ID of the mask
    experiment2ix::Dict{E, Int}     # experiment Id to index
end

function MaskedSetMosaic(mosaic::SetMosaic{T, S}, elmasks::AbstractMatrix{Bool},
                         setXexp_olaps::Matrix{MaskOverlap},
                         set_indices::Vector{Int},
                         experiment_ids::Union{AbstractVector{E}, Nothing} = nothing
) where {T, S, E}
    size(elmasks, 1) == nelements(mosaic) ||
        throw(ArgumentError("Elements mask length ($(size(elmasks, 1))) should match the number of elements ($(nelements(mosaic)))"))
    nmasks = size(elmasks, 2)
    (experiment_ids === nothing) || (nmasks == length(experiment_ids)) ||
        throw(ArgumentError("Number of masks ($nmasks) should match the number of mask IDs ($(length(masks_ids)))"))
    bit_elmasks = convert(BitMatrix, elmasks)
    bit_elunionmask = dropdims(any(bit_elmasks, dims=2), dims=2)
    @assert length(bit_elunionmask) == size(bit_elmasks, 1)
    MaskedSetMosaic(
        mosaic, bit_elmasks, bit_elunionmask,
        dropdims(sum(elmasks, dims=1), dims=1),
        setXexp_olaps,
        set_indices, SparseVector(nsets(mosaic), set_indices, collect(eachindex(set_indices))),
        collect(experiment_ids !== nothing ? experiment_ids : 1:nmasks),
        # FIXME use vector container if experiment_ids === nothing or if ids are integers
        Dict(id => ix for (ix, id) in enumerate(experiment_ids !== nothing ? experiment_ids : 1:nmasks)))
end

"""
    mask(mosaic::SetMosaic, elmasks::AbstractMatrix{Bool};
         [experiment_ids::Union{AbstractVector, AbstractSet, Nothing} = nothing],
         [min_nmasked=1], [max_setsize=nothing],
         [max_overlap_logpvalue=0.0]) -> MaskedSetMosaic

Construct [`MaskedSetMosaic`](@ref) from the [`SetMosaic`](@ref) and the collection of element masks.

## Arguments
* `min_nmasked`: the minimal number of masked elements in a set to include in the mosaic
* `max_setsize` (optional): ignore the annotation sets bigger than the specified size
* `max_overlap_logpvalue`: the threshold of Fisher's Exact Test log P-value of the overlap
   between the set and the mask for the inclusion of the set into the mosaic.
   *0* accepts all sets.
"""
function mask(mosaic::SetMosaic, elmasks::AbstractMatrix{Bool};
              experiment_ids::Union{AbstractVector, Nothing} = nothing,
              min_nmasked::Integer=1, max_setsize::Union{Integer, Nothing} = nothing,
              max_overlap_logpvalue::Number=0.0, # 0.0 would accept any overlap (as log(Fisher Exact Test P-value))
              max_min_overlap_logpvalue::Number=max_overlap_logpvalue
)
    size(elmasks, 1) == nelements(mosaic) ||
        throw(ArgumentError("Elements mask rows ($(size(elmasks, 1))) should match the number of elements ($(nelements(mosaic)))"))

    # calculate overlap sizes for each set and each mask
    nmasked_orgsets = nmasked_perset(mosaic, elmasks)

    # get the sets that pass the overlap filters
    setixs = Vector{Int}()
    ntotal = nelements(mosaic)
    ntotal_masked = dropdims(sum(elmasks, dims=1), dims=1)
    @inbounds for (setix, nmasked_permask) in enumerate(eachslice(nmasked_orgsets, dims=1))
        nset = setsize(mosaic, setix)
        min_logpvalue = 0.0
        for (expix, nmasked) in enumerate(nmasked_permask)
            (nmasked < min_nmasked) && continue # skip small overlap
            (max_setsize !== nothing) && (nset > max_setsize) && continue # skip very big and generic sets
            overlap_logpvalue = logpvalue(nmasked, nset, ntotal_masked[expix], ntotal)
            (overlap_logpvalue > max_overlap_logpvalue) && continue # skip non-signif overlaps
            min_logpvalue = min(min_logpvalue, overlap_logpvalue)
            if min_logpvalue <= max_min_overlap_logpvalue
                # set has passed the criteria
                push!(setixs, setix)
                break
            end
        end
    end

    # fill the overlaps matrix
    setXexp_olaps = Matrix{MaskOverlap}(undef, length(setixs), size(elmasks, 2))
    @inbounds for (expix, nmasked_perset) in enumerate(eachslice(nmasked_orgsets, dims=2))
        mask_olaps = view(setXexp_olaps, :, expix)
        masksize = ntotal_masked[expix]
        for (i, setix) in enumerate(setixs)
            nset = setsize(mosaic, setix)
            nmasked = nmasked_perset[setix]
            overlap_pvalue = logpvalue(nmasked, nset, masksize, ntotal)
            mask_olaps[i] = MaskOverlap(nmasked, nset - nmasked,
                                        overlap_pvalue <= max_overlap_logpvalue)
        end
    end

    return MaskedSetMosaic(mosaic, elmasks, setXexp_olaps, setixs, experiment_ids)
end

function mask(mosaic::SetMosaic{T}, elmasks::AbstractVector #= iterable with eltype()==Set{T} =#;
              kwargs...) where T
    @assert eltype(elmasks) === Set{T}
    mask(mosaic, Bool[in(e, elmask::Set{T}) for e in mosaic.ix2elm, elmask in elmasks]; kwargs...)
end

function mask(mosaic::SetMosaic{T}, elmasks::AbstractDict #= iterable with eltype()==Set{T} =#;
              kwargs...) where T
    @assert valtype(elmasks) === Set{T}
    mask(mosaic, Bool[in(e, elmask::Set{T}) for e in mosaic.ix2elm, elmask in values(elmasks)];
         experiment_ids = collect(keys(elmasks)), kwargs...)
end

nelements(mosaic::MaskedSetMosaic) = nelements(mosaic.original)
nsets(mosaic::MaskedSetMosaic) = size(mosaic.setXexp_olaps, 1) # only sets overlapping with masks
nsets(mosaic::MaskedSetMosaic, expix::Integer) = # number of sets overlapping with given mask
    sum(olap -> olap.used && olap.nmasked > 0,
        view(mosaic.setXexp_olaps, :, expix))
nmasks(mosaic::MaskedSetMosaic) = length(mosaic.total_masked)
nexperiments(mosaic::MaskedSetMosaic) = nmasks(mosaic)

_nmasked(mosaic::MaskedSetMosaic, expix::Int) = mosaic.total_masked[expix]
_nunmasked(mosaic::MaskedSetMosaic, expix::Int) = nelements(mosaic) - mosaic.total_masked[expix]

nmasked(mosaic::MaskedSetMosaic, mask::Any) = _nmasked(mosaic, mosaic.experiment2ix[mask])
nunmasked(mosaic::MaskedSetMosaic, mask::Any) = _nunmasked(mosaic, mosaic.experiment2ix[mask])

# FIXME swap exp and set args order
function overlap(mosaic::MaskedSetMosaic, exp::Any, set::Any)
    expix = mosaic.experiment2ix[exp] # throws if mask is not defined
    globsetix = mosaic.original.set2ix[set] # throws if set is not defined
    setix = mosaic.glob2loc_setix[globsetix]
    (setix == 0) && return missing # no overlap data for given set
    return mosaic.setXexp_olaps[setix, expix]
end

# FIXME swap mask and set args order
function logpvalue(mosaic::MaskedSetMosaic, exp::Any, set::Any)
    olap = overlap(mosaic, exp, set)
    ismissing(olap) && return missing
    return logpvalue(olap.nmasked, olap.nmasked + olap.nunmasked,
                     nmasked(mosaic, exp), nelements(mosaic))
end

function _setweight(mosaic::MaskedSetMosaic, setix::Integer, expix::Integer;
                    filter::Bool = true)
    olap = mosaic.setXexp_olaps[setix, expix]
    (ismissing(olap) || (filter && !olap.used)) && return missing
    return logpvalue(olap.nmasked, olap.nmasked + olap.nunmasked,
                     _nmasked(mosaic, expix), nelements(mosaic))
end

# copy everything, except the original mosaic (leave the reference to the same object)
Base.copy(mosaic::MaskedSetMosaic) =
    MaskedSetMosaic(mosaic.original, copy(mosaic.elmasks), copy(mosaic.elunionmask),
                    copy(mosaic.total_masked),
                    copy(mosaic.setXexp_olaps),
                    copy(mosaic.loc2glob_setix), copy(mosaic.glob2loc_setix),
                    copy(mosaic.ix2experiment), copy(mosaic.experiment2ix))


text:
using Compat.Test
include("../qcgan.jl")

@testset "quantum circuit gan - opdiff" begin
    N = 3
    target = rand_state(N)
    pairs_gen = [1=>2, 2=>3]
    pairs_dis = [1=>2, 2=>3, 3=>4]
    qcg = qcgan(target, 2, 2, pairs_gen, pairs_dis)
    op_expect = ()->expect(qcg.witness_op, apply!(copy(qcg.reg0), qcg.circuit))
    ggrad = opgrad(op_expect, qcg.grots)
    dgrad = opgrad(op_expect, qcg.drots)
    println(ggrad)
    @test isapprox(ggrad, num_gradient(op_expect, qcg.grots), atol=1e-4)
    @test isapprox(dgrad, num_gradient(op_expect, qcg.drots), atol=1e-4)
end

@testset "quantum circuit gan" begin
    N = 3
    target = rand_state(N)
    pairs_gen = [1=>2, 2=>3]
    pairs_dis = [1=>2, 2=>3, 3=>4]
    qcg = qcgan(target, 2, 2, pairs_gen, pairs_dis)
    ggrad, dgrad = gradient(qcg)
    nggrad= num_gradient(loss(qcg), qcg.grots)
    ndgrad= num_gradient(loss(qcg), qcg.drots)
    println(ndgrad)
    @test isapprox(ggrad, nggrad, atol=1e-3)
    @test isapprox(dgrad, ndgrad, atol=1e-3)
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
# try out some ideas for sequences
#   starting with the direct sequence
using Polylogarithms
using DataFrames, CSV
using PyPlot
using Printf
include("utilities.jl")
using SpecialFunctions
using Base.MathConstants
twoπ = 2.0*π

# chose a "hard value"
s = complex(1.5)
z = 10.0 + 10.0im
x = log(-z)

result1 = polylog(s, z; min_iterations=1000 )
result2 = Polylogarithms.polylog_series_4( s, z )

max_n = 10
results = zeros(Complex{Float64}, max_n)
rel_errors = zeros(Float64, max_n)
for k = 0:1:max_n-1
    results[k+1] = (-1.0)^k * (1.0 - 2.0^(1-2*k)) * twoπ^(2*k) * bernoulli(2*k,0.0) * x^(s-2*k) / ( SpecialFunctions.gamma(2*k+1) * SpecialFunctions.gamma(s+1-2*k) )
    # rel_errors[k] = abs( results[k] - Li) / abs(Li)
end

a = diff(results)
a = results

figure()
k = 0:1:max_n-1
semilogy(k, abs.(a), "-")
plot(k, abs.(real.(a)), "-.")
plot(k, abs.(imag.(a)),"--")
# plot(k_star, abs.(a[k_star]), "o")


text:
@testset "Problems" begin
  data2D = readgeotable(joinpath(datadir,"data2D.tsv"), coordnames=(:x,:y))
  data3D = readgeotable(joinpath(datadir,"data3D.tsv"))
  grid2D = RegularGrid(100,100)
  grid3D = RegularGrid(100,100,100)

  @testset "Estimation" begin
    # test basic problem interface
    problem3D = EstimationProblem(data3D, grid3D, :value)
    @test data(problem3D) == data3D
    @test domain(problem3D) == grid3D
    @test variables(problem3D) == Dict(:value => Float64)

    # problems with missing data have types inferred correctly
    img = Array{Union{Float64,Missing}}(rand(10,10))
    mdata = georef((var=img,))
    problem = EstimationProblem(mdata, grid2D, :var)
    @test variables(problem) == Dict(:var => Float64)

    # show methods
    problem2D = EstimationProblem(data2D, grid2D, :value)
    @test sprint(show, problem2D) == "2D EstimationProblem"
    @test sprint(show, MIME"text/plain"(), problem2D) == "2D EstimationProblem\n  data:      3 SpatialData{Float64,2}\n  domain:    100×100 RegularGrid{Float64,2}\n  variables: value (Float64)"

    if visualtests
      @plottest plot(problem2D,ms=2) joinpath(datadir,"estimation.png") !istravis
    end
  end

  @testset "Simulation" begin
    # test basic problem interface
    problem3D = SimulationProblem(data3D, grid3D, :value, 100)
    @test data(problem3D) == data3D
    @test domain(problem3D) == grid3D
    @test variables(problem3D) == Dict(:value => Float64)
    @test hasdata(problem3D)
    @test nreals(problem3D) == 100

    # problems with missing data have types inferred correctly
    img = Array{Union{Float64,Missing}}(rand(10,10))
    mdata = georef((var=img,))
    problem = SimulationProblem(mdata, grid2D, :var, 3)
    @test variables(problem) == Dict(:var => Float64)

    # specify type of variable explicitly
    problem = SimulationProblem(data3D, grid3D, :value => Float64, 100)
    @test variables(problem) == Dict(:value => Float64)

    # add variable not present in spatial data
    problem = SimulationProblem(data3D, grid3D, (:value => Float64, :other => Int), 100)
    @test variables(problem) == Dict(:value => Float64, :other => Int)
    @test isempty(datamap(problem)[:other])

    # infer type of variables in spatial data whenever possible
    problem = SimulationProblem(data3D, grid3D, (:value, :other => Int), 100)
    @test variables(problem) == Dict(:value => Float64, :other => Int)
    @test isempty(datamap(problem)[:other])

    # constructors without spatial data require variables with types
    problem = SimulationProblem(grid3D, :value => Float64, 100)
    @test variables(problem) == Dict(:value => Float64)
    @test isempty(datamap(problem)[:value])
    @test_throws MethodError SimulationProblem(grid3D, :value, 100)

    # show methods
    problem2D = SimulationProblem(data2D, grid2D, :value, 100)
    @test sprint(show, problem2D) == "2D SimulationProblem (conditional)"
    @test sprint(show, MIME"text/plain"(), problem2D) == "2D SimulationProblem (conditional)\n  data:      3 SpatialData{Float64,2}\n  domain:    100×100 RegularGrid{Float64,2}\n  variables: value (Float64)\n  N° reals:  100"

    if visualtests
      @plottest plot(problem2D,ms=2) joinpath(datadir,"simulation.png") !istravis
    end
  end

  @testset "Learning" begin
    Random.seed!(123)
    sdata = georef((x=rand(10),y=rand(10),z=rand(10)), 10rand(2,10))
    tdata = georef((x=rand(10,10),))
    rtask = RegressionTask(:x, :y)
    ctask = ClusteringTask(:x, :c)

    # test basic problem interface
    problem = LearningProblem(sdata, tdata, rtask)
    @test sourcedata(problem) == sdata
    @test targetdata(problem) == tdata
    @test task(problem) == rtask

    # show methods
    problem = LearningProblem(sdata, tdata, ctask)
    @test sprint(show, problem) == "2D LearningProblem"
    @test sprint(show, MIME"text/plain"(), problem) == "2D LearningProblem\n  source: 10 SpatialData{Float64,2}\n  target: 100 SpatialData{Float64,2}\n  task:   Clustering x → c"

    if visualtests
      @plottest plot(problem,ms=2) joinpath(datadir,"learning.png") !istravis
    end
  end
end


text:
using Triangulation_diagrams
using PyPlot

N = [ 0  1  2  4  4  5  3  3  2  1  5  6  6  2  6  4
      2  1  0  0  2  1  1  3  2  3  3  0  4  4  2  4 ]

T = [  1  3  3  5  5
       3  5 12 12 13
      14 14  5 13 14
       9  8  6 15 16
      10  9  7 11  8
       2  7  4  6 11 ]

E = [ 14 10  1  2  3  4 12 15 13 16
      10  1  2  3  4 12 15 13 16 14 ]

QN = 6

tri = Triangulation(N, T, E, QN)

figure(1)
draw_triangles(tri)
enumerate_vertices(tri)
enumerate_triangles(tri)
enumerate_bdry_edges(tri, 0.2)
mark_first_local_vertex(tri, 0.3)
axis("off")
savefig("quad_triangulation.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
module Test_pendulum_dae

#=
SavingCallback for DAEs: Not clear how to get du #547

https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/547

=#

using DifferentialEquations, Measurements, StaticArrays, PyPlot, Sundials

mutable struct Model{R<:Real}
    # parameters
    m::R  # mass
    L::R  # length
    d::R  # damping
    g::R  # gravity

	# variables
    ϕ::R  # angle
    ω::R  # der(angle)
    α::R  # der(der(angle))
	residue::R
    r::SVector{2,R}   # absolute position of end point

    der_x::Vector{R}

    function Model{R}(m::R, L::R, d::R, g::R) where {R<:Real}
        new(m, L, d, g, 0, 0, 0, 0, SVector{2,R}(0,0))
    end
end

Model(;m::R=0, L::R=0, d::R=0, g::R=0) where {R<:Real} = Model{R}(m,L,d,g)

"simplePendulum: Compute m::Model at actual time instant"
function pendulum!(m::Model, t, x, der_x)::Nothing
    m.ϕ = x[1]
    m.ω = x[2]
    m.α = (-m.m*m.g*m.L*sin(m.ϕ) - m.d*m.ω) / (m.m*m.L^2)
	m.residue = der_x[2] - m.α
    m.r = @SVector [m.L*sin(m.ϕ), -m.L*cos(m.ϕ)]
    return nothing
end

"SimplePendulum1: Called by integrator"
function pendulum1!(residue, der_x, x, m::Model, t)
    pendulum!(m, t, x, der_x)
    residue[1] = der_x[1] - m.ω
    residue[2] = m.residue
end

"SimplePendulum2: Called at every communication point"
function pendulum2!(x, t, integrator)
    m::Model = integrator.p
    # tmps = get_tmp_cache(integrator)
    # println("typeof(tmps) = ", typeof(tmps), ", length(tmps) = ", length(tmps))
    if t==0
        m.der_x = copy(integrator.du)   # Since IDA gives an error for integrator(t, Val{1]}) at the initial time instant
    else
        integrator(m.der_x, t, Val{1})
    end
    pendulum!(m, t, x, m.der_x)
    return deepcopy(m)
end

t_start = 0.0
t_inc   = 0.01
t_end   = 6.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

# Problem 1 with Float64
model1 = Model(m=1.0, L=1.0, d=0.2, g=9.81)
println("model1 = $model1")
x1₀     = [π/2, 0.0]

    # Initialize
    der_x1₀ = [0.0, 0.0]
    pendulum!(model1, 0.0, x1₀, der_x1₀)
    der_x1₀[1] = model1.ω
    der_x1₀[2] = model1.α

saved_values1=SavedValues(typeof(x1₀[1]), Model)
cb1 = SavingCallback(pendulum2!, saved_values1, saveat=tspan2)

prob1 = DAEProblem(pendulum1!, der_x1₀, x1₀, tspan, model1)
sol1 = solve(prob1, IDA(), reltol = 1e-6, saveat=t_inc, callback=cb1, save_timeseries=false)

#@time sol1 = solve(prob1, IDA(), reltol = 1e-6, saveat=t_inc)
#@time sol1 = solve(prob1, IDA(), reltol = 1e-6, saveat=t_inc)

# switch off storing of x, der_x in sol1 (but keep in saved_values1)
# sol1 = solve(prob1, IDA(), reltol = 1e-6, saveat=t_end, callback=cb1, save_timeseries=false)

t = saved_values1.t
m = saved_values1.saveval

figure(1)
clf()
plot(sol1.t, getindex.(sol1.u, 1), label="\$\\varphi\$")
#ϕ = getfield.(m, :ϕ)
#plot(t, ϕ, label="\$\\varphi\$")
grid(true)
legend()
title("Pendulum with Float64")

figure(2)
clf()

r = getfield.(m, :r)
r1 = getindex.(r,1)
r2 = getindex.(r,2)
plot(t, r1, label="r1")
plot(t, r2, label="r2")
grid(true)
legend()
title("Pendulum with Float64")

figure(3)
clf()
residue = getfield.(m, :residue)
plot(t, residue, label="residue")
grid(true)
legend()
title("Pendulum with Float64")

end


text:
# This file was generated, do not modify it. # hide
KNNClassifier = @load KNNClassifier

knnc = KNNClassifier(K=1)
classif = machine(knnc, X3, y)
fit!(classif, rows=train)
ŷ = predict_mode(classif, rows=test)
accuracy(ŷ, y[test]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using JUDI, SegyIO, LinearAlgebra, PyPlot, IterativeSolvers, JOLI

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v[:,Int(round(end/2)):end] .= 5f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
time = 2000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Source wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(time, dt, f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

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

# Write shots as segy files to disk
opt = Options()

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)

# Random weights (size of the model)
w = judiWeights(randn(Float32, model.n))

# Create operator for injecting the weights, multiplied by the provided wavelet(s)
Pw = judiLRWF(info, wavelet)

# Model observed data w/ extended source
lambda = 1f2
I = joDirac(info.n, DDT=Float32, RDT=Float32)
F = Pr*F*adjoint(Pw)
F̄ = [F; lambda*I]

# Simultaneous observed data
d_sim = F*w

# # Adjoint operation
w_adj = adjoint(F)*d_sim

# # LSQR
w_inv = 0f0 .* w
w_inv_no_damp = 0f0 .* w
lsqr!(w_inv, F̄, [d_sim; lambda*w]; maxiter=2, verbose=true, damp=1e2)
lsqr!(w_inv_no_damp, F, d_sim; maxiter=2, verbose=true, damp=1e2)

d_pred = F*w_inv;
d_pred_no_damp = F*w_inv_no_damp;

# Plot results
figure()
subplot(1,3,1)
imshow(d_sim.data[1], vmin=-5e2, vmax=5e2, cmap="gray"); title("Observed data")
subplot(1,3,2)
imshow(d_pred.data[1], vmin=-2e2, vmax=2e2, cmap="gray"); title("Predicted data")
subplot(1,3,3)
imshow(d_pred_no_damp.data[1], vmin=-2e2, vmax=2e2, cmap="gray"); title("Predicted data no damp")

figure()
subplot(2,2,1)
imshow(w.weights[1], vmin=-3, vmax=3, cmap="gray"); title("Weights")
subplot(2,2,2)
imshow(w_adj.weights[1], vmin=minimum(w_adj), vmax=maximum(w_adj), cmap="gray"); title("Adjoint")
subplot(2,2,3)
imshow(w_inv.weights[1], vmin=minimum(w_inv), vmax=maximum(w_inv), cmap="gray"); title("D-LSQR")
subplot(2,2,4)
imshow(w_inv_no_damp.weights[1], vmin=minimum(w_inv), vmax=maximum(w_inv), cmap="gray"); title("LSQR")

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using Printf: @printf
using MPI

function do_broadcast()
    comm = MPI.COMM_WORLD
    root = 0
    if MPI.Comm_rank(comm) == root
        println(repeat("-",78))
        println(" Running on $(MPI.Comm_size(comm)) processes")
        println(repeat("-",78))
    end

    MPI.Barrier(comm)

    N = 5

    if MPI.Comm_rank(comm) == root
        A = fill((1.0 + 2.0im), N)
    else
        A = Array{ComplexF64}(undef, N)
    end

    MPI.Bcast!(A, length(A), root, comm)

    @printf("[%02d] A:%s\n", MPI.Comm_rank(comm), A)

    if MPI.Comm_rank(comm) == root
        B = Dict("foo" => "bar")
    else
        B = nothing
    end

    B = MPI.bcast(B, root, comm)
    @printf("[%02d] B:%s\n", MPI.Comm_rank(comm), B)

    if MPI.Comm_rank(comm) == root
        f = x -> x^2 + 2x - 1
    else
        f = nothing
    end

    f = MPI.bcast(f, root, comm)
    @printf("[%02d] f(3):%d\n", MPI.Comm_rank(comm), f(3))
end

function main()
    MPI.Init()
    do_broadcast()
    MPI.Finalize()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DataAugmentation
using Images
using StaticArrays
using Test
using TestSetExtensions
using CoordinateTransformations

using DataAugmentation: Item, Transform, getrandstate, itemdata, setdata, ComposedProjectiveTransform,
    cropindices, makebounds, getprojection, boundsextrema, boundsranges, boundssize, offsetcropindices,
    CroppedProjectiveTransform, getbounds, project, project!, makebuffer, imagetotensor, imagetotensor!,
    normalize, normalize!, tensortoimage, denormalize, denormalize!
import DataAugmentation: apply, compose