text:
using ChaosTools
using Test

@testset "DyCA Tests" begin
    println("Running dyca() tests...")

    @testset "Eigenvalues of Roessler system" begin
        # Initial conditions from `dyca()` paper
        eigenthresold = collect(0.10:0.01:0.995)
        ds = Systems.roessler(rand(3); a = 0.15, b = 0.2, c = 10);
        ts = trajectory(ds, 1000.0, dt = 0.05);
        for thresold in eigenthresold
            eigenvalues,_ = dyca(Matrix(ts),thresold)
            @test sum([round(i) for i in eigenvalues]) == 2.0
        end
    end

    @testset "Eigenvalues of Embeded Roessler system" begin
        # Initial conditions from `dyca()` paper
        eigenthresold = collect(0.90:0.01:0.995)
        ds = Systems.roessler(rand(3); a = 0.15, b = 0.2, c = 10);
        ts = trajectory(ds, 1000.0, dt = 0.05);
        cols = columns(ts);
        for thresold in eigenthresold
            for i in 1:3
                Embedded_system = Matrix(embed(cols[i], 25, 1));
                eigenvals,_ = dyca(Embedded_system,thresold)
                # eigenvalues satisfy `dyca` condtiion and have imaginary part = 0.0
                wanted_eigenvalues = eigenvals[vec(0.999 .< abs.(eigenvals) .< 1.0) .& vec(imag(eigenvals) .== 0)]
                # due to floating point errors more than 2 eigenvalues might satisy the above condition
                # might be better to use `length(psum) in [1,2,3]`
                @test length(wanted_eigenvalues) in 0:4
            end
        end
    end
end


text:
using Soss

μdist = @model begin
    a ~ Normal()
    b ~ Normal()
    return a/b
end

σdist = @model begin
    x ~ Normal()
    return x^2
end

m = @model begin
    μ ~ μdist()
    σ ~ σdist()
    x ~ Normal(μ,σ) |> iid(2)
    return x
end

rand(m(),5)

simulate(m(), 5)

dynamicHMC(m() | (x = randn(2),))

ext: jl


text:
findnmax(x) = round(Int, 2 + x + 4x^(1 / 3))

struct Mie end

function scattering(m, x, ::Mie)

    nmax = findnmax(x)
    an, bn = mie_ab(m, x, nmax)

    return 2 / x^2 * sum((2n + 1) * (abs(an[n])^2 + abs(bn[n])^2) for n = 1:nmax)
end

absorption(m, x, s::Mie) = extinction(m, x, s) - scattering(m, x, s)

function extinction(m, x, ::Mie)

    nmax = findnmax(x)
    an, bn = mie_ab(m, x, nmax)

    return 2 / x^2 * sum((2n + 1) * real(an[n] + bn[n]) for n = 1:nmax)
end

function backscattering(m, x, ::Mie)

    nmax = findnmax(x)
    an, bn = mie_ab(m, x, nmax)

    return 1 / x^2 * abs(sum((2n + 1) * (-1)^n * (an[n] - bn[n]) for n = 1:nmax))^2
end

struct MieProblem
    m
    x
    nmax
    an
    bn

    function MieProblem(m, x)
        nmax = findnmax(x)
        an, bn = mie_ab(m, x, nmax)
        return new(m, x, nmax, an, bn)
    end
end

scattering(s::MieProblem) = 2 / s.x^2 * sum((2n + 1) * (abs(s.an[n])^2 + abs(s.bn[n])^2) for n = 1:s.nmax)
absorption(s::MieProblem) = extinction(s) - scattering(s)
extinction(s::MieProblem) = 2 / s.x^2 * sum((2n + 1) * real(s.an[n] + s.bn[n]) for n = 1:s.nmax)
backscattering(s::MieProblem) = 1 / s.x^2 * abs(sum((2n + 1) * (-1)^n * (s.an[n] - s.bn[n]) for n = 1:s.nmax))^2

@doc raw"""
    mie_ab(m, x, nmax)

Computes external field coefficients ``a_n`` and ``b_n`` based on inputs of `m` and ``x = \pi \, d_p / \lambda``.

# Arguments
- m: The complex refractive index with the convention ``m = n+ik``.
- x: The size parameter ``x = \pi \, d_p / \lambda``.

# Return values
- an, bn: Arrays of size ``n_max = 2 + x + 4x^{1/3}``

"""
function mie_ab(m, x, nmax)

    mx = m * x
    nmx = max(nmax, round(Int, abs(mx))) + 16

    n = 1:nmax

    nu = n .+ 0.5
    sx = sqrt(0.5 * π * x)

    px = sx * besselj.(nu, x)
    p1x = append!([sin(x)], px[1:end-1])

    chx = -sx * bessely.(nu, x)
    ch1x = append!([cos(x)], chx[1:end-1])

    gsx = px - 1im * chx
    gs1x = p1x - 1im * ch1x

    # B&H Equation 4.89
    Dn = zeros(ComplexF64, nmx)
    for i = nmx-1:-1:1
        Dn[i] = (i / mx) - (1 / (Dn[i+1] + i / mx))
    end

    D = Dn[2:nmax+1] # Dn(mx), drop terms beyond nMax

    da = D / m + n / x
    db = m * D + n / x

    an = @. (da * px - p1x) / (da * gsx - gs1x)
    bn = @. (db * px - p1x) / (db * gsx - gs1x)

    return an, bn
end

# spherical bessel functions
sbesselj(nu, x) = √(π / 2x) * besselj(nu + 0.5, x)
sbessely(nu, x) = √(π / 2x) * bessely(nu + 0.5, x)
shankelh1(nu, x) = sbesselj(nu, x) + 1im * sbessely(nu, x)

@doc raw"""
    mie_abcd(m, x, n::Int)

Calculate n-th order Mie coefficients a,b,c,d

# Arguments
- `m`: relative refractive index (``n_\mathrm{material} / n_\mathrm{environment}``)
- `x`: size parameter (wavenumber * radius)
- `n::Int`
"""
function mie_ab_n(m, x, n::Int)

    # Bessel/Hankel functions
    jx = sbesselj(n, x)
    jmx = sbesselj(n, m * x)
    hx = shankelh1(n, x)

    # derivative
    d1xjx = x * sbesselj(n - 1, x) - n * sbesselj(n, x)
    d1mxjmx = m * x * sbesselj(n - 1, m * x) - n * sbesselj(n, m * x)
    d1xhx = x * shankelh1(n - 1, x) - n * shankelh1(n, x)

    an = (m^2 * jmx * d1xjx - jx * d1mxjmx) / (m^2 * jmx * d1xhx - hx * d1mxjmx)
    bn = (jmx * d1xjx - jx * d1mxjmx) / (jmx * d1xhx - hx * d1mxjmx)

    # cn = (jx * d1xhx - hx * d1xjx) / (jmx * d1xhx - hx * d1mxjmx)
    # dn = (m * jx * d1xhx - m * hx * d1xjx) / (m^2 * jmx * d1xhx - hx * d1mxjmx)

    # return an, bn, cn, dn
    return an, bn
end

@doc raw"""
    mie_pi_tau(mu, nmax)

Calculates ``\pi_n`` and ``\tau_n``.

This function uses recurrence relations to calculate ``\pi_n`` and ``\tau_n``, beginning with ``\pi_0 = 1``, ``\pi_1 = 3 \mu`` (where ``mu`` is the cosine of the scattering angle), ``\tau_0 = \mu``, and ``\tau_1 = 3 cos(2 cos^{-1}(\mu))``:
```math
\pi_n = \frac{2n-1}{n-1} \mu \pi_{n-1} - \frac{n}{n-1} \pi_{n-2}
```
```math
\tau_n = n \mu \pi_n - (n+1) \pi_{n-1}
```

# Arguments
- mu : The cosine of the scattering angle.
- nmax: The number of elements to compute. Typically, ``n_{max} = floor(2+x+4x^{1/3})``, but can be given any integer.

# Return values
- p, t: The ``\pi_n`` and ``\tau_n`` arrays, of length `nmax`.
"""
function mie_pi_tau(mu, nmax)

    p = zeros(nmax)
    t = zeros(nmax)

    p[1] = 1
    p[2] = 3mu

    t[1] = mu
    t[2] = 3 * cos(2 * acos(mu))

    for n = 3:nmax
        p[n] = (2n - 1) / (n - 1) * mu * p[n-1] - n / (n - 1) * p[n-2]
        t[n] = n * mu * p[n] - (n + 1) * p[n-1]
    end

    return p, t
end

@doc raw"""
    mie_pi_tau(mu, n::Int)

See Sec. 4.3.1 of B&H ref of pi,tau: [D. Deirmendjian, “Electromagnetic Scattering on Spherical Polydispersions"]
"""
function mie_pi_tau_n(mu, n::Int)
    # See Sec. 4.3.1 of B&H
    # ref of pi,tau: [D. Deirmendjian, “Electromagnetic Scattering on Spherical Polydispersions"]

    if mu < -1 || mu > 1
        error("mu is cos(theta), -1 <= mu <= 1")
    end

    if n == 0
        pin = 0.0  # pi(0, mu) = 0
        taun = 0.0  # not used but tau(0, mu) = 0
    elseif n == 1
        pin = 1.0  # pi(1, mu) = 1
        taun = mu  # note that pi(0, mu) = 0
    elseif n > 1
        # recurrence, eq.(4.47) of B&H
        pin = (2n - 1) / (n - 1) * mu * mie_pi_tau_n(mu, n - 1)[1] - (n) / (n - 1) * mie_pi_tau_n(mu, n - 2)[1]
        taun = n * mu * pin - (n + 1) * mie_pi_tau_n(mu, n - 1)[1]
    else
        error("n should be larger than 0")
    end
    return pin, taun
end

@doc raw"""
    mie_scattering(m, x, [nmax])

Computes Mie efficencies *Q* and asymmetry parameter *g* of a single, homogeneous particle. Uses `mie_ab` to calculate ``a_n`` and ``b_n``, and then calculates *Q* via:
```math
Q_{ext} = \frac{2}{x^2} \sum \limits_{n=1}^{n_{max}} (2n+1) \: \text{Re} \left\{ a_n + b_n \right\}
```
```math
Q_{sca} = \frac{2}{x^2} \sum_{n=1}^{n_{max}} (2n+1) (|a_n|^2 + |b_n|^2)
```
```math
Q_{abs} = Q_{ext} - Q_{sca}
```
```math
Q_{back} = \frac{1}{x^2} \left| \sum\limits_{n=1}^{n_{max}} (2n+1) (-1)^n (a_n - b_n) \right|^2
```
```math
Q_{ratio} = \frac{Q_{back}}{Q_{sca}}
```
```math
g = \frac{4}{Q_{sca}x^2} \left[ \sum\limits_{n=1}^{n_{max}} \frac{n(n+2)}{n+1} \text{Re}\left\{a_n a_{n+1}^* + b_n b_{n+1}^*\right\} + \sum\limits_{n=1}^{n_{max}} \frac{2n+1}{n(n+1)} \text{Re}\left\{a_n b_n^*\right\}\right]
```
```math
Q_{pr} = Q_{ext} - g Q_{sca}
```
where asterisks denote the complex conjugates.

# Arguments
- m: The complex refractive index, with the convention *m = n+ik*.
- wavelength: The wavelength of incident light, in nanometers.
- diameter: The diameter of the particle, in nanometers.
- n_medium: The refractive index of the surrounding medium. This must be positive, nonzero, and real. Any imaginary part will be discarded.

# Return values
- qext, qsca, qabs, g, qpr, qback, qratio: The Mie efficencies described above.

# Examples

For example, compute the Mie efficencies of a particle ``300 \, [\mathrm{nm}]`` in diameter with ``m = 1.77+0.63i``, illuminated by ``\lambda = 375 \, [\mathrm{nm}]``:
```julia-repl
julia> mie_scattering(1.77+0.63im, 375, 300)

(2.858497199156411, 1.3149276685170943, 1.543569530639317, 0.2014551048135256)
```
"""
function mie_scattering(m, x, nmax = findnmax(x))

    if x ≈ 0
        return 0.0, 0.0, 0.0, 0.0
    end

    an, bn = mie_ab(m, x, nmax)

    # Note:
    # Q_i = \sigma_i / (\pi a^2)
    # where
    # - Q_i is the efficiency coefficients
    # - \sigma_{i} is the cross section of the respective process

    qsca = 2 / x^2 * sum((2n + 1) * (abs(an[n])^2 + abs(bn[n])^2) for n = 1:nmax)
    qext = 2 / x^2 * sum((2n + 1) * real(an[n] + bn[n]) for n = 1:nmax)
    qabs = qext - qsca
    qback = 1 / x^2 * abs(sum((2n + 1) * (-1)^n * (an[n] - bn[n]) for n = 1:nmax))^2

    # g = 4 / (qsca * x^2) * (sum(n * (n + 2) / (n + 1) * real(an[n]*an[n+1]'+bn[n]*bn[n+1]') for n = 1:nmax-1) + sum((2n + 1) / (n * (n + 1)) * real(an[n]*bn[n]')  for n = 1:nmax-1))
    # qpr = qext - qsca * g
    # qratio = qback / qsca

    return qext, qsca, qabs, qback
