text:
@doc """
    merge!(S::SeisData, U::SeisData[, prune_only=true])

Merge channels of two SeisData structures.

    merge!(S::SeisData[, prune_only=true])

"Flatten" a SeisData structure by merging channels with identical properties.

If `prune_only=true`, the only action taken is deletion of empty and duplicate
channels; `merge!(S, U, prune_only=true)` is identical to an in-place `S+U`.
""" merge!
function merge!(S::Y; v::Int64=KW.v, purge_only::Bool=false) where Y<:GphysData
  # Required preprocessing
  prune!(S)

  # Initialize variables
  ID    = getfield(S, :id)
  NAME  = getfield(S, :name)
  LOC   = getfield(S, :loc)
  FS    = getfield(S, :fs)
  GAIN  = getfield(S, :gain)
  RESP  = getfield(S, :resp)
  UNITS = getfield(S, :units)
  SRC   = getfield(S, :src)
  MISC  = getfield(S, :misc)
  NOTES = getfield(S, :notes)
  T     = getfield(S, :t)
  X     = getfield(S, :x)

  UID = unique(getfield(S, :id))
  cnt = 0
  N_warn = 0
  note_head = string(SeisIO.timestamp(), ": ")
  to_delete = Array{Int64,1}(undef, 0)
  while cnt < length(UID)
    cnt = cnt + 1
    id = getindex(UID, cnt)
    GRP = findall(S.id.==id)
    SUBGRPS = get_subgroups(LOC, FS, RESP, UNITS, GRP)

    for subgrp in SUBGRPS
      N = length(subgrp)
      if N > 1
        # Check for duplicates
        sort!(subgrp)
        u = trues(N)
        while N > 1
          t1 = getindex(T, getindex(subgrp, N))
          x1 = getindex(X, getindex(subgrp, N))

          # if a channel is already known to be a duplicate, skip it
          if getindex(u, N) == true
            j = N
            while j > 1
              j = j-1
              t2 = getindex(T, getindex(subgrp, j))
              x2 = getindex(X, getindex(subgrp, j))
              if t1 == t2 && x1 == x2
                setindex!(u, false, j)
              end
            end
          end

          # flag duplicates for deletion
          append!(to_delete, subgrp[u.==false])
          N = N-1
        end

        subgrp = subgrp[u.==true]
      end
      purge_only == true && continue

      N     = length(subgrp)
      i1    = getindex(subgrp, 1)
      fs    = getindex(FS, i1)
      if fs == 0.0
        Ω = merge_non_ts!(S, subgrp)
        append!(to_delete, subgrp[subgrp.!=Ω])
        continue
      end
      Δ     = round(Int64, sμ/fs)
      w_tmp = Array{Array{Int64,2},1}(undef, N)
      te    = Array{Int64,1}(undef, N)
      for i = 1:N
        m     = getindex(subgrp, i)
        w_m   = t_win(getindex(T, m), Δ)
        n_w   = size(w_m, 1)

        # Store: w_start, w_end, channel_number, window_number, x_start, x_end
        w = hcat(w_m, Array{Int64,2}(undef, n_w, 4))
        j  = 0
        ws = 0
        we = 0
        while j < n_w
          j       = j + 1
          we      = max(we, getindex(w_m, j, 2))
          setindex!(w, m, j, 3)                     # channel number
          setindex!(w, j, j, 4)                     # window number
          setindex!(w, ws+1, j, 5)                  # x_start
          ws      = ws + div(w[j,2]-w[j,1], Δ)+1    # x_end
          setindex!(w, ws, j, 6)
        end
        setindex!(te, we, i)
        setindex!(w_tmp, w, i)
      end
      W   = vcat(w_tmp...)
      ii  = sortperm(W[:,2], rev=true)
      W   = W[ii,:]

      # ====================================
      # Let Ω be the channel number in subgrp with the last end time
      Ω     = subgrp[argmax(te)]
      rest  = subgrp[subgrp.!=Ω]
      # ====================================

      # GAIN, MISC, NAME, NOTES, SRC ==========================================
      gain  = getindex(GAIN, Ω)
      notes = getindex(NOTES, Ω)
      misc  = getindex(MISC, Ω)
      for i in rest
        scalefac = gain / getindex(GAIN, i)
        if scalefac != 1.0
          rmul!(getindex(X, i), scalefac)
        end
        if getindex(SRC, i) != getindex(SRC, Ω)
          push!(notes, string(note_head, "+src: ", getindex(SRC, i)))
        end
        if getindex(NAME, i) != getindex(NAME, Ω)
          push!(notes, string(note_head, "alternate name: ", getindex(NAME, i)))
        end
        append!(notes, getindex(NOTES, i))
        merge!(misc, getindex(MISC, i))
      end
      sort!(notes)

      # EventTraceData extra fields ===========================================
      if Y == EventTraceData
        pha   = PhaseCat()
        az    = getindex(getfield(S, :az), Ω)
        baz   = getindex(getfield(S, :baz), Ω)
        dist  = getindex(getfield(S, :dist), Ω)
        for i in rest
          if az == 0.0
            θ = getindex(getfield(S, :az), i)
            (θ != 0.0) && (az = θ)
          end
          if baz == 0.0
            β = getindex(getfield(S, :baz), i)
            (β != 0.0) && (baz = β)
          end
          if dist == 0.0
            # Δ is already in use, so...
            d = getindex(getfield(S, :dist), i)
            (d != 0.0) && (dist = d)
          end
          merge!(pha, getindex(getfield(S, :pha), i))
        end
        # This guarantees that the phase catalog of Ω overwrites others
        merge!(pha, getindex(getfield(S, :pha), Ω))
        setindex!(getfield(S, :az),     az, Ω)
        setindex!(getfield(S, :baz),   baz, Ω)
        setindex!(getfield(S, :dist), dist, Ω)
        setindex!(getfield(S, :pha),   pha, Ω)
      end

      # T,X ===================================================================
      nW = size(W,1)
      (src, dest) = get_next_pair(W)
      while (src, dest) != (zeros(Int64, 7), zeros(Int64, 7))
        ts_i = src[1];  te_i = src[2];  p = src[3];  p_i = src[4];  os_p = src[5];  W_p = src[7]
        ts_j = dest[1]; te_j = dest[2]; q = dest[3]; q_i = dest[4]; os_q = dest[5]; W_q = dest[7]
        ts_max = max(ts_i, ts_j)
        te_min = min(te_i, te_j)
        nov = 1 + div(te_min - ts_max, Δ)
        Xq = getindex(X, q)

        # (1) determine the times and indices of overlap within each pair
        # a. determine sample times of overlap
        Ti = collect(ts_max:Δ:te_min)
        Tj = deepcopy(Ti)

        # b. get sample indices within each overlap window
        # i
        xsi_i = div(ts_max - ts_i, Δ) + os_p
        xei_i = xsi_i + nov - 1
        # j
        xsi_j = div(ts_max - ts_j, Δ) + os_q
        xei_j = xsi_j + nov - 1

        # (2) Extract sample windows
        Xi = getindex(getindex(X, p), xsi_i:xei_i)
        Xj = getindex(getindex(X, q), xsi_j:xei_j)
        lxp = length(getindex(X, p))
        lxq = length(getindex(X, q))

        # ================================================================
        # check for duplicate windows
        if (ts_i == ts_j) && (te_i == te_j) && (Xi == Xj)
          # delete time window
          W = W[setdiff(1:end, W_p), :]
        else
          # Check for misalignment:
          τ, χ, do_xtmerge, δj = check_alignment(Ti, Tj, Xi, Xj, Δ)
          if do_xtmerge
            xtmerge!(τ, χ, div(Δ,2))
          end
          if δj != 0
            xsi_i += δj
            xei_j -= δj
          end

          # (3) Merge X,T into S[q]
          deleteat!(Xq, xsi_j:xei_j)
          if xsi_j == 1
            prepend!(Xq, χ)
          else
            splice!(Xq, xsi_j:xsi_j-1, χ)
          end

          # (4) Adjust start, end indices of windows ≥ q_i in q
          # structure: w_start, w_end, channel_number, window_number, x_start, x_end
          nxq = length(Xq) - lxq
          i = 0
          while i < nW
            i += 1
            if W[i, 3] == q && W[i, 4] ≥ q_i
              W[i, 1] += nxq*Δ
              W[i, 2] += nxq*Δ
              W[i, 5] += nxq
              W[i, 6] += nxq
            end
          end

          #= if xsi_i ≤ os_p (which is always true, at this point in
          the control flow), we decrease W[W_p, 1:2] =#
          nxp = xei_i-xsi_i+1
          W[W_p, 1] -= δj*Δ
          W[W_p, 2] -= (nxp + δj)*Δ

          #= Control for when window P is emptied; the above two statements
          make this possible =#
          if (W[W_p, 2] < W[W_p, 1])
            W = W[setdiff(1:end, W_p), :]
          else
            W[W_p, 6] -= nxp
          end
        end
        # Sort by end time, to ensure we pick the window with latest end next
        k = sortperm(W[:,2], rev=true)
        W = W[k,:]
        nW = size(W,1)
        (src, dest) = get_next_pair(W)
      end
      # Repeat until no further merges are possible
      kk = sortperm(W[:,1])
      W = W[kk,:]

      #= At this point, we have nothing left that can be merged. So we're going
      to arrange T[subgrp] and X[subgrp] in windows using t_win =#

      n = size(W, 1)
      nx = broadcast(+, getindex(W, :, 6).-getindex(W, :,5), 1)
      X_Ω = Array{eltype(X[Ω]),1}(undef, sum(nx))
      xi = 1
      i = 0
      while i < n
        i   = i + 1
        p   = getindex(W, i, 3)
        lx  = getindex(nx, i)
        copyto!(X_Ω, xi, getindex(X, p), getindex(W, i, 5), lx)
        xi  = xi + lx
      end

      # Shrink W and eliminate windows with no actual gap between them
      m = trues(n)
      while n > 1
        if W[n-1, 2] + Δ == W[n,1]
          W[n-1, 2] = W[n,2]
          m[n] = false
        end
        n = n - 1
      end
      W = W[m, [1,2]]
      setindex!(T, w_time(W, Δ), Ω)
      setindex!(X, X_Ω, Ω)
      note!(S, Ω, string("merge!, combined channels ",
                          replace(repr(subgrp), ","=>""),
                          " (N = ", N, ") into :t, :x"))
      append!(to_delete, rest)
    end
    # Done with SUBGRPS
  end
  if N_warn > 0
    @info(string("Total warnings: ", N_warn))
  end
  deleteat!(S, to_delete)
  sort!(S)
  return nothing
end

# The following must be the same (or unset) to merge:
# :loc
# :resp
# :fs         cannot be unset; not a nullable array
# :units

# The following are easily dealt with:
# :gain       can be translated with no trouble
# :misc       can call merge! method for dictionaries
# :notes      append and sort
# :name       not important, can log extra names to :notes
# :src        use most recent, log extras to :notes


text:
#=
This script benchmarks the call to a length function as a condition bound
of a for loop verse declaring the iteration length prior to the for loop.

More effiecient to use the length function? Atleast that is what my machine is
suggesting.
=#
using BenchmarkTools

function forl(X)
    s = 0
    for i in 1:length(X)
        ((i % 2) > 0) ? (s+=1) : (s-=1)
    end
end

function ford(X)
    s = 0
    N = length(X)
    for i in 1:N
        ((i % 2) > 0) ? (s+=1) : (s-=1)
    end
end

x = rand(10^9)

bench_d = @benchmark $ford($x)
bench_l = @benchmark $forl($x)

println("\nDeclared ----------------------")
display(bench_d)

println("\n\nlength() Call ----------------------")
display(bench_l)

println()

