text:
"""
pdeF!.jl

This file contains everything you need to run the Ellptic PDE examples.
This includes the version with an explict sparse matrix Jacobian and
the fixed point formulations using the fish2d.jl preconditioner.

I've also parked the exact solution in here so you can do the grid refinement
study.

Look at pdeinit for the construction of the precomputed data. There is
a lot of it.

If you only want to run the examples, you should not have to look
at the code.
"""
### And now for the functions ...
"""
pdeF!(FV, u, pdata)

Residual using sparse matrix-vector multiplication
"""
function pdeF!(FV, u, pdata)
    D2 = pdata.D2
    CV = pdata.CV
    rhs = pdata.RHS
    p1 = pdata.jvect1
    #    FV .= D2 * u + 20.0 * u .* (CV * u) - rhs
    #    FV .= D2*u
    #    p1 .= CV*u
    mul!(FV, D2, u)
    mul!(p1, CV, u)
    p1 .*= 20.0
    p1 .*= u
    FV .+= p1
    FV .-= rhs
end

"""
pdeJ!(FP, F, u, pdata)

Sparse matrix Jacobian. The package does not do its own sparse
differencing. The Jacobian for this problem is easy enough to
compute analytically.
"""
function pdeJ!(FP, F, u, pdata)
    D2 = pdata.D2
    CV = pdata.CV
    CT = pdata.CT
    cu = CV * u
    #DC=spdiagm(0 => 20*cu); DU=spdiagm(0 => 20*u)
    DC = Diagonal(20 * cu)
    DU = Diagonal(20 * u)
    #
    # The easy way to compute the Jacobian is
    #FP .= D2 + DU*CV + DC
    # but you allocate yourself silly with that one.
    # So we preallocate room for DU*CV in CT and sum the terms for FP
    # one at a time. I have to use Diagonal instead of spdiagm if I want
    # mul! to work fast.
    #
    FP .= D2
    FP .+= DC
    mul!(CT, DU, CV)
    #CT .= CV; lmul!(DU,CT);
    FP .+= CT
    # I should be able to do mul!(FP,DU,CV), but it's 1000s of times slower.
end

"""
Jvec2d(v, FS, u, pdata)

Analytic Jacobian-vector product for PDE example
"""
function Jvec2d(v, FS, u, pdata)
    D2 = pdata.D2
    CV = pdata.CV
    CT = pdata.CT
    jvec = pdata.jvec
    p1 = pdata.jvect1
    #    jvec .= D2 * v
    #    p1 .= CV * u
    mul!(jvec, D2, v)
    mul!(p1, CV, u)
    p1 .*= 20.0
    p1 .*= v
    jvec .+= p1
    #    p1 .= CV * v
    mul!(p1, CV, v)
    p1 .*= 20.0
    p1 .*= u
    jvec .+= p1
    return jvec
end

"""
hardleft!(FV, u, pdata)
Convection-diffusion equation with left preconditioning hard-wired in

"""
function hardleft!(FV, u, pdata)
fdata=pdata.fdata
FV = pdeF!(FV,u,pdata)
FV .= Pfish2d(FV,fdata)
return FV
end

"""
hardleftFix!(FV, u, pdata)
Fixed point form of the left preconditioned nonlinear
convection-diffusion equation
"""
function hardleftFix!(FV, u, pdata)
FV = hardleft!(FV, u, pdata)
# G(u) = u - FV
axpby!(1.0, u, -1.0, FV)
return FV
end

"""
pdeinit(n)

collects the precomputed data for the elliptic pde example. This
includes

- the sparse matrix representation of the operators,
- the right side of the equation,
- the exact solution,
- the data that the fft-based fast Poisson solver (fish2d) needs
"""
function pdeinit(n)
    # Make the grids
    n2 = n * n
    h = 1.0 / (n + 1.0)
    x = collect(h:h:1.0-h)
    # collect the operators
    D2 = Lap2d(n)
    DX = Dx2d(n)
    DY = Dy2d(n)
    CV = (DX + DY)
    # I need a spare sparse matrix to save allocations in the Jacobian computation
    CT = copy(CV)
    # Exact solution and its derivatives
    uexact = solexact(x)
    dxe = dxexact(x)
    dye = dyexact(x)
    d2e = l2dexact(x)
    dxv = reshape(dxe, (n2,))
    dyv = reshape(dye, (n2,))
    d2v = reshape(d2e, (n2,))
    uv = reshape(uexact, (n2,))
    fdata = fishinit(n)
    # The right side of the equation
    RHS = d2v + 20.0 * uv .* (dxv + dyv)
    # preallocate a few vectors
    jvec = zeros(n2)
    jvect1 = zeros(n2)
    # Pack it and ship it.
    pdedata =
        (D2 = D2, CV = CV, CT = CT, RHS = RHS, jvec, jvect1, fdata = fdata, uexact = uexact)
end

"""
This collection of functions
builds u, u_x, u_y, and the negative Laplacian for the
example problem in the book. Here
u(x,y) = 10 x y (1-x)(1-y) exp(x^4.5)

which is the example from FA01.
"""

function w(x)
    w = 10.0 * x .* (1.0 .- x) .* exp.(x .^ (4.5))
end

function wx(x)
    wx = 4.5 * (x .^ (3.5)) .* w(x) + 10.0 * exp.(x .^ (4.5)) .* (1.0 .- 2.0 * x)
end

function wxx(x)
    wxx =
        (4.5 * 3.5) * (x .^ (2.5)) .* w(x) +
        4.5 * (x .^ (3.5)) .* wx(x) +
        +10.0 * 4.5 * (x .^ (3.5)) .* exp.(x .^ (4.5)) .* (1.0 .- 2.0 * x) +
        -20.0 * exp.(x .^ (4.5))
end

function v(x)
    v = x .* (1.0 .- x)
end

function vx(x)
    vx = 1.0 .- 2.0 * x
end

function vxx(x)
    vxx = -2.0 * ones(size(x))
end

function solexact(x)
    solexact = w(x) * v(x)'
end

function l2dexact(x)
    l2dexact = -(w(x) * vxx(x)') - (wxx(x) * v(x)')
end

function dxexact(x)
    dxexact = wx(x) * v(x)'
end

function dyexact(x)
    dxexact = w(x) * vx(x)'
end


text:
import Genie
import Logging, LoggingExtras
import Dates

function initialize_logging()
  date_format = "yyyy-mm-dd HH:MM:SS"

  logger =  if Genie.config.log_to_file
              isdir(Genie.config.path_log) || mkpath(Genie.config.path_log)
              LoggingExtras.TeeLogger(
                LoggingExtras.FileLogger(joinpath(Genie.config.path_log, "$(Genie.config.app_env)-$(Dates.today()).log"), always_flush = true, append = true),
                Logging.ConsoleLogger(stdout, Genie.config.log_level)
              )
            else
              Logging.ConsoleLogger(stdout, Genie.config.log_level)
            end

  timestamp_logger(logger) = LoggingExtras.TransformerLogger(logger) do log
    merge(log, (; message = "$(Dates.format(Dates.now(), date_format)) $(log.message)"))
  end

  LoggingExtras.TeeLogger(LoggingExtras.MinLevelLogger(logger, Genie.config.log_level)) |> timestamp_logger |> global_logger

  nothing
end

initialize_logging()

ext: jl


text:
import ..Interventions: intervene
import ..Util: mapf
export rid, ∥ᵈ
"Represents the parent variables of `child`"
struct Parents{T}
  child::T
end

"""
`rid(x, θ)`

Random Interventional Distribution of `x` given `θ`

# Arguments
`x` - random variable
`θ` - random variable that is causal parent of `x`

# Returns
random-variable valued random variable defined as:
  `ω -> intervene(x, θ => θ(ω))`
"""
rid(x, θ) = ω -> intervene(x, θ => θ(ω))

rid(x, θs::Tuple) = ω -> let θs_ = mapf(ω, θs)
                          intervene(x, θs, θs_)
                        end

@inline x ∥ᵈ θ = rid(x, θ)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
export TabularLearner

"""
    TabularLearner{S, T}

Use a `Dict{S,Vector{T}}` to store action probabilities.
"""
struct TabularLearner{S,T} <: AbstractLearner
    table::Dict{S,Vector{T}}
end

TabularLearner() = TabularLearner{Int,Float32}()
TabularLearner{S}() where {S} = TabularLearner{S,Float32}()
TabularLearner{S,T}() where {S,T} = TabularLearner(Dict{S,Vector{T}}())

(p::TabularLearner)(env::AbstractEnv) = p(ChanceStyle(env), env)

function (p::TabularLearner)(::ExplicitStochastic, env::AbstractEnv)
    if get_current_player(env) == get_chance_player(env)
        [a.prob for a::ActionProbPair in get_actions(env)]
    else
        p(DETERMINISTIC, env)  # treat it just like a normal one
    end
end

function (p::TabularLearner)(::RLBase.AbstractChanceStyle, env::AbstractEnv)
    get!(p.table, get_state(env)) do
        n = length(get_legal_actions(env))
        fill(1 / n, n)
    end
end

RLBase.update!(p::TabularLearner, experience::Pair) =
    p.table[first(experience)] = last(experience)


text:
struct FixedTimeZoneF <: TimeZone
    local_minus_utc::Int32
    function FixedTimeZoneF(secs::Int)
        if (-86_400 < secs) && (secs < 86_400)
           new(secs)
        else
            error("Invalid offset")
        end
    end
end

FixedTimeZoneF(s::Second) = FixedTimeZoneF(s.value)

function name(tz::FixedTimeZoneF)
    offset = tz.local_minus_utc

    # TODO: could use offset_string in utcoffsets.jl with some adaptation
    if offset < 0
        sig = '-'
        offset = -offset
    else
        sig = '+'
    end

    hour, rem = divrem(offset, 3600)
    minute, second = divrem(rem, 60)

    if hour == 0 && minute == 0 && second == 0
        name = "UTC"
    elseif second == 0
        name = @sprintf("UTC%c%02d:%02d", sig, hour, minute)
    else
        name = @sprintf("UTC%c%02d:%02d:%02d", sig, hour, minute, second)
    end
    return name
end


text:
using AffineArithmetic
using IntervalArithmetic
using Polynomials

using Test

@testset "Constructor" begin
    C = Affine(0.0)
    @test C == Affine(0.0, Float64[])

    C = Affine(1.0)
    @test C == Affine(1.0, Float64[])

    C = Affine(1.0, [3.0, 4.0])
    @test C.c == 1.0
    @test C.γ == [3.0, 4.0]

    reset_affine_index()
end

reset_affine_index()

@testset "Construction from intervals" begin

    X = 1..3
    X_a = Affine(X)
    @test X_a == Affine(2.0, [1.0])

    Y = 2..4
    Y_a = Affine(Y)
    @test Y_a == Affine(3.0, [0.0, 1.0])

    @test X_a + Y_a == Affine(5.0, [1.0, 1.0])
    @test X_a * Y_a == Affine(6.0, [3.0, 2.0, 1.0])
end

@testset "Small powers and range" begin

    X = 1..3
    Y = Affine(X)

    f(x) = x^2 - x + 1

    @test range(f(X)) == -1..9
    @test range(f(Y)) == -1..7   # affine is a little better

    g(x) = (x - 1)^3

    @test range(g(X)) == 0..8
    @test range(g(Y)) == -6..8   # affine is significantly worse
end

reset_affine_index()