end

@doc raw"""
    mie_scattering_n(m, x, nmax::Int=-1)

Calculate scattering efficiency of a sphere

# Arguments
- `m`: relative refractive index (``n_\mathrm{material} / n_\mathrm{environment}``)
- `x`: size parameter (wavenumber * radius)
- `namx`: maximun order of resonances being considered
      (nmax = 1: dipole, 2: dipole+quadrupole, ...,
       -1: sufficiently large value determined by x)

# Return values
- `Qsca`: Scattering efficiency (total up to n-th order resonances)
"""
function mie_scattering_n(m, x, nmax = findnmax(x))

    if x ≈ 0
        return 0.0
    end

    qsca, qext = 0.0, 0.0
    for n = 1:nmax
        an, bn = mie_ab_n(m, x, n)
        qsca += (2n + 1) * (abs(an)^2 + abs(bn)^2)
        qext += (2n + 1) * real(an + bn)
    end
    qsca *= 2 / x^2
    qext *= 2 / x^2

    qabs = qext - qsca

    return qext, qsca, qabs, nothing
end

@doc raw"""
    mie_S1_S2(m,x,mu)

Calculates ``S_1`` and ``S_2`` at ``\mu = cos(\theta)``, where ``\theta`` is the scattering angle.

Uses `mie_ab` to calculate `a_n` and `b_n`, and `mie_pi_tau` to calculate ``\pi_n`` and ``\tau_n``. ``S_1`` and ``S_2`` are calculated by:

```math
S_1 = \sum\limits_{n=1}^{n_{max}} \frac{2n+1}{n(n+1)} (a_n \pi_n + b_n \tau_n)
```
```math
S_2 = \sum\limits_{n=1}^{n_{max}} \frac{2n+1}{n(n+1)} (a_n \tau_n + b_n \pi_n)
```

# Parameters
- m: The complex refractive index with the convention *m = n+ik*.
- x: The size parameter ``x = \pi \, d_p / \lambda``.
- mu: The cosine of the scattering angle.

# Return values
- S1, S2: The ``S_1`` and ``S_2`` values.

"""
function mie_S1_S2(m, x, mu, nmax = findnmax(x))

    an, bn = mie_ab(m, x, nmax)
    pin, taun = mie_pi_tau(mu, nmax)

    S1 = sum(1:nmax) do n
        (2n + 1) / (n * (n + 1)) * (an[n] * pin[n] + bn[n] * taun[n])
    end

    S2 = sum(1:nmax) do n
        (2n + 1) / (n * (n + 1)) * (an[n] * taun[n] + bn[n] * pin[n])
    end

    return S1, S2
end

@doc raw"""
    mie_S1_S2(m, x, mu, n::Int)

calculate n-th order S1,S2 component of scattering matrix.
"""
function mie_S1_S2_n(m, x, mu, n)

    a, b = mie_ab_n(m, x, n)
    p, t = mie_pi_tau_n(mu, n)

    S1 = (2n + 1) / (n * (n + 1)) * (a * p + b * t)
    S2 = (2n + 1) / (n * (n + 1)) * (a * t + b * p)

    return S1, S2
end

@doc raw"""
    scattering_function(theta, m, x)

Creates arrays for plotting the angular scattering intensity functions in theta-space with parallel, perpendicular, and unpolarized light. Also includes an array of the angles for each step. This angle can be in either degrees, radians, or gradians for some reason. The angles can either be geometrical angle or the qR vector (see *Sorensen, M. Q-space analysis of scattering by particles: a review. J. Quant. Spectrosc. Radiat. Transfer 2013, 131, 3-12 <http://www.sciencedirect.com/science/article/pii/S0022407313000083>*). Uses `mie_S1_S2` to compute ``S_1`` and ``S_2``, then computes parallel, perpendicular, and unpolarized intensities by

```math
SL(\theta) = |S_1|^2
```
```math
SR(\theta) = |S_2|^2
```
```math
SU(\theta) = \frac{1}{2}(SR + SL)
```

# Arguments

- theta: An array of the angles used in calculations.
- m: The complex refractive index with the convention ``m = n+ik``.
- wavelength: The wavelength of incident light, in nanometers.
- diameter: The diameter of the particle, in nanometers.
- n_medium: The refractive index of the surrounding medium. This must be positive, nonzero, and real. Any imaginary part will be discarded.

# Return values

- SL: An array of the scattered intensity of left-polarized (perpendicular) light. Same size as the `theta` array.
- SR: An array of the scattered intensity of right-polarized (parallel) light. Same size as the `theta` array.
- SU: An array of the scattered intensity of unpolarized light, which is the average of SL and SR. Same size as the `theta` array.
"""
function scattering_function(theta, m, x)

    n = length(theta)

    SL = zeros(n)
    SR = zeros(n)
    SU = zeros(n)

    if x ≈ 0
        return SL, SR, SU
    end

    for i = 1:n
        mu = cos(theta[i])
        s1, s2 = mie_S1_S2(m, x, mu)
        SL[i] = abs(s1^2)
        SR[i] = abs(s2^2)
        SU[i] = (SR[i] + SL[i]) / 2
    end

    return SL, SR, SU
end

@doc raw"""
    scattering_function_n(theta, m, x)

"""
function scattering_function_n(theta, m, x)

    n = length(theta)

    SL = zeros(n)
    SR = zeros(n)
    SU = zeros(n)

    if x ≈ 0
        return SL, SR, SU
    end

    for i = 1:n
        mu = cos(theta[i])
        nmax = findnmax(x)
        s1, s2 = 0.0, 0.0
        for n = 1:nmax
            s1n, s2n = mie_S1_S2_n(m, x, mu, n)
            s1 += s1n
            s2 += s2n
        end

        SL[i] = abs(s1^2)
        SR[i] = abs(s2^2)
        SU[i] = (SR[i] + SL[i]) / 2
    end

    return SL, SR, SU
end

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
struct IfElseBoolVariate <: AbstractBoolVariate
    name::String
    ifboolvar::AbstractBoolVariate
    trueboolvar::AbstractBoolVariate
    falseboolvar::AbstractBoolVariate
end

Base.length(var::IfElseBoolVariate) = length(var.ifboolvar)