ext: jl

max_stars_repo_name: danphenderson/benchmarks

max_issues_repo_name: danphenderson/benchmarks

max_forks_repo_name: danphenderson/benchmarks


text:
function mol_volume_impl(mt::SingleΦT,model::SatLiquidModel,t)
    return mol_volume_impl(QuickStates.sat_t(),model,t)
end

function mol_volume_impl(mt::MultiΦT,model::SatLiquidModel,t,x)
    return mol_volume_impl(QuickStates.sat_tx(),model,t,x)
end

function mol_volume(model::SatLiquidModel,st::ThermodynamicState,unit=u"m^3/mol",mw=nothing)
    return mol_volume(state_type(st),model,st,unit)
end

function mass_volume(model::SatLiquidModel,st::ThermodynamicState,unit=u"m^3/kg",mw=nothing)
    return mol_volume(state_type(st),model,st,unit)
end

function total_volume(model::SatLiquidModel,st::ThermodynamicState,unit=u"m^3",mw=nothing)
    return total_volume(state_type(st),model,st,unit)
end

function mol_density(model::SatLiquidModel,st::ThermodynamicState,unit=u"mol/m^3",mw=nothing)
    return mol_density(state_type(st),model,st,unit)
end

function mass_density(model::SatLiquidModel,st::ThermodynamicState,unit=u"kg/m^3",mw=nothing)
    return mass_density(state_type(st),model,st,unit)
end