@testset "Matrix multiplication" begin
    A = 0.5 * [1 2; -1 1]

    X = Affine(-1..1)
    Y = Affine(-1..1)

    XX = [X, Y]
    @test XX == [Affine(0.0, [1.0]), Affine(0.0, [0.0, 1.0])]

    @test A * XX == [Affine(0.0, [0.5, 1.0]), Affine(0.0, [-0.5, 0.5])]
end

reset_affine_index()

@testset "Range of polynomial" begin
    # example from Rump 2015
    p = Poly([-3, 1])   # x - 3
    p2 = p^8

    x = 4 ± 1e-4
    y = Affine(x)

    @test (-70..70) ⊆ range(p2(x))
    @test range(p2(y)) ⊆ 0.998..1.002   # affine is extremely much better!
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

#=
Williams A, Barrus S, Morley R K, et al., 2005.
(https://dl.acm.org/doi/abs/10.1145/1198555.1198748)
=#
function intersection(f, r::Ray{Dim,T}, b::Box{Dim,T}) where {Dim,T}
  invdir = one(T) ./ direction(r)
  lo, up = coordinates.(extrema(b))
  orig = coordinates(origin(r))

  tmin = zero(T)
  tmax = typemax(T)

  # check for intersection with slabs along with each axis
  for i in 1:Dim
    imin = (lo[i] - orig[i]) * invdir[i]
    imax = (up[i] - orig[i]) * invdir[i]

    # swap variables if necessary
    invdir[i] < zero(T) && ((imin, imax) = (imax, imin))

    # the ray is on a face of the box, avoid NaN
    (isnan(imin) || isnan(imax)) && continue

    (tmin > imax || imin > tmax) && return @IT NoIntersection nothing f

    tmin = max(tmin, imin)
    tmax = min(tmax, imax)
  end

  tmin ≈ tmax && return @IT TouchingRayBox r(tmin) f

  return @IT CrossingRayBox Segment(r(tmin), r(tmax)) f
end


text:
"""
    cosine_vec(wv::WordVectors, wordvector, n=10 [;vocab=nothing])

Compute the cosine similarities and return best `n` positions and
calculated values between `wordvector` and the word vectors from `wv`.
A vocabulary mask `vocab` can be specified to consider only a subset
of word vectors.
"""
function cosine_vec(wv::WordVectors, wordvector, n=10; vocab=nothing)
    if vocab === nothing
        _vocab = wv.vocab
    else
        _vocab = intersect(vocab, wv.vocab)
    end
    idxs = [wv.vocab_hash[w] for w in _vocab]
    metrics = view(wv.vectors, :, idxs)' * wordvector
    topn_positions = sortperm(metrics, rev = true)[1:min(n, length(idxs))]
    return idxs[topn_positions], metrics[topn_positions]
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
# Helpful packages for working with images and factorizations
# using Pkg; Pkg.add("Images")
# using Pkg; Pkg.add("ImageMagick") # And this allows us to load JPEG-encoded images
using Images, LinearAlgebra, Interact

# ------------------------------------------------------------------------------------------
# ### Using a SVD to compress an image
#
# In this exercise, we'll use a singular value decomposition (SVD) to compress an image --
# so that we can store an image without keeping around "unnecessary" information.
#
# To start, let's define a singular value decomposition. In a SVD, we take a matrix $A$ and
# factorize it so that
#
# $$A = USV^T$$
#
# where matrices $U$ and $V$ are unitary and hold our singular vectors. Matrix $S$ is
# diagonal and stores our singular values in decreasing order from top/left to bottom/right.
#
# In Julia, our images are stored as arrays, so we can think of `yellowbanana` as a matrix
# ------------------------------------------------------------------------------------------

file = download("https://uploads6.wikiart.org/images/salvador-dali/the-persistence-of-memory-1931.jpg!Large.jpg")

img = load(file)

size(img)

img[24,24] # Each element in the array is a color

dump(img[24,24])

# ------------------------------------------------------------------------------------------
# We can extract each "channel" of red, green, and blue and view each independently:
# ------------------------------------------------------------------------------------------

channels = Float64.(channelview(img))
Gray.(channels[1, :, :])

# ------------------------------------------------------------------------------------------
# That means we can take the SVD of this image. So, we can store this picture of a banana as
# sets of singular vectors and singular values.
#
# **The reason this is important** is that we'll find that we do **not** need to keep track
# of *all* the singular vectors and *all* the singular values to store an image that still
# looks like a banana! This means we can choose to keep only the important information,
# throw away the rest, and thereby "compress" the image.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# If we don't throw away any data, we get back what we started with:
# ------------------------------------------------------------------------------------------

U, S, V = svd(channels[1,:,:])
Gray.(U * Diagonal(S) * V')

# ------------------------------------------------------------------------------------------
# But of course we're not doing any compression here — the sizes of U, S, and V are bigger
# than our original matrix! This is like the opposite of compression.  The key is that the
# values are stored in decreasing order so we can start throwing things away.
# ------------------------------------------------------------------------------------------

sum(length.((U, S, V)))

length(img)

Gray.(U[:, 1:25] * Diagonal(S[1:25]) * V[:, 1:25]')

sum(length.((U[:, 1:25], S[1:25], V[:, 1:25])))/length(img)

# ------------------------------------------------------------------------------------------
# Of course this is just one channel of the image. Let's put it all back together and see
# how we can compress the different channels to find an acceptable compression level.
# ------------------------------------------------------------------------------------------

function rank_approx(M, k)
    U, S, V = svd(M)

    M = U[:, 1:k] * Diagonal(S[1:k]) * V[:, 1:k]'

    M = min.(max.(M, 0.0), 1.)
end

n = 100
@manipulate for k1 in 1:n, k2 in 1:n, k3 in 1:n
    colorview(  RGB,
                rank_approx(channels[1,:,:], k1),
                rank_approx(channels[2,:,:], k2),
                rank_approx(channels[3,:,:], k3)
    )
end

# ------------------------------------------------------------------------------------------
# **So how can we use a SVD to determine what information in an image is really important?**
#
# The singular values tell us!
#
# If we have matrices $U$, $S$, and $V$ from our image, we can rebuild that image with the
# matrix product $USV^T$.
#
# Taking this matrix product is the same as adding together the outer products of each
# corresponding pair of vectors from $U$ and $V$, scaled by a singular value ($\sigma$) from
# $S$. In other words, for a (100 x 100) pixel image,
#
# $$A_{image} = USV^T = \sum_{i = 1}^{100} \sigma_i \mathbf{u_i}\mathbf{v_i'} $$
#
# Every outer product $u_i * v_i'$ creates a (100 x 100) matrix. Here we're summing together
# one hundred (100 x 100) matrices in order to create the original matrix $A_{image}$. The
# matrices at the beginning of the series -- those that are scaled by **large** singular
# values -- will be **much** more important in recreating the original matrix $A_{image}$.
#
# This means we can approximate $A_{image}$ as
#
# $$A_{image} \approx \sum_{i = 1}^{n} \sigma_i \mathbf{u_i}\mathbf{v_i'}$$
#
# where $n < 100$.
#
#
# #### Exercise
#
# Write a function called `compress_image`. Its input arguments should be an image and the
# factor by which you want to compress the image. A compressed grayscale image should
# display when `compress_image` is called.
#
# For example,
#
# ```julia
# compress_image("images/104_100.jpg", 33)
# ```
#
# will return a compressed image of a grayscale banana built using 3 singular values. (This
# image has 100 singular values, so use `fld(100, 33)` to determine how many singular values
# to keep. `fld` performs "floor" division.)
#
# *Hints*:
#
# * Perform the SVD on the `channelview` of a grayscale image.
# * In an empty input cell, execute `?svd` to find a function that wil perform an SVD for
# you.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia/long-version/050 Compressing an Image.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using ProjA
using Test

@testset "ProjA.jl" begin
    # Write your own tests here.
end


text:
import Base.==
using DataFrames

export
    Edge, Node, MolecularGraph, AbstractMolGraph, AbstractNode, getNumberOfEdges, getNumberOfNodes,
    deleteNode, newNode, collectPartnerEdges, deleteEdge, newEdge, collectPartnerNodes,
    breadthFirstSearch, collectNodes, collectEdges, eachEdge, eachNode

const DFR = DataFrameRow

"""
Interface which defines some getters, setters and constructors in conjunction
with [`AbstractMolGraph`](@ref).
"""
abstract type AbstractNode end

"""
Interface for graphs which describe some topology of biological systems.
"""
abstract type AbstractMolGraph end

"""
Type to represent Edges in [`MolecularGraph`](@ref).
"""
mutable struct Edge
    source_ ::Union{AbstractNode, DFR}
    target_ ::Union{AbstractNode, DFR}
    bond_   ::Union{Bond, DFR}
    Edge(source::Union{AbstractNode, DFR}, target::Union{AbstractNode, DFR},
            bond::Union{Bond, DFR}) = new(source, target, bond)
end
Base.show(io::IO, edge::Edge) = print(io, "E[$(edge.source_) | $(edge.target_)]")
(==)(x::Edge, y::Edge) = return (x === y || x.bond_ == y.bond_)

"""
Type to represent Nodes in [`MolecularGraph`](@ref).
"""
mutable struct Node <: AbstractNode
    adjacent_edges_ ::Vector{Edge}
    atom_           ::Union{Atom, DFR}
    Node(at::Union{Atom, DFR}) = new(Edge[], at)
end
Base.show(io::IO, node::AbstractNode) = print(io, "N[$(node.atom_.serial_)]")

#atoms_to_nodes_ relates atoms to nodes in this Graph and bodns_to_edges_ bonds to edges
 #
"""
Implements [`AbstractMolGraph`](@ref). Uses `Dicts` to manage the [`Node`](@ref)s
 and [`Edge`](@ref)s in it.
"""
mutable struct MolecularGraph <: AbstractMolGraph

    atoms_to_nodes_::Dict{Union{Atom, DFR},AbstractNode}
    bonds_to_edges_::Dict{Union{Bond, DFR},Edge}
    MolecularGraph() = new(Dict{Union{Atom, DFR},AbstractNode}(), Dict{Union{Bond, DFR},Edge}())
    MolecularGraph(NodeType::Type{typ}) where typ<:AbstractNode =
                            new(Dict{Union{Atom, DFR},NodeType}(), Dict{Union{Bond, DFR},Edge}())

end
Base.show(io::IO, gr::MolecularGraph) =
    print(io, "$(length(gr.atoms_to_nodes_)) nodes, $(length(gr.bonds_to_edges_)) edges")

"""
    getFirstNode(graph::MolecularGraph)
Returns the first Node in `graph`.
"""
getFirstNode(graph::MolecularGraph) = begin
    return first(values(graph.atoms_to_nodes_))
end

"""
    getFirstNode(graph::MolecularGraph)
Returns the first Edge in `graph`.
"""
getFirstEdge(graph::MolecularGraph) = begin
    return first(values(graph.bonds_to_edges_))
end

"""
    collectNodes(graph::MolecularGraph)
Returns a Dict-Values iterator over the nodes in the graph.
"""
collectNodes(graph::MolecularGraph) = begin
    return values(graph.atoms_to_nodes_)
end

"""
    eachPartnerNode(node::AbstractNode)
Gets all the nodes, which are connected to `node` via edges.
"""
eachPartnerNode(node::AbstractNode) = begin
    return [x.source_ == node ? x.target_ : x.source_ for x in values(node.adjacent_edges_)]
end

"""
    eachPartnerEdge(node::AbstractNode)
Returns a Dict-Values iterator over the edges adjacent to `node`.
"""
eachPartnerEdge(node::AbstractNode) = begin
    return values(node.adjacent_edges_)
end

"""
    eachEdge(graph::MolecularGraph)
Returns a Dict-Values iterator over the edges in `graph`.
"""
eachEdge(graph::MolecularGraph) = begin
    return values(graph.bonds_to_edges_)
end

"""
    eachNode(graph::MolecularGraph)
Returns a Dict-Values iterator over the nodes in `graph`.
"""
eachNode(graph::MolecularGraph) = begin
    return values(graph.atoms_to_nodes_)
end

"""
    countNodes(graph::MolecularGraph)
Returns the number of Nodes in `graph`.
"""
countNodes(graph::MolecularGraph) = begin
    return length(graph.atoms_to_nodes_)
end

"""
    countEdges(graph::MolecularGraph)
Returns the number of Nodes in `graph`.
"""
countEdges(graph::MolecularGraph) = begin
    return length(graph.bonds_to_edges_)
end

"""
    newNode(graph::MolecularGraph, at::Atom, NodeType::Type{typ} = Node) where typ<:AbstractNode
Creates and inserts a node of type `typ`, carrying `at` into `graph`.
"""
newNode(graph::MolecularGraph, at::Union{Atom,DFR}, NodeType::Type{typ} = Node) where typ<:AbstractNode = begin
    if haskey(graph.atoms_to_nodes_, at)
        return false
    end
    temp::NodeType = NodeType(at)
    graph.atoms_to_nodes_[at] = temp
    return temp
end

"""
    deleteEdge(graph::MolecularGraph, e::Edge)
Deletes `e` from all adjacency lists of `Node`s in `graph` and from `graph` itself.
"""
deleteEdge(graph::MolecularGraph, e::Edge) = begin
    source = e.source_
    target = e.target_
    delete!(graph.bonds_to_edges_,e.bond_)
    filter!(x -> !(x == e),source.adjacent_edges)
    filter!(x -> !(x == e),target.adjacent_edges)
end

"""
    deleteNode(graph::MolecularGraph,node::AbstractNode)
Deletes `node` from `graph` and deletes all edges adjacent to `node`.
"""
deleteNode(graph::MolecularGraph,node::AbstractNode) = begin
    for edge in node.adjacent_edges_
        deleteEdge(graph, edge)
    end
    delete!(graph.atoms_to_bodes_, node.atom_)
end

"""
    newEdge(graph::MolecularGraph, bond::Bond)
Creates an edge connecting the two `Node`s that represents the atoms that `bond` connects.
Has no effect if the respective nodes do not exist.
"""
newEdge(graph::MolecularGraph, bond::Union{Bond, DFR}) = begin
    if haskey(graph.bonds_to_edges_, bond)
        return false
    end

    if !haskey(graph.atoms_to_nodes_, getSource(bond)) || !haskey(graph.atoms_to_nodes_, getTarget(bond))
        return false
    end

    temp = Edge(graph.atoms_to_nodes_[getSource(bond)], graph.atoms_to_nodes_[getTarget(bond)], bond)
    graph.bonds_to_edges_[bond] = temp
    push!(graph.atoms_to_nodes_[getSource(bond)].adjacent_edges_, temp)
    push!(graph.atoms_to_nodes_[getTarget(bond)].adjacent_edges_, temp)
end

"""
    MolecularGraph(root::T, NodeType::Type{typ} = Node) where typ<:AbstractNode where T<:Union{System, Chain, Residue}
Constructs a `MolecularGraph` mirroring the structure of the tree rooted in `root`.
The nodes will be of type `typ`.
"""
MolecularGraph(root::Union{KernelInterface, DataFrameRow}, NodeType::Type{typ} = Node) where typ<:AbstractNode = begin
    graph = MolecularGraph(NodeType)
    for at in eachAtom(root)
        newNode(graph, at)
    end

    for bond in eachBond(root)
        newEdge(graph, bond)
    end
    return graph
end

"""
    MolecularGraph(dfs::DataFrameSystem) = MolecularGraph(dfs.systems[1,:system],
        NodeType::Type{typ} = Node) where typ<:AbstractNode where T<:Union{System, Chain, Residue}
Gets the first `System` in `dfs` and constructs a Graph from it.
"""

"""
    breadthFirstSearch(graph::MolecularGraph, root::Node, NodeType::Type{typ} = Node) where typ<:AbstractNode
Generic BFS algorithm, that returns the BFS of `graph` rooted in `root`. The resulting
graph will have nodes of type `typ`.
"""
breadthFirstSearch(graph::MolecularGraph, root::Node, NodeType::Type{typ} = Node) where typ<:AbstractNode = begin
    bfs = MolecularGraph(NodeType)
    q::Queue{NodeType} = Queue{NodeType}()
    visited_nodes::Set{NodeType} = Set{NodeType}()
    cur::NodeType = root
    newNode(bfs,root.atom_, NodeType)
    push!(visited_nodes, root)
    cur = root
    while true
        for edge in collectPartnerEdges(cur)
            neighbour::Node = edge.source_ == cur ? edge.target_ : edge.source_
            if !(neighbour in visited_nodes)
                push!(visited_nodes, neighbour)
                enqueue!(q, neighbour)

                newNode(bfs, neighbour.atom_, BackpointingNode)
                newEdge(bfs, edge.bond_)
            end
        end
        isempty(q) && break
        cur = dequeue!(q)
    end
    return bfs
end


text:
@testset "5.6 Inverse cosecant" begin
    include("5.6.1 u (a+b arccsc(c x))^n.jl")
    include("5.6.2 Inverse cosecant functions.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.6 Inverse cosecant/.jl


text:
using Distributions
import ClusteringGA: fitness

function fruspini()
    ruspini = [
        4	53;
        5	63;
        10	59;
        9	77;
        13	49;
        13	69;
        12	88;
        15	75;
        18	61;
        19	65;
        22	74;
        27	72;
        28	76;
        24	58;
        27	55;
        28	60;
        30	52;
        31	60;
        32	61;
        36	72;
        28	147;
        32	149;
        35	153;
        33	154;
        38	151;
        41	150;
        38	145;
        38	143;
        32	143;
        34	141;
        44	156;
        44	149;
        44	143;
        46	142;
        47	149;
        49	152;
        50	142;
        53	144;
        52	152;
        55	155;
        54	124;
        60	136;
        63	139;
        86	132;
        85	115;
        85	96;
        78	94;
        74	96;
        97	122;
        98	116;
        98	124;
        99	119;
        99	128;
        101	115;
        108	111;
        110	111;
        108	116;
        111	126;
        115	117;
        117	115;
        70	4;
        77	12;
        83	21;
        61	15;
        69	15;
        78	16;
        66	18;
        58	13;
        64	20;
        69	21;
        66	23;
        61	25;
        76	27;
        72	31;
        64	30
    ]
    return [ruspini[i, :] for i = 1:lastindex(ruspini, 1)]
end

function rand200()
    nd1 = MvNormal([0, 10], [1.7, 1.7])
    nd2 = MvNormal([20, 12], [0.7, 0.7])
    nd3 = MvNormal([10, 20], [1.0, 1.0])
    s = hcat(rand(nd1, 120), rand(nd2, 60), rand(nd3, 20))
    t = s'
    tobj = [t[i, :] for i = 1:lastindex(t, 1)]
    return tobj
end

@testset "Clustering Genetic Algorithms" begin
    d, r = cga(fruspini())
    @test counts(r) == [20, 23, 17, 15]
    @test 1.73 < fitness(assignments(r), d) < 1.74
    d, r = cga(rand200())
    @test counts(r) == [120, 60, 20]
    @test 1.75 < fitness(assignments(r), d) < 1.85
end


text:
# let julia handle python internally
ENV["PYTHON"] = ""

# add
package_names = [
    "JLD",
    "GridInterpolations",
    "PyCall",
    "PyPlot"
]
for name in package_names
    Pkg.add(name)
end

# clone
urls = [
    "https://github.com/sisl/AutomotiveDrivingModels.jl.git",
    "https://github.com/sisl/AutoViz.jl.git",
    "https://github.com/sisl/BayesNets.jl.git",
    "https://github.com/sisl/NGSIM.jl.git",
    "https://github.com/wulfebw/AutoRisk.jl.git",
]

packages = keys(Pkg.installed())
for url in urls
    try
        id1 = search(url, "https://github.com/")[end]
        offset = search(url[(id1[end]+1):end], "/")[end]
        package = url[(id1+offset+1): (search(url,".jl.git")[1]-1)]
        if !in(package, packages)
          Pkg.clone(url)
        else
          println("$(package) already exists. Not cloning.")
        end
    catch e
        println("Exception when cloning $(url): $(e)")
    end
end

Pkg.build("AutomotiveDrivingModels")
Pkg.build("AutoViz")
Pkg.build("BayesNets")


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# run configuration xml file
function runConfig(configFilename::String; doPrint::Bool = false)
	sim = initSim(configFilename; allowResim = true, allowWriteOutput = true, doPrint = doPrint)
	if !isempty(sim.reps)
		simulateReps!(sim)
	else
		openOutputFiles!(sim)
		simulate!(sim)
		closeOutputFiles!(sim)
		doPrint && printSimStats(sim)
	end
	if sim.writeOutput
		writeOutputFiles(sim)
		sim.writeOutput = sim.backup.writeOutput = false # need to switch off in order to re-run the simulation
	end
	return sim
end

# initialise simulation from input files
function initSim(configFilename::String;
	allowResim::Bool = false, createBackup::Bool = true, allowWriteOutput::Bool = false, doPrint::Bool = false)

	# read sim config xml file
	configFilename = configFilename |> interpolateString |> abspath
	global configFileDir = dirname(configFilename)
	rootElt = xmlFileRoot(configFilename)
	@assert(xName(rootElt) == "simConfig", string("xml root has incorrect name: ", xName(rootElt)))

	# for progress messages:
	t = 0.0
	printInitMsg(msg) = doPrint && (t = time(); print(msg))
	printInitTime() = doPrint && println(": ", round(time() - t, digits = 2), " seconds")

	##################
	# sim config

	printInitMsg("reading config file data")

	sim = Simulation()
	sim.configRootElt = rootElt

	# input
	sim.inputPath = joinPathIfNotAbs(configFileDir, eltContentInterpVal(rootElt, "inputPath")) # input path can be absolute, or relative to configFilename
	simFilesElt = findElt(rootElt, "simFiles")
	inputFiles = childrenNodeNames(simFilesElt)
	sim.inputFiles = Dict{String,File}()
	for inputFile in inputFiles
		file = File()
		file.path = joinPathIfNotAbs(sim.inputPath, eltContentInterpVal(simFilesElt, inputFile))
		file.name = basename(file.path)
		if inputFile != "rNetTravels" # do not need checksum of rNetTravels file
			file.checksum = fileChecksum(file.path)
		end
		sim.inputFiles[inputFile] = file
	end

	# output
	sim.writeOutput = allowWriteOutput && eltContentVal(rootElt, "writeOutput")
	sim.outputPath = joinPathIfNotAbs(configFileDir, eltContentInterpVal(rootElt, "outputPath")) # output path can be absolute, or relative to configFilename
	isdir(sim.outputPath) || mkdir(sim.outputPath)
	outputFilesElt = findElt(rootElt, "outputFiles")
	outputFiles = childrenNodeNames(outputFilesElt)
	sim.outputFiles = Dict{String,File}()
	for outputFile in outputFiles
		file = File()
		file.path = joinPathIfNotAbs(sim.outputPath, eltContentInterpVal(outputFilesElt, outputFile))
		file.name = basename(file.path)
		sim.outputFiles[outputFile] = file
	end

	# events file filter
	if containsElt(rootElt, "eventsFileFilter")
		eventsFiltered = eltContentVal(rootElt, "eventsFileFilter")
		@assert(isa(eventsFiltered, Vector{EventForm}))
		for eventForm in setdiff(instances(EventForm), eventsFiltered)
			sim.eventsFile.eventFilter[eventForm] = false
		end
	end

	printInitTime()

	##################
	# read simulation input files

	printInitMsg("reading input data")

	simFilePath(name::String) = sim.inputFiles[name].path
	hasInputFile(name::String) = haskey(sim.inputFiles, name)

	# read sim data
	sim.ambulances = readAmbsFile(simFilePath("ambulances"))
	sim.hospitals = readHospitalsFile(simFilePath("hospitals"))
	sim.stations = readStationsFile(simFilePath("stations"))

	# calls
	callGenConfig = nothing
	@assert(hasInputFile("calls") + hasInputFile("callGenConfig") == 1, "Need exactly one of these input files: calls, callGenConfig.")
	if hasInputFile("calls")
		(sim.calls, sim.startTime) = readCallsFile(simFilePath("calls"))
	elseif hasInputFile("callGenConfig")
		callGenConfig = readGenConfig(simFilePath("callGenConfig"))
		sim.startTime = callGenConfig.startTime
	else error()
	end
	sim.time = sim.startTime

	sim.numAmbs = length(sim.ambulances)
	sim.numCalls = length(sim.calls)
	sim.numHospitals = length(sim.hospitals)
	sim.numStations = length(sim.stations)

	# read network data
	sim.net = Network()
	net = sim.net # shorthand
	fGraph = net.fGraph # shorthand
	fGraph.nodes = readNodesFile(simFilePath("nodes"))
	(fGraph.arcs, arcTravelTimes) = readArcsFile(simFilePath("arcs"))

	# read rNetTravels from file, if saved
	rNetTravelsLoaded = NetTravel[]
	rNetTravelsFilename = ""
	if hasInputFile("rNetTravels")
		rNetTravelsFilename = simFilePath("rNetTravels")
		if isfile(rNetTravelsFilename)
			try
				rNetTravelsLoaded = readRNetTravelsFile(rNetTravelsFilename)
			catch
				@warn("failed to read rNetTravels file")
			end
		elseif !isdir(dirname(rNetTravelsFilename)) || basename(rNetTravelsFilename) == ""
			# rNetTravelsFilename is invalid
			rNetTravelsFilename = ""
		else
			# save net.rNetTravels to file once calculated
		end
	end

	# read misc
	sim.map = readMapFile(simFilePath("map"))
	map = sim.map # shorthand
	(sim.targetResponseDurations, sim.responseTravelPriorities) = readPrioritiesFile(simFilePath("priorities"))
	sim.travel = readTravelFile(simFilePath("travel"))
	# "demand" file can be slow to read, will not read here but read elsewhere when needed
	if hasInputFile("demandCoverage")
		sim.demandCoverage = readDemandCoverageFile(simFilePath("demandCoverage"))
	end
	if hasInputFile("mobilisationDelay")
		sim.mobilisationDelay = readMobilisationDelayFile(simFilePath("mobilisationDelay"))
	end

	printInitTime()

	##################
	# grid

	printInitMsg("placing nodes in grid")

	# hard-coded grid size
	# grid rects will be roughly square, with one node per square on average
	n = length(fGraph.nodes)
	xDist = map.xRange * map.xScale
	yDist = map.yRange * map.yScale
	nx = Int(ceil(sqrt(n * xDist / yDist)))
	ny = Int(ceil(sqrt(n * yDist / xDist)))

	sim.grid = Grid(map, nx, ny)
	grid = sim.grid # shorthand
	gridPlaceNodes!(map, grid, fGraph.nodes)
	printInitTime()

	doPrint && println("nodes: ", length(fGraph.nodes), ", grid size: ", nx, " x ", ny)

	##################
	# generate calls

	callSets = [Call[]] # init
	if hasInputFile("callGenConfig")
		printInitMsg("generating calls")

		numReps = sim.numReps = callGenConfig.numCallsFiles # shorthand
		@assert(numReps >= 1)
		callSets = [makeCalls(callGenConfig) for i = 1:numReps]

		if numReps > 1
			setSimReps!(sim, callSets)
		end

		sim.calls = callSets[1] # set sim.calls for first replication
		sim.numCalls = length(sim.calls)

		printInitTime()

		(numReps > 1) && doPrint && println("number of call sets: $numReps")
	end

	##################
	# network

	printInitMsg("initialising fGraph")
	initGraph!(fGraph)
	printInitTime()

	if any(arc -> isnan(arc.distance), fGraph.arcs)
		printInitMsg("calculating arc distances for arcs with NaN distance")
		setArcDistances!(fGraph, map)
		printInitTime()
	end

	printInitMsg("checking fGraph")
	checkGraph(fGraph, map)
	printInitTime()

	printInitMsg("initialising fNetTravels")
	initFNetTravels!(net, arcTravelTimes)
	printInitTime()

	printInitMsg("creating rGraph from fGraph")
	createRGraphFromFGraph!(net)
	printInitTime()
	doPrint && println("fNodes: ", length(net.fGraph.nodes), ", rNodes: ", length(net.rGraph.nodes))

	printInitMsg("checking rGraph")
	checkGraph(net.rGraph, map)
	printInitTime()

	if rNetTravelsLoaded != []
		doPrint && println("using data from rNetTravels file")
		try
			printInitMsg("creating rNetTravels from fNetTravels")
			createRNetTravelsFromFNetTravels!(net; rNetTravelsLoaded = rNetTravelsLoaded)
			printInitTime()
		catch
			doPrint && println()
			@warn("failed to use data from rNetTravels file")
			rNetTravelsLoaded = []
		end
	end
	if rNetTravelsLoaded == []
		printInitMsg("creating rNetTravels from fNetTravels, and shortest paths")
		createRNetTravelsFromFNetTravels!(net)
		printInitTime()
		if rNetTravelsFilename != ""
			printInitMsg("saving rNetTravels to file")
			writeRNetTravelsFile(rNetTravelsFilename, net.rNetTravels)
			printInitTime()
		end
	end

	##################
	# travel

	printInitMsg("initialising travel")

	travel = sim.travel # shorthand
	@assert(travel.setsStartTimes[1] <= sim.startTime)
	@assert(length(net.fNetTravels) == travel.numModes)
	for travelMode in travel.modes
		travelMode.fNetTravel = net.fNetTravels[travelMode.index]
		travelMode.rNetTravel = net.rNetTravels[travelMode.index]
	end

	printInitTime()

	##################
	# sim - ambulances, calls, hospitals, stations...

	printInitMsg("adding ambulances, calls, etc")

	# for each call, hospital, and station, find nearest node
	for calls in (sim.calls, callSets[2:end]...) # note that sim.calls should be callSets[1] (if callSets is not empty)
		for c in calls
			if c.nearestNodeIndex == nullIndex
				(c.nearestNodeIndex, c.nearestNodeDist) = findNearestNode(map, grid, fGraph.nodes, c.location)
			end
		end
	end
	for h in sim.hospitals
		(h.nearestNodeIndex, h.nearestNodeDist) = findNearestNode(map, grid, fGraph.nodes, h.location)
	end
	for s in sim.stations
		(s.nearestNodeIndex, s.nearestNodeDist) = findNearestNode(map, grid, fGraph.nodes, s.location)
	end

	# create event list
	# try to add events to eventList in reverse time order, to reduce sorting required
	sim.eventList = Vector{Event}()

	# add first call to event list
	addEvent!(sim.eventList, sim.calls[1])

	# create ambulance wake up events
	for a in sim.ambulances
		initAmbulance!(sim, a)
		# currently, this sets ambulances to wake up at start of sim, since wake up and sleep events are not in ambulances file yet
	end

	# init station stats
	for station in sim.stations
		station.numIdleAmbsTotalDuration = OffsetVector(zeros(Float,sim.numAmbs+1), 0:sim.numAmbs)
		station.currentNumIdleAmbsSetTime = sim.startTime
	end

	printInitTime()

	printInitMsg("storing times between fNodes and common locations")

	# for each station, find time to each node in fGraph for each travel mode, and vice versa (node to station)
	# for each node in fGraph and each travel mode, find nearest hospital
	# requires deterministic and static travel times

	commonFNodes = sort(unique(vcat([h.nearestNodeIndex for h in sim.hospitals], [s.nearestNodeIndex for s in sim.stations])))
	setCommonFNodes!(net, commonFNodes)

	# find the nearest hospital to travel to from each node in fGraph
	numFNodes = length(fGraph.nodes) # shorthand
	for fNetTravel in net.fNetTravels
		fNetTravel.fNodeNearestHospitalIndex = Vector{Int}(undef, numFNodes)
		travelModeIndex = fNetTravel.modeIndex # shorthand
		travelMode = travel.modes[travelModeIndex] # shorthand
		for node in fGraph.nodes
			# find nearest hospital to node
			minTime = Inf
			nearestHospitalIndex = nullIndex
			for hospital in sim.hospitals
				travelTime = shortestPathTravelTime(net, travelModeIndex, node.index, hospital.nearestNodeIndex)
				travelTime += offRoadTravelTime(travelMode, hospital.nearestNodeDist)
				if travelTime < minTime
					minTime = travelTime
					nearestHospitalIndex = hospital.index
				end
			end
			fNetTravel.fNodeNearestHospitalIndex[node.index] = nearestHospitalIndex
		end
	end

	printInitTime()

	##################
	# decision logic

	decisionElt = findElt(rootElt, "decision")
	sim.addCallToQueue! = eltContentVal(decisionElt, "callQueueing")
	sim.findAmbToDispatch! = eltContentVal(decisionElt, "dispatch")
	hasInputFile("redispatch") && (sim.redispatch = readRedispatchFile(simFilePath("redispatch")))

	# move up
	mud = sim.moveUpData # shorthand
	moveUpElt = findElt(decisionElt, "moveUp")
	moveUpModuleName = eltContent(moveUpElt, "module")

	if moveUpModuleName == "none"
		mud.useMoveUp = false
		mud.moveUpModule = nullMoveUpModule
		doPrint && println("not using move up")
	else
		printInitMsg("initialising move up")

		mud.useMoveUp = true

		if moveUpModuleName == "comp_table"
			mud.moveUpModule = compTableModule
			compTableElt = findElt(moveUpElt, "compTable")
			compTableFilename = joinPathIfNotAbs(sim.inputPath, eltContentInterpVal(compTableElt, "filename"))
			initCompTable!(sim, compTableFilename)

		elseif moveUpModuleName == "ddsm"
			mud.moveUpModule = ddsmModule
			ddsmElt = findElt(moveUpElt, "ddsm")
			initDdsm!(sim;
				coverFractionTargetT1 = eltContentVal(ddsmElt, "coverFractionTargetT1"),
				travelTimeCost = eltContentVal(ddsmElt, "travelTimeCost"),
				slackWeight = eltContentVal(ddsmElt, "slackWeight"),
				coverTimeDemandPriorities = eltContentVal(ddsmElt, "coverTimeDemandPriorities"),
				options = convert(Dict{Symbol,Any}, Dict(eltContentVal(ddsmElt, "options"))))

		elseif moveUpModuleName == "dmexclp"
			mud.moveUpModule = dmexclpModule
			dmexclpElt = findElt(moveUpElt, "dmexclp")
			initDmexclp!(sim; busyFraction = eltContentVal(dmexclpElt, "busyFraction"))

		elseif moveUpModuleName == "multi_comp_table"
			mud.moveUpModule = multiCompTableModule
			multiCompTableElt = findElt(moveUpElt, "multiCompTable")
			multiCompTableFilename = joinPathIfNotAbs(sim.inputPath, eltContentInterpVal(multiCompTableElt, "filename"))
			initMultiCompTable!(sim, multiCompTableFilename)

		elseif moveUpModuleName == "priority_list"
			mud.moveUpModule = priorityListModule
			priorityListElt = findElt(moveUpElt, "priorityList")
			priorityListFilename = joinPathIfNotAbs(sim.inputPath, eltContentInterpVal(priorityListElt, "filename"))
			initPriorityList!(sim, priorityListFilename)

		elseif moveUpModuleName == "zhang_ip"
			mud.moveUpModule = zhangIpModule
			zhangIpElt = findElt(moveUpElt, "zhangIp")
			zhangIpParamsFilename = joinPathIfNotAbs(sim.inputPath, eltContentInterpVal(zhangIpElt, "paramsFilename"))
			initZhangIp!(sim;
				paramsFilename = zhangIpParamsFilename)

		elseif moveUpModuleName == "temp0"
			mud.moveUpModule = temp0Module
			temp0Elt = findElt(moveUpElt, "temp0")
			initTemp0!(sim;
				busyFraction = eltContentVal(temp0Elt, "busyFraction"),
				travelTimeCost = eltContentVal(temp0Elt, "travelTimeCost"),
				maxIdleAmbTravelTime = eltContentVal(temp0Elt, "maxIdleAmbTravelTime"),
				maxNumNearestStations = eltContentVal(temp0Elt, "maxNumNearestStations"))

		elseif moveUpModuleName == "temp1"
			mud.moveUpModule = temp1Module
			temp1Elt = findElt(moveUpElt, "temp1")
			initTemp1!(sim;
				busyFraction = eltContentVal(temp1Elt, "busyFraction"),
				travelTimeCost = eltContentVal(temp1Elt, "travelTimeCost"),
				maxIdleAmbTravelTime = eltContentVal(temp1Elt, "maxIdleAmbTravelTime"),
				maxNumNearestStations = eltContentVal(temp1Elt, "maxNumNearestStations"))

		elseif moveUpModuleName == "temp2"
			mud.moveUpModule = temp2Module
			temp2Elt = findElt(moveUpElt, "temp2")
			initTemp2!(sim;
				busyFraction = eltContentVal(temp2Elt, "busyFraction"),
				travelTimeCost = eltContentVal(temp2Elt, "travelTimeCost"),
				maxIdleAmbTravelTime = eltContentVal(temp2Elt, "maxIdleAmbTravelTime"),
				maxNumNearestStations = eltContentVal(temp2Elt, "maxNumNearestStations"))
		else
			error("invalid move up module name given: ", moveUpModuleName)
		end

		printInitTime()

		doPrint && println("using move up module: ", moveUpModuleName)
	end

	##################
	# resimulation

	if allowResim
		sim.resim.use = eltContentVal(rootElt, "resim")
		if sim.resim.use
			printInitMsg("")
			initResim!(sim, doPrint = doPrint)
			doPrint && print("initialised resimulation")
			printInitTime()
		end
	end

	##################
	# statistics

	if hasInputFile("statsControl")
		printInitMsg("initialising statistics")
		stats = sim.stats # shorthand
		stats.doCapture = true
		dict = readStatsControlFile(simFilePath("statsControl"))
		for fname in (:periodDurationsIter, :warmUpDuration, :recordResponseDurationHist, :responseDurationHistBinWidth)
			setfield!(stats, fname, dict[string(fname)])
		end
		stats.nextCaptureTime = sim.startTime + (stats.warmUpDuration > 0 ? stats.warmUpDuration : first(stats.periodDurationsIter))
		printInitTime()
	end

	##################
	# misc

	sim.initialised = true # at this point, the simulation could be run

	# set sim replications
	if sim.numReps >= 1
		makeRepsRunnable!(sim) # need to do this last in sim init (but can do before backing up), as it copies fields from sim
	end

	##################
	# sim backup

	if createBackup
		printInitMsg("creating sim backup")
		backup!(sim) # for restarting sim
		printInitTime()
	end

	return sim
end


text:
# mbe_demo_pitl.jl
#
# This script runs the control law on the target processor.

# Include our custom models.
include("models.jl")

# We work inside a module because it's more flexible than the main workspace.
module MBEDemoPITL

#############
# Resources #
#############

# Use/import the modules we'll need.
using Eratosthenes # The simulation engine (setup, simulate, mc)
using ..MBEModels  # The MBEModels module included above (and "above" this module)
import YAML        # For loading the YAML file directly
import HDF5        # For loading the log file
using Plots        # For plotting our results

# Create a directory for the log file.
if !isdir("out"); mkdir("out"); end

########################
# Software-in-the-Loop #
########################

# Specify the name of the YAML file that defines our scenario.
yaml_file = joinpath(@__DIR__, "underactuated.yaml")

# Load the YAML file and change the model name to use the PITL version.
# Also, change the log name.
yaml = YAML.load_file(yaml_file)
yaml["vehicles"][1]["computers"][1]["software"][1]["model"] = "ReducedEffortUnderactuatedControllerPITL"
yaml["sim"]["log_file"] = "out/underactuated-pitl.h5"

# The parameters are now stored under constants/parameters, so move them. Also,
# set up the UDP parameters here.
yaml["vehicles"][1]["computers"][1]["software"][1]["constants"] = Dict(
    "parameters"=>yaml["vehicles"][1]["computers"][1]["software"][1]["constants"],
    "conn"=>Dict(
        "target_ip"=>ip"192.168.1.195",
        "send_port"=>2000,
        "receive_port"=>2001
    )
);

# Set up the senario and run it.
scenario = setup(yaml)
scenario = simulate(scenario)

#########
# Plots #
#########

# Read some data from the logs.
t, q_BI, ω_BI_B = HDF5.h5open(scenario.sim.log_file, "r") do logs
    (dropdims(read(logs, "/smallsat/body/state/time"), dims=1),
     read(logs, "/smallsat/body/state/data/q_BI"),
     read(logs, "/smallsat/body/state/data/ω_BI_B"))
end

# Choose a friendly plotting package.
plotlyjs()

# Define each plot that we'll need.
display(plot(t, transpose(q_BI),
             label  = ["q1" "q2" "q3" "q4"],
             xlabel = "Time (s)",
             ylabel = "Attitude Quaternion"))
display(plot(t, 180/π * transpose(ω_BI_B),
             label  = ["ω1" "ω2" "ω3"],
             xlabel = "Time (s)",
             ylabel = "Rotation Rate (deg/s)"))

# Compare to the baseline.
(q_BI_baseline, ω_BI_B_baseline) = HDF5.h5open("out/underactuated.h5", "r") do logs
    (read(logs, "/smallsat/body/state/data/q_BI"),
     read(logs, "/smallsat/body/state/data/ω_BI_B"))
end
display(plot(t, transpose(q_BI) - transpose(q_BI_baseline),
    label  = ["q1" "q2" "q3" "q4"],
    xlabel = "Time (s)",
    ylabel = "Quaternion Differences",
    title  = "Differences in Attitude Quaternions"))
display(plot(t, 180/π * transpose((ω_BI_B) - transpose(ω_BI_B_baseline)),
    label  = ["ω1" "ω2" "ω3"],
    xlabel = "Time (s)",
    ylabel = "Rate (deg/s)",
    title  = "Differences in Rotation Rate"))

end # module


text:
using Glutils
using ModernGL
using GLFW
using Images

window = GLFW.CreateWindow(800, 600, "Textures")
GLFW.MakeContextCurrent(window)

vao = vertexarray()
vbo = arraybuffer(Cfloat[
#    Positions    Texture Coords
     0.5,  0.5,     1.0, 1.0,
     0.5, -0.5,     1.0, 0.0,
    -0.5, -0.5,     0.0, 0.0,
    -0.5,  0.5,     0.0, 1.0,
], g"static draw")

ibo = indexbuffer(Cuint[
    0, 1, 3,
    1, 2, 3
], g"static draw")

bindbuffer(:Array, vbo)
vertexattrib(Cfloat, 4, @map {
    0 := (size = 2, offset = 0),
    1 := (size = 2, offset = 2),
})

glActiveTexture(g"texture0")

image   = Images.load(assets("lenna.png")) |> rotr90
texture = texture2D(image, size(image)..., @map {
    min_filter = g"linear",
    mag_filter = g"linear",
    wrap_s = g"clamp",
    wrap_t = g"clamp"
})

shader = sprogram(assets("texture.glsl"), file = true)
@uniforms u_Texture = shader
glClearColor(0.0, 0.0, 0.0, 1.0)

while true
    # Break the loop when the window closes
    GLFW.WindowShouldClose(window) && break
    # Wait for events
    GLFW.PollEvents()
    # Clear the screen
    glClear(g"color buffer bit")
    glUseProgram(shader)
    uniform!(u_Texture, 0)
    glDrawElements(g"triangles", 6, g"unsigned int", C_NULL)
    # Swap the buffers
    GLFW.SwapBuffers(window)
end

GLFW.DestroyWindow(window)


text:
import Base.show

type DEAResult
  eff::Float64
  sx::Array
  sy::Array
  peers::Array

  function DEAResult(eff::Float64 = NaN,sx::Array = [],sy::Array = [],peers::Array = [])
    new(eff,sx,sy,peers)
  end
end

function geteff(res::DEAResult)
  return res.eff
end

function geteff(res::Array{DEAResult})
  eff = Array{Float64}(length(res))
  for k in eachindex(res)
    eff[k] = geteff(res[k])
  end
  return eff
end

function getpeers(res::DEAResult)
  return res.peers
end

function getpeers(res::Array{DEAResult})
  peers = Array{Float64}(length(res),length(res))
  for k in eachindex(res)
    peers[k,:] = getpeers(res[k])
  end
  return peers
end

function getsx(res::DEAResult)
  return res.sx
end

function getsx(res::Array{DEAResult})
  sx = Array{Float64}(length(res),length(getsx(res[1])))
  for k in eachindex(res)
    sx[k,:] = getsx(res[k])
  end
  return sx
end

function getsy(res::DEAResult)
  return res.sy
end

function getsy(res::Array{DEAResult})
  sy = Array{Float64}(length(res),length(getsy(res[1])))
  for k in eachindex(res)
    sy[k,:] = getsy(res[k])
  end
  return sy
end

function Base.show(io::IO,res::DEAResult)
  write(io,"eff = $(geteff(res)), sx = $(getsx(res)), sy = $(getsy(res))\n, peers = $(getpeers(res))")
end

max_stars_repo_name: kepiej/NEPA


text:
"""
    basis(N, n)

Generate a basis state (a.k.a. Fock or number state) ket \$|n⟩\$, in a Hilbert space of size `N`. Note that the size of the Hilbert space must be *at least* `n+1`. The function `fock` is an alias for `basis`.

Returns a sparse vector.

# Example
```jldoctest
julia> ψ = basis(3,2)
3-d Ket{SparseVector{Float64,Int64},1} with dimensions 3
1.00∠0°|2⟩
```
"""
function basis(n::Integer,dims::Dims)
    N = prod(dims)
    N>n || throw(ArgumentError("basis level $n is too large for a $N-d space"))
    # Julia is 1-indexed!
    return Ket(SparseVector(N,[n+1],[1.0]),dims)
end

basis(N::Integer, n::Integer) = basis(n::Integer,(N,))

"""
    coherent(N, α, analytic=false)

Generate a coherent state ket \$|α⟩\$, in a Hilbert space of size `N`. To create a coherent density operator, use the `Operator` function: `Operator(coherent(N,n))`.

Two methods can be used for generating a coherent state: via application of a displacment operator on a ground state (the default), or analytically, with the formula

```math
|α⟩ = e^{-\\frac{|α|^2}{2}} \\sum_{n=0}^{N-1} \\frac{α^n}{\\sqrt{n!}} |n⟩.
```

While the operator method will return a normalized ket, the analytic method will not. Both methods converge as `N` gets larger. The analytic method is also much faster, especially for large `N`.

Returns a dense vector.

# Example
```jldoctest
julia> coherent(6,0.4+1im)
6-d Ket{Array{Complex{Float64},1},1} with dimensions 6
0.60∠68°|1⟩ + 0.56∠0°|0⟩ + 0.46∠136°|2⟩ + 0.29∠-155°|3⟩ + 0.15∠-87°|4⟩ +…
```
"""
function coherent(N::Integer, α::Number, analytic::Bool=false)
    if analytic let a = exp(-abs2(α)/2) # due to julia issue #15276
        x = [a*α^n/sqrtfact(n) for n = 0:N-1]
        return Ket(x,(N,))
    end else
        return Ket(data(displacementop(N,α))[:,1],(N,)) # first column of D(α)
    end
end

"""
    thermal(N, n)

Generate a thermal state density matrix \$ρ_n\$ with particle number `n`, in a Hilbert space of size `N`. A thermal state \$ρ_n\$ is a probabilistic mixture of basis states such that the expectation value of the number operator \$\\hat{n}\$ is `n`. Note that this is true only if \$N≫n\$. The returned density matrix is always normalized.

Returns a sparse matrix.

# Example
```jldoctest
julia> N=5; n=0.2;

julia> ρ = thermal(N,n)
5×5 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 5
 0.833441  0.0       0.0        0.0         0.0
 0.0       0.138907  0.0        0.0         0.0
 0.0       0.0       0.0231511  0.0         0.0
 0.0       0.0       0.0        0.00385852  0.0
 0.0       0.0       0.0        0.0         0.000643087

julia> expect(numberop(N),ρ)
0.19935691318327978
```
"""
function thermal(N::Integer, n::Real)
    β = log(1.0/n+1.0)
    rowval = collect(1:N)
    colptr = Vector{Int}(undef,N+1); colptr[1:N] = rowval; colptr[end] = N+1
    nzval  = normalize!([exp(-β*k) for k = 0:N-1],1)
    return Operator(SparseMatrixCSC(N,N,colptr,rowval,nzval),(N,))
end

"""
    maxmixed(N)

Generate a maximally mixed density matrix. The maximally mixed state is a mixture of basis states with uniform probability.

Returns a sparse matrix.

# Example
```jldoctest
julia> maxmixed(4)
4×4 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 4
 0.25  0.0   0.0   0.0
 0.0   0.25  0.0   0.0
 0.0   0.0   0.25  0.0
 0.0   0.0   0.0   0.25
```
"""
function maxmixed(N::Integer)
    rowval = collect(1:N)
    colptr = Vector{Int}(undef,N+1); colptr[1:N] = rowval; colptr[end] = N+1
    nzval  = Vector{Float64}(undef,N); fill!(nzval, 1/N)
    return Operator(SparseMatrixCSC(N,N,colptr,rowval,nzval),(N,))
end

"""
    maxentangled(n,N=2)

Generate a maximally entangled state between `n` `N`-d systems:

```math
|\\phi⟩=\\sum_{j=0}^{N-1}\\frac{1}{\\sqrt{N}}|j⟩^{⊗n}.
```

Tracing out all but one of the entangled systems results in a maximally mixed state.

# Example
```jldoctest
julia> ψ = maxentangled(3,4)
64-d Ket{SparseVector{Float64,Int64},3} with dimensions 4⊗4⊗4
0.50∠0°|0,0,0⟩ + 0.50∠0°|1,1,1⟩ + 0.50∠0°|2,2,2⟩ + 0.50∠0°|3,3,3⟩

julia> ptrace(ψ,(1,3))
4×4 Operator{Array{Float64,2},1} with dimensions 4
 0.25  0.0   0.0   0.0
 0.0   0.25  0.0   0.0
 0.0   0.0   0.25  0.0
 0.0   0.0   0.0   0.25
```
"""
function maxentangled(n::Int,N::Int=2)
    c = div(N^n-1,N-1)
    nzind = [1+m*c for m=0:N-1]
    nzval = fill(1/sqrt(N),N)
    return Ket(SparseVector(N^n,nzind,nzval),ntuple(_->N,Val(n)))
end

"""
    ket(state,dims=2)

Generate a state ket from a tuple of basis levels and a tuple of corresponding space dimensions. Note that each space dimension must be larger than the level by *at least* 1. If only an integer is passed to `dims`, all basis levels will have the same dimension.

Returns a sparse vector.

# Example
```jldoctest
julia> ψ = ket((3,0,1),(5,2,3))
30-d Ket{SparseVector{Float64,Int64},3} with dimensions 5⊗2⊗3
1.00∠0°|3,0,1⟩
```

See also: [`@qb_str`](@ref), for generating qubit states via a bitstring.
"""
ket(state::Dims{N},dim::Int=2) where {N} = ket(state,ntuple(_->dim,Val(N)))
function ket(state::Dims{N},dims::Dims{N}) where {N}
    @inbounds for i = 1:N
        dims[i]>state[i] || throw(ArgumentError("basis level $(state[i]) is too large for a $(dims[i])-d space"))
    end
    n = tensored_sub2ind(dims,state) # 0-based indexing in tindex, returns 1 indexed index
    return Ket(SparseVector(prod(dims),[n],[1.0]),dims)
end

"""
    qb"bitstring"

Generate a qubit state ket from the given bitstring. This string macro can also be invoked via the macro notation `@qb_str "bitstring"`.

Returns a sparse vector.

# Example
```jldoctest
julia> Ψ⁻ = normalize!(qb"01" - qb"10")
4-d Ket{SparseVector{Float64,Int64},2} with dimensions 2⊗2
0.71∠0°|0,1⟩ + 0.71∠180°|1,0⟩
```

See also: [`ket`](@ref), for generating arbitrary states.
"""
macro qb_str(state::AbstractString)
    quote
        dims = $(ntuple(_->2,length(state)))
        Ket(SparseVector(prod(dims),[$(parseb2(state)+1)],[1.0]),dims)
    end
end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module MPI

using Compat

using Libdl

@static if Compat.Sys.iswindows()
    const depfile = "win_mpiconstants.jl"
end

@static if Compat.Sys.isunix()
    const depfile = joinpath(dirname(@__FILE__), "..", "deps", "src", "compile-time.jl")
    isfile(depfile) || error("MPI not properly installed. Please run Pkg.build(\"MPI\")")
end

include(depfile)

include("mpi-base.jl")
include("cman.jl")

const mpitype_dict = Dict{DataType, Cint}()
const mpitype_dict_inverse = Dict{Cint, DataType}()

@static if Compat.Sys.isunix()
    for (jname, fname) in _mpi_functions
        Core.eval(MPI, :(const $jname = ($(QuoteNode(fname)),libmpi)))
    end
end

function __init__()
    # Note: older versions of OpenMPI (e.g. the version on Travis) do not
    # define MPI_CHAR and MPI_*INT*_T for Fortran, so we don't use them (yet).
    for (T,mpiT) in (Char => MPI_INTEGER4,   # => MPI_WCHAR, (note: wchar_t is 16 bits in Windows)
                     Int8 => MPI_INTEGER1,   # => MPI_INT8_T,
                     UInt8 => MPI_INTEGER1,  # => MPI_UINT8_T,
                     Int16 => MPI_INTEGER2,  # => MPI_INT16_T,
                     UInt16 => MPI_INTEGER2, # => MPI_UINT16_T,
                     Int32 => MPI_INTEGER4,  # => MPI_INT32_T,
                     UInt32 => MPI_INTEGER4, # => MPI_UINT32_T,
                     Int64 => MPI_INTEGER8,  # => MPI_INT64_T,
                     UInt64 => MPI_INTEGER8, # => MPI_UINT64_T,
                     Float32 => MPI_REAL4,
                     Float64 => MPI_REAL8,
                     ComplexF32 => MPI_COMPLEX8,
                     ComplexF64 => MPI_COMPLEX16)
        mpitype_dict[T] = mpiT
        mpitype_dict_inverse[mpiT] = T
    end
end

end


text:
#####
##### CVMix mixing depth model
#####

Base.@kwdef struct LMDMixingDepth{T} <: AbstractParameters
     CSL :: T = 0.1   # Surface layer fraction
     CRi :: T = 0.3   # Critical bulk Richardson number
     CKE :: T = 4.32  # Unresolved turbulence parameter
    CKE₀ :: T = 1e-11 # Minimum unresolved turbulence kinetic energy
end

function update_mixing_depth!(m::Model{K, NL, <:LMDMixingDepth}) where {K, NL}
    m.state.h  = mixing_depth(m)
    return nothing
end

bulk_richardson_number(m::AbstractModel, i) = KPP.bulk_richardson_number(
    m.solution.U, m.solution.V, m.solution.T, m.solution.S,
    m.state.Qb, m.mixingdepth.CKE, m.mixingdepth.CKE₀, m.mixingdepth.CSL, m.constants.g,
    m.constants.α, m.constants.β, i)

"""
    mixing_depth(model)

Calculate the mixing depth 'h' for `model`.
"""
function mixing_depth(m)
    ih₁ = m.grid.N + 1 # start at top.
    Ri₁ = bulk_richardson_number(m, ih₁) # should be 0.

    # Descend through grid until Ri rises above critical value
    while ih₁ > 1 && Ri₁ < m.mixingdepth.CRi
        ih₁ -= 1 # descend
        Ri₁ = bulk_richardson_number(m, ih₁)
    end

    # Edge cases:
    # 1. Mixing depth is 0:
    if ih₁ == m.grid.N + 1
        z★ = m.grid.zf[ih₁]

    # 2. Mixing depth is whole domain because Ri is always less than CRi:
    elseif ih₁ == 1 && Ri₁ < m.mixingdepth.CRi
        z★ = m.grid.zf[ih₁]

    # 3. Ri is infinite somewhere inside the domain.
    elseif !isfinite(Ri₁)
        z★ = m.grid.zc[ih₁]

    # Main case: mixing depth is in the interior.
    else # Ri₁ > CRi
        ΔRi = bulk_richardson_number(m, ih₁+1) - Ri₁ # <0 linearly interpolate to find h.
        # x = x₀ + Δx * (y-y₀) / Δy
        z★ = m.grid.zf[ih₁] + Δf(m.grid, ih₁) * (m.mixingdepth.CRi - Ri₁) / ΔRi
    end

    -z★ < 0 && @warn "mixing depth $(-z★) is negative"

    return -z★ # "depth" is negative height.
end

#####
##### ROMS mixing depth model
#####

Base.@kwdef struct ROMSMixingDepth{T} <: AbstractParameters
     CSL :: T = 0.1  # Surface layer fraction
     CRi :: T = 0.3  # Critical bulk Richardson number
     CKE :: T = 5.07 # Minimum unresolved turbulence kinetic energy
     CEk :: T = 0.0  # Turbulent Ekman depth parameter
end

h_criterion(::ROMSMixingDepth, grid) = FaceField(grid)

function update_mixing_depth!(m::Model{K, NL, <:ROMSMixingDepth}) where {K, NL}
    mixing_depth_criterion!(m.state.h_crit, m)
    m.state.h = mixing_depth(m)
    return nothing
end

h_weight(h, CSL, zf, i) = @inbounds -zf[i] / (CSL*h - zf[i])
h_weight(m, i) = h_weight(m.state.h, m.mixingdepth.CSL, m.grid.zf, i)

function h_kernel(U, V, T, S, CRi, CEk, g, α, β, f, i)
    @inbounds ∂z(U, i)^2 + ∂z(V, i)^2 - ∂B∂z(T, S, g, α, β, i)/CRi - CEk*f^2
end

h_kernel(m, i) = h_kernel(m.solution.U, m.solution.V, m.solution.T, m.solution.S,
                            m.mixingdepth.CRi, m.mixingdepth.CEk,
                            m.constants.g, m.constants.α, m.constants.β, m.constants.f, i)

function unresolved_kinetic_energy(m, i)
    @inbounds unresolved_kinetic_energy(-m.grid.zf[i],
        ∂B∂z(m.solution.T, m.solution.S, m.constants.g, m.constants.α, m.constants.β, i),
        m.state.Qb, m.mixingdepth.CKE, 0, m.constants.g, m.constants.α, m.constants.β)
end

"Calculate the mixing depth criterion function by integrating from z=0 downwards."
function mixing_depth_criterion!(h_crit, m)
    @inbounds h_crit[m.grid.N+1] = 0

    for i = m.grid.N:-1:1
        @inbounds h_crit[i] = h_crit[i+1] + h_weight(m, i) * h_kernel(m, i) * Δc(m.grid, i)
    end

    for i in eachindex(h_crit)
        @inbounds h_crit[i] -= unresolved_kinetic_energy(m, i) / m.grid.zf[i]
    end

    return nothing
end

linear_interp(y★, x₀, y₀, Δx, Δy) = x₀ + Δx * (y★ - y₀) / Δy

function mixing_depth(m::Model{K, NL, <:ROMSMixingDepth}) where {K, NL}
    ih₁ = findprev(x -> x<=0, m.state.h_crit.data, m.grid.N)
    @inbounds begin
        if ih₁ === nothing # Mixing depth is entire grid
            z★ = m.grid.zf[1]
        elseif ih₁ == m.grid.N # Mixing depth at surface?
            z★ = ifelse(m.state.h_crit[ih₁]==0, m.grid.zf[m.grid.N], m.grid.zf[m.grid.N+1])
        else # linearly interpolate
            # x = x₀ + Δx * (y-y₀) / Δy
            z★ = linear_interp(0, m.grid.zf[ih₁], m.state.h_crit[ih₁], Δf(m.grid, ih₁),
                                m.state.h_crit[ih₁+1] - m.state.h_crit[ih₁])
        end
    end

    return -z★
end


text:
###########################
### Auxiliary Functions ###
###########################

"""
    consword(w::Integer)

Create a word with n consecutive set bits as least significant.
"""
consword(n::T) where {T <: Integer} = (T(1) << n) - T(1)

#= Counts the number of bits in the integer w =#
nbits(w::T) where {T <: Integer} = 8*sizeof(T) - leading_zeros(w)
nbits(w::BigInt)                 = (w == 0) ? Int64(0) : ndigits(w, base=2)

#= Counting the number of set bits in the integer w =#
nsetbits(w::T) where {T <: Integer} = count_ones(w)

#= Counting the maximum number of consecutive set bits in an integer w =#
function nconsbits(w::T) where {T <: Integer}
    n = 0
    while w != Miss
        w &= (w << T(1))
        n += 1
    end # while
    n
end # function

#################
### Functions ###
#################

"""
    shift(w::Integer, c)

Left shifts bits once and adds new bit as LSB (according to character c).
Removes the MSB if character is a Hit.
"""
function shift(w::T, c) where {T <: Integer}
    if c == Miss
        return w << T(1)
    else
        return ((w << T(1)) | T(1)) & consword(T(nbits(w)))
    end
end

"""
    shift_rowhit(w::Integer, c, n::Integer)

Left shifts bits once and adds new bit as LSB (according to character c).
Returns the word n if the bits in nw coincide with the ones in n.
"""
function shift_rowhit(w::T, c, n::T) where {T <: Integer}
    #= Shift bits according to charcter c and default to n if last few char are equal to n =#
    if c == Miss
        w <<= T(1)
    else
        w = (w << T(1)) | T(1)
    end
    return (w & n == n) ? n : w
end


text:
using BenchmarkTools

const SUITE = BenchmarkGroup()

module BenchmarkFinSets
  include("FinSets.jl")
end

module BenchmarkGraphs
  include("Graphs.jl")
end

SUITE["FinSets"] = BenchmarkFinSets.SUITE
SUITE["Graphs"] = BenchmarkGraphs.SUITE


text:
using Test

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

a = fill(rank,size)
@test MPI.Alltoall(a, 1, comm) == collect(0:(size-1))
MPI.Finalize()


text:
export randprime

const RPLIM = 10^4

# not using keyword here is 5-10 percent faster in a quick check
@inline function _randprime(a,b,lim)
    n = zero(b)
    i = 1
    while i < lim
        n = rand(a:b)
        isprime(n) && break
        i += 1
    end
    i == lim ? zero(b) : n
end

randprime(a,b; lim::Int=RPLIM) = _randprime(a,b,lim)
randprime(b; lim::Int=RPLIM) = _randprime(convert(typeof(b),2),b,lim)

function randprime(a,b,dims...;lim::Int=RPLIM)
    arr = Array{typeof(a)}(dims...)
    @inbounds for i in 1:length(arr) arr[i] = _randprime(a,b,lim) end
    arr
end


text:
import Base: *, A_mul_B!
import Base: A_mul_Bt, A_mul_Bt!, A_mul_Bc, A_mul_Bc!
import Base: At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!
import Base: At_mul_Bt, At_mul_Bt!, Ac_mul_Bc, Ac_mul_Bc!
import Base.LinAlg.BLAS: gemm!
import Base.LinAlg: axpy!

A_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','N',one(T),A,B,zero(T),C)
(*){T}(A::KnetMatrix{T},B::KnetMatrix{T})=A_mul_B!(similar(A,(size(A,1),size(B,2))),A,B)

A_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','T',one(T),A,B,zero(T),C)
A_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(similar(A,(size(A,1),size(B,1))),A,B)
A_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(C,A,B)
A_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt(A,B)

At_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','N',one(T),A,B,zero(T),C)
At_mul_B{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(similar(A,(size(A,2),size(B,2))),A,B)
Ac_mul_B!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(C,A,B)
Ac_mul_B{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B(A,B)

At_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','T',one(T),A,B,zero(T),C)
At_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(similar(A,(size(A,2),size(B,2))),A,B)
Ac_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(C,A,B)
Ac_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt(A,B)

function gemm!{T}(transA::Char, transB::Char, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T})
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    size2(x,i)=(if ndims(x)<=2; size(x,i); elseif i==1; div(length(x),size(x,ndims(x))); elseif i==2; size(x,ndims(x)); else 1; end)
    if transA == 'N'
        m=size2(A,1); k=size2(A,2)
    else
        m=size2(A,2); k=size2(A,1)
    end
    if transB == 'N'
        n=size2(B,2); k==size2(B,1) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    else
        n=size2(B,1); k==size2(B,2) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    (m == size2(C,1) && n == size(C,2)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    lda = size2(A,1); ldb = size2(B,1); ldc = size2(C,1)
    if T<:Float32
        @cuda(cublas, cublasSgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle, transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    elseif T<:Float64
        @cuda(cublas, cublasDgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle, transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

function axpy!{T}(n::Integer, alpha::Number, x::KnetArray{T}, incx::Integer, y::KnetArray{T}, incy::Integer)
    length(x) == length(y) || throw(DimensionMismatch("$(map(size,(x,y)))"))
    alpha = T[alpha]
    if T<:Float32
        @cuda(cublas, cublasSaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle, n, alpha, x, incx, y, incy)
    elseif T<:Float64
        @cuda(cublas, cublasDaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle, n, alpha, x, incx, y, incy)
    else
        error("$T not supported")
    end
    return y
end

axpy!{T}(alpha::Number, x::KnetArray{T}, y::KnetArray{T})=axpy!(length(x),alpha,x,1,y,1)

export axpy!


text:
using NCDatasets
using Distributed
using SharedArrays
using Formatting

s_file  = ARGS[1]
d_file  = ARGS[2]
varname = ARGS[3]

wgt_file = "wgt_gx1v6_to_gx3v7.nc"

ds_w = Dataset(wgt_file, "r")
ds_s = Dataset(s_file, "r")
ds_d = Dataset(d_file, "c")

missing_value = 1e20

NN_idx = convert(Array{Integer}, nomissing(ds_w["NN_idx"][:], 0))
NNN = size(NN_idx)[1]

d_Nx = ds_w.dim["d_Nx"]
d_Ny = ds_w.dim["d_Ny"]
d_N = d_Nx * d_Ny

defDim(ds_d, "Nx", d_Nx)
defDim(ds_d, "Ny", d_Ny)
defDim(ds_d, "time", Inf)

d_var = defVar(ds_d, varname, Float64, ("Nx", "Ny"))
d_var.attrib["_FillValue"] = missing_value

# Write grid information
for (varname, vardata, dims) in (
    ("lat", ds_w["d_lat"][:], ("Nx", "Ny")),
    ("lon", ds_w["d_lon"][:], ("Nx", "Ny")),
)
    v = defVar(ds_d, varname, Float64, dims)
    v.attrib["_FillValue"] = missing_value
    v[:] = vardata
end

println("NNN: ", NNN)

d_data = zeros(Float64, d_Nx, d_Ny)

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

d_data .= 0.0
s_data = reshape(nomissing(ds_s[varname][:, :], NaN), :)

for i = 1 : length(d_data)
    if NN_idx[1, i] == 0
        d_data[i] = missing_value
        continue
    end

    for j = 1:NNN
        d_data[i] += s_data[NN_idx[j, i]]
    end
    d_data[i] /= NNN
end

d_var[:, :] = d_data

##############
#=

for t = 1 : ds_s.dim["time"]
    global d_data
    print("\rProcessing time: ", t)
    d_data .= 0.0
    s_data = reshape(nomissing(ds_s[varname][:, :, t], NaN), :)

    for i = 1 : length(d_data)
        if NN_idx[1, i] == 0
            d_data[i] = missing_value
            continue
        end

        for j = 1:NNN
            d_data[i] += s_data[NN_idx[j, i]]
        end
        d_data[i] /= NNN
    end

    d_var[:, :, t] = d_data
end

=#

println("done.")

close(ds_w)
close(ds_s)
close(ds_d)

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays
using AbstractFFTs

A = sprand(256,256,256,0.1)

P = plan_ifft(A)

@benchmark P*A

A = rand(256,256,256)

P = plan_ifft(A)

@benchmark P*A

filter_hash = fink_filter_hash(1, 8, nx=64, pc=1, wd=1)
filt_3d = fink_filter_bank_3dizer(filter_hash, 1, nz=64)

function fink_filter_bank_3dizer_AKS(hash, cz; nz=256)
    # -------- set parameters
    nx = convert(Int64,hash["npix"])
    n2d = size(hash["filt_value"])[1]
    dz = nz/2-1
    J = size(hash["j_value"])[1]
    L = size(hash["theta_value"])[1]

    im_scale = convert(Int8,log2(nz))
    # -------- number of bins in radial direction (size scales)
    K = (im_scale-2)*cz

    # -------- allocate output array of zeros
    J_L_K     = zeros(Int32, n2d*K,3)
    psi_index = zeros(Int32, J+1, L, K)
    k_value   = zeros(Float64, K)
    logr      = zeros(Float64,nz)
    filt_temp = zeros(Float64,nx,nx,nz)
    F_z       = zeros(Float64,nx,nx,nz)
    #F_p       = zeros(Float64,nx,nx,nz)

    filtind = fill(CartesianIndex{3}[], K*n2d)
    filtval = fill(Float64[], K*n2d)

    for z = 1:convert(Int64,dz+1)
        sz = mod(z+dz,nz)-dz-1
        z2 = sz^2
        logr[z] = 0.5*log2(max(1,z2))
    end

    @inbounds for k_ind = 1:K
        k = k_ind/cz
        k_value[k_ind] = k  # store for later
        krad  = im_scale-k-1
        Δk    = abs.(logr.-krad)
        kmask = (Δk .<= 1/cz)
        k_count = count(kmask)
        k_vals = findall(kmask)

        # -------- radial part
        #I have somehow dropped a factor of sqrt(2) which has been reinserted as a 0.5... fix my brain
        @views F_z[:,:,kmask].= reshape(1/sqrt(2).*cos.(Δk[kmask] .* (cz*π/2)),1,1,k_count)
        @inbounds for index = 1:n2d
            p_ind = hash["filt_index"][index]
            p_filt = hash["filt_value"][index]
            #F_p[p_ind,:] .= p_filt
            f_ind    = k_ind + (index-1)*K
            @views filt_tmp = p_filt.*F_z[p_ind,kmask]
            #temp_ind = findall((F_p .> 1E-13) .& (F_z .> 1E-13))
            #@views filt_tmp = F_p[temp_ind].*F_z[temp_ind]
            ind = findall(filt_tmp .> 1E-13)
            @views filtind[f_ind] = map(x->CartesianIndex(p_ind[x[1]][1],p_ind[x[1]][2],k_vals[x[2]]),ind)
            #filtind[f_ind] = temp_ind[ind]
            @views filtval[f_ind] = filt_tmp[ind]
            @views J_L_K[f_ind,:] = [hash["J_L"][index,1],hash["J_L"][index,2],k_ind-1]
            psi_index[hash["J_L"][index,1]+1,hash["J_L"][index,2]+1,k_ind] = f_ind
        end
    end

    # -------- metadata dictionary
    info=Dict()
    info["2d_npix"]         = hash["npix"]
    info["2d_j_value"]      = hash["j_value"]
    info["2d_theta_value"]  = hash["theta_value"]
    info["2d_psi_index"]    = hash["psi_index"]
    info["2d_phi_index"]    = hash["phi_index"]
    info["2d_J_L"]          = hash["J_L"]
    info["2d_pc"]           = hash["pc"]
    info["2d_wd"]           = hash["wd"]
    info["2d_fs_center_r"]  = hash["fs_center_r"]
    info["2d_filt_index"]   = hash["filt_index"]
    info["2d_filt_value"]   = hash["filt_value"]

    info["nz"]              = nz
    info["cz"]              = cz
    info["J_L_K"]           = psi_index
    info["psi_index"]       = psi_index
    info["k_value"]         = k_value
    info["filt_index"]      = filtind
    info["filt_value"]      = filtval

    return info
end

function DHC_compute_3d_AKS(image::Array{Float64,3}, filter_hash; FFTthreads=2)
    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny, Nz)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    S20 = zeros(Float64, Nf, Nf)
    S12 = zeros(Float64, Nf, Nf)
    S2  = zeros(Float64, Nf, Nf)  # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nz, Nf)   # this must be zeroed!
    im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, Nz, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny*Nz)
    norm_im ./= sqrt(Nx*Ny*Nz*S0[2])

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    #f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    #f_val2   = filter_hash2["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny, Nz)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P   = plan_ifft(im_fd_0)   # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval
                im_fdf_0_1[ind,f] = abs(zval)
            end
            S1[f] = S1tot/(Nx*Ny*Nz)  # image power ###why do we need to normalize again?
            @views im_rd_0_1[:,:,:,f] .= abs2.(P*zarr)
            zarr[f_i] .= 0
        end
    end
    append!(out_coeff, S1[:])

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    im_rd_0_1 .= sqrt.(im_rd_0_1)

    ## 2nd Order
    #Amat = reshape(im_fdf_0_1, Nx*Ny, Nf)
    #S12  = Amat' * Amat
    #Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
    #S20  = Amat' * Amat

    #append!(out_coeff, S20)
    #append!(out_coeff, S12)

    ## Traditional second order
    @inbounds for f1 = 1:Nf
        @views thisim = fft(im_rd_0_1[:,:,:,f1])  # Could try rfft here
        # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
        # Loop over f2 and do second-order convolution
        @inbounds for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i
            # sum im^2 = sum(|fft|^2/npix)
            @views S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny*Nz)
        end
    end
    append!(out_coeff, S2)

    return out_coeff
end

test_img = rand(64,64,64)
@benchmark DHC_compute_3d_AKS(test_img,filt_3d)

Profile.clear()
@profile DHC_compute_3d_AKS(test_img,filt_3d)

Juno.profiler()

function DHC_compute_AKS(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS12::Bool=false, doS20::Bool=false, norm=true, iso=false, FFTthreads=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS12        - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS12 S12 = zeros(Float64, Nf, Nf) end  # Fourier correlation
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS12 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if doS12 im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nf) end   # this must be zeroed!
    if anyrd im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf) end

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im ./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P = plan_ifft(im_fd_0) end  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval        # filter*image in Fourier domain
                if doS12 im_fdf_0_1[ind,f] = abs(zval) end
            end
            S1[f] = S1tot/(Nx*Ny)  # image power
            if anyrd
                im_rd_0_1[:,:,f] .= abs2.(P*zarr) end
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
            # Loop over f2 and do second-order convolution
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                # sum im^2 = sum(|fft|^2/npix)
                S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny)
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Fourier domain 2nd order
    if doS12
        Amat = reshape(im_fdf_0_1, Nx*Nx, Nf)
        S12  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S12[:] : S12[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Nx, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

test_img2 = rand(256,256)
filter_hash2 = fink_filter_hash(1, 8, nx=256, pc=1, wd=1)

Profile.clear()
@profile DHC_compute_AKS(test_img2,filter_hash2)
Juno.profiler()

@benchmark DHC_compute_AKS(test_img2,filter_hash2)

function DHC_compute_AKS3(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS12::Bool=false, doS20::Bool=false, norm=true, iso=false, FFTthreads=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS12        - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS12 S12 = zeros(Float64, Nf, Nf) end  # Fourier correlation
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS12 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if doS12 im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nf) end   # this must be zeroed!
    if anyrd im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf) end

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im ./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P = plan_ifft(im_fd_0) end  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    @inbounds for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            @inbounds for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval        # filter*image in Fourier domain
                if doS12 im_fdf_0_1[ind,f] = abs(zval) end
            end
            S1[f] = S1tot/(Nx*Ny)  # image power
            if anyrd
                @views im_rd_0_1[:,:,f] .= abs2.(P*zarr) end
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        P1 = plan_fft(im_rd_0_1[:,:,1])

        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        @inbounds for f1 = 1:Nf
            thisim = P1*im_rd_0_1[:,:,f1]  # Could try rfft here
            # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
            # Loop over f2 and do second-order convolution
            @inbounds for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                # sum im^2 = sum(|fft|^2/npix)
                S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny)
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Fourier domain 2nd order
    if doS12
        Amat = reshape(im_fdf_0_1, Nx*Nx, Nf)
        S12  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S12[:] : S12[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Nx, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

Profile.clear()
@profile DHC_compute_AKS3(test_img2,filter_hash2)
Juno.profiler()

@benchmark DHC_compute_AKS3(test_img2,filter_hash2)

out1 = DHC_compute_AKS3(test_img2,filter_hash2)
out2 = DHC_compute_AKS(test_img2,filter_hash2)

out1 == out2

##Hm

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
println("=== Tests for OED.jl  ===")
include("testOEDTypeA.jl")
include("testUnconstrainedOED.jl")
include("testEqConstrOED.jl")
include("testIneqConstrOED.jl")

include("testSolveKKT.jl")
include("testIPQP.jl")

include("testOEDTypeB.jl")
include("testRotation2D.jl")
include("testLinearInter.jl")
include("testGetLinearInterMatrixTranspose.jl")

include("testRunOEDTypeA.jl")
include("testRunOEDTypeB.jl")
println("=== [OED.jl] : all tests passed ===")


text:
# test basic manifold product behaviour

using ApproxManifoldProducts
using KernelDensityEstimate
using Test
using TensorCast
using Manifolds

## 3D

@testset "Basic 3D multiply and cross dimension covariance test..." begin

##

N = 100

pts1 = [[0.05*randn(2);0.75*randn()] for i in 1:N]
pts2 = [[0.05*randn(2);0.75*randn()] for i in 1:N]

P1 = manikde!(TranslationGroup(3),pts1)
P2 = manikde!(TranslationGroup(3),pts2)

P12 = P1*P2

@test typeof(P12._u0) <: Vector{Float64}

pts_ = getPoints(P12)

pts = AMP._pointsToMatrixCoords(P12.manifold, pts_)

@test 0.8*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[3,:]) .< 2.0)

##

M = SpecialEuclidean(2)
u0 = ProductRepr(zeros(2),[1 0; 0 1.0])
ϵ = identity_element(M, u0)

pts1 = [exp(M, ϵ, hat(M, ϵ, [0.05*randn(2);0.75*randn()])) for i in 1:N]
pts2 = [exp(M, ϵ, hat(M, ϵ, [0.05*randn(2);0.75*randn()])) for i in 1:N]

P1 = manikde!(M, pts1)
P2 = manikde!(M, pts2)

P12 = P1*P2

pts_ = getPoints(P12)

pts = AMP._pointsToMatrixCoords(P12.manifold, pts_)

@test 0.8*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[3,:]) .< 2.0)