function slice(boolvar::IfElseBoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    ifboolvar = boolvar.ifboolvar
    trueboolvar = boolvar.trueboolvar
    falseboolvar = boolvar.falseboolvar
    f = ((x, y, z) -> x ? y : z)
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip3(slice(ifboolvar, fromobs, toobs, slicelength), slice(trueboolvar, fromobs, toobs, slicelength),
                  slice(falseboolvar, fromobs, toobs, slicelength))
    mapslice3(f, slices, slicelength, Bool)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
@testset "Test sailing craft performance functions" begin
    path = ENV["HOME"]*"/.julia/v0.6/sail_route/src/data/first40_orgi.csv"
    twa, tws, perf = sail_route.load_file(path)
    @test twa[2] ≈ 60.0 atol=0.0
    @test tws[2] ≈ 8.0 atol=0.0
    @test perf[2, 2] ≈ 6.77 atol=0.0
    test_itp = sail_route.setup_interpolation(tws, twa, perf)
    @test sail_route.perf_interp(test_itp, 60.0, 8.0) ≈ 6.77 atol=0.0
end


text:
using Photon, LightXML

# This file contains a more elaborate example on training a model on real images.
# This example uses the VOC2007 image dataset. Rather then using segmentation,
# this example uses classification. Without a GPU, this training will be slow!
#
# You can read more about this dataset at:
#    http://host.robots.ox.ac.uk/pascal/VOC/voc2007/
#
# And you can download the dataset from:
#    http://host.robots.ox.ac.uk/pascal/VOC/voc2007/VOCtrainval_06-Nov-2007.tar

# Change location below to where you extracted the VOC2007 dataset
const VOC_DIR = "/home/peter/data/VOCdevkit/VOC2007/"

# Shouldn't be a need to change anything below this line
const IMAGEDIR = joinpath(VOC_DIR,"JPEGImages/")
const LABELDIR = joinpath(VOC_DIR,"Annotations/")
const LABELS = ["chair", "diningtable", "person", "boat", "pottedplant", "car", "sofa",
          "sheep", "cat", "aeroplane", "dog", "train", "bottle", "bird", "bicycle",
          "horse", "cow", "tvmonitor", "motorbike", "bus"]

"""
The labels are stored in XML files, so some parsing is required to retrieve
and transform them into onehot encodings. A single image can have multiple
classes (like it contains a car AND and bicycle).
"""
function get_labels(files::Vector{String}):: Vector{Array}
    labels = []
    for f in files
        bf = f[1:end-4]
        filename = joinpath(LABELDIR, bf * ".xml")

        # Get the class labels for this image. These functions are from
        # the excellent LightXML package
        xroot = root(parse_file(filename))
        objects = get_elements_by_tagname(xroot, "object")
        cls = content.(find_element.(objects, "name"))

        # Now onehot encode the labels
        oh = onehot(unique(cls), LABELS, getContext().dtype)
        push!(labels, oh)
    end
    labels
end

# Get all the image filenames
files = readdir(IMAGEDIR)
@info "Found $(length(files)) image files"

labels = get_labels(files)
@info "Parsed $(length(labels)) label files"

# Split into train and validation data
images = joinpath.(IMAGEDIR, files)
data_train, data_valid = ImageDataset(images, labels, resize=(200, 200)) |> Split()
data_train = data_train |> MiniBatch(16)
data_valid = data_valid |> MiniBatch(32, shuffle=false)

# Create a simple convolutional network
model = Sequential(
    Conv2D(16, 3, relu),
    MaxPool2D(),
    Conv2D(64, 3, relu),
    MaxPool2D(),
    Dense(128, relu),
    Dense(20),
    softmax
)

workout = Workout(model, CrossEntropyLoss())

fit!(workout, data_train, data_valid, epochs=50)
@info "Finished training for $(workout.epochs) epochs"

import Plots
plotmetrics(Plots, workout, [:loss, :val_loss])


text:
export DoubleLayerRotatedMW3D

type DoubleLayerRotatedMW3D{T} <: MaxwellOperator3D
    "im times the wavenumber"
    gamma::T
end

function quaddata(operator::DoubleLayerRotatedMW3D,
        local_test_basis::RTRefSpace,
        local_trial_basis::RTRefSpace,
        test_elements, trial_elements)

    test_quad_data  = quadpoints(local_test_basis,  test_elements,  (2,))
    trial_quad_data = quadpoints(local_trial_basis, trial_elements, (3,))

    return test_quad_data, trial_quad_data
end

function quadrule(operator::DoubleLayerRotatedMW3D,
        local_test_basis::RTRefSpace,
        local_trial_basis::RTRefSpace,
        test_id, test_element,
        trial_id, trial_element,
        quad_data)

    test_quad_rules  = quad_data[1]
    trial_quad_rules = quad_data[2]

    DoubleQuadStrategy(
        test_quad_rules[1,test_id],
        trial_quad_rules[1,trial_id]
    )
end

function integrand(op::DoubleLayerRotatedMW3D, kernel_vals, test_vals, test_nbd, trial_vals, trial_nbd)

    n = normal(test_nbd)
    g = test_vals[1]
    f = trial_vals[1]
    ∇G = kernel_vals.gradgreen

    return g ⋅ (n × (∇G × f))
end


text:
module TestModularityLouvain
using Test
using EcologicalNetworks

# Louvain modularity
A = [
  true true true false false false; true true true false false false;
  false false false true true true; false false false true true true
  ]
B = BipartiteNetwork(A)
m = louvain(B, collect(1:richness(B)))
@test m.Q ≈ 0.5
@test Qr(m) ≈ 1.0

end


text:
# quick test
using Caesar
using Neo4j

include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
addrdict["session"] = "SESSTURT45"  # "SESSROX"
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
addrdict["clearslamindb"] = "n"

sessionname = addrdict["session"]
backendset=0
reqbackendset = false
ready=1

cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
conn = cloudGraph.neo4j.connection

sn = length(sessionname) > 0 ? ":"*sessionname : ""
query = "match (n$(sn):POSE) where n.ready=$(ready) and exists(n.exVertexId)"
query = reqbackendset ? query*" and n.backendset=$(backendset)" : query
query = query*" return n.exVertexId, id(n)"
# query = query*" return n"

@show query

loadtx = transaction(conn)
cph = loadtx(query, submit=true)

ret = Array{Tuple{Int64,Int64},1}()
for data in cph.results[1]["data"]
  exvid, neoid = data["row"][1], data["row"][2]
  push!(ret, (exvid,neoid)  )
end

ret


text:
#
# Why is the julia version x10 slower
# than the javascript version?!
#

@enum DD left right
@enum DT round square curly floor ceil vert dbvert hidden
@enum OT text unary binary fname
@enum ST rm sf sc fr tt bb
# @enum UT null terminal styled stack fraction root delimt

PI = 'π'

struct Terminal
	op::OT
	ch::Char
end

struct Styled
	op::OT
	ch::Char
	tf::ST
	bf::Bool
	it::Bool
end

struct Stack{TGroup}
	sup::TGroup
	sub::TGroup
end

struct Fraction{TGroup}
	num::TGroup
	den::TGroup
end

struct Root{TGroup}
	rad::TGroup
	deg::TGroup
end

struct Delimit{TGroup}
	del::DT
	rdl::DT
	val::TGroup
	stk::Union{Stack,Nothing}
end

struct NullBox
end

struct TGroup
	data::Array{Union{Terminal, Styled, Stack{TGroup}, Fraction{TGroup}, Root{TGroup}, Delimit{TGroup}, NullBox}, 1}
end

TUnit = Union{Terminal,Styled,Stack{TGroup},Fraction{TGroup},Root{TGroup},Delimit{TGroup},NullBox}

struct ImageData
	data::Array{Char, 2}
	width::UInt32
	height::UInt32
end

struct Box
	render::Function
	width::UInt32
	height::UInt32
	depth::UInt32
	marleft::UInt32
	marright::UInt32
end

mutable struct BoxMargin
	marleft::UInt32
	marright::UInt32
end

function draw(canvas::ImageData, x::Number, y::Number, value::Char)
	canvas.data[begin + x, begin + y] = value
end

function compile_group(group::TGroup)
	len = length(group.data)
	if len === 1
		return compile_unit(group.data[begin])
	elseif len === 0
		return compile_unit(NullBox())
	end

	boxes = compile_unit.(group.data)
	margins = map(box -> BoxMargin(0, 0), group.data)

	for i in 0:len - 2
		margins[begin + i].marright = max(boxes[begin + i].marright, boxes[begin + i + 1].marleft)
		margins[begin + i + 1].marleft = 0
	end

	width = 0
	height = 0
	depth = 0

	for (box, margin) in zip(boxes, margins)
		width += margin.marleft + box.width + margin.marright
		height = max(height, box.height)
		depth = max(depth, box.depth)
	end

	return Box((canvas::ImageData, x::Number, y::Number) -> begin
		for (box, margin) in zip(boxes, margins)
			x += margin.marleft
			box.render(canvas, x, y)
			x += box.width
			x += margin.marright
		end
	end, width, height, depth, box[begin].marleft, box[end].marright)
end

function compile_stack(unit::Stack, hsp::UInt32, dsp::UInt32)
	sup = compile_group(unit.sup)
	sub = compile_group(unit.sub)

	return Box((canvas::ImageData, x::Number, y::Number) -> begin
		sup.render(canvas, x, y - sup.depth - hsp)
		sub.render(canvas, x, y + sub.height - dsp)
	end, max(sup.width, sub.width), sup.height + sup.depth + hsp, sub.height + sub.depth + dsp, 0, 0)
end

function compile_unit(unit::TUnit)
	type = typeof(unit)

	if type == NullBox
		return Box(() -> nothing, 0, 0, 0, 0, 0)

	elseif type == Terminal
		margin = unit.op == text || unit.op == unary ? 0 : 1

		return Box((canvas::ImageData, x::Number, y::Number) -> begin
			draw(canvas, x, y, unit.ch)
		end, 1, 1, 0, margin, margin)

	elseif type == Stack{TGroup}
		compile_stack(unit, 1, 0)

	elseif type == Fraction{TGroup}
		num = compile_group(unit.num)
		den = compile_group(unit.den)
		width = max(num.width, den.width)
		numoff = div(width - num.width, 2)
		denoff = div(width - den.width, 2)

		return Box((canvas::ImageData, x::Number, y::Number) -> begin
			num.render(canvas, x + 1 + numoff, y - num.depth - 1)
			den.render(canvas, x + 1 + denoff, y + den.height)

			for i in 0:width + 1
				draw(canvas, x + i, y, '─')
			end
		end, width + 2, num.height + num.depth + 1, den.height + den.depth, 1, 1)

	else
		error("Invalid unit type")
	end
end

const resbox = compile_unit(Fraction(TGroup([Terminal(text, 'a')]), TGroup([Terminal(text, 'b')])))
const canvaw = resbox.width
const canvah = resbox.height + resbox.depth
const cnvarr = Array{Char, 2}(undef, canvaw, canvah)
fill!(cnvarr, ' ')
const canvas = ImageData(cnvarr, canvaw, canvah)
resbox.render(canvas, 0, resbox.height - 1)
print("\n ")
for (i, value) in enumerate(canvas.data)
	print(value)

	if (i - 1) % canvaw == canvaw - 1
		print("\n ")
	end
end
print('\n')

max_stars_repo_name: esdmr/terminal-math

max_issues_repo_name: esdmr/terminal-math

max_forks_repo_name: esdmr/terminal-math


text:
# export FDSNevq, FDSNevt, FDSNsta
export FDSNsta

# =============================================================================
# No export

function fdsn_chp(chans::ChanOpts, v::Integer)
  # Parse channel config
  if isa(chans, String)
    C = parse_chstr(chans, ',', true, false)
  elseif isa(chans, Array{String, 1})
    C = parse_charr(chans, '.', true)
  else
    C = copy(chans)
  end
  minreq!(C)
  v > 1 && println(stdout, "Most compact request form = ", C)
  return C
end

"""
    S = FDSNsta(chans, KW)

Retrieve station/channel info for formatted parameter file (or string) `chans`
into an empty SeisData structure.

Standard keywords: rad, reg, src, to, v

Other keywords:
* msr: get MultiStage (full) responses?
* s: Start time
* t: Termination (end) time
* xf: Name of XML file to save station metadata

See also: `web_chanspec`, `parsetimewin`, `get_data!`, `SeisIO.KW`
"""
function FDSNsta( chans::ChanOpts="*";
                    msr::Bool              = false,          # MultiStageResp
                    rad::Array{Float64,1}  = KW.rad,         # Search radius
                    reg::Array{Float64,1}  = KW.reg,         # Search region
                      s::TimeSpec          = 0,              # Start
                    src::String            = KW.src,         # Source server
                      t::TimeSpec          = (-600),         # End or Length (s)
                     to::Int64             = KW.to,          # Read timeout (s)
                      v::Integer           = KW.v,           # Verbosity
                     xf::String            = "FDSNsta.xml"   # XML filename
                 )

  d0, d1 = parsetimewin(s, t)
  v > 0 && @info(tnote("Querying FDSN stations"))
  URL = string(fdsn_uhead(src), "station/1/query")
  BODY = "level=response\nformat=xml\n"
  wc = "*"

  # Add geographic search to BODY
  if !isempty(rad)
    BODY *= string( "latitude=", rad[1], "\n",
                    "longitude=", rad[2], "\n",
                    "minradius=", rad[3], "\n",
                    "maxradius=", rad[4], "\n")
  end
  if !isempty(reg)
    BODY *= string( "minlatitude=", reg[1], "\n",
                    "maxlatitude=", reg[2], "\n",
                    "minlongitude=", reg[3], "\n",
                    "maxlongitude=", reg[4], "\n" )
  end

  # Add channel search to BODY
  if chans == wc
    (isempty(reg) && isempty(rad)) && error("No query! Please specify a search radius, a rectangular search region, or some channels.")
    BODY *= string("* * * * ", d0, " ", d1, "\n")
  else
    C = fdsn_chp(chans, v)
    Nc = size(C,1)
    for i = 1:Nc
      str = ""
      for j = 1:4
        str *= (" " * (isempty(C[i,j]) ? wc : C[i,j]))
      end
      BODY *= string(str, " ", d0, " ", d1, "\n")
    end
  end
  if v > 1
    printstyled("request url:", color=:light_green)
    println(URL)
    printstyled("request body: \n", color=:light_green)
    println(BODY)
  end
  open(xf, "w") do io
    request("POST", URL, webhdr, BODY, response_stream=io)
  end

  # Build channel list
  v > 0 && @info(tnote("Building list of channels"))
  io = open(xf, "r")
  xsta = read(io, String)
  close(io)
  S = FDSN_sta_xml(xsta, msr, d0, d1, v)

  # ===================================================================
  # Logging
  note!(S, string( "+meta ¦ ", URL ))
  for i in 1:S.n
    id = split_id(S.id[i])
    if isempty(id[3])
      id[3] = "--"
    end
    note!(S, i, string("POST ¦ ", join(id, " "), " ", d0, " ", d1))
  end
  # ===================================================================

  return S
end

function FDSNget!(U::SeisData,
              chans::ChanOpts,
                 d0::String,
                 d1::String,
           autoname::Bool,
                fmt::String,
                msr::Bool,
                 nd::Real,
               opts::String,
                rad::Array{Float64,1},
                reg::Array{Float64,1},
                 si::Bool,
                src::String,
                 to::Int64,
                  v::Integer,
                  w::Bool,
                 xf::String,
                  y::Bool)

  parse_err = false
  n_badreq = 0
  wc = "*"
  fname = ""

  # (1) Time-space query for station info
  S = (if si
    FDSNsta(chans,
             msr = msr,
             rad = rad,
             reg = reg,
               s = d0,
             src = src,
               t = d1,
              to = to,
               v = v,
              xf = xf
            )
  else
    SeisData()
  end)

  # (1a) Can we autoname the file? True iff S.n == 1
  (S.n == 1) || (autoname = false)

  # (2) Build ID strings for data query
  ID_str = Array{String,1}(undef, S.n)
  if S.n > 0
    N_ch = S.n
    for i in 1:N_ch
      ID_mat = split(S.id[i], ".")
      ID_mat[isempty.(ID_mat)] .= wc
      ID_str[i] = join(ID_mat, " ")
    end
  else
    C = fdsn_chp(chans, v)[:,1:4]
    C[isempty.(C)] .= "*"
    N_ch = size(C,1)
    ID_str = [join(C[i, :], " ") for i in 1:N_ch]
  end
  if v > 1
    printstyled("data query strings:\n", color=:light_green)
    for i = 1:length(ID_str)
      println(stdout, ID_str[i])
    end
  end

  # (3) Data query
  v > 0 && @info(tnote("Data query begins"))
  URL = string(fdsn_uhead(src), "dataselect/1/query")
  if occursin("ncedc", URL) || occursin("scedc", URL)
    BODY = ""
    if fmt != "miniseed"
      @warn(string("format ", fmt, " ignored; server only allows miniseed."))
    end
  elseif occursin("ph5ws",URL)
    BODY = "reqtype=FDSN\n"
    if fmt ∉ ["mseed","miniseed","sac","segy1","segy2","geocsv","geocsv.tspair","geocsv.slist"]
      @warn(string("format ", fmt, " ignored; server only allows:\n" *
        "mseed\n" *
        "sac\n" *
        "segy1\n" *
        "geocsv\n" *
        "geocsv.tspair\n" *
        "geocsv.slist")
      )
    end
    fmt = fmt == "miniseed" ? "mseed" : fmt
    BODY *= "format=" * fmt * "\n"
    if !isempty(opts)
      OPTS = split(opts, "&")
      for opt in OPTS
        BODY *= string(opt, "\n")
      end
    end
  else
    BODY = "format=" * fmt * "\n"
    if !isempty(opts)
      OPTS = split(opts, "&")
      for opt in OPTS
        BODY *= string(opt, "\n")
      end
    end
  end

  # Set the data source
  fill!(S.src, URL)

  # Create variables for query
  ts = tstr2int(d0)
  ti = round(Int64, nd*86400000000)
  te = tstr2int(d1)
  t1 = deepcopy(ts)
  rn = 0
  while t1 < te
    rn += 1
    os = rn > 1 ? 1 : 0
    t1 = min(ts + ti, te)
    s_str = int2tstr(ts + os)
    t_str = int2tstr(t1)
    qtail = string(" ", s_str, " ", t_str, "\n")
    QUERY = identity(BODY)
    for i = 1:N_ch
      QUERY *= ID_str[i]*qtail
    end
    if v > 1
      printstyled(string("request url: ", URL, "\n"), color=:light_green)
      printstyled(string("request body: \n", QUERY), color=:light_green)
    end

    # Request via "POST"
    if w
      if fmt == "miniseed"
        ext = "mseed"
      else
        ext = fmt
      end

      # Generate filename
      yy = s_str[1:4]
      mm = s_str[6:7]
      dd = s_str[9:10]
      HH = s_str[12:13]
      MM = s_str[15:16]
      SS = s_str[18:19]
      nn = lpad(div(parse(Int64, s_str[21:26]), 1000), 3, '0')
      jj = lpad(md2j(yy, mm, dd), 3, '0')
      if autoname
        fname = join([yy, jj, HH, MM, SS, nn, S.id[1], "R", ext], '.')
      else
        fname = join([yy, jj, HH, MM, SS, nn, "FDSNWS", src, ext], '.')
      end
      safe_isfile(fname) && @warn(string("File ", fname, " contains an identical request; overwriting."))
      open(fname, "w") do io
        request("POST", URL, webhdr, QUERY, readtimeout=to, response_stream=io)
      end
      io = open(fname, "r")
      parsable = true
    else
      (R, parsable) = get_http_post(URL, QUERY, to)
      io = IOBuffer(R)
    end
    (v > 1) && println("parsable = ", parsable)

    # Parse data (if we can)
    if parsable
      if fmt == "mseed" || fmt == "miniseed"
        parsemseed!(S, io, KW.nx_add, KW.nx_add, true, v)
      elseif fmt == "geocsv" || fmt == "geocsv.tspair"
        read_geocsv_tspair!(S, io)
      elseif fmt == "geocsv.slist"
        read_geocsv_slist!(S, io)
      else
        parse_err = true
        n_badreq += 1
        push!(S, SeisChannel(id = string("XX.FMT..", lpad(n_badreq, 3, "0")),
                             misc = Dict{String,Any}(  "url" => URL,
                                                      "body" => QUERY,
                                                       "raw" => read(io))))
        note!(S, S.n, "unparseable format; raw bytes in :misc[\"raw\"]")
      end
    else
      # Should only happen with an error message (parsable to String) in io
      parse_err = true
      n_badreq += 1
      push!(S, SeisChannel(id = string("XX.FAIL..", lpad(n_badreq, 3, "0")),
                           misc = Dict{String,Any}( "url" => URL,
                                                   "body" => QUERY,
                                                    "msg" => String(read(io)))))
      note!(S, S.n, "request failed; response in :misc[\"msg\"]")
    end
    close(io)
    ts += ti
  end

  # ===================================================================
  # Logging
  note!(S, string( "+source ¦ ", URL ))
  for i in 1:S.n
    id = split_id(S.id[i])
    if isempty(id[3])
      id[3] = "--"
    end
    note!(S, i, string("POST ¦ ", join(id, " "), " ", d0, " ", d1))
    if w
      wstr = string(timestamp(), " ¦ write ¦ get_data(\"FDSN\" ... w=true) ¦ wrote raw download to file ", fname)
      for i in 1:S.n
        push!(S.notes[i], wstr)
      end
    end
  end
  # ===================================================================

  append!(U,S)
  # Done!
  v > 0 && @info(string(timestamp(), ": done FDSNget query."))
  v > 1 && @info(string("n_badreq = ", n_badreq))
  return parse_err
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function geometric_mean(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "geometricMean", jdouble, (Vector{jdouble},), arg0)
end

function geometric_mean(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "geometricMean", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function max(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "max", jdouble, (Vector{jdouble},), arg0)
end

function max(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "max", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function mean(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "mean", jdouble, (Vector{jdouble},), arg0)
end

function mean(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "mean", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function mean_difference(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(StatUtils, "meanDifference", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function min(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "min", jdouble, (Vector{jdouble},), arg0)
end

function min(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "min", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function mode(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "mode", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function mode(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "mode", Vector{jdouble}, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function normalize(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "normalize", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function percentile(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(StatUtils, "percentile", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function percentile(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint, arg3::jdouble)
    return jcall(StatUtils, "percentile", jdouble, (Vector{jdouble}, jint, jint, jdouble), arg0, arg1, arg2, arg3)
end

function population_variance(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "populationVariance", jdouble, (Vector{jdouble},), arg0)
end

function population_variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(StatUtils, "populationVariance", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function population_variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jdouble, arg2::jint, arg3::jint)
    return jcall(StatUtils, "populationVariance", jdouble, (Vector{jdouble}, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function population_variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "populationVariance", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function product(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "product", jdouble, (Vector{jdouble},), arg0)
end

function product(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "product", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function sum(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "sum", jdouble, (Vector{jdouble},), arg0)
end

function sum(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "sum", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function sum_difference(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(StatUtils, "sumDifference", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function sum_log(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "sumLog", jdouble, (Vector{jdouble},), arg0)
end

function sum_log(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "sumLog", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function sum_sq(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "sumSq", jdouble, (Vector{jdouble},), arg0)
end

function sum_sq(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "sumSq", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function variance(::Type{StatUtils}, arg0::Vector{jdouble})
    return jcall(StatUtils, "variance", jdouble, (Vector{jdouble},), arg0)
end

function variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(StatUtils, "variance", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jdouble, arg2::jint, arg3::jint)
    return jcall(StatUtils, "variance", jdouble, (Vector{jdouble}, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end

function variance(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(StatUtils, "variance", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function variance_difference(::Type{StatUtils}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(StatUtils, "varianceDifference", jdouble, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end


text:
@doc raw"""
	spectral_moment(order::Int, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator; nodes::Int=31, control_freqs::Vector{Float64}=[1e-3, 1e-1, 1.0, 10.0, 100.0, 300.0] ) where {S<:Real,T<:Real}

Compute spectral moment of a specified order.

Evaluates the expression:
```math
	m_k = 2\int_{0}^{\infty} \left(2\pi f\right)^k |H(f;f_n,\zeta_n)|^2 |A(f)|^2 df
```
where ``k`` is the order of the moment.

Integration is performed using Gauss-Legendre integration using `nodes` nodes and weights.
The integration domain is partitioned over the `control_freqs` as well as two inserted frequencies at `f_n/1.5` and `f_n*1.5` in order to ensure good approximation of the integral around the `sdof` resonant frequency.

See also: [`spectral_moments`](@ref)
"""
function spectral_moment(order::Int, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator; nodes::Int=31, control_freqs::Vector{Float64}=[1e-3, 1e-1, 1.0, 10.0, 100.0, 300.0] ) where {S<:Real,T<:Real}
	# pre-allocate for squared fourier amplitude spectrum
	if S <: Dual
		Af2 = Vector{S}(undef, nodes)
		Hf2 = Vector{S}(undef, nodes)
		int_sum = zero(S)
	else
		U = get_parametric_type(fas)
		Af2 = Vector{U}(undef, nodes)
		Hf2 = Vector{U}(undef, nodes)
		int_sum = zero(U)
	end
	# partition the integration domain to make sure the integral captures the key change in the integrand
	f_n = sdof.f_n
	# note that we start slightly above zero to avoid a numerical issue with the frequency dependent Q(f) gradients
	fidLO = findall(control_freqs .< f_n/1.5)
	fidHI = findall(control_freqs .> f_n*1.5)
	# perform the integration with a logarithmic transformation
	lnflims = log.([ control_freqs[fidLO]; f_n/1.5; f_n*1.5; control_freqs[fidHI] ])

	# compute the Gauss Legendre nodes and weights
	xi, wi = gausslegendre(nodes)

	for i in 2:length(lnflims)
		@inbounds dfac = (lnflims[i]-lnflims[i-1])/2
		@inbounds pfac = (lnflims[i]+lnflims[i-1])/2
		lnfi = @. dfac * xi + pfac
		fi = exp.(lnfi)
		squared_transfer!(Hf2, fi, sdof)
		squared_fourier_spectrum!(Af2, fi, m, r_ps, fas)
		# note that the integrand here is scaled by exp(lnfi)=fi for the logarithmic transformation of the integrand
		Yf2 = @. (2π * fi)^order * Hf2 * Af2 * fi
		# weighted combination of amplitudes with Gauss-Legendre weights
		int_sum += dfac * dot( wi, Yf2 )
	end
    return 2.0 * int_sum
end

@doc raw"""
	spectral_moments(order::Vector{Int}, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator; nodes::Int=31, control_freqs::Vector{Float64}=[1e-3, 1e-1, 1.0, 10.0, 100.0, 300.0] ) where {S<:Real,T<:Real}

Compute a vector of spectral moments for the specified `order`.

Evaluates the expression:
```math
	m_k = 2\int_{0}^{\infty} \left(2\pi f\right)^k |H(f;f_n,\zeta_n)|^2 |A(f)|^2 df
```
for each order, where ``k`` is the order of the moment.

Integration is performed using Gauss-Legendre integration using `nodes` nodes and weights.
The integration domain is partitioned over the `control_freqs` as well as two inserted frequencies at `f_n/1.5` and `f_n*1.5` in order to ensure good approximation of the integral around the `sdof` resonant frequency.

See also: [`spectral_moment`](@ref), [`spectral_moments_gk`](@ref)
"""
function spectral_moments(order::Vector{Int}, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator; nodes::Int=31, control_freqs::Vector{Float64}=[1e-3, 1e-1, 1.0, 10.0, 100.0, 300.0] ) where {S<:Real,T<:Real}
	# pre-allocate for squared fourier amplitude spectrum
	if S <: Dual
		Af2 = Vector{S}(undef, nodes)
		Hf2 = Vector{S}(undef, nodes)
		int_sumi = zeros(S, length(order))
	else
		U = get_parametric_type(fas)
		Af2 = Vector{U}(undef, nodes)
		Hf2 = Vector{U}(undef, nodes)
		int_sumi = zeros(U, length(order))
	end
	# partition the integration domain to make sure the integral captures the key change in the integrand
	f_n = sdof.f_n
 	# note that we start slightly above zero to avoid a numerical issue with the frequency dependent Q(f) gradients
	fidLO = findall(control_freqs .< f_n/1.5)
	fidHI = findall(control_freqs .> f_n*1.5)
	# perform the integration with a logarithmic transformation
	lnflims = log.([ control_freqs[fidLO]; f_n/1.5; f_n*1.5; control_freqs[fidHI] ])

	# compute the Gauss Legendre nodes and weights
	xi, wi = gausslegendre(nodes)

	# make sure the orders are listed as increasing for the following loop approach
	sort!(order)
	dorder = diff(order)

	for i in 2:length(lnflims)
		@inbounds dfac = (lnflims[i]-lnflims[i-1])/2
		@inbounds pfac = (lnflims[i]+lnflims[i-1])/2
		lnfi = @. dfac * xi + pfac
		fi = exp.(lnfi)
		squared_transfer!(Hf2, fi, sdof)
		squared_fourier_spectrum!(Af2, fi, m, r_ps, fas)
		# compute default zeroth order integrand amplitude
		# note that the integrand here is scaled by exp(lnfi)=fi for the logarithmic transformation of the integrand
		Yf2 = @. Hf2 * Af2 * fi

		for (idx, o) in enumerate(order)
			if idx == 1
				Yf2 .*= (2π * fi).^o
			else
				@inbounds Yf2 .*= (2π * fi).^(dorder[idx-1])
			end
			# weighted combination of amplitudes with Gauss-Legendre weights
			@inbounds int_sumi[idx] += dfac * dot( wi, Yf2 )
		end
	end
    return 2.0 * int_sumi
end

@doc raw"""
	spectral_moments_gk(order::Vector{Int}, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator) where {S<:Real,T<:Real}

Compute a vector of spectral moments for the specified `order` using Gauss-Kronrod integration from the `QuadGK.jl` package.

Evaluates the expression:
```math
	m_k = 2\int_{0}^{\infty} \left(2\pi f\right)^k |H(f;f_n,\zeta_n)|^2 |A(f)|^2 df
```
for each order, where ``k`` is the order of the moment.

Integration is performed using adaptive Gauss-Kronrod integration with the domain split over two intervals from ``[0,f_n]`` and ``[f_n,\infty]`` to ensure that the resonant peak is not missed.

Note that due to the default tolerances, the moments computed by this method are more accurate than those from `spectral_moments` using the Gauss-Legendre approximation. However, this method is also significantly slower, and cannot be used within an automatic differentiation environment.

See also: [`spectral_moment`](@ref), [`spectral_moments`](@ref)
"""
function spectral_moments_gk(order::Vector{Int}, m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator) where {S<:Real,T<:Real}
	int_sumi = zeros(length(order))
	for (i, o) in enumerate(order)
		moment_integrand(f) = squared_transfer(f, sdof) * squared_fourier_spectral_ordinate(f, m, r_ps, fas) * (2π*f)^o
		int_sumi[i] = quadgk(moment_integrand, 0.0, sdof.f_n, Inf)[1]
	end
    return 2.0 * int_sumi
end

@doc raw"""
	zeros_extrema_frequencies(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator) where {S<:Real,T<:Real}

Defines the frequencies of extrema and zero-crossings using moments ``m_0``, ``m_2`` and ``m_4``. Returns a tuple of ``(f_z,f_e)``.

The frequency of zero crossings is:
```math
f_{z} = \frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}}
```
The frequency of extrema is:
```math
f_{e} = \frac{1}{2\pi}\sqrt{\frac{m_4}{m_2}}
```

See also: [`zeros_extrema_numbers`](@ref)
"""
function zeros_extrema_frequencies(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator) where {S<:Real,T<:Real}
	mi = spectral_moments([0, 2, 4], m, r_ps, fas, sdof)
	m_0 = mi[1]
	m_2 = mi[2]
	m_4 = mi[3]
	fzero = sqrt( m_2 / m_0 ) / (2π)
	fextrema = sqrt( m_4 / m_2 ) / (2π)
	return fzero, fextrema
end

@doc raw"""
	zeros_extrema_numbers(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Real,T<:Real}

Defines the numbers of extrema and zero-crossings using moments ``m_0``, ``m_2`` and ``m_4``. Returns a tuple of ``(2f_z D_{ex}, 2f_e D_{ex})``.

The frequency of zero crossings is:
```math
n_{z} = \frac{D_{ex}}{\pi}\sqrt{\frac{m_2}{m_0}}
```
The frequency of extrema is:
```math
n_{e} = \frac{D_{ex}}{\pi}\sqrt{\frac{m_4}{m_2}}
```
In both cases, the excitaion duration, ``D_{ex}`` is computed from the `excitation_duration` function.

See also: [`zeros_extrema_frequencies`](@ref), [`excitation_duration`](@ref)
"""
function zeros_extrema_numbers(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Real,T<:Real}
    fz, fe = zeros_extrema_frequencies(m, r_ps, fas, sdof)
    Dex = excitation_duration(m, r_ps, fas, rvt)
    return 2fz*Dex, 2fe*Dex
end

@doc raw"""
	rvt_response_spectral_ordinate(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Real,T<:Real}

Response spectral ordinate (units of ``g``) for the specified scenario.

The spectral ordinate is computed using the expression:
```math
S_a = \psi \sqrt{\frac{m_0}{D_{rms}}}
```
where ``\psi`` is the peak factor computed from `peak_factor`, ``m_0`` is the zeroth order spectral moment from `spectral_moment`, and ``D_{rms}`` is the RMS duration computed from `rms_duration`.

See also: [`rvt_response_spectral_ordinate`](@ref), [`rvt_response_spectrum`](@ref), [`rvt_response_spectrum!`](@ref)
"""
function rvt_response_spectral_ordinate(m::S, r_ps::T, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Real,T<:Real}
	# get the duration metrics
	Drms, Dex, Dratio = rms_duration(m, r_ps, fas, sdof, rvt)
	# get the rms response
	m_0 = spectral_moment(0, m, r_ps, fas, sdof)
	y_rms = sqrt( m_0 / Drms )
	# get the peak factor
	pf = peak_factor(m, r_ps, Dex, m_0, fas, sdof, rvt)
	# response spectral value (in g)
	Sa = pf * y_rms / 9.80665
	return Sa
end

@doc raw"""
	rvt_response_spectral_ordinate(period::U, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Float64}

Response spectral ordinate (units of ``g``) for the specified scenario.

The spectral ordinate is computed using the expression:
```math
S_a = \psi \sqrt{\frac{m_0}{D_{rms}}}
```
where ``\psi`` is the peak factor computed from `peak_factor`, ``m_0`` is the zeroth order spectral moment from `spectral_moment`, and ``D_{rms}`` is the RMS duration computed from `rms_duration`.

See also: [`rvt_response_spectral_ordinate`](@ref), [`rvt_response_spectrum`](@ref), [`rvt_response_spectrum!`](@ref)
"""
function rvt_response_spectral_ordinate(period::U, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Float64}
  	# create a sdof instance
  	sdof = Oscillator(1.0/period)
	return rvt_response_spectral_ordinate(m, r_ps, fas, sdof, rvt)
end

@doc raw"""
	rvt_response_spectrum(period::Vector{U}, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Float64}

Response spectrum (units of ``g``) for the vector of periods `period` and the specified scenario.

Each spectral ordinate is computed using the expression:
```math
S_a = \psi \sqrt{\frac{m_0}{D_{rms}}}
```
where ``\psi`` is the peak factor computed from `peak_factor`, ``m_0`` is the zeroth order spectral moment from `spectral_moment`, and ``D_{rms}`` is the RMS duration computed from `rms_duration`. The various terms are all functions of the oscillator period.

See also: [`rvt_response_spectral_ordinate`](@ref), [`rvt_response_spectrum!`](@ref)
"""
function rvt_response_spectrum(period::Vector{U}, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Float64}
	if S <: Dual
		Sa = zeros(S,length(period))
	else
		V = get_parametric_type(fas)
		Sa = zeros(V,length(period))
	end
	for i in 1:length(period)
		@inbounds Sa[i] = rvt_response_spectral_ordinate(period[i], m, r_ps, fas, rvt)
  	end
  	return Sa
end

@doc raw"""
	rvt_response_spectrum!(Sa::Vector{U}, period::Vector{V}, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}

In-place response spectrum (units of ``g``) for the vector of periods `period` and the specified scenario.

Each spectral ordinate is computed using the expression:
```math
S_a = \psi \sqrt{\frac{m_0}{D_{rms}}}
```
where ``\psi`` is the peak factor computed from `peak_factor`, ``m_0`` is the zeroth order spectral moment from `spectral_moment`, and ``D_{rms}`` is the RMS duration computed from `rms_duration`. The various terms are all functions of the oscillator period.

See also: [`rvt_response_spectral_ordinate`](@ref), [`rvt_response_spectrum!`](@ref)
"""
function rvt_response_spectrum!(Sa::Vector{U}, period::Vector{V}, m::S, r_ps::T, fas::FourierParameters, rvt::RandomVibrationParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}
	for i in 1:length(period)
		@inbounds Sa[i] = rvt_response_spectral_ordinate(period[i], m, r_ps, fas, rvt)
  	end
end


text:
struct DFP{T1} <: QuasiNewton{T1}
   approx::T1
end
# function update!(scheme::DFP, B::InverseApprox, s, y)
#    B.A = B.A + s*s'/dot(s, y) - B.A*y*y'*B.A/(y'*B.A*y)
# end
# function update!(scheme::DFP, B::DirectApprox, s, y)
#    B.A = (I - y*s'/dot(y, s))*B.A*(I - s*y'/dot(y, s)) + y*y'/dot(y, s)
# end

function update(H, s, y, scheme::DFP{<:InverseApprox})
    ρ = inv(dot(s, y))
    if ρ > 1e6
        H = H + ρ*s*s' - H*(y*y')*H/(y'*H*y)
    end
    H
end
function update(B, s, y, scheme::DFP{<:DirectApprox})
    ρ = inv(dot(s, y))

    if ρ > 1e6
        C = (I - ρ*y*s')
        B = C*B*C' + ρ*y*y'
    else
        if isa(B, UniformScaling)
            B = B + zero(eltype(s))* y * s'
        end
    end
    B
end
function update!(H, s, y, scheme::DFP{<:InverseApprox})
    ρ = inv(dot(s, y))
    if ρ > 1e6
        H .+= ρ*s*s' - H*(y*y')*H/(y'*H*y)
    end
    H
end
function update!(B, s, y, scheme::DFP{<:DirectApprox})
    ρ = inv(dot(s, y))
    # so right now, we just skip the update if dsy is zero
    # but we might do something else here
    if ρ > 1e6
        C = (I - ρ*y*s')
        B .= C*B*C' + ρ*y*y'
    end
    B
end
update!(A::UniformScaling, s, y , scheme::DFP{<:InverseApprox}) = update(A, s, y, scheme)
update!(A::UniformScaling, s, y , scheme::DFP{<:DirectApprox}) = update(A, s, y, scheme)


text:
#=

reveal=true anstelle von show=true ?
reveal(p)

2D Arrays in 1D plot wie in plts: Label als array, color als array.
Das geht auch einfach in Makie

pass=(...): kwargs für plotter
=#

initialize!(p, ::Type{PlotsType})=nothing

function reveal(p::GridVisualizer,::Type{PlotsType})
    Plots=p.Plotter
    subplots=Plots.Plot[]

    for i=1:size(p.subplots,1)
        for j=1:size(p.subplots,2)
            ctx=p.subplots[i,j]
            if haskey(ctx,:ax)
                push!(subplots,ctx[:ax])
                delete!(ctx,:ax)
            else
                push!(subplots,Plots.plot(legend=false,grid=false,border=:none))
            end
        end
    end
    p=Plots.plot(subplots...,layout=p.context[:layout],size=p.context[:resolution])
    Plots.gui(p)
    p
end

function reveal(ctx::SubVisualizer,TP::Type{PlotsType})
    if ctx[:show]||ctx[:reveal]
        reveal(ctx[:GridVisualizer],TP)
    end
end

function save(fname,scene,Plots,::Type{PlotsType})
   Plots.savefig(scene,fname)
end

function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}},grid)
    Plots=ctx[:Plotter]

    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    p=ctx[:ax]

    cellregions=grid[CellRegions]
    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]
    ncellregions=grid[NumCellRegions]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbfaceregions=grid[NumBFaceRegions]

    xmin=minimum(coord)
    xmax=maximum(coord)
    h=(xmax-xmin)/20.0

    cmap=region_cmap(ncellregions)
    for icell=1:num_cells(grid)
        x1=coord[1,cellnodes[1,icell]]
        x2=coord[1,cellnodes[2,icell]]
        Plots.plot!(p,[x1,x1],[-h,h],linewidth=0.5,color=:black,label="")
        Plots.plot!(p,[x2,x2],[-h,h],linewidth=0.5,color=:black,label="")
        Plots.plot!(p,[x1,x2],[0,0],linewidth=3.0,color=cmap[cellregions[icell]],label="")
    end

    cmap=bregion_cmap(nbfaceregions)
    for ibface=1:num_bfaces(grid)
        if bfaceregions[ibface]>0
            x1=coord[1,bfacenodes[1,ibface]]
            Plots.plot!(p,[x1,x1],[-2*h,2*h],linewidth=3.0,color=cmap[bfaceregions[ibface]],label="")
        end
    end
    reveal(ctx,TP)
end

function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}},grid)
    Plots=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    p=ctx[:ax]
    cellregions=grid[CellRegions]
    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]
    ncellregions=grid[NumCellRegions]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbfaceregions=grid[NumBFaceRegions]

    cmap=region_cmap(ncellregions)
    for icell=1:num_cells(grid)
        inode1=cellnodes[1,icell]
        inode2=cellnodes[2,icell]
        inode3=cellnodes[3,icell]
        # https://github.com/JuliaPlots/Plots.jl/issues/605
        tri=Plots.Shape([coord[1,inode1],coord[1,inode2], coord[1,inode3]],[coord[2,inode1],coord[2,inode2],coord[2,inode3]])
        Plots.plot!(p,tri,color=cmap[cellregions[icell]],label="")
    end
    for icell=1:num_cells(grid)
        inode1=cellnodes[1,icell]
        inode2=cellnodes[2,icell]
        inode3=cellnodes[3,icell]
        Plots.plot!(p, [coord[1,inode1],coord[1,inode2]],[coord[2,inode1],coord[2,inode2]]  ,linewidth=0.5,color=:black,label="")
        Plots.plot!(p, [coord[1,inode1],coord[1,inode3]],[coord[2,inode1],coord[2,inode3]]  ,linewidth=0.5,color=:black,label="")
        Plots.plot!(p, [coord[1,inode2],coord[1,inode3]],[coord[2,inode2],coord[2,inode3]]  ,linewidth=0.5,color=:black,label="")
    end

    cmap=bregion_cmap(nbfaceregions)
    for ibface=1:num_bfaces(grid)
        inode1=bfacenodes[1,ibface]
        inode2=bfacenodes[2,ibface]
        Plots.plot!(p,[coord[1,inode1],coord[1,inode2]],[coord[2,inode1],coord[2,inode2]]  ,linewidth=5,
                    color=cmap[bfaceregions[ibface]],
                    label="")
    end
    reveal(ctx,TP)
end

function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}},grid, func)

    legpos=Dict(
    :none => :none,
    :best => :best,
    :lt => :topleft,
    :ct => :topcenter,
    :rt => :topright,
    :lc => :centerleft,
    :rc => :centerright,
    :lb => :bottomleft,
    :cb => :bottomcenter,
    :rb => :bottomright
)

    Plots=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    p=ctx[:ax]

    coord=grid[Coordinates]
    xmin=coord[1,1]
    xmax=coord[1,end]
    ymin=func[1]
    ymax=func[end]
    xlimits=ctx[:xlimits]
    ylimits=ctx[:flimits]
    if xlimits[1]<xlimits[2]
        xmin=xlimits[1]
        xmax=xlimits[2]
    end
    if ylimits[1]<ylimits[2]
        ymin=ylimits[1]
        ymax=ylimits[2]
    end

    Plots.plot!(p,[xmin,xmax], [ymin,ymax], seriestype = :scatter,
                makersize=0,
                markercolor=:white,
                markerstrokecolor=:white,
                label="")

    color=ctx[:color]
    if ctx[:cellwise]
        cellnodes=grid[CellNodes]
        for icell=1:num_cells(grid)
            i1=cellnodes[1,icell]
            i2=cellnodes[2,icell]
            x1=coord[1,i1]
            x2=coord[1,i2]
            if icell==1 && ctx[:label] !=" "
                Plots.plot!(p,[x1,x2],[func[i1],func[i2]],linecolor=Plots.RGB(color...),label=ctx[:label])
            else
                Plots.plot!(p,[x1,x2],[func[i1],func[i2]],linecolor=Plots.RGB(color...),label="")
            end
        end
    else
        markevery=ctx[:markevery]
        markershape=ctx[:markershape]
        X=vec(grid[Coordinates])
        if markershape==:none
            Plots.plot!(p,X,func,
                        linecolor=Plots.RGB(color),
                        linewidth=ctx[:linewidth],
                        linestyle=ctx[:linestyle],
                        legend=legpos[ctx[:legend]],
                        label=ctx[:label])
        else
            #Trick plots to use markers
            Plots.plot!(p,X,func,
                        linecolor = Plots.RGB(color),
                        linewidth=ctx[:linewidth],
                        linestyle=ctx[:linestyle],
                        label="")
            Plots.plot!(p,[X[1]], [func[1]],
                        markershape = markershape,
                        label = ctx[:label],
                        markersize = ctx[:markersize] ,
                        linecolor = Plots.RGB(color),
                        linewidth=ctx[:linewidth],
                        legend=legpos[ctx[:legend]],
                        linestyle=ctx[:linestyle],
                        markercolor = Plots.RGB(color))
            @views Plots.plot!(p,X[1:markevery:end],func[1:markevery:end],
                               markercolor=Plots.RGB(color),
                               label="",
                               linecolor = :white,
                               markershape = markershape ,
                               markersize = ctx[:markersize] ,
                               lines=false)
        end
    end
    reveal(ctx,TP)
end

"""
$(SIGNATURES)
Return rectangular grid data + function to be splatted into Plots calls
"""
function rectdata(grid,U)
    if dim_grid(grid)==1 && haskey(grid,XCoordinates)
        return grid[XCoordinates],U
    end
    if dim_grid(grid)==2 && haskey(grid,XCoordinates) && haskey(grid,YCoordinates)
        X=grid[XCoordinates]
        Y=grid[YCoordinates]
        return X,Y,transpose(reshape(U,length(X),length(Y)))
    end
    nothing
end

function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}},grid, func)
    rdata=rectdata(grid,func)
    if rdata==nothing
        return nothing
    end
    Plots=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    p=ctx[:ax]
    umin=minimum(func)
    umax=maximum(func)
    flimits=ctx[:flimits]
    if flimits[1]<flimits[2]
        umin=flimits[1]
        umax=flimits[2]
    end
    levels=collect(umin:(umax-umin)/(ctx[:isolines]-1):umax)
    Plots.contourf!(p,rdata...,aspect_ratio=ctx[:aspect],fill=ctx[:colormap],c=:black,levels=levels)
    reveal(ctx,TP)
end

function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{3}}, grid)
    Plots=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    reveal(ctx,TP)
end

function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{3}}, grid,func)
    Plots=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=Plots.plot(title=ctx[:title])
    end
    reveal(ctx,TP)
end


text:
# # Sensitivity Analysis of SVM

#md # [![](https://img.shields.io/badge/GitHub-100000?style=for-the-badge&logo=github&logoColor=white)](@__REPO_ROOT_URL__/examples/sensitivity-analysis-svm.jl)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/sensitivity-analysis-svm.ipynb)

# This notebook illustrates sensitivity analysis of data points in a [Support Vector Machine](https://en.wikipedia.org/wiki/Support-vector_machine) (inspired from [@matbesancon](http://github.com/matbesancon)'s [SimpleSVMs](http://github.com/matbesancon/SimpleSVMs.jl).)

# For reference, Section 10.1 of https://online.stat.psu.edu/stat508/book/export/html/792 gives an intuitive explanation of what it means to have a sensitive hyperplane or data point. The general form of the SVM training problem is given below (without regularization):

# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & \sum_{i=1}^{N} \xi_{i} \\
# \mbox{s.t.} & \xi_{i} \ge 0 \quad i=1..N  \\
#             & y_{i} (w^T X_{i} + b) \ge 1 - \xi[i]\\
# \end{array}
# \end{split}
# ```
# where
# - `X`, `y` are the `N` data points
# - `ξ` is the soft-margin loss.

# ## Define and solve the SVM

# Import the libraries.

using SCS, DiffOpt, LinearAlgebra, JuMP
import Random, Plots

# Construct separable, non-trivial data points.

N = 100
D = 2
Random.seed!(62)
X = vcat(randn(N ÷ 2, D), randn(N ÷ 2, D) .+ [4.0, 1.5]')
y = append!(ones(N ÷ 2), -ones(N ÷ 2));

# Let's initialize a special model that can understand sensitivities

model = Model(() -> diff_optimizer(SCS.Optimizer))

# Add the variables

@variable(model, l[1:N])
@variable(model, w[1:D])
@variable(model, b);

# Add the constraints.

@constraint(
    model,
    1.0 * l ∈ MOI.Nonnegatives(N),
)
@constraint(
    model,
    cons,
    y .* (X * w .+ b) + l .- 1 ∈ MOI.Nonnegatives(N),
);

# Define the objective and solve

@objective(
    model,
    Min,
    sum(l),
)

optimize!(model) # solve

# We can visualize the separating hyperplane.

loss = objective_value(model)
wv = value.(w)
bv = value(b)

svm_x = [0.0, 3.0]
svm_y = (-bv .- wv[1] * svm_x )/wv[2]

p = Plots.scatter(X[:,1], X[:,2], color = [yi > 0 ? :red : :blue for yi in y], label = "")
Plots.yaxis!(p, (-2, 4.5))
Plots.plot!(p, svm_x, svm_y, label = "loss = $(round(loss, digits=2))", width=3)

# # Experiments
# Now that we've solved the SVM, we can compute the sensitivity of optimal values -- the separating hyperplane in our case -- with respect to perturbations of the problem data -- the data points -- using DiffOpt. For illustration, we've explored two questions:

# - How does a change in labels of the data points (`y=1` to `y=-1`, and vice versa) affect the position of the hyperplane? This is achieved by finding the gradient of `w`, `b` with respect to `y[i]`, the classification label of the ith data point.
# - How does a change in coordinates of the data points, `X`, affects the position of the hyperplane? This is achieved by finding gradient of `w`, `b` with respect to `X[i]`, 2D coordinates of the data points.

# Note that finding the optimal SVM can be modelled as a conic optimization problem:

# ```math
# \begin{align*}
# & \min_{x \in \mathbb{R}^n} & c^T x \\
# & \text{s.t.}               & A x + s = b  \\
# &                           & b \in \mathbb{R}^m  \\
# &                           & s \in \mathcal{K}
# \end{align*}
# ```

# where
# ```math
# \begin{align*}
# c &= [l_1 - 1, l_2 -1, ... l_N -1, 0, 0, ... 0 \text{(D+1 times)}] \\\\

# A &=
# \begin{bmatrix}
#  -l_1 &    0 & ... &    0 &            0 & ... & 0 & 0  \\
#     0 & -l_2 & ... &    0 &            0 & ... & 0 & 0  \\
#     : &    : & ... &    : &            0 & ... & 0 & 0  \\
#     0 &    0 & ... & -l_N &            0 & ... & 0 & 0  \\
#     0 &    0 & ... &    0 & -y_1 X_{1,1} & ... & -y_1 X_{1,N} & -y_1  \\
#     0 &    0 & ... &    0 & -y_2 X_{2,1} & ... & -y_1 X_{2,N} & -y_2  \\
#     : &    : & ... &    : &           :  & ... &          :   & :   \\
#     0 &    0 & ... &    0 & -y_N X_{N,1} & ... & -y_N X_{N,N} & -y_N  \\
# \end{bmatrix} \\\\

# b &= [0, 0, ... 0 \text{(N times)}, l_1 - 1, l_2 -1, ... l_N -1] \\\\

# \mathcal{K} &= \text{Set of Nonnegative cones}
# \end{align*}
# ```

# ## Experiment 1: Gradient of hyperplane wrt the data point labels

# Construct perturbations in data point labels `y` without changing the data point coordinates `X`.

∇ = Float64[]
dy = zeros(N);

# begin differentiating
for Xi in 1:N
    dy[Xi] = 1.0  # set

    MOI.set(
        model,
        DiffOpt.ForwardInConstraint(),
        cons,
        MOI.Utilities.vectorize(dy .* MOI.SingleVariable(b)),
    )

    DiffOpt.forward(model)

    dw = MOI.get.(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        w
    )
    db = MOI.get(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        b
    )
    push!(∇, norm(dw) + norm(db))

    dy[Xi] = 0.0  # reset the change made above
end
normalize!(∇);

# Visualize point sensitivities with respect to separating hyperplane. Note that the gradients are normalized.

p2 = Plots.scatter(
    X[:,1], X[:,2],
    color = [yi > 0 ? :red : :blue for yi in y], label = "",
    markersize = ∇ * 20,
)
Plots.yaxis!(p2, (-2, 4.5))
Plots.plot!(p2, svm_x, svm_y, label = "loss = $(round(loss, digits=2))", width=3)

# ## Experiment 2: Gradient of hyperplane wrt the data point coordinates

# Similar to previous example, construct perturbations in data points coordinates `X`.

∇ = Float64[]
dX = zeros(N, D);

# begin differentiating
for Xi in 1:N
    dX[Xi, :] = ones(D)  # set

    for i in 1:D
        MOI.set(
            model,
            DiffOpt.ForwardInConstraint(),
            cons,
            MOI.Utilities.vectorize(dX[:,i] .* MOI.SingleVariable(w[i])),
        )
    end

    DiffOpt.forward(model)

    dw = MOI.get.(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        w
    )
    db = MOI.get(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        b
    )
    push!(∇, norm(dw) + norm(db))

    dX[Xi, :] = zeros(D)  # reset the change made ago
end
normalize!(∇);

# We can visualize point sensitivity with respect to the separating hyperplane. Note that the gradients are normalized.

p3 = Plots.scatter(
    X[:,1], X[:,2],
    color = [yi > 0 ? :red : :blue for yi in y], label = "",
    markersize = ∇ * 20,
)
Plots.yaxis!(p3, (-2, 4.5))
Plots.plot!(p3, svm_x, svm_y, label = "loss = $(round(loss, digits=2))", width=3)


text:
function a3fn(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            n_samples::Int) where T1 <: AbstractFloat
    bins = calculate_area_bins(verts, triangles)
    angles = Vector{Float64}(undef, n_samples)
    @inbounds for i in 1:n_samples
        points = [point_sampling(verts, triangles, bins) for _ in 1:3]
        v1 = points[2] - points[1]
        v2 = points[3] - points[1]
        angles[i] = angle_between_vector(v1, v2)
    end
    return angles
end
function a3fn(mesh::Mesh, n_samples::Int)
    verts, triangles = decompose_vertex_and_face(mesh)
    return a3fn(verts, triangles, n_samples)
end
function d1fn(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            n_samples::Int) where T1 <: AbstractFloat
    bins = calculate_area_bins(verts, triangles)
    centroid_point = centroid(verts)
    distances = Vector{Float64}(undef, n_samples)
    @inbounds for i in 1:n_samples
        point = point_sampling(verts, triangles, bins)
        distances[i] = euclidean(centroid_point, point)
    end
    return distances
end
function d1fn(mesh::Mesh, n_samples::Int)
    verts, triangles = decompose_vertex_and_face(mesh)
    return d1fn(verts, triangles, n_samples)
end
function d2fn(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            n_samples::Int) where T1 <: AbstractFloat
    bins = calculate_area_bins(verts, triangles)
    distances = Vector{Float64}(undef, n_samples)
    @inbounds for i in 1:n_samples
        points = [point_sampling(verts, triangles, bins) for _ in 1:2]
        distances[i] = euclidean(points[1], points[2])
    end
    return distances
end
function d2fn(mesh::Mesh, n_samples::Int)
    verts, triangles = decompose_vertex_and_face(mesh)
    return d2fn(verts, triangles, n_samples)
end
function d3fn(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            n_samples::Int) where T1 <: AbstractFloat
    bins = calculate_area_bins(verts, triangles)
    sqrt_areas = Vector{Float64}(undef, n_samples)
    @inbounds for i in 1:n_samples
        points = [point_sampling(verts, triangles, bins) for _ in 1:3]
        area = triangle_area(points[1], points[2], points[3])
        sqrt_areas[i] = sqrt(area)
    end
    return sqrt_areas
end
function d3fn(mesh::Mesh, n_samples::Int)
    verts, triangles = decompose_vertex_and_face(mesh)
    return d3fn(verts, triangles, n_samples)
end
function d4fn(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            n_samples::Int) where T1 <: AbstractFloat
    bins = calculate_area_bins(verts, triangles)
    cbrt_vols = Vector{Float64}(undef, n_samples)
    @inbounds for i in 1:n_samples
        points = [point_sampling(verts, triangles, bins) for _ in 1:4]
        vol = tetrahedron_volume(points[1], points[2], points[3], points[4])
        cbrt_vols[i] = cbrt(vol)
    end
    return cbrt_vols
end
function d4fn(mesh::Mesh, n_samples::Int)
    verts, triangles = decompose_vertex_and_face(mesh)
    return d4fn(verts, triangles, n_samples)
end


text:
__precompile__()

module ERFA

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("ERFA is not properly installed. Please run Pkg.build(\"ERFA\")")
end

include("erfa_common.jl")
include("deprecated.jl")

include("a.jl")
include("b.jl")
include("c.jl")
include("d.jl")
include("e.jl")
include("f.jl")
include("g.jl")
include("h.jl")
include("i.jl")
include("j.jl")
include("l.jl")
include("n.jl")
include("o.jl")
include("p.jl")
include("r.jl")
include("s.jl")
include("t.jl")
include("u.jl")
include("x.jl")

function ASTROM(pmt, eb::AbstractArray, eh::AbstractArray, em, v::AbstractArray, bm1, bpn::AbstractArray, along, phi, xpl, ypl, sphi, cphi, diurab, eral, refa, refb)
    ASTROM(pmt,
              (eb[1], eb[2], eb[3]),
              (eh[1], eh[2], eh[3]),
              em,
              (v[1], v[2], v[3]),
              bm1,
              (bpn[1], bpn[2], bpn[3], bpn[4], bpn[5], bpn[6], bpn[7], bpn[8], bpn[9]),
              along,
              phi,
              xpl,
              ypl,
              sphi,
              cphi,
              diurab,
              eral,
              refa,
              refb)
end

function LDBODY(bm, dl, pv::AbstractArray)
    LDBODY(bm, dl, (pv...,))
end

end # module


text:
using Test, Pkg
using MPI

if haskey(Pkg.installed(), "CuArrays")
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0
isroot = rank == root

counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
ref = collect(Iterators.flatten([fill(r, counts[r+1]) for r = 0:size-1]))

for T in Base.uniontypes(MPI.MPIDatatype)
    A = isroot ? ArrayType{T}(ref) : ArrayType{T}(undef, 0)

    # Allocating
    B = MPI.Scatterv(A, counts, root, comm)
    @test B isa ArrayType{T}
    @test B == ArrayType{T}(fill(rank, counts[rank+1]))

    # Non Allocating
    B = ArrayType{T}(undef, counts[rank+1])
    MPI.Scatterv!(A, B, counts, root, comm)
    @test B == ArrayType{T}(fill(rank, counts[rank+1]))

    # IN_PLACE
    B = isroot ? copy(A) : ArrayType{T}(undef, counts[rank+1])
    MPI.Scatterv_in_place!(B, counts, root, comm)
    if isroot
        @test B == A
    else
        @test B == ArrayType{T}(fill(rank, counts[rank+1]))
    end
end

MPI.Finalize()
@test MPI.Finalized()


text:
using Test;
using Combinatorics;
# LEARNINGS:
# - for run2: using Combinatorics + powerset(nums) ... is the right logic but WAAAAY too slow

function run(fname)
    differences = Dict()
    # initialize
    differences[1] = 0
    differences[2] = 0
    differences[3] = 0

    f = open(fname, "r")
    lines = readlines(f)
    nums = map((x) -> parse(Int, x), lines)
    sort!(nums)
    println(nums)

    # START edge case (wall has joltage of 0)
    diff = nums[1] - 0
    differences[diff] += 1

    for i in collect(2:length(nums))
        diff = nums[i] - nums[i - 1]
        differences[diff] += 1
    end

    # END edge case (there's always one last difference of 3)
    differences[3] += 1

    ones = differences[1]
    threes = differences[3]
    result = ones * threes
    println("[", fname, "] Ones = ", ones, " Threes = ", threes, " ... result = ", result)
    return result
end

function run2(fname)
    f = open(fname, "r")
    lines = readlines(f)
    nums = map((x) -> parse(Int, x), lines)
    sort!(nums)

    valid_count = 0
    min_val = 0
    max_val = nums[length(nums)] + 3

    nums = [min_val, nums..., max_val]

    # NEW ALGO:
    # - find all jumps of 3.
    # - partition the list into sublists, separated by those jumps.
    all_sets = []
    current_set = [nums[1]]
    for i in collect(2:length(nums))
        diff = nums[i] - nums[i - 1]
        if diff == 3
            push!(all_sets, current_set)
            current_set = [nums[i]]
        else
            push!(current_set, nums[i])
        end
    end
    # add last set
    push!(all_sets, current_set)

    result = 1
    for s in all_sets
        # - run algo on those smaller sublists
        # - multiply results
        vs = valid_subsets(s)
        println("S = ", s, " .. valid subsets = ", length(vs))
        for vvs in vs
            println("\t", vvs)
        end
        result *= length(vs)
    end

    println("[", fname, "] [part II] result = ", result)
    return result
end

function valid_subsets(s)
    if length(s) < 3
        return 1
    end

    out = []
    total = 0
    inner_items = powerset(s[2:length(s) - 1])
    for ii in inner_items
        vs = [first(s), ii..., last(s)]
        if is_valid(vs)
            push!(out, vs)
        end
    end
    return out
end

function is_valid(l)
    for i in 2:length(l)
        if (l[i] - l[i - 1]) > 3
            return false
        end
    end
    return true
end

@test run("10ex1.txt") == 7 * 5
@test run("10ex2.txt") == 22 * 10
run("10.txt")

@test run2("10ex1.txt") == 8
@test run2("10ex2.txt") == 19208
run2("10.txt")


text:
#Do the Proximal Policy Optimization update:
function do_PPO_update(all_traj::Traj_Data, θa::Any, θv::Any, params::Params, opt_para_a, opt_para_v)
    for k = 1: all_traj.n_traj
        ĝ_a = grad_PPO(θa,params,all_traj.traj[k])

        ĝ_a2 = grad_PPO2(θa,params,all_traj.traj[k])

        figure(2)
        clf()
        plot(net_to_vec(ĝ_a))
        plot(net_to_vec(ĝ_a2))
        sleep(1)

        ĝ_v = grad_value(θv,params,all_traj.traj[k])
        Knet.update!(θa, ĝ_a, opt_para_a)
        Knet.update!(θv, ĝ_v, opt_para_v)
    end
    return θa
end

#Constructed to surrogate Objective to be optimized
@everywhere function loss_PPO(θa::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    eps    = params.eps
    σ      = params.σ

    #Init
    L_CLIP  = 0
    s_t = traj.s_t[1,:]

    #Hope: improved peroformance by going backwards...
    for i_t = 1:l_traj
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t,:]
        r_t = traj.r_t[i_t]
        A_t = traj.A_t[i_t]
        p_a_t_old = traj.p_a_t[i_t,:]

        #Only for Segway
        s_mod = s_t
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_t[1:2]
        s_mod[3]     = sin(s_t[3])
        s_mod[4]     = cos(s_t[3])
        s_mod[5:end] = s_t[4:end]
        #-------------------------------=#

        out_t = mlp_nlayer_policy_norm(θa,s_mod,params) .+ 0.5 #Evaluate the Network
        p_a_t = get_π_prob(out_t,σ,a_t)

        ratio = prod(p_a_t ./ p_a_t_old) #Get the Ratio! prod(p_a_t./p_a_t_old)

        #Do the conditioned Clipping! -Hard-CLIP!
        clip_ratio = clamp(ratio, 1-eps, 1+eps)

        #Take the maximum value of both
        loc_L = max(ratio*A_t, clip_ratio*A_t) #Take the max of both!

        L_CLIP += loc_L
    end

    #Return the final Surrogate Objective
    return L_CLIP
end

@everywhere grad_PPO = grad(loss_PPO,1)

#Constructed to surrogate Objective to be optimized
@everywhere function loss_PPO2(θa::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    eps    = params.eps
    σ      = params.σ

    #Init
    L_CLIP  = 0
    s_t = traj.s_t[1,:]

    #Hope: improved peroformance by going backwards...
    for i_t = 1:l_traj
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t,:]
        r_t = traj.r_t[i_t]
        A_t = traj.A_t[i_t]
        p_a_t_old = traj.p_a_t[i_t,:]

        #Only for Segway
        s_mod = s_t
        #=---------------------------------
        s_mod        = zeros(8)
        s_mod[1:2]   = s_t[1:2]
        s_mod[3]     = sin(s_t[3])
        s_mod[4]     = cos(s_t[3])
        s_mod[5:end] = s_t[4:end]
        #-------------------------------=#

        out_t = mlp_nlayer_policy_norm(θa,s_mod,params) #Evaluate the Network
        p_a_t = get_π_prob(out_t,σ,a_t)

        ratio = prod(p_a_t ./ p_a_t_old) #Get the Ratio! prod(p_a_t./p_a_t_old)

        #Do the conditioned Clipping! -Soft-CLIP!
        clip_ratio = 1 + eps*tanh.((ratio-1) / eps)

        #Take the maximum value of both
        loc_L = max(ratio*A_t, clip_ratio*A_t) #Take the max of both!

        L_CLIP += loc_L
    end

    #Return the final Surrogate Objective
    return L_CLIP
end

@everywhere grad_PPO2 = grad(loss_PPO2,1)

#Parallel PPO update
#-------------------------------------------------------------------------------
"""
Function enabling parallel execution of the episodes by setting up the channels
"""
function parallel_PPO(params::Params,θa::Any,θv::Any,all_traj::Traj_Data,opt_para_a, opt_para_v)
        #Take all partial traj and shuffle their order!
        n_traj = all_traj.n_traj
        shuffle_arr = shuffle(1:n_traj)

        active_workers = workers()
        n_workers = length(active_workers)
        ĝ_a = copy(θa)
        #ĝ_v = copy(θv)

        # Create Channels storing results and starting values
        input_channel  = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks and traj
        output_channel = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_gradients, p, params, input_channel, output_channel)
            end
        end

        # Start optimizations
        n_parallel = Int(ceil(n_traj / n_workers))
        prog = Progress(n_traj,1)

        for i_parallel = 1 : n_parallel

            #1. Send the work tasks:
            for i_w = 1 : n_workers
                i_true = ((i_parallel-1)*n_workers + i_w) #get true order
                if i_true <= n_traj
                    i_traj = shuffle_arr[i_true] #shuffle order
                    traj_local = all_traj.traj[i_traj]
                    @async put!(input_channel, (traj_local, θa, θv))
                end
            end

            #2. Take the results
            ĝ_a *= 0
            #ĝ_v *= 0
            for i_w = 1: n_workers
                i_true = ((i_parallel-1) * n_workers + i_w)
                if i_true <= n_traj
                    #Take from Channel
                    gradients = take!(output_channel)
                    #Reorder
                    ĝ_a += gradients[1]
                    #ĝ_v += gradients[2]
                end
            end

            #Update the Networks
            Knet.update!(θa, ĝ_a, opt_para_a)
        end

        # Tell workers that its done
        NaN_mat = zeros(length(params.batch_size))
        fill!(NaN_mat,NaN)
        ending_traj = init_traj(params,params.batch_size)
        ending_traj.r_t = NaN_mat

        for j = 1 : n_workers
                @async put!(input_channel, (ending_traj, θa, θv))
        end

        #Close Channels
        sleep(0.1)
        close(input_channel)
        close(output_channel)

        #Return the updated Networks
        return (θa, θv)
end

# Function running on the worker
@everywhere function distributed_gradients(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel)

    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            #Take the new Input
            input = take!(input_channel)

            #Reorder data:
            traj_local = input[1]
            θa         = input[2]
            θv         = input[3]

            #Check if u reached final
            if (sum(isnan, traj_local.r_t) > 0)
                      #gc()
                      break
            end

            #DO THE WORK:
            #-----------------------------------------------------------------------
            #Calculate the Policy Gradient using PPO
            ĝ_a_loc = grad_PPO(θa,params,traj_local)

            #Calculate the Value Function gradient using TD(1)
            ĝ_v_loc = 0#grad_value(θv,params,traj_local)
            #-----------------------------------------------------------------------

            put!(output_channel, (ĝ_a_loc, ĝ_v_loc))

        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            break
        end
    end

    return nothing
end

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

#Constructed to surrogate Objective to be optimized
@everywhere function loss_PPO_combined(θ::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = traj.len
    eps    = params.eps
    σ      = params.σ

    #Init
    L_CLIP  = 0
    s_t = traj.s_t[1,:]

    #Hope: improved peroformance by going backwards...
    for i_t = 1:l_traj

        #1. Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t,:]
        r_t = traj.r_t[i_t]
        A_t = traj.A_t[i_t]
        R_t = traj.R_t[i_t]
        p_a_t_old = traj.p_a_t[i_t,:]
        val_t_old = traj.val_t[i_t,:]

        #2. Evaluate the Network
        out_t, val_t = mlp_nlayer(θ,s_t,params)

        #3. Get Policy update
        p_a_t = get_π_prob(out_t,σ,a_t) #Calculate Action Probability

        ratio = prod(p_a_t) / prod(p_a_t_old) #Get the Ratio!

        #Do the conditioned Clipping!
        clip_ratio = 1 + eps*tanh.((ratio-1) / eps)

        #Local Actor Error
        loc_AE = max(ratio*A_t, clip_ratio*A_t) #Take the max of both!

        #4. Get Value Update
        loc_VE = (val_t - R_t)^2

        L_CLIP += loc_AE + loc_VE
    end

    #Return the final Surrogate Objective
    return L_CLIP
end

@everywhere grad_PPO_combined = grad(loss_PPO_combined,1)

#Parallel PPO update
#-------------------------------------------------------------------------------
"""
Function enabling parallel execution of the episodes by setting up the channels
"""
function parallel_PPO_combined(params::Params,θ::Any,all_traj::Traj_Data,opt_para)
        n_traj = all_traj.n_traj
        active_workers = workers()
        n_workers = length(active_workers)

        # Create Channels storing results and starting values
        input_channel  = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks and traj
        output_channel = RemoteChannel(()->Channel{Any}(length(workers())+1)) # will contain networks

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_gradients_combined, p, params, input_channel, output_channel)
            end
        end

        # Start optimizations
        n_parallel = Int(ceil(n_traj / n_workers))
        prog = Progress(n_traj,1)

        for i_parallel = 1 : n_parallel

            #1. Send the work tasks:
            for i_w = 1 : n_workers
                i_traj = ((i_parallel-1)*n_workers + i_w)
                if i_traj <= n_traj
                    traj_local = all_traj.traj[i_traj]
                    @async put!(input_channel, (traj_local, θ))
                end
            end

            #2. Take the results
            for i_w = 1: n_workers
                i_traj = ((i_parallel-1)*n_workers + i_w)
                if i_traj <= n_traj
                    #Take from Channel
                    ĝ = take!(output_channel)
                    #Update the Networks
                    Knet.update!(θ, ĝ, opt_para)
                end
            end
        end

        # Tell workers that its done
        NaN_mat = zeros(length(params.batch_size))
        fill!(NaN_mat,NaN)
        ending_traj = init_traj(params,params.batch_size)
        ending_traj.r_t = NaN_mat

        for j = 1 : n_workers
                @async put!(input_channel, (ending_traj, θ))
        end

        #Close Channels
        sleep(0.1)
        close(input_channel)
        close(output_channel)

        #Return the updated Networks
        return θ
end

# Function running on the worker
@everywhere function distributed_gradients_combined(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel)

    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            #Take the new Input
            input = take!(input_channel)

            #Reorder data:
            traj_local = input[1]
            θ          = input[2]

            #Check if u reached final
            if (sum(isnan, traj_local.r_t) > 0)
                      break
            end

            #DO THE WORK:
            #-----------------------------------------------------------------------
            #Calculate the Policy Gradient using PPO
            ĝ = grad_PPO_combined(θ,params,traj_local)
            #-----------------------------------------------------------------------

            put!(output_channel, ĝ)

        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            break
        end
    end

    return nothing
end

#Parallel PPO update
#-------------------------------------------------------------------------------
"""
Function enabling parallel execution of the episodes by setting up the channels
"""
function parallel_PPO2(params::Params,θa::Any,θv::Any,all_traj::Traj_Data,opt_para_a, opt_para_v)
        #Take all partial traj and shuffle their order!
        n_traj = all_traj.n_traj
        shuffle_arr = shuffle(1:n_traj)

        active_workers = workers()
        n_workers = length(active_workers)
        ĝ_a = copy(θa)
        #ĝ_v = copy(θv)

        # Create Channels storing results and starting values
        input_channel  = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks and traj
        output_channel = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain networks

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_gradients2, p, params, input_channel, output_channel)
            end
        end

        # Start optimizations
        n_parallel = Int(ceil(n_traj / n_workers))
        prog = Progress(n_traj,1)

        for i_parallel = 1 : n_parallel

            #1. Send the work tasks:
            for i_w = 1 : n_workers
                i_true = ((i_parallel-1)*n_workers + i_w) #get true order
                if i_true <= n_traj
                    i_traj = shuffle_arr[i_true] #shuffle order
                    traj_local = all_traj.traj[i_traj]
                    @async put!(input_channel, (traj_local, θa, θv))
                end
            end

            #2. Take the results
            ĝ_a *= 0
            #ĝ_v *= 0
            for i_w = 1: n_workers
                i_true = ((i_parallel-1) * n_workers + i_w)
                if i_true <= n_traj
                    #Take from Channel
                    gradients = take!(output_channel)
                    #Reorder
                    ĝ_a += gradients[1]
                    #ĝ_v += gradients[2]
                end
            end

            #Update the Networks
            Knet.update!(θa, ĝ_a, opt_para_a)
        end

        # Tell workers that its done
        NaN_mat = zeros(length(params.batch_size))
        fill!(NaN_mat,NaN)
        ending_traj = init_traj(params,params.batch_size)
        ending_traj.r_t = NaN_mat

        for j = 1 : n_workers
                @async put!(input_channel, (ending_traj, θa, θv))
        end

        #Close Channels
        sleep(0.1)
        close(input_channel)
        close(output_channel)

        #Return the updated Networks
        return (θa, θv)
end

# Function running on the worker
@everywhere function distributed_gradients2(params::Params, input_channel::RemoteChannel, output_channel::RemoteChannel)

    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            #Take the new Input
            input = take!(input_channel)

            #Reorder data:
            traj_local = input[1]
            θa         = input[2]
            θv         = input[3]

            #Check if u reached final
            if (sum(isnan, traj_local.r_t) > 0)
                      #gc()
                      break
            end

            #DO THE WORK:
            #-----------------------------------------------------------------------
            #Calculate the Policy Gradient using PPO
            ĝ_a_loc = grad_PPO2(θa,params,traj_local)

            #Calculate the Value Function gradient using TD(1)
            ĝ_v_loc = 0#grad_value(θv,params,traj_local)
            #-----------------------------------------------------------------------

            put!(output_channel, (ĝ_a_loc, ĝ_v_loc))

        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            break
        end
    end

    return nothing
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
module GravMag

using LinearAlgebra

# Type declarations
struct Point{T}
    x::T
    y::T
    z::T
end

abstract type Datum end

struct ScalarDatum <: Datum
    p::Point
    d::Float64
    datatype::String
    unit::String
end

ScalarDatum(p, d) = ScalarDatum(p, d, "dT", "nT")
struct VectorDatum <: Datum
	p::Point
	d::Array{Float64, 1}
	datatype::String
	unit::String
end

VectorDatum(p, d) = VectorDatum(p, d, "B", "nT")

include("gbox.jl")
include("mprism.jl")
include("dipole.jl")
include("tools.jl")

const mbox = mprism

export Point, Datum, ScalarDatum, VectorDatum
export gbox, mbox
export mprism, dipole, dircos, mprism_, mprism!

end # module

max_stars_repo_name: ruboerner/GravMag

max_issues_repo_name: ruboerner/GravMag

max_forks_repo_name: ruboerner/GravMag