function mol_volume(mt::SingleSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    res = mol_volume_impl(mt,model,t)
    return convert_unit(u"m^3/mol",unit,res)
end

function mol_volume(mt::MultiSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st,nothing,mw)
    res = mol_volume_impl(mt,model,t,x)
    return convert_unit(u"m^3/mol",unit,res)
end

function mass_volume(mt::SingleSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    res = mol_volume_impl(mt,model,t)
    return convert_unit(u"m^3/kg",unit,res)
end

function mass_volume(mt::MultiSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st,nothing,mw)
    res = mol_volume_impl(mt,model,t,x)/molar_mass(FromState(),st,u"kg/mol",mw)
    return convert_unit(u"m^3/kg",unit,)
end

function total_volume(mt::SingleSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    res = mol_volume_impl(mt,model,t)
    return convert_unit(u"m^3",unit,res)
end

function total_volume(mt::MultiSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st,nothing,mw)
    res = mol_volume_impl(mt,model,t,x)*moles(FromState(),st,u"mol",mw)
    return convert_unit(u"m^3",unit,)
end

function mol_density(mt::SingleSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    res = mol_volume_impl(mt,model,t)
    return convert_unit(u"mol/m^3",unit,one(res)/res)
end

function mol_density(mt::MultiSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st,nothing,mw)
    res = mol_volume_impl(mt,model,t,x)
    return convert_unit(u"mol/m^3",unit,one(res)/res)
end

function mass_density(mt::SingleSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    res = mol_volume_impl(mt,model,t)
    return convert_unit(u"kg/m^3",unit,one(res)/res)
end

function mass_density(mt::MultiSatT,model::SatLiquidModel,st::ThermodynamicState,unit,mw)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st,nothing,mw)
    res = mol_volume_impl(mt,model,t,x)/molar_mass(FromState(),st,u"kg/mol",mw)
    return convert_unit(u"kg/m^3",unit,one(res)/res)
end


text:
"""
    PencilArray(pencil::Pencil, data::AbstractArray{T,N})

Create array wrapper with pencil decomposition information.

The array dimensions and element type must be consistent with those of the given
pencil.

!!! note "Index permutations"

    If the `Pencil` has an associated index permutation, then `data` must have
    its dimensions permuted accordingly (in *memory* order).

    Unlike `data`, the resulting `PencilArray` should be accessed with
    unpermuted indices (in *logical* order).

    ##### Example

    Suppose `pencil` has local dimensions `(10, 20, 30)` before permutation, and
    has an asociated permutation `(2, 3, 1)`.
    Then:
    ```julia
    data = zeros(20, 30, 10)       # parent array (with dimensions in memory order)

    u = PencilArray(pencil, data)  # wrapper with dimensions (10, 20, 30)
    @assert size(u) === (10, 20, 30)

    u[15, 25, 5]          # BoundsError (15 > 10 and 25 > 20)
    u[5, 15, 25]          # correct
    parent(u)[15, 25, 5]  # correct

    ```

!!! note "Extra dimensions"

    The data array can have one or more extra dimensions to the right (slow
    indices), which are not affected by index permutations.

    ##### Example

    ```julia
    dims = (20, 30, 10)
    PencilArray(pencil, zeros(dims...))        # works (scalar)
    PencilArray(pencil, zeros(dims..., 3))     # works (3-component vector)
    PencilArray(pencil, zeros(dims..., 4, 3))  # works (4×3 tensor)
    PencilArray(pencil, zeros(3, dims...))     # fails
    ```

---

    PencilArray{T}(undef, pencil::Pencil, [extra_dims...])

Allocate an uninitialised `PencilArray` that can hold data in the local pencil.

Extra dimensions, for instance representing vector components, can be specified.
These dimensions are added to the rightmost (slowest) indices of the resulting
array.

# Example
Suppose `pencil` has local dimensions `(20, 10, 30)`. Then:
```julia
PencilArray{Float64}(undef, pencil)        # array dimensions are (20, 10, 30)
PencilArray{Float64}(undef, pencil, 4, 3)  # array dimensions are (20, 10, 30, 4, 3)
```

More examples:

```jldoctest
julia> pen = Pencil((20, 10, 12), MPI.COMM_WORLD);

julia> u = PencilArray{Float64}(undef, pen);

julia> summary(u)
"20×10×12 PencilArray{Float64, 3}(::Pencil{3, 2, NoPermutation})"

julia> PencilArray{Float64}(undef, pen, 4, 3) |> summary
"20×10×12×4×3 PencilArray{Float64, 5}(::Pencil{3, 2, NoPermutation})"

```
"""
struct PencilArray{
        T,
        N,
        A <: AbstractArray{T,N},
        Np,  # number of "spatial" dimensions (i.e. dimensions of the Pencil)
        E,   # number of "extra" dimensions (= N - Np)
        P <: Pencil,
    } <: AbstractArray{T,N}
    pencil   :: P
    data     :: A
    space_dims :: Dims{Np}  # spatial dimensions in *logical* order
    extra_dims :: Dims{E}

    function PencilArray(pencil::Pencil{Np, Mp} where {Np, Mp},
                         data::AbstractArray{T, N}) where {T, N}
        P = typeof(pencil)
        A = typeof(data)
        Np = ndims(pencil)
        E = N - Np
        size_data = size(data)

        geom_dims = ntuple(n -> size_data[n], Np)  # = size_data[1:Np]
        extra_dims = ntuple(n -> size_data[Np + n], E)  # = size_data[Np+1:N]

        dims_local = size_local(pencil, MemoryOrder())

        if geom_dims !== dims_local
            throw(DimensionMismatch(
                "array has incorrect dimensions: $(size_data). " *
                "Local dimensions of pencil: $(dims_local)."))
        end

        space_dims = permutation(pencil) \ geom_dims  # undo permutation

        new{T, N, A, Np, E, P}(pencil, data, space_dims, extra_dims)
    end
end

function PencilArray{T}(init, pencil::Pencil, extra_dims::Vararg{Integer}) where {T}
    dims = (size_local(pencil, MemoryOrder())..., extra_dims...)
    PencilArray(pencil, Array{T}(init, dims))
end

pencil_type(::Type{PencilArray{T,N,A,M,E,P}}) where {T,N,A,M,E,P} = P

# This is called by `summary`.
function Base.showarg(io::IO, u::PencilArray, toplevel)
    toplevel || print(io, "::")
    print(io, nameof(typeof(u)), '{', eltype(u), ", ", ndims(u), '}')
    if toplevel
        print(io, '(')
        Base.showarg(io, pencil(u), false)
        print(io, ')')
    end
    nothing
end

"""
    PencilArrayCollection

`UnionAll` type describing a collection of [`PencilArray`](@ref)s.

Such a collection can be a tuple or an array of `PencilArray`s.

Collections are **by assumption** homogeneous: each array has the same
properties, and in particular, is associated to the same [`Pencil`](@ref)
configuration.

For convenience, certain operations defined for `PencilArray` are also defined
for `PencilArrayCollection`, and return the same value as for a single
`PencilArray`.
Some examples are [`pencil`](@ref), [`range_local`](@ref) and
[`get_comm`](@ref).

Also note that functions from `Base`, such as `size`, `ndims` and `eltype`, are **not**
overloaded for `PencilArrayCollection`, since they already have a definition
for tuples and arrays (and redefining them would be type piracy...).
"""
const PencilArrayCollection =
    Union{Tuple{Vararg{A}}, AbstractArray{A}} where {A <: PencilArray}

collection_size(x::Tuple{Vararg{<:PencilArray}}) = (length(x), )
collection_size(x::AbstractArray{<:PencilArray}) = size(x)
collection_size(::PencilArray) = ()

# This is convenient for iterating over one or more PencilArrays.
# A single PencilArray is treated as a "collection" of one array.
collection(x::PencilArrayCollection) = x
collection(x::PencilArray) = (x, )

const MaybePencilArrayCollection = Union{PencilArray, PencilArrayCollection}

function _apply(f::Function, x::PencilArrayCollection, args...; kwargs...)
    a = first(x)
    if !all(b -> pencil(a) === pencil(b), x)
        throw(ArgumentError("PencilArrayCollection is not homogeneous"))
    end
    f(a, args...; kwargs...)
end

"""
    size(x::PencilArray)

Return local dimensions of a `PencilArray` in logical order.

Same as `size_local(x, LogicalOrder())` (see [`size_local`](@ref)).
"""
Base.size(x::PencilArray) = (x.space_dims..., extra_dims(x)...)

"""
    size_local(x::PencilArray, [order = LogicalOrder()])
    size_local(x::PencilArrayCollection, [order = LogicalOrder()])

Local dimensions of the data held by the `PencilArray`.

If `order = LogicalOrder()`, this is the same as `size(x)`.

See also [`size_local(::Pencil)`](@ref).
"""
size_local(x::MaybePencilArrayCollection, args...; kwargs...) =
    (size_local(pencil(x), args...; kwargs...)..., extra_dims(x)...)

Base.axes(x::PencilArray) = permutation(x) \ axes(parent(x))

"""
    similar(x::PencilArray, [element_type=eltype(x)], [dims])

Returns an array similar to `x`.

The actual type of the returned array depends on whether `dims` is passed:

- if `dims` is *not* passed, then a `PencilArray` of same dimensions of `x` is
  returned.

- otherwise, an array similar to that wrapped by `x` (typically a regular
  `Array`) is returned, with the chosen dimensions.

# Examples

```jldoctest
julia> pen = Pencil((20, 10, 12), MPI.COMM_WORLD);

julia> u = PencilArray{Float64}(undef, pen);

julia> similar(u) |> summary
"20×10×12 PencilArray{Float64, 3}(::Pencil{3, 2, NoPermutation})"

julia> similar(u, ComplexF32) |> summary
"20×10×12 PencilArray{ComplexF32, 3}(::Pencil{3, 2, NoPermutation})"

julia> similar(u, (4, 3, 8)) |> summary
"4×3×8 Array{Float64, 3}"

julia> similar(u, (4, 3)) |> summary
"4×3 Matrix{Float64}"

julia> similar(u, ComplexF32) |> summary
"20×10×12 PencilArray{ComplexF32, 3}(::Pencil{3, 2, NoPermutation})"

julia> similar(u, ComplexF32, (4, 3)) |> summary
"4×3 Matrix{ComplexF32}"
```
"""
function Base.similar(x::PencilArray, ::Type{S}) where {S}
    dims_perm = permutation(x) * size(x)
    PencilArray(x.pencil, similar(x.data, S, dims_perm))
end

Base.similar(x::PencilArray, ::Type{S}, dims::Dims) where {S} =
    similar(parent(x), S, dims)

# Use same index style as the parent array.
Base.IndexStyle(::Type{<:PencilArray{T,N,A}} where {T,N}) where {A} =
    IndexStyle(A)

# Overload Base._sub2ind for converting from Cartesian to linear index.
@inline function Base._sub2ind(x::PencilArray, I...)
    # _sub2ind(axes(x), I...)  <- default implementation for AbstractArray
    J = permutation(x) * I
    Base._sub2ind(parent(x), J...)
end

# Linear indexing
@propagate_inbounds @inline Base.getindex(x::PencilArray, i::Integer) =
    x.data[i]

@propagate_inbounds @inline Base.setindex!(x::PencilArray, v, i::Integer) =
    x.data[i] = v

# Cartesian indexing: assume input indices are unpermuted, and permute them.
# (This is similar to the implementation of PermutedDimsArray.)
@propagate_inbounds @inline Base.getindex(
        x::PencilArray{T,N}, I::Vararg{Int,N}) where {T,N} =
    x.data[_genperm(x, I)...]

@propagate_inbounds @inline Base.setindex!(
        x::PencilArray{T,N}, v, I::Vararg{Int,N}) where {T,N} =
    x.data[_genperm(x, I)...] = v

@inline function _genperm(x::PencilArray{T,N}, I::NTuple{N,Int}) where {T,N}
    # Split "spatial" and "extra" indices.
    M = ndims_space(x)
    E = ndims_extra(x)
    @assert M + E === N
    J = ntuple(n -> I[n], Val(M))
    K = ntuple(n -> I[M + n], Val(E))
    perm = permutation(x)
    ((perm * J)..., K...)
end

@inline _genperm(x::PencilArray, I::CartesianIndex) =
    CartesianIndex(_genperm(x, Tuple(I)))

"""
    pencil(x::PencilArray)

Return decomposition configuration associated to a `PencilArray`.
"""
pencil(x::PencilArray) = x.pencil
pencil(x::PencilArrayCollection) = _apply(pencil, x)

"""
    parent(x::PencilArray)

Return array wrapped by a `PencilArray`.
"""
Base.parent(x::PencilArray) = x.data

# This enables aliasing detection (e.g. using Base.mightalias) on PencilArrays.
Base.dataids(x::PencilArray) = Base.dataids(parent(x))

# This is based on strides(::PermutedDimsArray)
function Base.strides(x::PencilArray)
    s = strides(parent(x))
    permutation(x) * s
end

"""
    pointer(x::PencilArray)

Return pointer to the start of the underlying data.

Use with caution: this may not make a lot of sense if the underlying data is not
contiguous or strided (e.g. if the `PencilArray` is wrapping a non-strided
`SubArray`).
"""
Base.pointer(x::PencilArray) = pointer(parent(x))

"""
    ndims_extra(::Type{<:PencilArray})
    ndims_extra(x::PencilArray)
    ndims_extra(x::PencilArrayCollection)

Number of "extra" dimensions associated to `PencilArray`.

These are the dimensions that are not associated to the domain geometry.
For instance, they may correspond to vector or tensor components.

These dimensions correspond to the rightmost indices of the array.

The total number of dimensions of a `PencilArray` is given by:

    ndims(x) == ndims_space(x) + ndims_extra(x)

"""
ndims_extra(x::MaybePencilArrayCollection) = length(extra_dims(x))
ndims_extra(::Type{<:PencilArray{T,N,A,M,E}}) where {T,N,A,M,E} = E

"""
    ndims_space(x::PencilArray)
    ndims_space(x::PencilArrayCollection)

Number of dimensions associated to the domain geometry.

These dimensions correspond to the leftmost indices of the array.

The total number of dimensions of a `PencilArray` is given by:

    ndims(x) == ndims_space(x) + ndims_extra(x)

"""
ndims_space(x::PencilArray) = ndims(x) - ndims_extra(x)
ndims_space(x::PencilArrayCollection) = _apply(ndims_space, x)

"""
    extra_dims(x::PencilArray)
    extra_dims(x::PencilArrayCollection)

Return tuple with size of "extra" dimensions of `PencilArray`.
"""
extra_dims(x::PencilArray) = x.extra_dims
extra_dims(x::PencilArrayCollection) = _apply(extra_dims, x)

"""
    size_global(x::PencilArray, [order = LogicalOrder()])
    size_global(x::PencilArrayCollection, [order = LogicalOrder()])

Global dimensions associated to the given array.

By default, the logical dimensions of the dataset are returned.

See also [`size_global(::Pencil)`](@ref).
"""
size_global(x::MaybePencilArrayCollection, args...; kw...) =
    (size_global(pencil(x), args...; kw...)..., extra_dims(x)...)

"""
    sizeof_global(x::PencilArray)
    sizeof_global(x::PencilArrayCollection)

Global size of array in bytes.
"""
sizeof_global(x::PencilArray) = prod(size_global(x)) * sizeof(eltype(x))
sizeof_global(x::PencilArrayCollection) = sum(sizeof_global, x)

"""
    range_local(x::PencilArray, [order = LogicalOrder()])
    range_local(x::PencilArrayCollection, [order = LogicalOrder()])

Local data range held by the `PencilArray`.

By default the dimensions are returned in logical order.
"""
range_local(x::MaybePencilArrayCollection, args...; kw...) =
    (range_local(pencil(x), args...; kw...)..., map(Base.OneTo, extra_dims(x))...)

"""
    range_remote(x::PencilArray, coords, [order = LogicalOrder()])
    range_remote(x::PencilArrayCollection, coords, [order = LogicalOrder()])

Get data range held by the `PencilArray` in a given MPI process.

The location of the MPI process in the topology is determined by the `coords`
argument, which can be given as a linear or Cartesian index.

See [`range_remote(::Pencil, ...)`](@ref range_remote(::Pencil, ::Integer,
::LogicalOrder)) variant for details.
"""
range_remote(x::MaybePencilArrayCollection, args...) =
    (range_remote(pencil(x), args...)..., map(Base.OneTo, extra_dims(x))...)

"""
    get_comm(x::PencilArray)
    get_comm(x::PencilArrayCollection)

Get MPI communicator associated to a pencil-distributed array.
"""
get_comm(x::MaybePencilArrayCollection) = get_comm(pencil(x))

"""
    permutation(::Type{<:PencilArray})
    permutation(x::PencilArray)
    permutation(x::PencilArrayCollection)

Get index permutation associated to the given `PencilArray`.

Returns `NoPermutation()` if there is no associated permutation.
"""
function permutation end

function permutation(::Type{A}) where {A <: PencilArray}
    P = pencil_type(A)
    perm = permutation(P)
    E = ndims_extra(A)
    append(perm, Val(E))
end

permutation(x::PencilArray) = permutation(typeof(x))
permutation(x::PencilArrayCollection) = _apply(permutation, x)

"""
    topology(x::PencilArray)
    topology(x::PencilArrayCollection)

Get [`MPITopology`](@ref) associated to a `PencilArray`.
"""
topology(x::MaybePencilArrayCollection) = topology(pencil(x))

Base.zero(x::PencilArray) = fill!(similar(x), zero(eltype(x)))

function Base.fill!(A::PencilArray, x)
    fill!(parent(A), x)
    A
end


text:
using Test #, Random

@testset "Testing solution to Exercise 1" begin

println("All tests for exercise 1 are in the notebook.")

end; # Exercise 1


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# A* shortest-path algorithm

"""
    struct AStar <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [A* search algorithm](http://en.wikipedia.org/wiki/A%2A_search_algorithm).
An optional `heuristic` function may be supplied. If missing, the heuristic is set to
`n -> 0`.

### Implementation Notes
`AStar` supports the following shortest-path functionality:
- non-negative distance matrices / weights
- single destination
"""
struct AStar{F<:Function} <: ShortestPathAlgorithm
    heuristic::F
end

AStar(T::Type=Float64) = AStar(n -> zero(T))

struct AStarResult{T, U<:Integer} <: ShortestPathResult
    path::Vector{U}
    dist::T
end

function reconstruct_path!(total_path, # a vector to be filled with the shortest path
    came_from, # a vector holding the parent of each node in the A* exploration
    end_idx, # the end vertex
    g) # the graph

    curr_idx = end_idx
    while came_from[curr_idx] != curr_idx
        pushfirst!(total_path, came_from[curr_idx])
        curr_idx = came_from[curr_idx]
    end
    push!(total_path, end_idx)
end

function a_star_impl!(g, # the graph
    goal, # the end vertex
    open_set, # an initialized heap containing the active vertices
    closed_set, # an (initialized) color-map to indicate status of vertices
    g_score, # a vector holding g scores for each node
    f_score, # a vector holding f scores for each node
    came_from, # a vector holding the parent of each node in the A* exploration
    distmx,
    heuristic)

    T = eltype(g)
    total_path = Vector{T}()

    @inbounds while !isempty(open_set)
        current = dequeue!(open_set)

        if current == goal
            reconstruct_path!(total_path, came_from, current, g)
            return total_path
        end

        closed_set[current] = true

        for neighbor in outneighbors(g, current)
            closed_set[neighbor] && continue

            tentative_g_score = g_score[current] + distmx[current, neighbor]

            if tentative_g_score < g_score[neighbor]
                g_score[neighbor] = tentative_g_score
                priority = tentative_g_score + heuristic(neighbor)
                enqueue!(open_set, neighbor, priority)
                came_from[neighbor] = current
            end
        end
    end
    return total_path
end

function calc_dist(path, distmx)
    T = eltype(distmx)
    isempty(path) && return typemax(T)
    rest_of_path = copy(path)
    dist = zero(T)
    u = pop!(rest_of_path)
    while !isempty(rest_of_path)
        v = pop!(rest_of_path)
        dist += distmx[u, v]
        u = v
    end
    return dist
end

function shortest_paths(g::AbstractGraph, s::Integer, t::Integer, distmx::AbstractMatrix, alg::AStar)
    T = eltype(distmx)

    # if we do checkbounds here, we can use @inbounds in a_star_impl!
    checkbounds(distmx, Base.OneTo(nv(g)), Base.OneTo(nv(g)))

    open_set = PriorityQueue{Integer, T}()
    enqueue!(open_set, s, 0)

    closed_set = zeros(Bool, nv(g))

    g_score = fill(Inf, nv(g))
    g_score[s] = 0

    f_score = fill(Inf, nv(g))
    f_score[s] = alg.heuristic(s)

    came_from = -ones(Integer, nv(g))
    came_from[s] = s

    path = a_star_impl!(g, t, open_set, closed_set, g_score, f_score, came_from, distmx, alg.heuristic)
    return AStarResult(path, calc_dist(path, distmx))
end

shortest_paths(g::AbstractGraph, s::Integer, t::Integer, alg::AStar) = shortest_paths(g, s, t, weights(g), alg)
paths(s::AStarResult) = [s.path]
paths(s::AStarResult, v::Integer) = throw(ArgumentError("AStar produces at most one path."))
dists(s::AStarResult) = [[s.dist]]
dists(s::AStarResult, v::Integer) = throw(ArgumentError("AStar produces at most one path."))


text:
using Dash
using PlotlyBase
app = dash()
app.layout = html_div() do
    dcc_graph(id = "graph",
        figure = Plot(scatter(x=1:10, y = 1:10))
    ),
    html_button("draw", id = "draw"),
    html_div("", id = "status")
end

callback!(app, Output("graph", "figure"), Output("status", "children"), Input("draw", "n_clicks")) do nclicks
    plot = isnothing(nclicks) ?
           no_update() :
           Plot([scatter(x=1:10, y = 1:10), scatter(x=1:10, y = 1:2:20)])
    status = isnothing(nclicks) ? "first" : "second"
    return (plot, status)
end

run_server(app)


text:
using BinaryBuilder, BinaryProvider

# Collection of sources required to build MySQL
sources = Dict(
    Linux(:i686, :glibc) => "https://dev.mysql.com/get/Downloads/MySQL-8.0/mysql-8.0.15-linux-glibc2.12-i686.tar.xz"=>"b5a18de4e0b8c9209286d887bf187b8e7396e43d4b367870ca870ed95302fc7e",
    Linux(:x86_64, :glibc) => "https://dev.mysql.com/get/Downloads/MySQL-8.0/mysql-8.0.15-linux-glibc2.12-x86_64.tar.xz"=>"f3f1fd7d720883a8a16fe8ca3cb78150ad2f4008d251ce8ac0a2c676e2cf1e1f",
    MacOS() => "https://dev.mysql.com/get/Downloads/MySQL-8.0/mysql-8.0.15-macos10.14-x86_64.tar.gz"=>"f6b1313e89b549947fa774e160a31cf051742110f7f27beadcdc0b4ebea7baa9",
    FreeBSD(:x86_64) => "https://dev.mysql.com/get/Downloads/MySQL-8.0/mysql-8.0.15-freebsd11-x86_64.tar.gz"=>"6099b7fc5444c183d0e1ca8973b32429c58060548c15a2056ed2d81269184a39",
    Windows(:x86_64) => "https://dev.mysql.com/get/Downloads/MySQL-8.0/mysql-8.0.15-winx64.zip"=>"16bf0ad985dd862c81f030438e7907872c2989fbb67548177f5853803d3fb803",
)

product_hashes = Dict()

mkpath(joinpath(pwd(), "products"))

for (platform, (url, hash)) in sources
    println("downloading $url")
    if platform == Windows(:x86_64)
        @assert BinaryProvider.download_verify(url, hash, joinpath(pwd(), "mysql.zip"); force=true, verbose=true)
        println("unzipping $(joinpath(pwd(), "mysql.zip"))")
        success(`unzip $(joinpath(pwd(), "mysql.zip"))`)
        root = joinpath(pwd(), splitext(basename(url))[1])
    else
        @assert BinaryProvider.download_verify_unpack(url, hash, pwd(); force=true, ignore_existence=true, verbose=true)
        root = joinpath(pwd(), splitext(splitext(basename(url))[1])[1])
    end
    println("removing uneeded files...")
    rm(joinpath(root, "bin"); force=true, recursive=true)
    rm(joinpath(root, "docs"); force=true, recursive=true)
    rm(joinpath(root, "share"); force=true, recursive=true)
    filepath = joinpath(pwd(), "products", string("MySQL.$(triplet(platform)).tar.gz"))
    println("packaging $root into $filepath...")
    @show readdir(root)
    BinaryProvider.package(root, filepath)
    product_hashes[triplet(platform)] = open(filepath) do file
        BinaryBuilder.bytes2hex(BinaryBuilder.sha256(file))
    end
    rm(root; force=true, recursive=true)
end

# The products that we will ensure are always built
products(prefix) = Product[
    LibraryProduct(prefix, "libmysqlclient", :libmysql)
]

# If we're only reconstructing a build.jl file on Travis, grab the information and do it
if !haskey(ENV, "TRAVIS_REPO_SLUG") || !haskey(ENV, "TRAVIS_TAG")
    error("Must provide repository name and tag through Travis-style environment variables!")
end
repo_name = ENV["TRAVIS_REPO_SLUG"]
tag_name = ENV["TRAVIS_TAG"]
bin_path = "https://github.com/$(repo_name)/releases/download/$(tag_name)"
dummy_prefix = Prefix(pwd())
print_buildjl(pwd(), "MySQL", v"8.0.15", products(dummy_prefix), product_hashes, bin_path)

println("Writing out the following reconstructed build.jl:")
print_buildjl(stdout, products(dummy_prefix), product_hashes, bin_path)


text:
#
# Convenience simulators for estimator types
#

"""
    SingleStateSimulator(sys::AbstractSystem, obs::AbstractObserver,
                         state::AbstractUncertainState)

Single-state simulator type that mimics the structure of `Estimator` types. An
instance of SingleStateSimulator can be used to simulate a sequence of truth
states and measurements in the presence of system and observer noise.

In general, the best way to create a simulator is through the `make_simulator`
method.
"""
immutable SingleStateSimulator{T,S<:AbstractAbsoluteState{T}}
    sys::AbstractSystem{T}
    obs::AbstractObserver{T}
    true_state::S
end
function SingleStateSimulator(sys::AbstractSystem, obs::AbstractObserver,
                              state::AbstractUncertainState)
    return SingleStateSimulator(sys, obs, sample(state))
end

"""
    MultiStateSimulator(simulator_bank::Vector{SingleStateSimulator{T,S}})

Multi-state simulator type that contains a bank of `SingleStateSimulator`. This
simulator type mimics the structure of a `MultiTargetFilter` and can be used to
simulate a sequence of truth states and measurements in the presence of system
and observer noise for a collection of systems.

In general, the best way to create a simulator is through the `make_simulator`
method.
"""
immutable MultiStateSimulator{T,S<:AbstractAbsoluteState{T}}
    simulator_bank::Vector{SingleStateSimulator{T,S}}
end

"""
    make_simulator(est::Estimator)

Convenience function for simulator creation. Automatically chooses the
appropriate simulator type for the input estimator.
"""
function make_simulator(est::Estimator)
    return SingleStateSimulator(est.sys, est.obs, est.estimate)
end
function make_simulator(mtf::MultiTargetFilter)
    simulator_bank = [make_simulator(filter) for filter in mtf.filter_bank]
    return MultiStateSimulator(simulator_bank)
end

"""
    simulate(state::AbstractState, sys::AbstractSystem, obs::AbstractObserver, t::Real)

Simulate a combined state prediction and observation with noise. Returns two
absolute states. The first represents the "truth" state and the second is a
simulated measurement of that state.
"""
function simulate(state::AbstractState, sys::AbstractSystem,
                  obs::AbstractObserver, t::Real)
   simulated_state = simulate(state, sys, t)
   simulated_measurement = simulate(simulated_state, obs)
   return simulated_state, simulated_measurement
end

"""
    simulate(ssm:SingleStateSimulator, t::Real)

Simulate a state and measurement for a given time using the internal model
and truth state. Advances the simulator's internal truth state to the given
time. Returns a tuple of (true_state, measurement).
"""
function simulate(ssm::SingleStateSimulator, t::Real)
    true_state, measurement = simulate(ssm.true_state, ssm.sys, ssm.obs, t)
    ssm.true_state .= true_state
    return true_state, measurement
end

"""
    simulate(msm:MultiStateSimulator, t::Real)

Simulate a state and measurement for a given time using a bank of internal
models and truth states. Returns a vector of (true_state, measurement) tuples.
"""
function simulate(msm::MultiStateSimulator, t::Real)
    sim_tuples = [simulate(ssm.true_state, ssm.sys, ssm.obs, t) for
                  ssm in msm.simulator_bank]
    return sim_tuples
end


text:
export closure, small_generating_set

################################################################################
#
#  Find identity with respect to operation
#
################################################################################

# It is assumed that the elements have finite order with respect to op.
function find_identity(S, op)
  @assert length(S) > 0
  g = S[1]
  h = g
  while true
    hh = op(h, g)
    if hh == g
      return h
    end
    h = hh
  end
end

################################################################################
#
#  Small generating set
#
################################################################################

function small_generating_set(G, op)
  i = find_identity(G, op)
  return small_generating_set(G, op, i)
end

function small_generating_set(G, op, id)
  orderG = length(G)

  if length(G) == 1
    return G
  end

  function non_trivial_randelem()
    x = rand(G)
    while (x == id)
      x = rand(G)
    end
    return x
  end

  firsttry = 10
  secondtry = 20
  # First try one element
  for i in 1:firsttry
    gen = non_trivial_randelem()
    if length(closure([gen], op, id)) == orderG
      return [gen]
    end
  end

  for i in 1:secondtry
    gens = [non_trivial_randelem(), non_trivial_randelem()]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end

  # Now use that unconditionally log_2(|G|) elements generate G

  b = ceil(Int, log(2, orderG))
  @assert orderG <= 2^b

  j = 0
  while true
    if j > 2^20
      error("Something wrong with generator search")
    end
    j = j + 1
    gens = [non_trivial_randelem() for i in 1:b]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end
end

################################################################################
#
#  Computing closure under group operation
#
################################################################################

# It is assumed that S is nonempty and that the group generated by S under op
# is finite.
function closure(S, op)
  i = find_identity(S, op)
  return closure(S, op, i)
end

function closure(S, op, id)
  if length(S) == 0
    return [id]
  elseif length(S) == 1
    return _closing_under_one_generator(S[1], op, id)
  else
    return _closing_under_generators_dimino(S, op, id)
  end
end

function _closing_under_generators_naive(S, op, id)
  list = push!(copy(S), id)
  stable = false
  while !stable
    stable = true
    for g in list
      for s in S
        m = op(g, s)
        if !(m in list)
          push!(list, m)
          stable = false
        end
      end
    end
  end
  return list
end

function _closing_under_one_generator(x, op, id)
  elements = [x]
  y = x
  while !(y == id)
    y = op(y, x)
    push!(elements, y)
  end
  return elements
end

function _closing_under_generators_dimino(S, op, id)

  t = length(S)
  order = 1
  elements = [id]
  g = S[1]

  while g != id
    order = order +1
    push!(elements, g)
    g = op(g, S[1])
  end

  for i in 2:t
    if !(S[i] in elements)
      previous_order = order
      order = order + 1
      push!(elements, S[i])
      for j in 2:previous_order
        order = order + 1
        push!(elements, op(elements[j], S[i]))
      end

      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          s = S[k]
          elt = op(elements[rep_pos], s)
          if !(elt in elements)
            order = order + 1
            push!(elements, elt)
            for j in 2:previous_order
              order = order + 1
              push!(elements, op(elements[j], elt))
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  return elements
end

################################################################################
#
#  Multiplication table
#
################################################################################

# Construct multiplication table of G under op
function _multiplication_table(G, op)
  l = length(G)
  z = Matrix{Int}(l, l)
  for i in 1:l
    for j in 1:l
      p = op(G[i], G[j])
      for k in 1:l
        if p == G[k]
          z[i, j] = k
          break
        end
      end
    end
  end
  return z
end

################################################################################
#
#  Generic group given by multiplication table
#
################################################################################

# The elements are just 1..n and i * j = G.mult_table[i, j]
mutable struct GrpGen
  identity::Int
  order::Int
  mult_table::Array{Int, 2}
  gens::Array{Int, 1}

  function GrpGen(M::Array{Int, 2})
    z = new()
    z.mult_table = M
    z.identity = find_identity(z)
    return z
  end
end

struct GrpGenElem
  group::GrpGen
  i::Int
end

################################################################################
#
#  Compute generic group from anything
#
################################################################################

function generic_group(G, op)
  return GrpGen(_multiplication_table(G, op))
end

################################################################################
#
#  Parent
#
################################################################################

function parent(g::GrpGenElem)
  return g.group
end

################################################################################
#
#  Construct the ith element
#
################################################################################

function getindex(G::GrpGen, i::Int)
  return GrpGenElem(G, i)
end

function ==(g::GrpGenElem, h::GrpGenElem)
  return parent(g) == parent(h) && g.i == h.i
end

################################################################################
#
#  Order
#
################################################################################

function order(G::GrpGen)
  return size(G.mult_table, 1)
end

################################################################################
#
#  Order of an element
#
################################################################################

function order(g::GrpGenElem)
  k = 2
  h = g * g
  while h != g
    h = g * h
    k = k + 1
  end
  return k - 1
end

################################################################################
#
#  Identity
#
################################################################################

function find_identity(G::GrpGen)
  return _find_identity(G.mult_table)
end

function _find_identity(m::Array{Int, 2})
  return find_identity([1], (i, j) -> m[i, j])
end

function id(G::GrpGen)
  return GrpGenElem(G, G.identity)
end

################################################################################
#
#  Multiplication
#
################################################################################

function *(g::GrpGenElem, h::GrpGenElem)
  G = parent(g)
  return GrpGenElem(G, G.mult_table[g.i, h.i])
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(g::GrpGenElem)
  G = parent(g)
  m = G.mult_table
  i = g.i
  l = size(m, 1)
  ide = id(G)
  for j in 1:l
    if m[i, j] == G.identity
      return GrpGenElem(G, j)
    end
  end
  error("Not possible")
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, G::GrpGen)
  print(io, "Generic group with multiplication table\n")
  println(io, G.mult_table)
end

function Base.show(io::IO, g::GrpGenElem)
  print(io, "Element of generic group ($(g.i))\n")
end

################################################################################
#
#  Is abelian?
#
################################################################################

function isabelian(G::GrpGen)
  return defines_abelian_group(G.mult_table)
end

function defines_abelian_group(m)
  return issymmetric(m)
end

################################################################################
#
#  Isomorphism for specific groups
#
################################################################################

function isisomorphic_to_16T7(G::GrpGen)
  ordershape = [ (1, 1), (2, 3), (4, 12) ]
  # Use algorithm of Tarjahn (or something like that)
  # 16T7 has the following presentation:
  #
  #  $.1^4 = Id($)
  #  $.2^4 = Id($)
  #  $.3^2 = Id($)
  #  $.2^-1 * $.1^2 * $.2^-1 = Id($)
  #  $.1^-1 * $.2^-1 * $.1 * $.2^-1 = Id($)
  #  $.1^-1 * $.3 * $.1 * $.3 = Id($)
  #  $.2^-1 * $.3 * $.2 * $.3 = Id($)

  relations = [[(2, -1), (1, 2), (2, -1)],
               [(1, -1), (2, -1), (1, 1), (2, -1)],
               [(1, -1), (3, 1), (1, 1), (3, 1)],
               [(2, -1), (3, 1), (2, 1), (3, 1)]]

  l = order(G)

  if l != 16
    return false
  end

  if isabelian(G)
    return false
  end

  elements_by_orders = Dict{Int, Array{GrpGenElem, 1}}()

  for i in 1:l
    g = G[i]
    o = order(g)
    if haskey(elements_by_orders, o)
      push!(elements_by_orders[o], g)
    else
      elements_by_orders[o] = [g]
    end
  end

  for (o, no) in ordershape
    if !haskey(elements_by_orders, o)
      return false
    else
      elts = elements_by_orders[o]
      if length(elts) != no
        return false
      end
    end
  end

  ide = id(G)

  img_of_gens = Iterators.product(elements_by_orders[4],
                                  elements_by_orders[4],
                                  elements_by_orders[2])

  for (g1, g2, g3) in img_of_gens
    g2inv = inv(g2)
    g1inv = inv(g1)
    #  $.2^-1 * $.1^2 * $.2^-1 = Id($)
    z = g2inv
    z = z * g1
    z = z * g1
    z = z * g2inv
    if z != ide
      continue
    end
    #  $.1^-1 * $.2^-1 * $.1 * $.2^-1 = Id($)
    z = g1inv
    z = z * g2inv
    z = z * g1
    z = z * g2inv
    if z != ide
      continue
    end
    #  $.1^-1 * $.3 * $.1 * $.3 = Id($)
    z = g1inv
    z = z * g3
    z = z * g1
    z = z * g3
    if z != ide
      continue
    end
    #  $.2^-1 * $.3 * $.2 * $.3 = Id($)
    z = g2inv
    z = z * g3
    z = z * g2
    z = z * g3
    if z != ide
      continue
    end

    cl = closure([g1, g2, g3], *, ide)

    if length(cl) < 16
      continue
    end
    return true
  end
  return false
end

function isisomorphic_to_8T5(G::GrpGen)
  l = order(G)

  if l != 8
    return false
  end

  # First check if it is abelian or not.
  if isabelian(G)
    return false
  end

  # Now G is D_8 (aka D_4 or 8T4) or Q_8
  # But D_8 has 5 elements of order 2, where as Q_8 has only 1 element of order 2.

  z = 0
  for i in 1:l
    g = G[i]
    o = order(g)
    if o == 2
      z += 1
    end
    if z > 1
      return false
    end
  end
  @assert z == 1
  return true
end


text:
include(joinpath(@__DIR__, "affine_arithmetic.jl"))
include(joinpath(@__DIR__, "enumeration.jl"))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


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

# Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015.
create_basis_and_eval(
    :Wedge6,
    "6 node linear prismatic/wedge element",
    [
     (0.0, 0.0, -1.0), # N1
     (1.0, 0.0, -1.0), # N2
     (0.0, 1.0, -1.0), # N3
     (0.0, 0.0,  1.0), # N4
     (1.0, 0.0,  1.0), # N5
     (0.0, 1.0,  1.0), # N6
    ],
    [
     :(1/2 * (1-w) * (1-u-v)), # N1
     :(1/2 * (1-w) * u), # N2
     :(1/2 * (1-w) * v), # N3
     :(1/2 * (1+w) * (1-u-v)), # N4
     :(1/2 * (1+w) * u), # N5
     :(1/2 * (1+w) * v), # N6
    ],
   )

# Basis functions are from ABAQUS theory manual
create_basis_and_eval(
    :Wedge15,
    "15 node quadratic prismatic/wedge element",
    [
     (0.0, 0.0, -1.0), # N1
     (1.0, 0.0, -1.0), # N2
     (0.0, 1.0, -1.0), # N3
     (0.0, 0.0,  1.0), # N4
     (1.0, 0.0,  1.0), # N5
     (0.0, 1.0,  1.0), # N6
     (0.5, 0.0, -1.0), # N7
     (0.5, 0.5, -1.0), # N8
     (0.0, 0.5, -1.0), # N9
     (0.5, 0.0,  1.0), # N10
     (0.5, 0.5,  1.0), # N11
     (0.0, 0.5,  1.0), # N12
     (0.0, 0.0,  0.0), # N13
     (1.0, 0.0,  0.0), # N14
     (0.0, 1.0,  0.0), # N15
    ],
    [
     :(u^2*w^2 - u^2*w + 2*u*v*w^2 - 2*u*v*w - 3*u*w^2/2 + 3*u*w/2 + v^2*w^2 - v^2*w - 3*v*w^2/2 + 3*v*w/2 + w^2/2 - w/2),
     :(u^2*w^2 - u^2*w - u*w^2/2 + u*w/2),
     :(v^2*w^2 - v^2*w - v*w^2/2 + v*w/2),
     :(u^2*w^2 + u^2*w + 2*u*v*w^2 + 2*u*v*w - 3*u*w^2/2 - 3*u*w/2 + v^2*w^2 + v^2*w - 3*v*w^2/2 - 3*v*w/2 + w^2/2 + w/2),
     :(u^2*w^2 + u^2*w - u*w^2/2 - u*w/2),
     :(v^2*w^2 + v^2*w - v*w^2/2 - v*w/2),
     :(-2*u^2*w^2 + 2*u^2*w - 2*u*v*w^2 + 2*u*v*w + 2*u*w^2 - 2*u*w),
     :(2*u*v*w^2 - 2*u*v*w),
     :(-2*u*v*w^2 + 2*u*v*w - 2*v^2*w^2 + 2*v^2*w + 2*v*w^2 - 2*v*w),
     :(2*u^2*w^2 - 2*u^2*w - 4*u^2 + 2*u*v*w^2 - 2*u*v*w - 4*u*v - 2*u*w^2 + 2*u*w + 4*u),
     :(-2*u*v*w^2 + 2*u*v*w + 4*u*v),
     :(2*u*v*w^2 - 2*u*v*w - 4*u*v + 2*v^2*w^2 - 2*v^2*w - 4*v^2 - 2*v*w^2 + 2*v*w + 4*v),
     :(-2*u^2*w^2 + 2*u^2 - 4*u*v*w^2 + 4*u*v + 3*u*w^2 - 3*u - 2*v^2*w^2 + 2*v^2 + 3*v*w^2 - 3*v - w^2 + 1),
     :(-2*u^2*w^2 + 2*u^2 + u*w^2 - u),
     :(-2*v^2*w^2 + 2*v^2 + v*w^2 - v),
    ],
   )


text:
export PlainLogicCircuit, PlainLogicLeafNode, PlainLogicInnerNode,
    PlainLiteralNode, PlainConstantNode, PlainTrueNode, PlainFalseNode,
    Plain⋀Node, Plain⋁Node

#####################
# Plain logic nodes without additional fields
#####################

"""
Root of the plain logic circuit node hierarchy
"""
abstract type PlainLogicCircuit <: LogicCircuit end

"""
A plain logical leaf node
"""
abstract type PlainLogicLeafNode <: PlainLogicCircuit end

"""
A plain logical inner node
"""
abstract type PlainLogicInnerNode <: PlainLogicCircuit end

"""
A plain logical literal leaf node, representing the positive or negative literal of its variable
"""
mutable struct PlainLiteralNode <: PlainLogicLeafNode
    literal::Lit
    data
    counter::UInt32
    PlainLiteralNode(l) = new(l, nothing, 0)
end

"""
A plain logical constant leaf node, representing true or false
"""
abstract type PlainConstantNode <: PlainLogicInnerNode end

"""
Plain constant true node
"""
mutable struct PlainTrueNode <: PlainConstantNode
    data
    counter::UInt32
    PlainTrueNode() = new(nothing, 0)
end

"""
Plain constant false node
"""
mutable struct PlainFalseNode <: PlainConstantNode
    data
    counter::UInt32
    PlainFalseNode() = new(nothing, 0)
end

"""
A plain logical conjunction node (And node)
"""
mutable struct Plain⋀Node <: PlainLogicInnerNode
    children::Vector{PlainLogicCircuit}
    data
    counter::UInt32
    Plain⋀Node(c) = new(c, nothing, 0)
end

"""
A plain logical disjunction node (Or node)
"""
mutable struct Plain⋁Node <: PlainLogicInnerNode
    children::Vector{PlainLogicCircuit}
    data
    counter::UInt32
    Plain⋁Node(c) = new(c, nothing, 0)
end

#####################
# traits
#####################

@inline GateType(::Type{<:PlainLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:PlainConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Plain⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Plain⋁Node}) = ⋁Gate()

#####################
# methods
#####################

"Get the logical constant in a given constant leaf node"
@inline constant(n::PlainTrueNode)::Bool = true
@inline constant(n::PlainFalseNode)::Bool = false

"Get the children of a given inner node"
@inline children(n::PlainLogicInnerNode) = n.children

function conjoin(arguments::Vector{<:PlainLogicCircuit};
                 reuse=nothing)
    @assert length(arguments) > 0
    reuse isa Plain⋀Node && children(reuse) == arguments && return reuse
    return Plain⋀Node(arguments)
end

function disjoin(arguments::Vector{<:PlainLogicCircuit};
                    reuse=nothing)
    @assert length(arguments) > 0
    reuse isa Plain⋁Node && children(reuse) == arguments && return reuse
    return Plain⋁Node(arguments)
end

negate(a::PlainLiteralNode) = compile(PlainLiteralNode, -a.literal)

# claim `PlainLogicCircuit` as the default `LogicCircuit` implementation
compile(::Type{LogicCircuit}, args...) =
    compile(PlainLogicCircuit, args...)

compile(::Type{<:PlainLogicCircuit}, b::Bool) =
    b ? PlainTrueNode() : PlainFalseNode()

compile(::Type{<:PlainLogicCircuit}, l::Lit) =
    PlainLiteralNode(l)

function compile(::Type{<:PlainLogicCircuit}, circuit::LogicCircuit)
    f_con(n) = compile(PlainLogicCircuit, constant(n))
    f_lit(n) = compile(PlainLogicCircuit, literal(n))
    f_a(_, cns) = conjoin(cns)
    f_o(_, cns) = disjoin(cns)
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, PlainLogicCircuit)
end

fully_factorized_circuit(::Type{LogicCircuit}, n::Int) =
    fully_factorized_circuit(PlainLogicCircuit, n)

function fully_factorized_circuit(::Type{<:PlainLogicCircuit}, n::Int)
    ors = map(1:n) do v
        v = Var(v)
        pos = compile(PlainLogicCircuit, var2lit(v))
        neg = compile(PlainLogicCircuit, -var2lit(v))
        pos | neg
    end
    and = conjoin(ors)
    disjoin([and]) # see logistic circuits bias term
end


text:
##############################################################################
##
## Type with stored qr
##
##############################################################################
type DenseQRAllocatedSolver{Tqr <: StridedMatrix, Tu <: AbstractVector} <: AbstractAllocatedSolver
    qr::Tqr
    u::Tu
    function DenseQRAllocatedSolver(qr, u)
        length(u) == size(qr, 1) || throw(DimensionMismatch("u must have length size(J, 1)"))
        new(qr, u)
    end
end

function DenseQRAllocatedSolver{Tqr, Tu <: AbstractVector}(qr::Tqr, u::Tu)
    DenseQRAllocatedSolver{Tqr, Tu}(qr, u)
end

##############################################################################
##
## solve J'J \ J'y by QR
##
##############################################################################

function AbstractAllocatedSolver{Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg, solver::QR)
    return DenseQRAllocatedSolver(similar(nls.J), _zeros(nls.y))
end

function A_ldiv_B!(x::AbstractVector, J::StridedMatrix,  y::AbstractVector, A::DenseQRAllocatedSolver)
    u, qr = A.u, A.qr
    copy!(qr, J)
    copy!(u, y)
    A_ldiv_B!(qrfact!(qr, Val{true}), u)
    for i in 1:length(x)
        x[i] = u[i]
    end
    return x, 1
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y by QR
##
##############################################################################

function AbstractAllocatedSolver{Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt, solver::QR)
    qr = zeros(eltype(nls.J), length(nls.y) + length(nls.x), length(nls.x))
    u = zeros(length(nls.y) + length(nls.x))
    return DenseQRAllocatedSolver(qr, u)
end

function A_ldiv_B!(x::AbstractVector, J::StridedMatrix, y::AbstractVector,
                damp::AbstractVector, A::DenseQRAllocatedSolver, verbose::Bool = false)
    u, qr = A.u, A.qr

    # transform dammp
    length(u) ==  length(y) + length(x) || throw(DimensionMismatch("length(u) should equal length(x) + length(y)"))

    # update qr as |J; diagm(damp)|
    fill!(qr, zero(eltype(qr)))
    for j in 1:size(J, 2)
        for i in 1:size(J, 1)
            qr[i, j] = J[i, j]
        end
    end

    leny = length(y)
    for i in 1:length(damp)
        qr[leny + i, i] = sqrt(damp[i])
    end

    # update u as |J; 0|
    fill!(u, zero(eltype(u)))
    for i in 1:length(y)
        u[i] = y[i]
    end
    if verbose
        @show mean(u)
        sleep(1)
    end

    # solve
    A_ldiv_B!(qrfact!(qr, Val{true}), u)
    if verbose
        @show mean(u)
        sleep(1)
    end
    for i in 1:length(x)
        x[i] = u[i]
    end
    return x, 1
end


text:
using CausalForest
using Test

@testset "causal forest non centered test" begin

    indices,X,T,Y = load_data("causal")
    Xtest = rand(Float64, (1000, size(X, 2)))
    Xtest[:, 1] =  LinRange(0, 1, 1000)
    true_effect = Xtest[:, 1].>0.5

    cf1 = build_forest(false, false, false, Y, T, X, false, -1, 100)
    pred1 = apply_forest(cf1, Xtest)
    comparison1 = isapprox.(true_effect, pred1, atol=0.3)
    @test (sum(comparison1)/1000) > 0.7

    cf2 = build_forest(false, false, true, Y, T, X, false, -1, 100)
    pred2 = apply_forest(cf2, Xtest)
    comparison2 = isapprox.(true_effect, pred2, atol=0.3)
    @test (sum(comparison2)/1000) > 0.7

    cf3 = build_forest(false, true, false, Y, T, X, false, -1, 100)
    pred3 = apply_forest(cf3, Xtest)
    comparison3 = isapprox.(true_effect, pred3, atol=0.3)
    @test (sum(comparison3)/1000) > 0.7

    cf4 = build_forest(false, true, true, Y, T, X, false, -1, 100)
    pred4 = apply_forest(cf4, Xtest)
    comparison4 = isapprox.(true_effect, pred4, atol=0.3)
    @test (sum(comparison4)/1000) > 0.7

end

@testset "causal forest centered test" begin

    indices,X,T,Y = load_data("causal")
    Xtest = rand(Float64, (1000, size(X, 2)))
    Xtest[:, 1] =  LinRange(0, 1, 1000)
    true_effect = Xtest[:, 1].>0.5

    cf1 = build_forest(true, false, false, Y, T, X, false, -1, 100)
    pred1 = apply_forest(cf1, Xtest)
    comparison1 = isapprox.(true_effect, pred1, atol=0.4)
    @test (sum(comparison1)/1000) > 0.6

    cf2 = build_forest(true, false, true, Y, T, X, false, -1, 100)
    pred2 = apply_forest(cf2, Xtest)
    comparison2 = isapprox.(true_effect, pred2, atol=0.3)
    @test (sum(comparison2)/1000) > 0.7

    cf3 = build_forest(true, true, false, Y, T, X, false, -1, 100)
    pred3 = apply_forest(cf3, Xtest)
    comparison3 = isapprox.(true_effect, pred3, atol=0.3)
    @test (sum(comparison3)/1000) > 0.7

    cf4 = build_forest(true, true, true, Y, T, X, false, -1, 100)
    pred4 = apply_forest(cf4, Xtest)
    comparison4 = isapprox.(true_effect, pred4, atol=0.3)
    @test (sum(comparison4)/1000) > 0.7

end


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
# Author: Michael Broxton (broxton@stanford.edu)
#
# Copyright (C) 2015 Stanford University.
# All rights reserved.
#
# __END_LICENSE__

#__precompile__()   # uncomment to enable precompilation of the full module
module WaveOptics
using FFTW
import ThreeDeconv

# --------
# Includes

# Core
include("core/fftfreq.jl")
include("core/fftpad.jl")
include("core/meshgrid.jl")

# PointSpreadFunction
include("psf/coordinates.jl")
include("psf/recipe.jl")
include("psf/utility.jl")
include("psf/transform.jl")
include("psf/apsf.jl")
include("psf/debye.jl")
include("psf/defocus.jl")

end # WaveOptics module


text:
# ---
# title: 873. Length of Longest Fibonacci Subsequence
# id: problem873
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/length-of-longest-fibonacci-subsequence/description/>
# hidden: true
# ---
#
# A sequence `X_1, X_2, ..., X_n` is _fibonacci-like_ if:
#
#   * `n >= 3`
#   * `X_i + X_{i+1} = X_{i+2}` for all `i + 2 <= n`
#
# Given a **strictly increasing**  array `A` of positive integers forming a
# sequence, find the **length** of the longest fibonacci-like subsequence of
# `A`.  If one does not exist, return 0.
#
# ( _Recall that a subsequence is derived from another sequence`A` by deleting
# any number of elements (including none) from `A`, without changing the order
# of the remaining elements.  For example, `[3, 5, 8]` is a subsequence of `[3,
# 4, 5, 6, 7, 8]`._)
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,4,5,6,7,8]
#     Output: 5
#     Explanation: The longest subsequence that is fibonacci-like: [1,2,3,5,8].
#
#
# **Example 2:**
#
#
#
#     Input: [1,3,7,11,12,14,18]
#     Output: 3
#     **Explanation** :
#     The longest subsequence that is fibonacci-like:
#     [1,11,12], [3,11,14] or [7,11,18].
#
#
#
#
# **Note:**
#
#   * `3 <= A.length <= 1000`
#   * `1 <= A[0] < A[1] < ... < A[A.length - 1] <= 10^9`
#   * _(The time limit has been reduced by 50% for submissions in Java, C, and C++.)_
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Emulates the legacy `use_response_type=false` response behavior using the `AWS.Response`
function legacy_response(
    request::AWS.Request, response::AWS.Response; return_headers=nothing
)
    response_dict_type = something(request.response_dict_type, LittleDict)
    return_headers = something(return_headers, false)

    # When a user defined I/O stream is passed in use the actual `HTTP.Response` body
    # instead of the `AWS.Response` body which requires the I/O stream to be seekable.
    body = request.response_stream !== nothing ? response.response.body : response.body

    # The stored service name is always lowercase and may not match the module name
    # specified by the user. We'll assume that the typical casing used is titlecase.
    alt_service = "@service $(titlecase(request.service)) use_response_type=true"

    # When a user specifies a `response_dict_type` we'll update the deprecations to show how
    # to use this type.
    # Note: Using overly terse function name to stick within line length. A more descriptive
    # function name would be `response_dict_type_str`.
    tstr = if request.response_dict_type !== nothing
        str -> "$(request.response_dict_type)($str)"
    else
        identity
    end

    # For HEAD request, return headers...
    if request.request_method == "HEAD"
        Base.depwarn(
            "Using \"HEAD\" in AWS requests to return headers is deprecated, " *
            "use `$alt_service` to return an `AWS.Response` allowing for " *
            "header access via " *
            "`$(tstr("response.headers"))`.",
            :legacy_response,
        )

        return response_dict_type(response.headers)
    end

    # Return response stream if requested...
    if something(request.return_stream, false)
        Base.depwarn(
            "The keyword `return_stream` is deprecated, " *
            "use `$alt_service` instead to return an `AWS.Response` allowing for " *
            "streaming access via `response.io`.",
            :legacy_response,
        )

        # When `return_stream` was `true` the backends would create different I/O types to
        # use. We'll replicate that behavior.
        if request.response_stream === nothing
            io = if request.backend isa HTTPBackend
                Base.BufferStream()
            else
                IOBuffer()
            end

            write(io, _rewind(read, response.io))
            request.response_stream = io
        end

        # Emulate HTTP 0.9.14 behavior of always closing the passed in stream. Doing this is
        # particularly important for `Base.BufferStream`.
        if request.backend isa HTTPBackend
            close(request.response_stream)
        end

        return request.response_stream
    elseif request.response_stream !== nothing
        # Emulate HTTP 0.9.14 behavior of always closing the passed in stream and the
        # `read_body` behavior of closing all non-`IOBuffer` streams when `return_stream` is
        # not `true`. Doing this is particularly important for `Base.BufferStream`.
        if request.backend isa HTTPBackend || !(request.response_stream isa IOBuffer)
            close(request.response_stream)
        end
    end

    # Return raw data if requested...
    if something(request.return_raw, false)
        if return_headers
            Base.depwarn(
                "The keywords `return_raw` and `return_headers` are deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw data and header access via " *
                "`response.body`/`response.header`.",
                :legacy_response,
            )
        else
            Base.depwarn(
                "The keyword `return_raw` is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw data access via " *
                "`response.body`.",
                :legacy_response,
            )
        end

        return (return_headers ? (body, response.headers) : body)
    end

    # Parse response data according to mimetype...
    mime = HTTP.header(response.response, "Content-Type", "")

    if isempty(mime)
        if length(body) > 5 && body[1:5] == b"<?xml"
            mime = "text/xml"
        end
    end

    body_str = String(copy(body))

    if occursin(r"/xml", mime)
        if return_headers
            Base.depwarn(
                "The keyword `return_headers` is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "parsed XML and header access via " *
                "`$(tstr("parse(response)"))`/`$(tstr("response.header"))` respectively",
                :legacy_response,
            )
        else
            Base.depwarn(
                "Returning the parsed AWS response is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "parsed XML access via " *
                "`$(tstr("parse(response.body)"))`.",
                :legacy_response,
            )
        end

        xml_dict_type = response_dict_type{Union{Symbol,String},Any}
        xml = parse_xml(body_str)
        root = XMLDict.root(xml.x)

        return (
            if return_headers
                (xml_dict(root, xml_dict_type), response_dict_type(response.headers))
            else
                xml_dict(root, xml_dict_type)
            end
        )
    elseif occursin(r"/x-amz-json-1.[01]$", mime) || endswith(mime, "json")
        if return_headers
            Base.depwarn(
                "The keyword `return_headers` is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "parsed JSON and header access via " *
                "`$(tstr("parse(response)"))`/`$(tstr("response.header"))` respectively",
                :legacy_response,
            )
        else
            Base.depwarn(
                "Returning the parsed AWS response is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "parsed JSON access via " *
                "`$(tstr("parse(response.body)"))`.",
                :legacy_response,
            )
        end

        info = isempty(body) ? nothing : JSON.parse(body_str; dicttype=response_dict_type)
        return (return_headers ? (info, response_dict_type(response.headers)) : info)
    elseif startswith(mime, "text/")
        if return_headers
            Base.depwarn(
                "The keyword `return_headers` is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw string data and header access via " *
                "`String(response.body)`/`$(tstr("response.header"))` respectively.",
                :legacy_response,
            )
        else
            Base.depwarn(
                "Returning the raw AWS response body is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw string data access via " *
                "`String(response.body)`.",
                :legacy_response,
            )
        end

        return (
            return_headers ? (body_str, response_dict_type(response.headers)) : body_str
        )
    else
        if return_headers
            Base.depwarn(
                "The keyword `return_headers` is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw data and header access via " *
                "`response.body`/`response.header` respectively.",
                :legacy_response,
            )
        else
            Base.depwarn(
                "Returning the raw AWS response body is deprecated, " *
                "use `$alt_service` instead to return an `AWS.Response` allowing for " *
                "raw data access via " *
                "`response.body`.",
                :legacy_response,
            )
        end

        return (return_headers ? (body, response.headers) : body)
    end
end


text:
using Orekit
using Test

import Dates

Orekit.init()

@testset "Orekit" begin
    @testset "Time" begin
        TAI = get_tai()
        dt = Dates.DateTime(2020, 1, 1, 0, 0, 12, 500)
        ep1 = AbsoluteDate(2020, 1, 1, 0, 0, 12.5, TAI)
        ep2 = AbsoluteDate(dt, TAI)
        ep3 = AbsoluteDate(string(dt), TAI)
        dt1 = Dates.DateTime(ep1, TAI)
        dt2 = Dates.DateTime(ep2, TAI)
        dt3 = Dates.DateTime(ep3, TAI)
        @test dt == dt1
        @test dt == dt2
        @test dt == dt3
    end
end


text:
function parsefields(line::AbstractString, delim, quotes, escape, trimwhitespace::Bool)
    fields = SubString{String}[]
    isquoted = falses(0)
    badbreak = false
    eol = lastindex(line)
    fieldstart = firstindex(line)
    # delimrange = search(line, delim, fieldstart)
    delimrange = something(findnext(isequal(delim), line, fieldstart), 0)
    delimstart, nextfieldstart = first(delimrange), nextind(line, last(delimrange))
    done = false
    while !done
        fieldend = delimstart == 0 ? eol : prevind(line, delimstart)
        f = SubString(line, fieldstart, fieldend)
        field, quoted, badbreak = checkfield(f, quotes, escape, trimwhitespace)
        if badbreak
            if delimstart < 1
                return fields, isquoted, badbreak
            else
                delimrange = something(findnext(isequal(delim), line, nextfieldstart), 0)
                # delimrange = search(line, delim, nextfieldstart)
                delimstart, nextfieldstart = first(delimrange), nextind(line, last(delimrange))
                continue
            end
        end
        push!(fields, field)
        push!(isquoted, quoted)
        if nextfieldstart <= fieldstart
            done = true
        else
            fieldstart = nextfieldstart
        end
        # delimrange = search(line, delim, fieldstart)
        delimrange = something(findnext(isequal(delim), line, fieldstart), 0)
        delimstart, nextfieldstart = first(delimrange), nextind(line, last(delimrange))
    end
    return fields, isquoted, badbreak
end

function parsefields(line::AbstractString, delim, quotes::Missing, escape::Missing, trimwhitespace::Bool)
    fields = split(line, delim)
    isquoted = falses(length(fields))
    if trimwhitespace
        for (i, f) in enumerate(fields)
            fields[i] = strip(fields[i])
        end
    end
    return fields, isquoted, false
end


text:
"""
A collection of abstract types representing problems and methods for their solution
"""
module AbstractProblems

export Problem, Algorithm, IterativeAlgorithm, solve

"""
Abstract type representing a problem to solve
"""
abstract type Problem end

"""
Abstract type containing any algorithm
"""
abstract type Algorithm end

"""
Iterative methods form an important class of algorithms with iteratively adjust solution.

Concrete types of the `IterativeAlgorithm` should contain the initial value, tolerance and maximum number of iterations and
are used more for convenience. These can be obtained by fucntions `initial`, `tolerance`, `maxit` fuctions. If applied the
abstract types, these fucntions return `missing` and can trigger using of the default values.

In addition, the concrete types contain instructions on whether or not to keeep the history of the convergence and
some snapshots of the inner state of an iterator. These values are given by functions `keephistory` and `snapshots`.

"""
abstract type IterativeAlgorithm <: Algorithm end

initial(alg::IterativeAlgorithm) = missing
tolerance(alg::IterativeAlgorithm) = missing
maxit(alg::IterativeAlgorithm) = missing
keephistory(alg::IterativeAlgorithm) = false
snapshots(alg::IterativeAlgorithm) = Int64[]

"""
    solve(p::Problem,x⁰,alg::Algorithm)

    solve(p::Problem,alg::IterativeAlgorithm; x⁰, ϵ, maxit, keephistory, snapshshots)

Solve problem `p`, using method `alg`. For iterative algorithms the arguments may be specified separately.
Optionally keep the error history and the iteration snapshots.

"""
solve(p::Problem, alg::Algorithm) = error("Don't know how to solve ", typeof(p), " with method ", typeof(alg))

#unpack the values of iterative algorithm
solve(p::Problem, alg::IterativeAlgorithm; x⁰ = initial(alg), ϵ =tolerance(alg), maxit = maxit(alg), keephistory = keephistory(alg), snapshots = snapshots(alg)) = solve(p, alg, x⁰, ϵ, maxit, keephistory, snapshots)
# solve(p::Problem, alg::IterativeAlgorithm) = solve(p, alg, initial(alg), tolerance(alg), maxit(alg), keephistory(alg), snapshots(alg))

# Now proceed with the unpacked default or specified by the user values
function solve(p::Problem, alg::IterativeAlgorithm,  args...)
    error("Don't know how to solve ", typeof(p), " with method ", typeof(alg))
end

end

# these types and fucntions will be modified by including file, so we import them.
import .AbstractProblems: Problem, Algorithm, IterativeAlgorithm, solve, initial, tolerance, maxit, keephistory, snapshots


text:
# convert_a3d

# Function that allows federation by setting the `output_format`  in the Stanmodel().

function convert_a3d(a3d_array, cnames, ::Val{:mcmcchains})
  pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
  p = filter(p -> !(p in  pi), cnames)

  MCMCChains.Chains(a3d_array,
    cnames,
    Dict(
      :parameters => p,
      :internals => pi
    )
  )
end


text:
@testset "square" begin
    b1 = [0, 0, 5, 10]
    b2 = [0, 0, 10, 5]

    @test square(b1) == [-2.5, 0, 7.5, 10]
    @test square(b2) == [0, -2.5, 10, 7.5]
end


text:
function get_value(::Type{NewcombOperators}, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(NewcombOperators, "getValue", jdouble, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end


text:
// Assigning double to int variable.

int main () {
 int x = 1.0;
 return 0 ;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
function dynmodel!(xt::AbstractVector, ut::AbstractVector, params::AbstractVector)
    ## Parameters and constants
    dt = 0.01
    # motor
    k_motor = exp(params[1])
    # k1_left = exp(params[1]) # Motor constant: current -> torque
    # k1_right = exp(params[2]) # Motor constant: current -> torque
    # segway
    Mb = 1.76 # kg # Masse des Körpers
    Mw = 0.147 # kg # Masse der Räder
    R = 0.07 # m #Radius der Räder
    cz = 0.077 # m #Schwerpunkt des Körpers Höhe über Radachse
    b = 0.1985/2 # m #Half the distance between the wheels (OO_w in paper)
    Ixx = (0.166^2+0.21^2)*Mb/12 + Mb*(cz-0.07)^2 # kg m^2 #Inertia of the body
    Iyy = (0.072^2+0.21^2)*Mb/12 + Mb*(cz-0.07)^2 # kg m^2 #Inertia of the body
    Izz = (0.072^2+0.166^2)*Mb/12 # kg m^2 #Inertia of the body
    Iwa = Mw*R^2/2 # kg m^2 #Inertia of the wheels about their axis
    Iwd = Mw*b^2 # kg m^2 #Inertia of the wheels about their axis
    seg = Segway_unmutable(Mb, Mw, R, cz, b, Ixx, Iyy, Izz, Iwa, Iwd)

    # friction
    fl_1 = exp(params[2]) # left wheel
    fr_1 = exp(params[3]) # right wheel
    fl_2 = exp(params[4]) # left wheel
    fr_2 = exp(params[5]) # right wheel

    ## Evaluate dynamics
    # get wheel rotation speeds
    v_left = xt[6]-xt[7]*b
    v_right = xt[6]+xt[7]*b
    omega_right = v_right/R - xt[5]
    omega_left = v_left/R - xt[5]
    # apply simplified motor model (omega taken out in favor of viscose friction)
    tau_right = k_motor*ut[2]
    tau_left = k_motor*ut[1]
    # apply simplified friction model
    tau_right = frictionmodel(tau_right, omega_right, fr_1, fr_2)
    tau_left = frictionmodel(tau_left, omega_left, fl_1, fl_2)
    # apply segway model
    tau = [tau_right, tau_left]
    segway_rk4!(xt, tau, dt, seg)

    return xt
end

@inline function frictionmodel(tau, omega, param1, param2)
    # Viskoser Reibterm
    dtau_v = 0.0

    # Haftreibung (geglättet)
    dtau_h = tau # Case tau < Fs and no rotation
    if (abs(omega) < 1e-3) && (param2 < tau) # Case tau > Fs and no rotation
        dtau_h = sign(tau)*param2
    elseif (abs(omega) > 1e-3)
        dtau_h = sign(omega)*param2
        dtau_v = param1*omega
    end

    return tau-dtau_v-dtau_h
end

# Does not overwrite x
function dynmodel(x::AbstractVector, u::AbstractVector, params::AbstractVector)
    x2 = copy(x)
    dynmodel!(x2, u, params)
    return x2
end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function segway_cost(xt::AbstractVector, u::AbstractVector)
    cost = (xt[4]^2 + 0.1*xt[5]^2
            + 10*xt[2]^2 + 10*xt[1]^2 # or squared here...
            + 0.1* xt[7]^2
            + 0.01*(u'u)) # punish for falling over

    reward = 0.0 #reward for moving along the x-axis

    if xt[4]^2 > 2 #check if lager than pi/2
        terminal = 1
        cost *= 1.2
    else
        terminal = 0
    end

    if sum(isnan, xt) > 0
        terminal = 1
    end

    return cost-reward, terminal
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using DataStructures: CircularBuffer
include("AlgorithmicEnv.jl")

mutable struct ReversedAdditionEnv <: GridAlgorithmicEnv
    MIN_REWARD_SHORTFALL_FOR_PROMOTION::Float32
    base::Int
    episode_total_reward::Float32
    reward_shortfalls::CircularBuffer{Float32}
    charmap::Array{Char, 1}
    min_length::Int
    action_space::TupleSpace
    observation_space::Discrete

    MOVEMENTS::NTuple
    READ_HEAD_START::Array{Int8, 1}

    rows::Int

    target
    input_data
    time
    read_head_position
    write_head_position
    last_action
    last_reward
end

ReversedAdditionEnv(base::Int=3; rows::Int=2) =
    ReversedAdditionEnv(-1f0, # MIN_REWARD_SHORTFALL_FOR_PROMOTION
            base, # last
            0f0, # episode_total_reward
            CircularBuffer{Float32}(10),  # reward_shortfalls
            push!(['0'+i for i=0:base-1], ' '),  # charmap
            2,  # starting min_length
            TupleSpace([Discrete(2), Discrete(2), Discrete(base)]),  # action_space
            Discrete(base+1),  # observation_space
            (:left, :right, :up, :down), [1, 1],  # MOVEMENTS and READ_HEAD_START
            rows,
            Int8[], Int8[], 0, [1, 1], [1, 1], nothing, 0.0)

function target_from_input_data!(env::ReversedAdditionEnv)
    curry = 0
    env.target = []
    for digits in env.input_data
        total = sum(digits) + curry
        push!(env.target, total % env.base)
        curry = div(total, env.base)
    end
    curry > 0 && (push!(env.target, curry))
end

time_limit(env::ReversedAdditionEnv) = len(env.input_data) * 4


text:
# ------------ Choice sites ------------ #

@inline function (ctx::GenerateContext)(call::typeof(trace),
                                        addr::T,
                                        d::Distribution{K}) where {T <: Address, K}
    visit!(ctx, addr)
    if has_value(ctx.target, addr)
        s = getindex(ctx.target, addr)
        score = logpdf(d, s)
        add_choice!(ctx, addr, score, s)
        increment!(ctx, score)
    else
        s = rand(d)
        add_choice!(ctx, addr, logpdf(d, s), s)
    end
    return s
end

# ------------ Learnable ------------ #

@inline function (ctx::GenerateContext)(fn::typeof(learnable), addr::Address)
    visit!(ctx, addr)
    haskey(ctx.params, addr) && return getindex(ctx.params, addr)
    error("(learnable): parameter not provided at address $addr.")
end

# ------------ Black box call sites ------------ #

@inline function (ctx::GenerateContext)(c::typeof(trace),
                                        addr::T,
                                        call::Function,
                                        args...) where T <: Address
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ss = get_sub(ctx.target, addr)
    ret, cl, w = generate(ss, ps, call, args...)
    add_call!(ctx, addr, cl)
    increment!(ctx, w)
    return ret
end

@inline function (ctx::GenerateContext)(c::typeof(trace),
                                        addr::T,
                                        call::G,
                                        args...) where {G <: GenerativeFunction,
                                                        T <: Address}
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ss = get_sub(ctx.target, addr)
    tr, w = generate(call, args, JChoiceMap(ss))
    ret = get_retval(tr)
    add_call!(ctx, addr, DynamicCallSite(get_choices(tr), get_score(tr), get_gen_fn(tr), get_args(tr), ret))
    increment!(ctx, w)
    return ret
end

# ------------ Convenience ------------ #

function generate(opt::J, target::L, params, fn::Function, args...) where {J <: CompilationOptions, L <: AddressMap}
    ctx = Generate(opt, DynamicTrace(), target, params)
    ret = ctx(fn, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.weight
end

function generate(target::L, params, fn::Function, args...) where L <: AddressMap
    ctx = Generate(DefaultPipeline(), DynamicTrace(), target, params)
    ret = ctx(fn, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.weight
end

function generate(target::L, fn::Function, args...) where L <: AddressMap
    ctx = Generate(DefaultPipeline(), DynamicTrace(), target, Empty())
    ret = ctx(fn, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.weight
end

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct GenerateContext{T <: AddressMap, K <: AddressMap, P <: Parameters} <: ExecutionContext
     tr::T
     target::K
     weight::Float64
     score::Float64
     visited::Visitor
     params::P
end
```
`GenerateContext` is used to generate traces, as well as record and accumulate likelihood weights given observations at addressed randomness.

Inner constructors:
```julia
GenerateContext(tr::T, target::K) where {T <: AddressMap, K <: AddressMap} = new{T, K}(tr, target, 0.0, Visitor(), Parameters())
GenerateContext(tr::T, target::K, params::P) where {T <: AddressMap, K <: AddressMap, P <: Parameters} = new{T, K, P}(tr, target, 0.0, Visitor(), params)
```
Outer constructors:
```julia
Generate(target::AddressMap) = GenerateContext(AddressMap(), target)
Generate(target::AddressMap, params) = GenerateContext(AddressMap(), target, params)
Generate(tr::AddressMap, target::AddressMap) = GenerateContext(tr, target)
```
""", GenerateContext)

@doc(
"""
```julia
ret, cl, w = generate(target::L, fn::Function, args...; params = Parameters()) where L <: AddressMap
ret, cs, w = generate(target::L, fn::typeof(rand), d::Distribution{K}; params = Parameters()) where {L <: AddressMap, K}
ret, v_cl, w = generate(target::L, fn::typeof(markov), call::Function, len::Int, args...; params = Parameters()) where L <: AddressMap
ret, v_cl, w = generate(target::L, fn::typeof(plate), call::Function, args::Vector; params = Parameters()) where L <: AddressMap
ret, v_cl, w = generate(target::L, fn::typeof(plate), d::Distribution{K}, len::Int; params = Parameters()) where {L <: AddressMap, K}
```

`generate` provides an API to the `GenerateContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, a `RecordSite` instance specialized to the call, and the score/weight `w` computed with respect to the constraints `target`.
""", generate)


text:
module PyTorchTraining

using Flux
using FluxTraining, FluxTraining.Events, FluxTraining.Phases
using FluxTraining: runstep, runepoch
using InlineTest
using PyCall
using Markdown
import Zygote

include("tensor.jl")
include("model.jl")
include("optimisers.jl")
include("loss.jl")

include("callback.jl")

const torch = PyNULL()

function __init__()
    copy!(torch, pyimport("torch"))
end

export PyTorchBackend, totensor, fromtensor

end


text:
using Pseudoseq, Pseudoseq.Sequencing, Documenter, Literate

Documenter.post_status(; type = "pending", repo = "github.com/bioinfologics/Pseudoseq.jl.git")

SEQ_EXAMPLES = [joinpath("examples/sequencing", f) for f in ("pe-example.jl",
                                                             "se-example.jl",
                                                             "tg-example.jl")]

#BAG_EXAMPLES = [joinpath("examples/build-a-genome", f) for f in ("build-a-yeast.jl",)]

OUTPUT = joinpath(@__DIR__, "src/man")
SEQ_OUTPUT = joinpath(OUTPUT, "sequencing/examples")
#BAG_OUTPUT = joinpath(OUTPUT, "build-a-genome/examples")
mkdir(SEQ_OUTPUT)
#mkdir(BAG_OUTPUT)

cp("examples/sequencing/ecoli-ref.fasta", "docs/src/man/sequencing/examples/ecoli-ref.fasta")
#cp("examples/build-a-genome/yeast-chr1.fasta", "docs/src/man/build-a-genome/examples/yeast-chr1.fasta")

for ex in SEQ_EXAMPLES
    Literate.markdown(ex, SEQ_OUTPUT)
    Literate.notebook(ex, SEQ_OUTPUT)
    Literate.script(ex, SEQ_OUTPUT)
end

#for ex in BAG_EXAMPLES
#    Literate.markdown(ex, BAG_OUTPUT)
#    Literate.notebook(ex, BAG_OUTPUT)
#    Literate.script(ex, BAG_OUTPUT)
#end

makedocs(
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "GITHUB_ACTIONS")
    ),
    modules = [Pseudoseq, Sequencing],
    sitename = "Pseudoseq.jl",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "PuzzleMaker" => "man/puzzlemaker/concepts.md",
            "Sequencing" => [
                "Core concepts & workflow" => "man/sequencing/concepts.md",
                "Examples" => [
                    "Paired end reads" => "man/sequencing/examples/pe-example.md",
                    "Long single end reads" => "man/sequencing/examples/se-example.md",
                    "Tagged paired end reads" => "man/sequencing/examples/tg-example.md"
                ]
            ]
        ],
        "API" => [
            "Sequencing" => [
                "Molecule Pool" => "api/pool.md",
                "Reads" => "api/reads.md",
            ]
        ]
    ],
    authors = "Ben J. Ward."
)

# Scripts will generate some read file output FASTQs, we don't really want to
# keep these, and their size might trip up CI when we upload docs to gh-pages,
# so we find all files ending in .fastq, and delete them.
run(`find . -name "*.fastq" -type f -delete`)

deploydocs(
    repo = "github.com/bioinfologics/Pseudoseq.jl.git",
    push_preview = true
)


text:
U(T) = Union{Missing, T}

function add_row!(t, named_tuple)
	push!(rows(t), named_tuple)
    table(t, pkey = t.pkey, copy = false)
    return
end

macro check_not_missing(heading, args...)
    expr = Expr(:block)
    for arg in args
        push!(expr.args, :(!(($arg)[1] isa Missing) || (warn("$($heading) $(($arg)[2]) is missing and must be set."); return )))
    end
    esc(expr)
end

macro check_type(T, heading, args...)
    expr = Expr(:block)
    for arg in args
        push!(expr.args, :(($arg)[1] isa $T || (warn("$($heading) $(($arg)[2]) must be of type $($T)."); return )))
    end
    esc(expr)
end

newid(t) = reduce(max, t, select = :ID, 0) + 1

macro check_valid_id(t, field, id, heading)
    return esc(:($id ∈ columns($t).$field || (warn("$($heading) id does not exist."); return )))
end

function DIR(repository)
    d = Base.find_in_path(repository)
    gitdir = joinpath(d[1:search(d, "\\src\\").start-1], ".git")
    gitdir
end

function HEAD(repository)
    gitdir = DIR(repository)
    commit = strip(readstring(`git --git-dir $gitdir rev-parse HEAD`))
    commit
end

function STASH(repository)
    gitdir = DIR(repository)
    run(`git --git-dir $gitdir stash`)
    return
end

function CHECKOUT(repository, commit)
    gitdir = DIR(repository)
    run(`git --git-dir $gitdir checkout $commit`)
    return
end

function editvalue!(t, f, i, v)
    t.columns.columns[f][i] = v
    return
end