# using KernelDensityEstimatePlotting
# using Gadfly
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE([P1;P2;P12], c=["red";"blue";"magenta"],levels=1) |> PDF("/tmp/test.pdf")
# run(`evince /tmp/test.pdf`)

##

end

#


text:
Causes.notimplemented(
    f,
    args::Tuple,
    kwargs::Union{NamedTuple,Iterators.Pairs} = NamedTuple(),
) = Err(NotImplementedError(f, args, kwargs))

Causes.empty(container) = Err(TryExperimental.EmptyError(container))

struct EmptyError <: TryExperimental.EmptyError
    container::Any
end
# TODO: check if it is better to not "type-erase"

TryExperimental.EmptyError(container) = EmptyError(container)

Causes.closed(container) = Err(TryExperimental.ClosedError(container))

struct ClosedError <: TryExperimental.ClosedError
    container::Any
end
# TODO: check if it is better to not "type-erase"

TryExperimental.ClosedError(container) = ClosedError(container)


text:
function delta2Jl!(
    delta :: Matrix, AField :: Matrix, BField :: Matrix, Lap :: Matrix,
    Out :: Matrix, decay :: Real)

    d1, d2 = size(Out)
    for j in 1:d2
      for i in 1:d1
        Out[i,j] = delta[i,j] / (1 + AField[i,j]) - decay * BField[i,j] + Lap[i,j]
      end
    end
end