text:
"""
    DMU()
---
# Description
Prepare data and driver files for DMU(https://dmu.ghpc.au.dk/).
One can prepare a sample dataset with `ABG.simu_qp(h2=0.8)`.
Since there is no pedigree info in toy dataset, that part is ignored here.

# Parameters

Below are for DMU1, a must-run procedure of DMU.

## ANALYSE
- `task` = 1, variance components estimation using REML
- `method` = 31, AI (average information), (updates ⊄ par space) && (AI + EM)
- `scaling` = 0, no data scaling before computation
- `test_prt` = 0, standard

## DATA
- `FMT` = ASCII
- (#`int`, #`real`, `miss`) = (1, 1, -99)
- `fn` = project.dat, file name

## VARIABLE
- ID phenotype

## MODEL
- 1          # 1 trait(s)
- 0          # One line for each traits.  Usually zero
- 1 2 2 2 1  #
- 1 1
- 0
- 0

## VAR_STR
- `r_factor`: 1, structure number
- `type`: GREL,
- `FMT`: ASCII, I will try BINARY here
- `fn`: toy.giv

## PRIOR
-
1 1 1 139666
2 1 1 224499

## DMUAI
10
1.0d-7
1.0d-6
1
0
"""
function DMU(;
             dir = abg_rst,
             project = "toy",
             desc = ["Analysis on toy data"],
             giv = "toy.giv")
    subtitle("Prepare $project data and driver files for DMU")

    item("The driver file")
    driver = joinpath(dir, "toy.drv")

    drvctt = ["\$COMMENT",
              join(desc, '\n'),
              "",
              "\$ANALYSE 1 31 0 0", # task=1, var
              "",
              "\$DATA ASCII (2, 2, -99) $project.dat",
              "",
              "\$VARIABLE",
              "ID phenotype",
              "",
              "\$VAR_STR 1 GREL BINARY $giv",
              "",
              "\$PRIOR",
              "",
              "\$DMUAI",
              "",
              ]
    write(driver, join(drvctt, '\n'))
    open(joinpath(dir, "toy.dat"), "w") do io
        id = 1
        for line in eachline(joinpath(dir, "toy.ph"))
            println(io, id, ' ', line)
            id += 1
        end
    end
    done()
end

#= sample driver for DMU1
$COMMENT
Dynamisk hy-definisjon
Analyse av kg melk, 1-3 laktasjon, gjentak-dyremodell random HY gen grup

$ANALYSE 1 31 0 0

$DATA ASCII (2,2,-99) fat.dat

$VARIABLE
# 1 2
id  line
#1 2
yld wg

$MODEL
1
0
1 2 2 2 1
1 1
0
0

$VAR_STR 1 GREL ASCII /home/theom/data/ulf/4k_qtl/10000genots/chr1/mcmc2/melb/dmu_/g.giv

$PRIOR
1 1 1 139666
2 1 1 224499

$DMUAI
10
1.0d-7
1.0d-6
1
0
=#


text:
import Base.show
import DataStructures
import Random

struct Node{InputType, MetricReturnType}
    index::Int
    data::InputType
    radius::MetricReturnType
    min_dist::MetricReturnType
    max_dist::MetricReturnType
    n_data::Int
    left_child::Union{Node{InputType, MetricReturnType}, Nothing}
    right_child::Union{Node{InputType, MetricReturnType}, Nothing}
end

function Base.show(io::IO, n::Node)
    print("Node $(typeof(n).parameters[1]): $(n.data), index $(n.index) radius $(n.radius)")
end

"""
    VPTree(data::Vector{InputType}, metric; threaded=nothing)

Construct Vantage Point Tree with a vector of `data` and given a callable `metric`.
`threaded` uses threading is only avaible in Julia 1.3+ to parallelize construction of the Tree.
When not explicitly set, is set to true when the necessary conditions are met.

## Example:

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
vptree = VPTree(data, Levenshtein())
query = "blau"
radius = 2
data[find(vptree, query, radius)]
# 2-element Array{String,1}:
#  "bla"
#  "blub"
n_neighbors = 3
data[find_nearest(vptree, query, n_neighbors)]
# 3-element Array{String,1}:
#  "baube"
#  "blub"
#  "bla"
```
"""
struct VPTree{InputType, MetricReturnType}
    data::Vector{InputType}
    metric
    root::Node{InputType, MetricReturnType}
    threaded::Bool
    function VPTree(data::Vector{InputType}, metric; threaded=nothing) where InputType
        threaded = _check_threaded(threaded)
        @assert length(data) > 0 "Input data must contain at least one point."
        @assert !isempty(methods(metric)) "`metric` must be callable, was: $(metric)"
        MetricReturnType = typeof(metric(data[1], data[1]))
        indexed_data = Random.shuffle!(collect(enumerate(data)))
        root = threaded ? _construct_tree_rec_threaded!(indexed_data, metric, MetricReturnType) : _construct_tree_rec!(indexed_data, metric, MetricReturnType)
        new{InputType, MetricReturnType}(data, metric, root, threaded)
    end
end

function _check_threaded(threaded)
    if !isnothing(threaded) && threaded == true
        if VERSION < v"1.3-DEV"
            @warn "incompatible julia version for `threaded=true`: $VERSION, requires >= v\"1.3\", setting `threaded=false`"
            threaded = false
        elseif Threads.nthreads() == 1
            @warn "`threaded = true`, but `Threads.nthreads() == 1`, setting `threaded=false`"
            threaded = false
        end
    end
    if isnothing(threaded)
        threaded = VERSION >= v"1.3-DEV" && Threads.nthreads() > 1
    end
    threaded
end

function VPTree(data::Vector{T}, metric::Function, MetricReturnType) where T
    VPTree(data, metric)
end

const SMALL_DATA = 1000

@deprecate VPTree(data::Vector, metric::Function, MetricReturnType::DataType) VPTree(data::Vector, metric)

function _construct_tree_rec_threaded!(data::AbstractVector{Tuple{Int, InputType}}, metric, MetricReturnType) where InputType
    if isempty(data)
        return nothing
    end
    n_data = length(data)
    if n_data == 1
        return Node(data[1][1], data[1][2], zero(MetricReturnType), zero(MetricReturnType), zero(MetricReturnType), n_data, nothing, nothing)
    end
    vantage_point = data[1]
    rest = view(data, 2:length(data))
    distances = [metric(d[2], vantage_point[2]) for d in rest]
    i_middle = div(n_data - 1, 2) + 1
    distance_order = sortperm(distances, alg=PartialQuickSort(i_middle))
    permute!(rest, distance_order)

    min_dist, max_dist = extrema(distances)
    radius = distances[distance_order[i_middle]]
    left_rest = view(rest, 1:i_middle)

    should_spawn = length(rest) > SMALL_DATA
    left_node = if should_spawn
        Threads.@spawn _construct_tree_rec_threaded!(left_rest, metric, MetricReturnType)
    else
        _construct_tree_rec!(left_rest, metric, MetricReturnType)
    end

    right_rest = view(rest, i_middle + 1:length(rest))
    right_node = _construct_tree_rec_threaded!(right_rest, metric, MetricReturnType)

    if should_spawn
        Node{InputType, MetricReturnType}(vantage_point[1], vantage_point[2], radius,  min_dist, max_dist, n_data, fetch(left_node), right_node)
    else
        Node{InputType, MetricReturnType}(vantage_point[1], vantage_point[2], radius,  min_dist, max_dist, n_data, left_node, right_node)
    end
end

function _construct_tree_rec!(data::AbstractVector{Tuple{Int, InputType}}, metric, MetricReturnType) where InputType
    if isempty(data)
        return nothing
    end
    n_data = length(data)
    if n_data == 1
        return Node(data[1][1], data[1][2], zero(MetricReturnType), zero(MetricReturnType), zero(MetricReturnType), n_data, nothing, nothing)
    end
    vantage_point = data[1]
    rest = view(data, 2:length(data))
    distances = [metric(d[2], vantage_point[2]) for d in rest]
    i_middle = div(n_data - 1, 2) + 1
    distance_order = sortperm(distances, alg=PartialQuickSort(i_middle))
    permute!(rest, distance_order)

    left_rest = view(rest, 1:i_middle)

    left_node = _construct_tree_rec!(left_rest, metric, MetricReturnType)

    right_rest = view(rest, i_middle + 1:length(rest))
    right_node = _construct_tree_rec!(right_rest, metric, MetricReturnType)

    min_dist, max_dist = extrema(distances)
    radius = distances[distance_order[i_middle]]

    Node{InputType, MetricReturnType}(vantage_point[1], vantage_point[2], radius,  min_dist, max_dist, n_data, left_node, right_node)
end

"""
    find(vptree::VPTree{InputType, MetricReturnType}, query::InputType, radius::MetricReturnType, skip=nothing)::Vector{Int}

Find all items in `vptree` within `radius` of `query` with respect to the metric defined in the VPTree.
Returns indices into VPTree.data. The optional `skip` argument is a function `f(::Int)::Bool` which can
be used to omit points in the tree from the search based on their index.

## Example

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
metric = (a, b) -> evaluate(Levenshtein(),a,b)
vptree = VPTree(data, metric)
query = "blau"
radius = 2
data[find(vptree, query, radius)]
# 2-element Array{String,1}:
#  "bla"
#  "blub"
```
"""
function find(vptree::VPTree{InputType, MetricReturnType}, query::InputType, radius::MetricReturnType, skip=nothing)::Vector{Int} where {InputType, MetricReturnType}
    if vptree.threaded
        results_threads = [Vector{Int}() for i in 1:Threads.nthreads()]
        _find_threaded(vptree.root, query, radius, results_threads, vptree.metric, skip)
        return reduce(vcat, results_threads)
    else
        results = Vector{Int}()
        _find(vptree.root, query, radius, results, vptree.metric, skip)
        return results
    end
end

function _find(vantage_point, query, radius, results, metric, skip)
    distance = metric(vantage_point.data, query)
    if distance <= radius && (isnothing(skip) || !skip(vantage_point.index))
        push!(results, vantage_point.index)
    end
    if distance - radius <= vantage_point.radius && !isnothing(vantage_point.left_child) && distance + radius >= vantage_point.min_dist
        _find(vantage_point.left_child, query, radius, results, metric, skip)
    end
    if distance + radius >= vantage_point.radius && !isnothing(vantage_point.right_child) && distance - radius <= vantage_point.max_dist
        _find(vantage_point.right_child, query, radius, results, metric, skip)
    end
end

function _find_threaded(vantage_point, query, radius, results, metric, skip)
    distance = metric(vantage_point.data, query)
    if distance <= radius && (isnothing(skip) || !skip(vantage_point.index))
        push!(results[Threads.threadid()], vantage_point.index)
    end
    goleft = distance + radius >= vantage_point.radius && !isnothing(vantage_point.right_child) && distance - radius <= vantage_point.max_dist
    if distance - radius <= vantage_point.radius && !isnothing(vantage_point.left_child) && distance + radius >= vantage_point.min_dist
        if goleft && vantage_point.n_data > SMALL_DATA
            r = Threads.@spawn _find_threaded(vantage_point.left_child, query, radius, results, metric, skip)
        else
            _find_threaded(vantage_point.left_child, query, radius, results, metric, skip)
        end
    end
    if goleft
        _find_threaded(vantage_point.right_child, query, radius, results, metric, skip)
    end
    if @isdefined r
        wait(r)
    end
end

"""
    find_nearest(vptree::VPTree{InputType, MetricReturnType}, query::InputType, n_neighbors::Int, skip=nothing)::Vector{Int}

Find `n_neighbors` items in `vptree` closest to `query` with respect to the metric defined in the VPTree.
Returns indices into VPTree.data. The optional `skip` argument is a function `f(::Int)::Bool` which can
be used to omit points in the tree from the search based on their index.

## Example:

```julia
using VPTrees
using StringDistances

data = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
metric = (a, b) -> evaluate(Levenshtein(),a,b)
vptree = VPTree(data, metric)
query = "blau"
n_neighbors = 3
data[find_nearest(vptree, query, n_neighbors)]
# 3-element Array{String,1}:
#  "baube"
#  "blub"
#  "bla"
```
"""
function find_nearest(vptree::VPTree{InputType, MetricReturnType}, query::InputType, n_neighbors::Int, skip=nothing)::Vector{Int} where {InputType, MetricReturnType}
    @assert n_neighbors > 0 "Can't search for fewer than 1 neighbors"
    candidates = DataStructures.BinaryMaxHeap{Tuple{MetricReturnType, Int}}()
    _find_nearest(vptree.root, query, n_neighbors, candidates, vptree.metric, skip)
    [t[2] for t in candidates.valtree]
end

function _find_nearest(vantage_point, query, n_neighbors, candidates, metric, skip)
    distance = metric(vantage_point.data, query)
    radius = length(candidates) < n_neighbors ? typemax(typeof(vantage_point.radius)) : first(candidates)[1]
    if distance < radius && (isnothing(skip) || !skip(vantage_point.index))
        push!(candidates, (distance, vantage_point.index))
        if length(candidates) > n_neighbors
            pop!(candidates)
        end
    end
    # Switch radius to one side to prevent overflow.
    if distance < vantage_point.radius
        # Switch radius to one side to prevent overflow.
        if distance - vantage_point.radius <= radius && !isnothing(vantage_point.left_child) && distance - vantage_point.min_dist >= -radius
            _find_nearest(vantage_point.left_child, query, n_neighbors, candidates, metric, skip)
        end
        if distance - vantage_point.radius > -radius && !isnothing(vantage_point.right_child) && distance - vantage_point.max_dist <= radius
            _find_nearest(vantage_point.right_child, query, n_neighbors, candidates, metric, skip)
        end
    else
        if distance - vantage_point.radius > -radius && !isnothing(vantage_point.right_child) && distance - vantage_point.max_dist <= radius
            _find_nearest(vantage_point.right_child, query, n_neighbors, candidates, metric, skip)
        end
        if distance - vantage_point.radius <= radius && !isnothing(vantage_point.left_child) && distance - vantage_point.min_dist >= -radius
            _find_nearest(vantage_point.left_child, query, n_neighbors, candidates, metric, skip)
        end
    end
end

ext: jl


text:
@testset "H" begin
    @testset "hx2dp" begin
        @test hx2dp("2A^3") == 672.0
        @test_throws SpiceError hx2dp("2A^3000")
        @test_throws SpiceError hx2dp("-2A^3000")
        @test_throws SpiceError hx2dp("foo")
    end
    @testset "hrmint" begin
        xvals = [-1.0, 0.0, 3.0, 5.0]
        yvals = [6.0, 3.0, 5.0, 0.0, 2210.0, 5115.0, 78180.0, 109395.0]
        answer, deriv = hrmint(xvals, yvals, 2)
        @test answer ≈ 141.0
        @test deriv  ≈ 456.0
        @test_throws ArgumentError hrmint(xvals, yvals[1:4], 2)
        @test_throws SpiceError hrmint([0.0, 0.0], yvals[1:4], 2)
    end
    @testset "halfpi" begin
        @test SPICE._halfpi() ≈ π/2
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SiloBuilder"
version = v"0.1.0"

# Collection of sources required to build SiloBuilder
sources = [
    "https://wci.llnl.gov/content/assets/docs/simulation/computer-codes/silo/silo-4.10.2/silo-4.10.2-bsd.tar.gz" =>
    "4b901dfc1eb4656e83419a6fde15a2f6c6a31df84edfad7f1dc296e01b20140e",

]

# Bash recipe for building across all platforms
script = raw"""
set -e

cd $WORKSPACE/srcdir/silo-4.10.2-bsd/
./configure --disable-silex --disable-browser --prefix=$prefix --enable-shared --host=$target

# We need to do some patching to cross-compile for mingw32
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
CURRENTFILE=src/silo/silo_win32_compatibility.h
cp ${CURRENTFILE} ${CURRENTFILE}.orig
sed -e 's/#include <sys\types.h>/#include <sys\/types.h>/' -e 's/#include <sys\stat.h>/#include <sys\/stat.h>/' ${CURRENTFILE}.orig > ${CURRENTFILE}

CURRENTFILE=SiloWindows/include/config.h
cp ${CURRENTFILE} ${CURRENTFILE}.orig

# We don't have FPCLASS
sed 's/#define HAVE_FPCLASS 1/#define HAVE_FPCLASS 0/' ${CURRENTFILE}.orig > ${CURRENTFILE}

CURRENTFILE=tools/silock/silock.c
cp ${CURRENTFILE} ${CURRENTFILE}.orig

# Generate silock source file patch (we don't build silock)
cat > $WORKSPACE/silock.c.patch << 'END'
--- tools/silock/silock.c.orig
+++ tools/silock/silock.c
@@ -111,26 +111,7 @@
    static char lastDir[1024], lastVar[1024];
    char errMsg[128];

-   /* try to produce a useful error message regarding the kind of NaN */
-#ifdef HAVE_FPCLASS
-   {  fpclass_t theClass = fpclass(value);
-      switch (theClass)
-      {
-      case FP_SNAN:    strcpy(errMsg,"signaling NaN"); break;
-      case FP_QNAN:    strcpy(errMsg,"quiet NaN"); break;
-      case FP_NINF:    strcpy(errMsg,"negative infinity"); break;
-      case FP_PINF:    strcpy(errMsg,"positive infinity"); break;
-      case FP_NDENORM: strcpy(errMsg,"negative denormalized non-zero"); break;
-      case FP_PDENORM: strcpy(errMsg,"positive denormalized non-zero"); break;
-      case FP_NZERO:   strcpy(errMsg,"negative zero"); break;
-      case FP_PZERO:   strcpy(errMsg,"positive zero"); break;
-      case FP_NNORM:   strcpy(errMsg,"negative normalized non-zero"); break;
-      case FP_PNORM:   strcpy(errMsg,"positive normalized non-zero"); break;
-      }
-   }
-#else
    strcpy(errMsg, "unkown NaN");
-#endif

    if (!disableVerbose)
    {
END

# Patch silock source file
patch -l ${CURRENTFILE}.orig $WORKSPACE/silock.c.patch -o ${CURRENTFILE}
fi

# Generate Silock Makefile patch (we don't build silock)
CURRENTFILE=tools/silock/Makefile
cp ${CURRENTFILE} ${CURRENTFILE}.orig
cat > $WORKSPACE/silock-Makefile.patch << 'END'
--- tools/silock/Makefile.orig
+++ tools/silock/Makefile
@@ -401,7 +401,6 @@
        rm -f $$list
 silock$(EXEEXT): $(silock_OBJECTS) $(silock_DEPENDENCIES)
        @rm -f silock$(EXEEXT)
-       $(CXXLINK) $(silock_OBJECTS) $(silock_LDADD) $(LIBS)

 mostlyclean-compile:
        -rm -f *.$(OBJEXT)
END

# Patch silock Makefile
patch -l ${CURRENTFILE}.orig $WORKSPACE/silock-Makefile.patch -o ${CURRENTFILE}

# Finally, finish making library
make VERBOSE=1
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    MacOS(:x86_64)
]

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

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jgoldfar/ClpBuilder

max_issues_repo_name: jgoldfar/ClpBuilder

max_forks_repo_name: jgoldfar/ClpBuilder


text:
module MultipleTestingDoc

using MultipleTesting
using Lexicon

cd( joinpath(Pkg.dir("MultipleTesting"), "doc") )

config = Config(md_permalink = false,
                md_subheader = :category,
                include_internal = false,
                metadata_order = Symbol[],
                mathjax = true)
index = save("api.md", MultipleTesting, config);
save("index.md", Index([index]));

## notebooks
nb_files = ["pi0-estimation"]
cd( joinpath(Pkg.dir("MultipleTesting"), "doc", "notebooks") )

for nb in nb_files
    run(`jupyter nbconvert --to markdown "$nb.ipynb"`)
end

end


text:
module RetryRequest

import ..HTTP
using ..Sockets
using ..IOExtras
using ..MessageRequest
using ..Messages
import ..@debug, ..DEBUG_LEVEL, ..sprintcompact

export retrylayer

"""
    retrylayer(req) -> HTTP.Response

Retry the request if it throws a recoverable exception.

`Base.retry` and `Base.ExponentialBackOff` implement a randomised exponentially
increasing delay is introduced between attempts to avoid exacerbating network
congestion.

Methods of `isrecoverable(e)` define which exception types lead to a retry.
e.g. `HTTP.IOError`, `Sockets.DNSError`, `Base.EOFError` and `HTTP.StatusError`
(if status is ``5xx`).
"""
function retrylayer(handler)
    return function(req::Request; retry::Bool=true, retries::Int=4, retry_non_idempotent::Bool=false, kw...)
        if !retry || retries == 0
            # no retry
            return handler(req; kw...)
        end
        retry_request = Base.retry(handler,
            delays=ExponentialBackOff(n = retries),
            check=(s, ex)->begin
                retry = isrecoverable(ex, req, retry_non_idempotent, get(req.context, :retrycount, 0))
                if retry
                    @debug 1 "🔄  Retry $ex: $(sprintcompact(req))"
                    reset!(req.response)
                else
                    @debug 1 "🚷  No Retry: $(no_retry_reason(ex, req))"
                end
                return s, retry
            end)

        return retry_request(req; kw...)
    end
end

isrecoverable(e) = false
isrecoverable(e::IOError) = true
isrecoverable(e::Sockets.DNSError) = true
isrecoverable(e::HTTP.StatusError) = e.status == 403 || # Forbidden
                                     e.status == 408 || # Timeout
                                     e.status >= 500    # Server Error

isrecoverable(e, req, retry_non_idempotent, retrycount) =
    isrecoverable(e) &&
    isbytes(req.body) &&
    isbytes(req.response.body) &&
    (retry_non_idempotent || retrycount == 0 || isidempotent(req))
    # "MUST NOT automatically retry a request with a non-idempotent method"
    # https://tools.ietf.org/html/rfc7230#section-6.3.1

function no_retry_reason(ex, req)
    buf = IOBuffer()
    show(IOContext(buf, :compact => true), req)
    print(buf, ", ",
        ex isa HTTP.StatusError ? "HTTP $(ex.status): " :
        !isrecoverable(ex) ?  "$ex not recoverable, " : "",
        !isbytes(req.body) ? "request streamed, " : "",
        !isbytes(req.response.body) ? "response streamed, " : "",
        !isidempotent(req) ? "$(req.method) non-idempotent" : "")
    return String(take!(buf))
end

end # module RetryRequest


text:
using ArnoldiMethod, LinearAlgebra
using BenchmarkTools

function bench_schur(n = 30)
    H = triu(rand(n, n), -1)
    Q = Matrix{Float64}(I, n, n)
    a = @benchmark ArnoldiMethod.local_schurfact!(HH, 1, $n, QQ) setup = (HH = copy($H); QQ = copy($Q))
    b = @benchmark eigvals(HH) setup = (HH = copy($H))
    a, b
end


text:
"""
    surface(cmd0::String="", arg1=[]; kwargs...)

Reads randomly-spaced (x,y,z) triples and produces a binary grid file of gridded values z(x,y) by solving:

		(1 - T) * L (L (z)) + T * L (z) = 0

Full option list at [`surface`](http://gmt.soest.hawaii.edu/doc/latest/surface.html)

Parameters
----------

- $(GMT.opt_R)
- **I** : **inc** : -- Str or Number --

	*x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#i)
- **A** : **aspect_ratio** : -- Number --

    Aspect ratio. If desired, grid anisotropy can be added to the equations.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#a)
- **C** : **convergence** : -- Number --

	Convergence limit. Iteration is assumed to have converged when the maximum absolute change in any
	grid value is less than convergence_limit.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#c)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = surface(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#g)
- **Ll** : **lower** : -- Str or Number --

	Impose limits on the output solution. lower sets the lower bound. lower can be the name of a grid
	file with lower bound values, a fixed value, d to set to minimum input value,
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#l)
- **Lu** : **upper** : -- Str or Number --

    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#l)
- **N** : **max_iterations** : -- Number --

	Number of iterations. Iteration will cease when convergence_limit is reached or when number of
	iterations reaches max_iterations.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#n)
- **Q** : **suggest** : -- Bool or [] --

    Suggest grid dimensions which have a highly composite greatest common factor.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#q)
- **S** : **search_radius** : -- Number or Str --

    Sets the resolution of the projected grid that will be created.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#s)
- **T** : **tension** : -- Number or Str --

    Tension factor[s]. These must be between 0 and 1.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#t)
- $(GMT.opt_V)
- **Z** : **over_relaxation** : -- Str or GMTgrid --

    Over-relaxation factor. This parameter is used to accelerate the convergence; it is a number between 1 and 2.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#z)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function surface(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("surface", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_a(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :aspect_ratio])
	cmd = add_opt(cmd, 'C', d, [:C :convergence])
	cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, "Ll", d, [:Ll :lower])
	cmd = add_opt(cmd, "Lu", d, [:Ll :upper])
	cmd = add_opt(cmd, 'N', d, [:N :max_iterations])
	cmd = add_opt(cmd, 'Q', d, [:Q :suggest])
	cmd = add_opt(cmd, 'S', d, [:S :search_radius])
	cmd = add_opt(cmd, 'T', d, [:T :tension])
	cmd = add_opt(cmd, 'Z', d, [:Z :over_relaxation])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "surface", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
surface(arg1=[], cmd0::String=""; kw...) = surface(cmd0, arg1; kw...)


text:
using Bio.Seq
using BioBenchmarks

function searchall(seq, pat, k)
    query = ApproximateSearchQuery(pat)
    count = 0
    from = 0
    while (from = approxsearchindex(seq, query, k, from + 1)) > 0
        count += 1
    end
    return count
end

chrom = first(open(datafile("chr1.fa"), FASTA)).seq
@benchmark searchall(chrom, dna"ATTACGTGACGT", 2)


text:
"""
Provides a definition for a reverse mutation operator.
"""
type ReverseMutationDefinition <: MutationDefinition
  stage::AbstractString
  rate::Float

  ReverseMutationDefinition(stage::AbstractString, rate::Float) = new(stage, rate)
end

type ReverseMutation <: Mutation
  stage::AbstractString
  rate::Float
  representation::Representation

  ReverseMutation(stage::AbstractString, rate::Float, rep::Representation) =
    new(stage, rate, rep)
end

"""
Composes a reverse mutation operator from its definition.
"""
function compose!(d::ReverseMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  ReverseMutation(d.stage, d.rate, sp.stages[d.stage].representation)
end

"""
Reverses the order of all genes between two randomly selected positions, i and j,
along the chromosome.

**Parameters:**

* `stage::AbstractString`, the name of the developmental stage that this
  operator should be applied to. Defaults to the genotype if no stage is
  specified.
* `rate::Float`, the probability that this operator will perform a reverse upon
   invocation.
"""
reverse() = reverse(0.5)
reverse(rate::Float) = reverse("", rate)
reverse(stage::AbstractString) = reverse(stage, 0.5)
reverse(stage::AbstractString, rate::Float) =
  ReverseMutationDefinition(stage, rate)

"""
Possibly performs a reverse mutation operation on a provided chromosome,
with a probability given by the mutation rate of the underlying operator. If the
mutation is performed, the order of all genes within a randomly selected
sub-sequence of the given chromosome are reversed.
"""
function mutate!{T}(o::SingleSwapMutation, input::Vector{T})
  rand() < o.rate && return input

  # Perform the reverse between two randomly selected points.

  return input
end


text:
struct BFloat16 <: AbstractFloat
    val::UInt16
end

const MAGIC = Float32(1 + 2^16 * eps(Float32) / 2)
Base.convert(::Type{Float32}, x::BFloat16) = reinterpret(Float32, UInt32(x.val) << 16)
Base.convert(::Type{BFloat16}, x::Float32) = BFloat16(UInt16(reinterpret(UInt32, x * MAGIC) >> 16))
Base.convert(::Type{BFloat16}, x::Float64) = convert(BFloat16, convert(Float32, x))

BFloat16(x::Float32) = convert(BFloat16, x)

Base.zero(::Type{BFloat16}) = BFloat16(zero(UInt16))
dnnl_type(::Type{BFloat16}) = Lib.dnnl_bf16

Base.promote_rule(::Type{BFloat16}, ::Type{Float32}) = Float32

# This definition is a little weird, but should work just fine for now.
Base.:+(a::BFloat16, b::BFloat16) = convert(Float32, a) + convert(Float32, b)
Base.:-(a::BFloat16, b::BFloat16) = convert(Float32, a) - convert(Float32, b)
#Random.randn(::Type{BFloat16}) = convert(BFloat16, randn(Float32))

Base.eps(x::BFloat16) = 2^16 * convert(Float32, x)
Base.eps(::Type{BFloat16}) = 2^16 * eps(Float32)
LinearAlgebra.norm(x::BFloat16) = LinearAlgebra.norm(convert(Float32, x))
Base.abs2(x::BFloat16) = Base.abs2(convert(Float32, x))

#####
##### Extended Memory
#####

# Keeping the weights as BF16 has nice performance implications, but doesn't tend to
# converge as quickly as Float32 training.
#
# This struct will mirror the wrapped BFloat16 memory with another BFloat16 that will hold
# the lower mantissa bits for better numerical properties.
struct Mirrored{A <: AbstractArray{BFloat16}, B <: AbstractArray{UInt16}}
    base::A
    mantissa::B
end

function Mirrored(base::AbstractArray{BFloat16})
    mantissa = similar(base, UInt16)
    return Mirrored(base, mantissa)
end

function Base.parent(x::Mirrored)
    return Mirrored(parent(x.base), parent(x.mantissa))
end

Base.eachindex(x::Mirrored) = eachindex(x.base, b.mantissa)

Base.@propagate_inbounds function Base.getindex(x::Mirrored, i::Int)
    (hi, lo) = (x.base[i], x.mantissa[i])
    return reinterpret(Float32, UInt32(hi.val) << 16 + lo)
end

Base.@propagate_inbounds function Base.setindex!(x::Mirrored, v::Float32, i::Int)
    vint = reinterpret(UInt32, v)
    hi = UInt16(vint >> 16)
    lo = UInt16(vint & typemax(UInt16))
    (x.base[i], x.mantissa[i]) = (hi, lo)
end


text:
#################### No-U-Turn Sampler ####################

#################### Types and Constructors ####################

mutable struct NUTSTune <: SamplerTune
  logfgrad::Union{Function, Missing}
  adapt::Bool
  alpha::Float64
  epsilon::Float64
  epsilonbar::Float64
  gamma::Float64
  Hbar::Float64
  kappa::Float64
  m::Int
  mu::Float64
  nalpha::Int
  t0::Float64
  target::Float64
  tree_depth::Int

  NUTSTune() = new()

  function NUTSTune(x::Vector, epsilon::Real, logfgrad::Union{Function, Missing};
                    target::Real=0.6, tree_depth::Int=10)
    new(logfgrad, false, 0.0, epsilon, 1.0, 0.05, 0.0, 0.75, 0, NaN, 0, 10.0,
        target, tree_depth)
  end
end

NUTSTune(x::Vector{Float64}, logfgrad::Function, ::NullFunction; args...) =
  NUTSTune(x, nutsepsilon(x, logfgrad); args...)

NUTSTune(x::Vector{Float64}, logfgrad::Function; args...) =
  NUTSTune(x, nutsepsilon(x, logfgrad), logfgrad; args...)

NUTSTune(x::Vector, epsilon::Real; args...) =
  NUTSTune(x, epsilon, missing; args...)

const NUTSVariate = SamplerVariate{NUTSTune}

#################### Sampler Constructor ####################
"""
    NUTS(params::ElementOrVector{Symbol}; dtype::Symbol=:forward, args...)

Construct a `Sampler` object for NUTS sampling, with the algorithm’s step size
parameter adaptively tuned during burn-in iterations. Parameters are assumed
to be continuous, but may be constrained or unconstrained.

Returns a `Sampler{NUTSTune}` type object.

* `params`: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic `unlist()` function.

* `dtype` : differentiation for gradient calculations. Options are
    * `:central` : central differencing
    * `:forward` : forward differencing.

* `args...`: additional keyword arguments to be passed to the `NUTSVariate` constructor.
"""
function NUTS(params::ElementOrVector{Symbol}; dtype::Symbol=:forward, args...)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
    f = x -> logpdfgrad!(block, x, dtype)
    v = SamplerVariate(block, f, NullFunction(); args...)
    sample!(v, f, adapt=model.iter <= model.burnin)
    relist(block, v)
  end
  Sampler(params, samplerfx, NUTSTune())
end

#################### Sampling Functions ####################

sample!(v::NUTSVariate; args...) = sample!(v, v.tune.logfgrad; args...)
"""
    sample!(v::NUTSVariate, logfgrad::Function; adapt::Bool=false)

Draw one sample from a target distribution using the NUTS sampler. Parameters
are assumed to be continuous and unconstrained.

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::NUTSVariate, logfgrad::Function; adapt::Bool=false)
  tune = v.tune
  setadapt!(v, adapt)
  if tune.adapt
    tune.m += 1
    nuts_sub!(v, tune.epsilon, logfgrad)
    p = 1.0 / (tune.m + tune.t0)
    tune.Hbar = (1.0 - p) * tune.Hbar +
                p * (tune.target - tune.alpha / tune.nalpha)
    tune.epsilon = exp(tune.mu - sqrt(tune.m) * tune.Hbar / tune.gamma)
    p = tune.m^-tune.kappa
    tune.epsilonbar = exp(p * log(tune.epsilon) +
                          (1.0 - p) * log(tune.epsilonbar))
  else
    if (tune.m > 0) tune.epsilon = tune.epsilonbar end
    nuts_sub!(v, tune.epsilon, logfgrad)
  end
  v
end

function setadapt!(v::NUTSVariate, adapt::Bool)
  tune = v.tune
  if adapt && !tune.adapt
    tune.m = 0
    tune.mu = log(10.0 * tune.epsilon)
  end
  tune.adapt = adapt
  v
end

function nuts_sub!(v::NUTSVariate, epsilon::Real, logfgrad::Function)
  n = length(v)
  x, r, logf, grad = leapfrog(v.value, randn(n), zeros(n), 0.0, logfgrad)
  logp0 = logf - 0.5 * dot(r)
  logu0 = logp0 + log(rand())
  xminus = xplus = x
  rminus = rplus = r
  gradminus = gradplus = grad
  j = 0
  n = 1
  s = true
  while s || j < v.tune.tree_depth
    pm = 2 * (rand() > 0.5) - 1
    if pm == -1
      xminus, rminus, gradminus, _, _, _, xprime, nprime, sprime, alpha,
        nalpha = buildtree(xminus, rminus, gradminus, pm, j, epsilon, logfgrad,
                           logp0, logu0)
    else
      _, _, _, xplus, rplus, gradplus, xprime, nprime, sprime, alpha, nalpha =
        buildtree(xplus, rplus, gradplus, pm, j, epsilon, logfgrad, logp0,
                  logu0)
    end
    if sprime && rand() < nprime / n
      v[:] = xprime
    end
    j += 1
    n += nprime
    s = sprime && nouturn(xminus, xplus, rminus, rplus)
    v.tune.alpha, v.tune.nalpha = alpha, nalpha
  end
  v
end

function leapfrog(x::Vector{Float64}, r::Vector{Float64}, grad::Vector{Float64},
                  epsilon::Real, logfgrad::Function)
  r += (0.5 * epsilon) * grad
  x += epsilon * r
  logf, grad = logfgrad(x)
  r += (0.5 * epsilon) * grad
  x, r, logf, grad
end

function buildtree(x::Vector{Float64}, r::Vector{Float64},
                   grad::Vector{Float64}, pm::Integer, j::Integer,
                   epsilon::Real, logfgrad::Function, logp0::Real, logu0::Real)
  if j == 0
    xprime, rprime, logfprime, gradprime = leapfrog(x, r, grad, pm * epsilon,
                                                    logfgrad)
    logpprime = logfprime - 0.5 * dot(rprime)
    nprime = Int(logu0 < logpprime)
    sprime = logu0 < logpprime + 1000.0
    xminus = xplus = xprime
    rminus = rplus = rprime
    gradminus = gradplus = gradprime
    alphaprime = min(1.0, exp(logpprime - logp0))
    nalphaprime = 1
  else
    xminus, rminus, gradminus, xplus, rplus, gradplus, xprime, nprime, sprime,
      alphaprime, nalphaprime = buildtree(x, r, grad, pm, j - 1, epsilon,
                                          logfgrad, logp0, logu0)
    if sprime
      if pm == -1
        xminus, rminus, gradminus, _, _, _, xprime2, nprime2, sprime2,
          alphaprime2, nalphaprime2 = buildtree(xminus, rminus, gradminus, pm,
                                                j - 1, epsilon, logfgrad, logp0,
                                                logu0)
      else
        _, _, _, xplus, rplus, gradplus, xprime2, nprime2, sprime2,
          alphaprime2, nalphaprime2 = buildtree(xplus, rplus, gradplus, pm,
                                                j - 1, epsilon, logfgrad, logp0,
                                                logu0)
      end
      if rand() < nprime2 / (nprime + nprime2)
        xprime = xprime2
      end
      nprime += nprime2
      sprime = sprime2 && nouturn(xminus, xplus, rminus, rplus)
      alphaprime += alphaprime2
      nalphaprime += nalphaprime2
    end
  end
  xminus, rminus, gradminus, xplus, rplus, gradplus, xprime, nprime, sprime,
    alphaprime, nalphaprime
end

function nouturn(xminus::Vector{Float64}, xplus::Vector{Float64},
                 rminus::Vector{Float64}, rplus::Vector{Float64})
  xdiff = xplus - xminus
  dot(xdiff, rminus) >= 0 && dot(xdiff, rplus) >= 0
end

#################### Auxilliary Functions ####################

function nutsepsilon(x::Vector{Float64}, logfgrad::Function)
  n = length(x)
  _, r0, logf0, grad0 = leapfrog(x, randn(n), zeros(n), 0.0, logfgrad)
  epsilon = 1.0
  _, rprime, logfprime, gradprime = leapfrog(x, r0, grad0, epsilon, logfgrad)
  prob = exp(logfprime - logf0 - 0.5 * (dot(rprime) - dot(r0)))
  pm = 2 * (prob > 0.5) - 1
  while prob^pm > 0.5^pm
    epsilon *= 2.0^pm
    _, rprime, logfprime, _ = leapfrog(x, r0, grad0, epsilon, logfgrad)
    prob = exp(logfprime - logf0 - 0.5 * (dot(rprime) - dot(r0)))
  end
  epsilon
end


text:
__precompile__()

module Torch

export THArray

include("api.jl")
include("extras.jl")

function __init__()
  Libdl.dlopen(joinpath(@__DIR__, "..", "deps", "libTH.dylib"))
end

end # module


text:
using Turing, Test
import Random

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@stage_testset "dynamichmc" "dynamichmc.jl" begin
    import DynamicHMC
    Random.seed!(100)

    @test DynamicPPL.alg_str(Sampler(DynamicNUTS(), gdemo_default)) == "DynamicNUTS"

    chn = sample(gdemo_default, DynamicNUTS(), 10_000)
    check_gdemo(chn)

    chn2 = sample(gdemo_default, Gibbs(PG(15, :s), DynamicNUTS(:m)), 10_000)
    check_gdemo(chn2)

    chn3 = sample(gdemo_default, Gibbs(DynamicNUTS(:s), ESS(:m)), 10_000)
    check_gdemo(chn3)
end


text:
## Case study from batch run
ENV["GKS_ENCODING"] = "utf-8"
using EstimateSDE
using Serialization
using Pipe: @pipe
using Plots
using DataFrames
using Statistics
using StatsPlots
using DataStructures

## Load settings and aggregated results
results_dir = settings_resultsdir((name="experiment01",))
S = deserialize(joinpath(results_dir,"settings.jld"))
df_eval = deserialize(joinpath(results_dir,"model-evaluation.jld")).df
df_data = deserialize(joinpath(results_dir,"data.jld")).df
df_eval = innerjoin(df_eval,df_data,on=:index)

df_kld_mean = deserialize(joinpath(results_dir,"aggregate.jld")).df

report_dir = joinpath(results_dir,"results")

## Select case
df = filter(x->x.proclabel==:ICMLA_C&&x.method=="MLE",df_eval)
index = df.index[end]
df = filter(x->x.index==index,df)
display(select(df,[:proclabel,:N,:engine,:kld]))

## Data
row_data = filter(x->x.index==index,df_data) |> first
data = row_data.data
t,y = (getfield.(data,x) for x in (:t,:y))
plot(t,y,ms=2,markershape=:square,xlabel="t",ylabel="y",legend=false)
plot!(title=string("Process ",row_data.proclabel," Tav = ","N = ",length(data)))

## Plot traces
Tkld = S.mle_trace.Tkld
proc = testprocess(row_data.proclabel)
duration = row_data.N*row_data.Tav
t_kld = (0:Tkld:duration)

traces = map(eachrow(df)) do x
  x.engine=>x.mle_trace
end |> OrderedDict

for (k,trace) in traces
  plts = reporttrace(data,trace,Tkld;proc=proc,t_kld=t_kld)
  plot!(plts.fit_order_wide,title=k)
end

## Estimation rerun
# ENV["JULIA_DEBUG"] = EstimateSDE
# r = filter(x->x.engine==ElimQ,df) |> first
# S.mle_trace
# spec = (engine=r.engine,init_actual=r.init_actual)
# trace = batch_mle_trace_single(S.mle_trace,spec,row_data,throw_errors=false)
# plts = reporttrace(data,trace,Tkld,proc=proc,t_kld=t_kld)
# plot!(plts.fit_order_wide,title="ElimQ rerun")
# nothing


text:
# These rule provide basic promotion and conversion rules to used Dual numbers with the rest of numbers.

function Base.promote_rule(::Type{ForwardDiff.Dual{T}}, ::Type{ForwardDiff.Dual{R}}) where {T,R}
    ForwardDiff.Dual{promote_type(T,R)}
  end

function Base.promote_rule(::Type{ForwardDiff.Dual{T}}, ::Type{R}) where {T,R}
    ForwardDiff.Dual{promote_type(T,R)}
end

function Base.convert(::Type{ForwardDiff.Dual{T}}, x::ForwardDiff.Dual) where T
    ForwardDiff.Dual(convert(T, x.x), convert(T, x.ϵ))
  end

function Base.convert(::Type{ForwardDiff.Dual{T}}, x::Real) where T
    ForwardDiff.Dual(convert(T, x), zero(T))
end


text:
export @timedone, fasthash, cache, dictcache
export loadedmodules, reloadmodules

macro timedone(a, ex)
    quote
        local t = time()
        println($a, " ... ")
        flush(stdout)
        r = $ex
        local elapsed = round((time()-t)*100)/100
        println("    done, $elapsed sec!")
        flush(stdout)
        r
    end
end

# OPTIM remove collect:
fasthash(a::Array{T}) where {T<:Number} = bytes2hex(eltype(a)!=UInt8 ? sha256(collect(reinterpret(UInt8, vec(a)))) : sha256(a))

fasthash(a::Number) = fasthash(Any[[a], "__julianumber"])
fasthash(a::Char) = fasthash(Any[string(a), "__juliachar"])
fasthash(a::AbstractString) = bytes2hex(sha256(a))
fasthash(a::Symbol) = fasthash(Any[string(a),"__juliasymbol"])
fasthash(a::Tuple) = fasthash(Any[a..., "__juliatuple"])
fasthash(f::Function) = fasthash(Any[string(f), "__juliafunction"])
fasthash(a::UnitRange) = fasthash(Any[a.start, a.stop, "__juliaunitrange"])
fasthash(a::AbstractRange) = fasthash(Any[a.start, a.step, a.stop, "__juliarange"])
fasthash(a::StepRangeLen) = fasthash(Any[a.ref, a.step, a.len, a.offset, "__juliafloatrange"])
fasthash(a::Bool) = fasthash("$(a)_juliabool")
fasthash(a) = fasthash(repr(a))

function fasthash(a::Array)
    if length(a)>100000 && isa(a, DenseArray) && eltype(a)<:Number
        fasthash(Any[a[1:911:end],sum(a)])
    else
        r = Base.map(fasthash, a)
        fasthash(join(r))
    end
end
function fasthash(a::Dict)
    K = collect(keys(a))
    ind = sortperm(Base.map(string, K))
    return fasthash([(k, fasthash(a[k])) for k in K[ind]])
end
fasthash(f::Function, args, version) = fasthash(Any[string(f), args..., version])

cache(f::FD.Callable, f2::FD.Callable, args...) = error("cache(f,f2,...) not supported")
function cache(f::FD.Callable, args...; version = "0 none set", kargs...)
    mkpath("zzz")
    h = hash(fasthash([f; args; version; kargs]))
    filename = "zzz/$(string(f))-$h.jls"

    prefix = ""
    if existsfile(filename)
        try
            return @timedone "reloading $(string(f))" open(deserialize, filename)
        catch
        end
        println("  reloading failed")
        prefix = "re"
    end

    @timedone "$(prefix)computing $(string(f))" begin
        r = f(args...; kargs...)
        open(s->serialize(s,r),filename,"w")
        r
    end
end
cache(a,f::FD.Callable,args...; kargs...) = cache(f, Any[a, args...]...; kargs...)

function dictcache(f, args...; version = "0 none set", filepath = "cache.dictfile")
    dictopen(filepath) do cache
        h = fhash(f, args, version)
        key = tuple(string(f), version, h)
        #@show "before" key
        if haskey(cache, key...)
             @timedone "reloading $(string(f))" cache[key...]
        else
            @timedone "computing $(string(f))" begin
              r = f(args...)
              cache[key...] = r
              #@show keys(cache, key[1:end-1])
              @assert in(key[end], keys(cache, key[1:end-1]...))
              r
            end
        end
        #@show "after" key keys(cache, string(f), version)
    end
end

function loadedmodules()
    a = names(Main)
    r = Any[]
    for x in a
        #try
            if typeof(eval(Main,x))==Module && !in(x,Any[:Main, :Core])
                push!(r, fullname(eval(Main,x))[1])
            end
        #catch e
        #    Base.display_error(e)
        #end
    end
    r
end

function reloadmodules()
    for x in vcat(:__anon__, sort(loadedmodules())) # HACK to make sure HDF5 is loaded before JLD
        if in(x, Any[:MAT_HDF5, :MAT_v5, :IPythonDisplay, :IJulia, :ReactiveContexts])
            try
				a = "@everywhere using $x"
				println("trying: $a")
                eval(Main, parse(a))
            catch e
				showerror(e)
            end
        end
    end
end


text:
module MIMEBundles

export MIMEBundle, mimerepr, stringmime

using JSON

include("./mimerepr.jl")
include("./mimes.jl")

struct MIMEBundle
    data::Dict{MIME, Union{String, JSONText}}
    metadata::Dict{Any, Any}
end

function MIMEBundle(obj)
    data = Dict{MIME, Any}()
    metadata = Dict()
    for mime in MIMEOrVec[DEFAULT_MIMES..., extramimes(obj)...]
        mime, result = mimerepr(mime, obj)
        if result !== nothing
            if isjsonmime(mime)
                result = JSONText(result)
            end
            data[mime] = result
        end
    end
    return MIMEBundle(data, metadata)
end

end # module


text:
#####
##### Get counters and timing numbers
#####

const DATAFILE = joinpath(@__DIR__, "..", "data.jls")

# A lot of this is copied directly from the `CountOnMe` project repo.
# TODO: Can we figure out an API for this to roll it into CounterTools for making
# these kinds of measurements a lot easier??
getdb() = ispath(DATAFILE) ? ExperimentsDB.load(DATAFILE) : ExperimentsDB.Database()

# Alow for debugging
@static if DEBUG
    save(x) = nothing
else
   save(x) = ExperimentsDB.save(DATAFILE, x)
end

# By default, we run on Socket 1
default_cpuset() = 25:48
default_socket() = CounterTools.IndexZero(1)
default_cha_sample_cpu() = CounterTools.IndexZero(72)
daemon_port() = 2000

const CoreSelectRegister = CounterTools.CoreSelectRegister
const UncoreSelectRegister = CounterTools.UncoreSelectRegister

function dictpush!(d, x, v::T) where {T}
    vec = get!(d, x, T[])
    push!(vec, v)
    return nothing
end

#####
##### Parameters for the benchmarks we are running.
#####

Base.@kwdef struct BenchmarkParameters
    # Dimensions of the internal arrays
    arraydims::Tuple
    totalsize::Int
    flushpercent::Float32
    eviction_policy::String
    # sampletime for measurements
    sampletime::TimePeriod
    # Brief description of the benchmark
    description::String
    kernel::String
    iterations::Int
    # number of threads
    num_threads::Int
    mode::String
end

# Convert this parameter to a dict.
todict(x::BenchmarkParameters) = SortedDict(f => getproperty(x, f) for f in fieldnames(typeof(x)))

function imcmeasurements(params, data::Dict{String,Any}; applyfilter = true)
    eventlist = [
        "CAS Count Read"    => UncoreSelectRegister(; event = 0x04, umask = 0x03),
        "CAS Count Write"   => UncoreSelectRegister(; event = 0x04, umask = 0x0C),
        "PMM Read"          => UncoreSelectRegister(; event = 0xE3),
        "PMM Write"         => UncoreSelectRegister(; event = 0xE7),
        #"Tag Hit"           => UncoreSelectRegister(; event = 0xD3, umask = 0x01),
        #"Tag Miss Clean"    => UncoreSelectRegister(; event = 0xD3, umask = 0x02),
        #"Tag Miss Dirty"    => UncoreSelectRegister(; event = 0xD4, umask = 0x04),
        #"Uncore Clocks"     => UncoreSelectRegister(; event = 0x00, umask = 0x00),
    ]

    # Optionally filter
    if applyfilter
        filter!(x -> !haskey(data, makekey(x)), eventlist)
    end

    # wrap this up in an iterators.partition to only yield a maximum of 4 values at a time.
    # only 4 core counters are programmable at a time.
    numcounters = 4
    socket = default_socket()
    return map(Iterators.partition(eventlist, numcounters)) do x
        keys = makekey.(x)
        events = last.(x)
        return keys => MattDaemon.@measurement CounterTools.IMCMonitor(Tuple(events), socket)
    end
end

#####
##### Running and Recording Logic
#####

# Use a `maybeiterate` function to
# - Iterate if the iterate if the iterator isn't done.
# - Do nothing if it IS done.
maybeiterate(itr, state::Tuple) = iterate(itr, last(state))
maybeiterate(itr, ::Nothing) = nothing

maybefirst(state::Tuple) = first(state)
maybefirst(::Nothing) = nothing

_first(x::AbstractString) = x
_first(x::Pair) = first(x)

makekey(x...) = join(_first.(x), " - ")

function runner(@nospecialize(f), params::BenchmarkParameters; force = false)
    # Now, grab the database.
    # See if we have entries for this.
    db = getdb()
    paramsdict = todict(params)

    printstyled(stdout, "Current Parameters\n"; color = :cyan)
    for (k,v) in paramsdict
        println("    ", k, " => ", v)
    end

    entry = ExperimentsDB.findonly(db, paramsdict)
    if isnothing(entry)
        println("Entry Does not yet exist in Database")
        data = Dict{String,Any}()
        db[paramsdict] = data
    else
        println("Found Entry in Database!")
        data = entry[]
    end

    # This may look overly compilcated, but that's because it's designed to be
    # quickly augmented with `core` and `cha` measurements if desired.
    imc = imcmeasurements(params, data; applyfilter = !force)
    iters = (imc,)
    states = iterate.(iters)
    order = (:imc,)
    timestamp = MattDaemon.@measurement SystemSnoop.Timestamp()

    # Iterate until all states are exhausted.
    firstiter = true
    while !all(isnothing, states)
        # Run the function once to trigger compilation
        if firstiter
            f()
            firstiter = false
        end

        measurements = (pretime = timestamp,)
        # Create a measurements tuple.
        #
        # We begin by unpacking the iterator tuples to get actual values from these iterators.
        values = maybefirst.(states)

        for (i, v) in enumerate(values)
            # If this iterator is empty, don't apply it to any measurement.
            isnothing(v) && continue

            # Create a Named tuple using this symbol
            symbol = order[i]
            keys = first(v)
            measurement = last(v)

            nt = NamedTuple{(symbol,)}((measurement,))
            measurements = merge(measurements, nt)

            # Print out a helpful message
            name = uppercase(String(symbol))
            printstyled(stdout, "Running $name Measurements: "; color = :green)
            println(stdout, keys)
        end

        # Apply one final timestamp
        measurements = merge(measurements, (posttime = timestamp,))

        # Construct a payload and launch the function.
        payload = MattDaemon.ServerPayload(params.sampletime, measurements)
        @time run_data, _, runtime = MattDaemon.run(
            f,
            payload,
            daemon_port();
            sleeptime = Second(5)
        )

        @show runtime
        dictpush!(data, "Runtimes", runtime)
        # Print some statistics on measurement times.
        measurement_times = Dates.value.(Dates.Millisecond.(run_data.posttime .- run_data.pretime))
        printstyled(stdout, "Measurement Times\n"; color = :yellow)
        println("    Minimum: ", minimum(measurement_times), " (ms)")
        println("    Maximum: ", maximum(measurement_times), " (ms)")
        println("    Average: ", mean(measurement_times), " (ms)")

        #####
        ##### Post processing
        #####

        for (i,v) in enumerate(values)
            isnothing(v) && continue

            symbol = order[i]
            keys = first(v)

            counter_values = CounterTools.aggregate.(diff(getproperty(run_data, symbol)))
            for (i, m) in enumerate(keys)
                data[makekey(m)] = getindex.(counter_values, i)
            end
        end

        states = maybeiterate.(iters, states)
        save(db)
    end
    return nothing
end

#####
##### The actual entry points
#####

_evict_string(::CachedArrays.Select) = "Select"
_evict_string(::CachedArrays.RandomPolicy) = "Random"

function tests_1d(totalsize, arraysize, flushpercent, eviction_policy)
    println("Booting")
    # if in 2LM mode, all of the temporary arrays should go into local memory
    if IS_2LM
        # Round up for headroom
        localsize = round(Int, totalsize * 1.5)
        remotesize = 2^24
    # Size local mamory to be a little smaller than the Local cache to allow for other
    # program activities that don't belong to our heap.
    else
        localsize = DEBUG ? totalsize >> 1 : 180_000_000_000
        remotesize = DEBUG ? 2 * totalsize : 1_000_000_000_000
    end

    println("Making Manager")
    manager = CachedArrays.CacheManager(
        ENV["JULIA_PMEM_PATH"];
        localsize = localsize,
        remotesize = remotesize,
        # Setup tunables
        policy = eviction_policy,
        flushpercent = convert(Float32, flushpercent),
        gc_before_evict = false,
    )
    println("Manager Made")

    arrays = alloc_1d(manager, totalsize, arraysize)
    println("Allocated Arrays")
    mode = IS_2LM ?  "2LM" : "1LM"

    # If we have a small total size, do more iterations.
    if totalsize < 180_000_000_000
        iterations = 4
    else
        iterations = 1
    end

    #####
    ##### Tests to Run
    #####

    # Sequential Copy
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),
        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Sequential Copy",
        kernel = string(sequential_copy),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> sequential_copy(arrays; iterations = iterations)
    runner(f, params)

    # Sequential Accumulate
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),
        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Sequential Accumulate",
        kernel = string(sequential_accumulation),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> sequential_accumulation(arrays; iterations = iterations)
    runner(f, params)

    # Sequential zero
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),
        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Sequential Zero",
        kernel = string(sequential_zero),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> sequential_zero(arrays; iterations = iterations)
    runner(f, params)

    # Sequential Bump
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),
        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Sequential Bump",
        kernel = string(sequential_bump),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> sequential_bump(arrays)
    runner(f, params)
    return nothing
end

function alloc_tests(totalsize, arraysize, flushpercent, eviction_policy)
    # if in 2LM mode, all of the temporary arrays should go into local memory
    if IS_2LM
        # Round up for headroom
        localsize = round(Int, totalsize * 1.5)
        remotesize = 4096
    # Size local mamory to be a little smaller than the Local cache to allow for other
    # program activities that don't belong to our heap.
    else
        localsize = DEBUG ? totalsize >> 1 : 180_000_000_000
        remotesize = DEBUG ? 2 * totalsize : 500_000_000_000
    end

    manager = CachedArrays.CacheManager(
        ENV["JULIA_PMEM_PATH"];
        localsize = localsize,
        remotesize = remotesize,
        # Setup tunables
        policy = eviction_policy,
        flushpercent = convert(Float32, flushpercent),
        gc_before_evict = true,
    )

    mode = IS_2LM ?  "2LM" : "1LM"

    # If we have a small total size, do more iterations.
    if totalsize < 180_000_000_000
        iterations = 4
    else
        iterations = 1
    end

    numelements = div(arraysize, sizeof(Float32))

    # Sequential zero
    params = BenchmarkParameters(;
        arraydims = (numelements,),
        totalsize = totalsize,
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),
        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Alloc and Free",
        kernel = string(alloc_and_dealloc),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> alloc_and_dealloc(manager, totalsize, arraysize)

    # Run once to warm up.
    @show f()
    runner(f, params)
end

function tests_2d(totalsize, arraysize, flushpercent, eviction_policy)
    if IS_2LM
        localsize = round(Int, 1.5 * totalsize)
        remotesize = 2^24
    else
        localsize = DEBUG ? totalsize >> 1 : 180_000_000_000
        remotesize = DEBUG ? 2 * totalsize : 1_000_000_000_000
    end

    manager = CachedArrays.CacheManager(
        ENV["JULIA_PMEM_PATH"];
        localsize = localsize,
        remotesize = remotesize,
        # Setup tunables
        policy = eviction_policy,
        flushpercent = convert(Float32, flushpercent),
        gc_before_evict = true,
    )

    # Setup and get arrays
    arrays = alloc_2d(manager, totalsize, arraysize)
    mode = IS_2LM ?  "2LM" : "1LM"

    #####
    ##### Tests to Run
    #####

    iterations = 1

    # Matrix Mult Forward
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),

        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Matrix MarixMultForward",
        kernel = string(matrix_mult_forward),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> matrix_mult_forward(arrays; iterations = iterations)
    runner(f, params; force = true)

    # Matrix Mult Forward and Backward
    params = BenchmarkParameters(;
        arraydims = size(first(arrays)),
        totalsize = sum(sizeof, arrays),
        flushpercent = flushpercent,
        eviction_policy = _evict_string(eviction_policy),

        # Sample twice per second
        sampletime = Millisecond(500),
        description = "Matrix Mult Forward and Back",
        kernel = string(matrix_mult_forward_and_back),
        iterations = iterations,
        mode = mode,
        num_threads = length(default_cpuset())
    )

    f = () -> matrix_mult_forward_and_back(arrays; iterations = iterations)
    runner(f, params; force = true)
end


text:
# Required packages:
using MATLAB
import IterativeSolvers
using ArgParse
using Iterators

immutable SolveInfo
    x :: Vector
    st_res :: Float64
    residuals :: Vector
    end_res :: Float64
    cyclecomplexity :: Float64
    iters :: Int
    setup_time :: Float64
    solve_time :: Float64
end

"""
Work per digit of accuracy

wda(cycle complexity, change in residual, # of cycles)
"""
wda(cc, dr, cycles) = -1 / log10(ecf(cc, dr, cycles))

"""
Effective convergence factor

ecf(cycle complexity, change in residual, # of cycles)
"""
ecf(cc, dr, cycles) = dr ^ (1/(cc*cycles))

# relaxation is one of sgs, gs, jacobi
# cycleIndex is 1.5 for W, 1 for V
# aggregation is one of affinity-mex, affinity-energy-mex
# minres is true or false
function solve_lamg_matlab(L, b, relaxation, cycle_index, aggregation, minres, tol=1e-8)
    put_variable(:L, L)
    put_variable(:b, b)
    eval_string("""tic;
                   solver = Solvers.newSolver('lamg', 'relaxType', '$relaxation', 'cycleIndex', $cycle_index, 'aggregationUpdate', '$aggregation', 'minRes', $minres);
                   setup = solver.setup('laplacian', L)
                   tsetup = toc;
                   tic;
                   [x, ~, res, details] = solver.solve(setup, b, 'errorReductionTol', $tol);
                   tsolve = toc;
                   """)
    cc = mat"setup.cycleComplexity"
    residuals = mat"details.errorNormHistory"
    iters = mat"length(details.errorNormHistory)"
    tsetup = get_variable(:tsetup)
    tsolve = get_variable(:tsolve)
    x = get_variable(:x)

    SolveInfo(x, norm(b), residuals, norm(b - L*x), cc, iters, tsetup, tsolve)
end

function solve_cmg_matlab(L, b, tol=1e-8)
    put_variable(:L, L)
    put_variable(:b, b)
    eval_string("""tic;
                   solver = Solvers.newSolver('cmg');
                   setup = solver.setup('laplacian', L);
                   tsetup = toc;
                   tic;
                   [x, ~, res, details] = solver.solve(setup, b, 'errorReductionTol', $tol);
                   tsolve = toc;
                   """)
    cc = mat"setup.cycleComplexity"
    residuals = mat"details.errorNormHistory"
    iters = mat"length(details.errorNormHistory)"
    tsetup = get_variable(:tsetup)
    tsolve = get_variable(:tsolve)
    x = get_variable(:x)

    SolveInfo(x, norm(b), residuals, norm(b - L*x), cc, iters, tsetup, tsolve)
end

function solve_jlamg(L, b, tol=1e-8)
    h = build_hierarchy(L)
    x, ci = solve_cg(h, zeros(b), b)
    SolveInfo(x, norm(b), ci.residuals, norm(b - L * x), ci.iteration_complexity, ci.iterations, 0, 0)
end

function solve_augtree(L, b, tol=1e-8)
    f = augTreeLapPrecon(L)
    fact_nnz = nnz(f.subSolver.contents.solvers[1]) # I assume this is the correct way to measure work
    (x, ch) = IterativeSolvers.cg(L, b, Pl=f, tol=tol, maxiter=500, log=true)

    SolveInfo(x, norm(b), ch.data[:resnorm], norm(b - L*x), fact_nnz/nnz(L) + 1, ch.iters)
end

function test_solvers(L, dump_file)
    # construct rhs
    b = rand(size(L,1))
    b = b - mean(b)

    relaxation = ["sgs", "gs", "jacobi"]
    cycle_index = [1.0, 1.5]
    aggregation = [#="affinity",=# "affinity-energy-mex"]
    minres = [true, false]

    #= lamg = solve_lamg_matlab(L, b) =#
    #= cmg = solve_cmg_matlab(L, b) =#
    #= jlamg = solve_jlamg(L, b) =#
    #= augtree = solve_augtree(L, b) =#

    println("    method iters       ∆r    work    wda    setup    solve")
    for (n,f) in [ ("lamg", solve_lamg_matlab)
                 ]
        #
        for (relax, ci, agg, mr) in product(relaxation, cycle_index, aggregation, minres)
            s = f(L, b, relax, ci, agg, mr)
            @printf "%s %f %s %s\n" relax ci agg mr
            @printf "%10s %5d %5.2e %7.2f %6.2f %5.2e %5.2e\n" n s.iters s.end_res/s.st_res s.cyclecomplexity*s.iters wda(s.cyclecomplexity, s.end_res/s.st_res, s.iters) s.setup_time s.solve_time

            if !isa(dump_file, Void)
                s = @sprintf "%s,%d,%e,%f,%f,%e,%e,%s,%f,%s,%s\n" n s.iters s.end_res/s.st_res s.cyclecomplexity*s.iters wda(s.cyclecomplexity, s.end_res/s.st_res, s.iters) s.setup_time s.solve_time relax ci agg mr
                write(dump_file, s)
            end
        end
    end
end

# replace this with your matrix loading code
using MAT
L=matread(expanduser(ARGS[1]))["L"]

s = ArgParseSettings()
@add_arg_table s begin
    "--dump"
        arg_type = String
        help = "Dump outputs to csv file"
    "MAT"
        help = "Graph Laplacian file"
        required = true
end

parsed_args = parse_args(ARGS, s)

dump_file = nothing
if !isa(parsed_args["dump"], Void)
    dump_file = open(parsed_args["dump"], "w")
    write(dump_file, "solver,iters,rel_error,work,wda,setup_time,solve_time,relaxation,cycle_index,aggregation,minres\n")
end

# L should be a graph lacplacian (not an adjacency matrix)
test_solvers(L, dump_file)

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
# Unit tests for judiRHS and judiWavefield (without PDE solves)
# Philipp Witte (pwitte.slim@gmail.com)
# May 2018
#

using JUDI.TimeModeling, SegyIO, Test, LinearAlgebra

# Example structures

example_info(; n=(120,100), nsrc=2, ntComp=1000) = Info(prod(n), nsrc, ntComp)
example_model(; n=(120,100), d=(10f0, 10f0), o=(0f0, 0f0), m=randn(Float32, n)) = Model(n, d, o, m)

function example_rec_geometry(; nsrc=2, nrec=120)
    xrec = range(50f0, stop=1150f0, length=nrec)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nrec)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

function example_src_geometry(; nsrc=2)
    xrec = range(100f0, stop=1000f0, length=nsrc)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nsrc)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

########################################################### judiRHS ####################################################

@testset "judiRHS Unit test" begin

    # Constructor
    nsrc = 2
    info = example_info(nsrc=nsrc)
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    data = Array{Array}(undef, nsrc)
    for j=1:nsrc
        data[j] = randn(Float32, rec_geometry.nt[j], length(rec_geometry.xloc[j]))
    end
    rhs = judiRHS(info, rec_geometry, data)

    @test isequal(typeof(rhs), judiRHS{Float32})
    @test isequal(rhs.geometry, rec_geometry)

    # conj, transpose, adjoint
    @test isequal(size(rhs), size(conj(rhs)))
    @test isequal(reverse(size(rhs)), size(transpose(rhs)))
    @test isequal(reverse(size(rhs)), size(adjoint(rhs)))

    # +, -
    info = example_info(nsrc=nsrc)
    rec_geometry = example_rec_geometry(nsrc=nsrc)
    src_geometry = example_src_geometry(nsrc=nsrc)
    data1 = Array{Array}(undef, nsrc)
    data2 = Array{Array}(undef, nsrc)
    for j=1:nsrc
        data1[j] = randn(Float32, rec_geometry.nt[j], length(rec_geometry.xloc[j]))
        data2[j] = randn(Float32, src_geometry.nt[j], length(src_geometry.xloc[j]))
    end
    rhs1 = judiRHS(info, rec_geometry, data1)
    rhs2 = judiRHS(info, src_geometry, data2)

    rhs_sum = rhs1 + rhs2
    rhs_sub = rhs1 - rhs2

    @test isequal(size(rhs_sum), size(rhs1))
    @test isequal(size(rhs_sub), size(rhs1))

    @test isequal(length(rhs_sum.geometry.xloc[1]), length(rhs1.geometry.xloc[1]) + length(rhs2.geometry.xloc[1]))
    @test isequal(length(rhs_sub.geometry.xloc[1]), length(rhs1.geometry.xloc[1]) + length(rhs2.geometry.xloc[1]))

    @test isequal(size(rhs_sum.data[1])[2], size(rhs1.data[1])[2] + size(rhs2.data[1])[2])
    @test isequal(size(rhs_sub.data[1])[2], size(rhs1.data[1])[2] + size(rhs2.data[1])[2])

    # get index
    rhs_sub = rhs[1]
    @test isequal(rhs_sub.info.nsrc, 1)
    @test isequal(typeof(rhs_sub.geometry), GeometryIC)
    @test isequal(typeof(rhs.data), Array{Array, 1})
    @test isequal(length(rhs_sub), Int(length(rhs)/2))

    rhs_sub = rhs[1:2]
    @test isequal(rhs_sub.info.nsrc, 2)
    @test isequal(typeof(rhs_sub.geometry), GeometryIC)
    @test isequal(typeof(rhs.data), Array{Array, 1})
    @test isequal(length(rhs_sub), length(rhs))

end

######################################################### judiWavefield ################################################


text:
@testset "Test API datatypes" begin
    test_function(x) = x

    # Test type Observable:
    obs1 = @inferred Observable(test_function)
    @test obs1.func == test_function
    @test obs1.min == -Inf
    @test obs1.max == Inf

    obs2 = @inferred Observable(test_function, min=-5, max=5)
    @test obs2.func == test_function
    @test obs2.min == -5
    @test obs2.max == 5

    # Test type Measurement:
    meas1 = @inferred Measurement(test_function, 50.0, uncertainties=(unc1=0.5, unc2=0.3), active=true)
    meas2 = @inferred Measurement(obs1, 50.0, uncertainties=(unc1=0.5, unc2=0.3), active=true)
    @test meas1 == meas2
    @test keys(meas1.uncertainties) == (:unc1, :unc2)

    #Test type MeasurementDistribution:
    measdist1 = MeasurementDistribution([obs1, obs1], [50.0, 49.9],
                uncertainties=(unc1=[0.5, 0.4], unc2=[0.3, 0.29]), active=[true, false])
    measdist2 = MeasurementDistribution([obs1, obs1], [50.0, 49.9],
                uncertainties=(unc1=[0.5, 0.4], unc2=[0.3, 0.29]), active=true)
    @test measdist1.active == [true, false]
    @test measdist2.active == [true, true]

    #Test type Correlation:
    corr1 = Correlation([1 0.5; 0.5 1])
    corr2 = Correlation([1 0.5; 0.5 1], active=false)
    @test corr1.active == true
    @test corr2.active == false
    @test corr1.matrix ==  [1 0.5; 0.5 1]

    #Test type NuisanceCorrelation:
    nui1 = @inferred NuisanceCorrelation(:stat, :meas1, :meas2, Normal(0, 1))
    @test nui1.unc_key == :stat
    @test nui1.meas1 == :meas1
    @test nui1.meas2 == :meas2
    @test nui1.prior == Normal(0, 1)

    nui2 = NuisanceCorrelation(:stat, :meas1, :meas2, -1..1)
    @test nui2.prior == Uniform(-1, 1)

    nui3 = NuisanceCorrelation(:stat, :meas1, :meas2, 0)
    @test nui3.prior == BAT.ConstValueDist(0)
end


text:
function isotropic_turbulence(spectrum, p)

    # --> Get the simulation parameters.
    @unpack Lx, Ly, nx, ny = p

    # --> Wavenumbers.
    α = fftfreq(nx, nx/Lx) * 2π
    β = fftfreq(ny, ny/Ly) * 2π

    # --> Random Fourier modes with unit energy.
    ω = randn(ny, nx)
    ω .-= mean(ω)

    ω = fft(ω)
    ω ./= abs.(ω)

    # --> Apply the isotropic enstrophy spectrum
    for j = 1:nx, i = 1:ny
        # --> Wavenumber
        k = √(α[j]^2 + β[i]^2)

        # --> Rescale the Fourier mode.
        ω[i, j] = i == j == 1 ? zero(ω[i, j]) : ω[i, j] * √spectrum(k)
    end

    # --> Zero-out the last wavenumber.
    ω[:, nx÷2+1] .= 0
    ω[ny÷2+1, :] .= 0

    # --> Scale the vorticity field to the desired Reynolds number.
    ω = scale_by_re(ω, p)

    return ω
end

function scale_by_re(ω, p)

    # --> Simulation parameters.
    @unpack Lx, Ly, nx, ny = p

    # --> Rescale the vorticity field to enforce the desired Re.
    u, v = compute_velocity(ω, p)

    E = mean(u.^2 + v.^2)

    ω /= √E * 2π

    return ω
end

function isolated_vortex(x, y ; x₀=0.0, y₀=0.0, Γ=1.0, r=1.0)

    # --> Center the mesh around the (x₀, y₀) position.
    Lx = maximum(abs.(x))
    Ly = maximum(abs.(y))
    x, y = collect(x .- x₀), collect(y .- y₀)

    # --> Enforce periodicity.
    x[x .> Lx] .= x[x .> Lx] .- 2Lx
    x[x .< -Lx] .= x[x .< -Lx] .+ 2Lx

    y[y .> Ly] .= y[y .> Ly] .- 2Ly
    y[y .< -Ly] .= y[y .< -Ly] .+ 2Ly

    # --> Isolated vortex vorticity field.
    r² = (x').^2 .+ y.^2
    ω = Γ/(π*r^2) .* exp.(-r²./r^2)

    # --> Map to spectral space.
    ω = fft(ω)

    return ω
end

function dipole(p)

    # --> Unpack parameters.
    @unpack Lx, Ly, nx, ny = p

    # --> Get mesh.
    x, y = mesh(p)

    # --> Dipole.
    ω = isolated_vortex(x, y ; Γ=1.0, r=0.5, y₀=0.5)
    ω += isolated_vortex(x, y ; Γ=-1.0, r=0.5, y₀=-0.5)

    # --> Remove the mean.
    ω[1, 1] = 0.0

    # --> Scale vorticity field to the desired Reynolds number.
    ω = scale_by_re(ω, p)

    return ω
end

function taylor_green(p)

    # --> Construct the mesh.
    x, y = mesh(p)

    # --> Taylor-Green vortices.
    ω = -2 .* cos.(x') .* cos.(y)

    # --> Map to spectral space.
    ω = fft(ω)

    return ω
end


text:
using ModelingToolkit
using LightGraphs
using DiffEqBase
using Test
using UnPack

# Define some variables
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)
lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)

lowered_eqs = [D(xˍt) ~ T*x,
               D(yˍt) ~ T*y - g,
               D(x) ~ xˍt,
               D(y) ~ yˍt,
               0 ~ x^2 + y^2 - L^2,]
@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g]) == lowered_sys
@test isequal(equations(lowered_sys), lowered_eqs)

# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
       D(y) ~ z,
       D(w) ~ T*x,
       D(z) ~ T*y - g,
       0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)

pendulum = initialize_system_structure(pendulum)
sss = structure(pendulum)
@unpack graph, fullvars, varassoc = sss
@test StructuralTransformations.matching(sss, varassoc .== 0) == map(x -> x == 0 ? StructuralTransformations.UNASSIGNED : x, [0, 0, 0, 0, 1, 2, 3, 4, 0])

sys, assign, eqassoc = StructuralTransformations.pantelides!(pendulum)
sss = structure(sys)
@unpack graph, fullvars, varassoc = sss
scc = StructuralTransformations.find_scc(graph, assign)
@test sort(sort.(scc)) == [
                           [1],
                           [2],
                           [3, 4, 7, 8, 9],
                           [5],
                           [6],
                          ]

@test graph.fadjlist == sort.([
 [5, 3],               # 1
 [6, 4],               # 2
 [7, 9, 1],            # 3
 [8, 9, 2],            # 4
 [2, 1],               # 5
 [2, 1, 6, 5],         # 6
 [5, 3, 10, 7],        # 7
 [6, 4, 11, 8],        # 8
 [2, 1, 6, 5, 11, 10], # 9
])
#                  [1, 2, 3, 4,   5,   6,  7,  8, 9,   10,   11]
#                  [x, y, w, z, xˍt, yˍt, w', z', T, xˍt', yˍt']
@test varassoc  == [5, 6, 7, 8,  10,  11,    0,    0, 0,      0,      0]
#1: D(x) ~ w
#2: D(y) ~ z
#3: D(w) ~ T*x
#4: D(z) ~ T*y - g
#5: 0 ~ x^2 + y^2 - L^2
# ----
#6: D(eq:5) -> 0 ~ 2xx'+ 2yy'
#7: D(eq:1) -> D(D(x)) ~ D(w) -> D(xˍt) ~ D(w) -> D(xˍt) ~ T*x
#8: D(eq:2) -> D(D(y)) ~ D(z) -> D(y_t) ~ T*y - g
#9: D(eq:6) -> 0 ~ 2xx'' + 2x'x' + 2yy'' + 2y'y'
#                 [1, 2, 3, 4, 5, 6, 7, 8, 9]
@test eqassoc == [7, 8, 0, 0, 6, 9, 0, 0, 0]

using ModelingToolkit
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
D = Differential(t)
idx1_pendulum = [D(x) ~ w,
                 D(y) ~ z,
                 #0 ~ x^2 + y^2 - L^2,
                 D(w) ~ T*x,
                 D(z) ~ T*y - g,
                 # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0,
                 # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y))
                 #0 ~ 2x*w + 2y*z,
                 # D(D(x)) ~ D(w) and substitute the rhs
                 D(xˍt) ~ T*x,
                 # D(D(y)) ~ D(z) and substitute the rhs
                 D(yˍt) ~ T*y - g,
                 # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and
                 # substitute the rhs
                 0 ~ 2x*(T*x) + 2*xˍt*xˍt + 2y*(T*y - g) + 2*yˍt*yˍt]
idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum)

using OrdinaryDiffEq
using LinearAlgebra
prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
        #  [x, y, w, z, xˍt, yˍt, T]
           [1, 0, 0, 0,   0,   0, 0.0],# 0, 0, 0, 0],
           (0, 100.0),
           [1, 9.8],
           mass_matrix=calculate_massmatrix(first_order_idx1_pendulum))
sol = solve(prob, Rodas5());
#plot(sol, vars=(1, 2))

new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))

prob_auto = ODEProblem(ODEFunction(new_sys),
        #  [xˍt, yˍt, x, y,   T]
           [0,     0, 1, 0, 0.0],# 0, 0, 0, 0],
           (0, 100.0),
           [1, 9.8],
           mass_matrix=calculate_massmatrix(new_sys))
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(3, 4))

# Define some variables
@parameters t L g
@variables x(t) y(t) T(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)

# Turn into a first order differential equation system
first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)

# Perform index reduction to get an Index 1 DAE
new_sys = dae_index_lowering(first_order_sys)

u0 = [
  D(x)    => 0.0,
  D(y)    => 0.0,
  x       => 1.0,
  y       => 0.0,
  T       => 0.0
]

p = [
    L => 1.0,
    g => 9.8
]

prob_auto = ODEProblem(new_sys,u0,(0.0,100.0),p)
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(D(x), y))

###
### More BLT/SCC tests
###

# Test Tarjan (1972) Fig. 3
g = [
     [2],
     [3,8],
     [4,7],
     [5],
     [3,6],
     Int[],
     [4,6],
     [1,7],
    ]
graph = StructuralTransformations.BipartiteGraph(8, 8)
for (eq, vars) in enumerate(g), var in vars
    add_edge!(graph, eq, var)
end
scc = StructuralTransformations.find_scc(graph)
@test scc == [
              [6],
              [3, 4, 5, 7],
              [1, 2, 8],
             ]


text:
#Notice that before run this script you should firstly use
#Pkg.add(Pkg.PackageSpec(;name="AbstractMCMC", version="0.5.3"))
#to downgrade your AbstractMCMC package so it'll be compatible
#with current tutuorial from Turing.ml, but noticed this can be
#a dangerous action as julia will automatically downgrade all
#dependents of AbstractMCMC

import AbstractMCMC
using Distributions
using Random
using MCMCChains
using StatsPlots

# Define a sampler type.
struct MetropolisHastings{T, D} <: AbstractMCMC.AbstractSampler
    init_θ::T
    proposal::D
end

# Default constructors.
MetropolisHastings(init_θ::Real) = MetropolisHastings(init_θ, Normal(0,1))
MetropolisHastings(init_θ::Vector{<:Real}) = MetropolisHastings(init_θ, MvNormal(length(init_θ),1))

# Define a model type. Stores the log density function.
struct DensityModel{F<:Function} <: AbstractMCMC.AbstractModel
    ℓπ::F
end

# Create a very basic Transition type, only stores the
# parameter draws and the log probability of the draw.
struct Transition{T, L}
    θ::T
    lp::L
end

# Store the new draw and its log density.
Transition(model::DensityModel, θ) = Transition(θ, ℓπ(model, θ))

# Define the first step! function, which is called at the
# beginning of sampling. Return the initial parameter used
# to define the sampler.
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    N::Integer,
    ::Nothing;
    kwargs...
)
    return Transition(model, spl.init_θ)
end

# Define a function that makes a basic proposal depending on a univariate
# parameterization or a multivariate parameterization.
propose(spl::MetropolisHastings, model::DensityModel, θ::Real) =
    Transition(model, θ + rand(spl.proposal))
propose(spl::MetropolisHastings, model::DensityModel, θ::Vector{<:Real}) =
    Transition(model, θ + rand(spl.proposal))
propose(spl::MetropolisHastings, model::DensityModel, t::Transition) =
    propose(spl, model, t.θ)

# Calculates the probability `q(θ|θcond)`, using the proposal distribution `spl.proposal`.
q(spl::MetropolisHastings, θ::Real, θcond::Real) = logpdf(spl.proposal, θ - θcond)
q(spl::MetropolisHastings, θ::Vector{<:Real}, θcond::Vector{<:Real}) =
    logpdf(spl.proposal, θ - θcond)
q(spl::MetropolisHastings, t1::Transition, t2::Transition) = q(spl, t1.θ, t2.θ)

# Calculate the density of the model given some parameterization.
ℓπ(model::DensityModel, θ) = model.ℓπ(θ)
ℓπ(model::DensityModel, t::Transition) = t.lp

# Define the other step function. Returns a Transition containing
# either a new proposal (if accepted) or the previous proposal
# (if not accepted).
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    ::Integer,
    θ_prev::Transition;
    kwargs...
)
    # Generate a new proposal.
    θ = propose(spl, model, θ_prev)

    # Calculate the log acceptance probability.
    α = ℓπ(model, θ) - ℓπ(model, θ_prev) + q(spl, θ_prev, θ) - q(spl, θ, θ_prev)

    # Decide whether to return the previous θ or the new one.
    if log(rand(rng)) < min(α, 0.0)
        return θ
    else
        return θ_prev
    end
end

# A basic chains constructor that works with the Transition struct we defined.
function AbstractMCMC.bundle_samples(
    rng::AbstractRNG,
    ℓ::DensityModel,
    s::MetropolisHastings,
    N::Integer,
    ts::Vector{<:Transition},
    chain_type::Type{Any};
    param_names=missing,
    kwargs...
)
    # Turn all the transitions into a vector-of-vectors.
    vals = copy(reduce(hcat,[vcat(t.θ, t.lp) for t in ts])')

    # Check if we received any parameter names.
    if ismissing(param_names)
        param_names = ["Parameter $i" for i in 1:(length(first(vals))-1)]
    end

    # Add the log density field to the parameter names.
    push!(param_names, "lp")

    # Bundle everything up and return a Chains struct.
    return Chains(vals, param_names, (internals=["lp"],))
end

# Generate a set of data from the posterior we want to estimate.
data = rand(Normal(5, 3), 30)

# Define the components of a basic model.
insupport(θ) = θ[2] >= 0
dist(θ) = Normal(θ[1], θ[2])
density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf

# Construct a DensityModel.
model = DensityModel(density)

# Set up our sampler with initial parameters.
spl = MetropolisHastings([0.0, 0.0])

# Sample from the posterior.
chain = sample(model, spl, 100000; param_names=["μ", "σ"])

summarystats(chain)
plot(chain; size=(840, 600))
savefig("Tutorial-AbstarctMCMC.png")


text:
#fi = "test_22.jl"
read_all_lines(fi) = readlines(fi)
function read_all_lines(fis::Vector)
    lins = String[]
    for fi in fis
        push!(lins, readlines(fi)...)
    end
    lins
end

function line_no(fi, no)
    read_all_lines(fi)[no]
end
function words(fi)
    wd = String[]
    for line in read_all_lines(fi)
        for w in split(line, r"\W")
            push!(wd, w)
        end
    end
    wd
end
#words(fi)
function word_lastline(fi)
    wdi = Dict{String, Int}()
    for (n, line) in enumerate(read_all_lines(fi))
        words = split(line, r"\W")
        for w in words
            push!(wdi, w => n)
        end
    end
    wdi
end
#word_lastline(fi)

function word_count(fi)
    wdi = Dict{String, Int}()
    for word in keys(word_lastline(fi))
        push!(wdi, word => get(wdi, word, 0) + 1)
    end
    wdi
end
#word_count(fi)

function singlewords_whichline(fi)
    wdi = Dict{String, Int}()
    w_ll = word_lastline(fi)
    for (wd, co) in word_count(fi)
        if co == 1
            push!(wdi, wd => w_ll[wd])
        end
    end
    wdi
end

#singlewords_whichline(fi)
function first_let_line(fi)
    lins = readlines(fi)
    n = length(lins)
    crite(l) = startswith(strip(l), "let") || startswith(strip(l), "#let")
    findfirst(crite, lins)
end
function referred_files(fi)
    reffis = String[]
    for (no, line) in enumerate(readlines(fi))
        if startswith(strip(line), "include(")
            sp, lp = findall("\"", line)
            push!(reffis, line[sp[1] + 1:lp[1] - 1])
        end
    end
    reffis
end

"""
    loopfiles(;mustinclude = "")

For each file, print a list of candidates for unecessary import, ie occuring in the first ten lines of each file and only once"
"""
function loopfiles(;mustinclude = "")
    excludelist = ["MechanicalSketch", "__DIR__", "using", "import", "Test"]
    for fi in filter(s->endswith(s, ".jl"), readdir())
        if contains(fi, mustinclude)
            referredfiles = referred_files(fi)
            fis = [fi; referredfiles...]
            letatline = first_let_line(fi)
            printstyled(stdout, "\n", fi, ":\n"; color=:yellow)
            wd_li = filter(singlewords_whichline(fis)) do (wd, lin)
                lin < letatline &&
                    length(wd) > 1 &&
                        !isnumeric(first(wd)) &&
                            !isdefined(Base, Symbol(wd)) &&
                                wd ∉ excludelist
            end
            byline = collect(keys(sort(wd_li, byvalue= true)))
            w_ll = word_lastline(fis)
            println("Single occurence words:")
            plno = 0
            lstr = ""
            for wd in byline
                lno = wd_li[wd]
                if lno > plno || wd == last(byline)
                    rstr = " at line $plno"
                    plno > 0 && println(lpad(lstr, 70), rstr)
                    lstr = wd
                    plno = lno
                else
                    lstr *= ", " * wd
                end
            end
        end
    end
end
println("""Example: loopfiles(mustinclude = "22")""")


text:
# Preprocessing for the big text files

const SOS = "<s>"
const EOS = "</s>"
const UNK = "<unk>"
const LLIMIT = 3
const ULIMIT = 300

type Data
    word_to_index::Dict{AbstractString, Int}
    index_to_word::Vector{AbstractString}
    sequences::Array
end

function Data(datafile; word_to_index=nothing, vocabfile=nothing)

    vocab_exists = (word_to_index != nothing)
    if !vocab_exists
        word_to_index = Dict{AbstractString, Int}(SOS=>1, EOS=>2, UNK=>3)
        if vocabfile != nothing
            info("Working with provided vocabfile: $vocabfile")
            V = vocab_from_file(vocabfile)
        end
    end

    stream = open(datafile)
    sequences = Array[]
    for line in eachline(stream)
        words = Int32[]
        push!(words, word_to_index[SOS])
        for word in split(line)
            if !vocab_exists && vocabfile != nothing && !(word in V)
                word = UNK
            end
            if vocab_exists
                index = get(word_to_index, word, word_to_index[UNK])
            else
                index = get!(word_to_index, word, 1+length(word_to_index))
            end
            push!(words, index)
        end
        push!(words, word_to_index[EOS])
        (length(words) < ULIMIT) && (length(words)>LLIMIT) && push!(sequences, words)
    end
    close(stream)
    vocabsize = length(word_to_index)
    index_to_word = Array(AbstractString, vocabsize)
    for (word, index) in word_to_index
        index_to_word[index] = word
    end
    sort!(sequences, by=x->length(x), rev=true)
    Data(word_to_index, index_to_word, sequences)
end

""" Creates a set that contains all the words in that file, vocab file given as each vocab in a single line
    sorted_counted represents pure create_vocab.sh output
"""
function vocab_from_file(vocabfile::AbstractString)
    V = Set{AbstractString}()
    open(vocabfile) do file
        for line in eachline(file)
            words= split(line)
            @assert(length(words) == 1, "The vocabulary file seems broken")
            push!(V, words[1])
        end
    end
    return V
end

""" To be able to understand the sequence structure of the words """
function create_metadata(d::Data)
    sequences = d.sequences
    metadata = Dict{Int32, Int32}()
    for item in sequences
        if length(item) in keys(metadata)
            metadata[length(item)] += 1
        else
            metadata[length(item)] = 1
        end
    end
    return metadata
end

function batch_info(metadata::Dict{Int32, Int32}, batchsize::Int)
    c1 = 0
    c2 = 0
    info("There are $(length(metadata)) different sequence lengths")
    for item in keys(metadata)
        if metadata[item]< batchsize
            c1 += 1
        end
        if div(batchsize,2) < metadata[item] < batchsize
            c2 += 1
        end
    end
    info("There are $c1 sequences less than batchsize")
    info("There are $c2 sequences bigger than half batchsize and less than batchsize")
end

function create_offset(d::Data)
    offsets = Int32[]
    s = d.sequences
    curr = prev = length(s[1])
    push!(offsets, 1)
    for i=2:length(s)
        curr = length(s[i])
        if curr != prev
            push!(offsets, i)
            prev = curr
        end
    end
    (offsets[end] != Int32(length(s))) && push!(offsets, length(s))
    return offsets
end

function modify_offset(of::Array{Int32, 1}, batchsize::Int)
    result = Int32[]
    for i = 1:length(of) - 1 # TODO take care of the last entry
        nit = of[i]
        push!(result, nit)
        nit += batchsize
        while (nit < of[i+1])
            push!(result, nit)
            nit += batchsize
        end
    end
    push!(result, of[end])
    return result
end

function create_minibatch(sequences::Array{Array,1}, batchsize::Int, st::Int, en::Int, vocabsize::Int)
    st = Int64(st)
    seqlen = length(sequences[st])
    data = [ falses(batchsize, vocabsize) for i=1:seqlen]
    sentences = sequences[st:en]
    noi = length(sentences)
    for cursor=1:seqlen
        for row=1:noi
            index = sentences[row][cursor]
            data[cursor][row, index] = 1
        end
    end
    return data
end

function ibuild_sentence(index_to_word::Array{AbstractString,1}, sequence::Array{BitArray{2}, 1}, kth::Int)
    sentence = Any[]
    for i=1:length(sequence)
        z = find(x->x==true, sequence[i][kth, :])
        append!(sentence, index_to_word[z])
    end
    return sentence
end

function next_seq(sequences::Array{Array, 1}, mof::Array{Int32, 1}, index::Int, vocabsize::Int, batchsize::Int)
    st = Int64(mof[index])
    en = mof[index+1]-1
    seq = create_minibatch(sequences, batchsize, st, en, vocabsize)
    return seq
end

# sample usage:
# @time d = Data("../bilstm-in-Knet8/data/nounCompound/google_sets/noun_google_data";vocabfile ="../bilstm-in-Knet8/data/nounCompound/google_sets/google_noun10k.vocab");

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
using Base: @kwdef

_checktype(type) = (type != :reg && type != :wire) &&
    error("Cannot create net with type $type (use :reg or :wire)")
_checkclass(class) = (class != :input &&
                      class != :output &&
                      class != :internal &&
                      class != :constant &&
                      class != :parameter) &&
    error("Cannot create net with class $class (use :input, :output, :internal, :constant, or :parameter)")
_checksize(size) = any(size .< 1) &&
    error("Cannot create net with size $size (size[i] must be >= 1)")

@kwdef struct Net
    name::String
    type::Symbol = :wire
    class::Symbol = :internal
    signed::Bool = false
    size::Tuple{Int, Int}

    function Net(name, type, class, signed, size)
        _checktype(type)
        _checkclass(class)
        _checksize(size)

        new(name, type, class, signed, size)
    end
end

isreg(x::Net) = (x.type == :reg)
iswire(x::Net) = !isreg(x)

isinput(x::Net) = (x.class == :input)
isoutput(x::Net) = (x.class == :output)
isinternal(x::Net) = (x.class == :internal)
isconstant(x::Net) = (x.class == :constant)
isparameter(x::Net) = (x.class == :parameter)

issigned(x::Net) = x.signed

const Netlist = Vector{Net}

inputs(n::Netlist) = filter(isinput, n)
outputs(n::Netlist) = filter(isoutput, n)

find(n::Netlist, x::String) = findfirst(λ -> λ.name == x, n)
contains(n::Netlist, x::String) = !isnothing(find(n, x))

function getsize(n::Netlist, x::String)
    i = find(n, x)
    isnothing(i) && error("Cannot get size for net $x since it does not exist in netlist.")

    return n[i].size
end
getsize(n::Netlist, x::Vector{String}) = map(λ -> getsize(n, λ), x)

delete!(n::Netlist, x::String) = deleteat!(n, find(n, x))

function update!(n::Netlist, x::Net)
    i = findfirst(λ -> λ.name == x.name, n)
    if isnothing(i)
        push!(n, x)
    else
        n[i] = x
    end

    return n
end

function setreg!(n::Netlist, x::String)
    i = find(n, x)
    isnothing(i) && error("Cannot set net $x as register because it does not exist in netlist.")
    update!(n, Net(name = x, type = :reg, class = n[i].class, signed = n[i].signed, size = n[i].size))

    return n
end

function setwire!(n::Netlist, x::String)
    i = find(n, x)
    isnothing(i) && error("Cannot set net $x as wire because it does not exist in netlist.")
    update!(n, Net(name = x, type = :wire, class = n[i].class, signed = n[i].signed, size = n[i].size))

    return n
end

function setsigned!(n::Netlist, x::String, signed::Bool)
    i = find(n, x)
    isnothing(i) && error("Cannot set net $x as signed = $signed because it does not exist in netlist.")
    update!(n, Net(name = x, type = n[i].type, class = n[i].class, signed = signed, size = n[i].size))

    return n
end


text:
export Batched

"""
Batched operation in Julia.

This module wraps some useful batched operation with a plain
for-loop on CPU. All the functions in this module are defined
with gradients in YAAD.
"""
module Batched
import LinearAlgebra
import LinearAlgebra: BLAS

# last dim is batch dim
for (gemm, elty) in
        ((:dgemm_,:Float64),
         (:sgemm_,:Float32),
         (:zgemm_,:ComplexF64),
         (:cgemm_,:ComplexF32))
    @eval begin
        function gemm!(transA::AbstractChar, transB::AbstractChar, alpha::($elty), A::AbstractArray{$elty, 3}, B::AbstractArray{$elty, 3}, beta::($elty), C::AbstractArray{$elty, 3})
            @assert !BLAS.has_offset_axes(A, B, C)
            @assert size(A, 3) == size(B, 3) == size(C, 3) "batch size mismatch"
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            BLAS.chkstride1(A)
            BLAS.chkstride1(B)
            BLAS.chkstride1(C)

            ptrA = Base.unsafe_convert(Ptr{$elty}, A)
            ptrB = Base.unsafe_convert(Ptr{$elty}, B)
            ptrC = Base.unsafe_convert(Ptr{$elty}, C)

            for k in 1:size(A, 3)
                ccall((LinearAlgebra.BLAS.@blasfunc($gemm), BLAS.libblas), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt},
                     Ref{BLAS.BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BLAS.BlasInt},
                     Ptr{$elty}, Ref{BLAS.BlasInt}, Ref{$elty}, Ptr{$elty},
                     Ref{BLAS.BlasInt}),
                     transA, transB, m, n,
                     ka, alpha, ptrA, max(1,stride(A,2)),
                     ptrB, max(1,stride(B,2)), beta, ptrC,
                     max(1,stride(C,2)))

                ptrA += size(A, 1) * size(A, 2) * sizeof($elty)
                ptrB += size(B, 1) * size(B, 2) * sizeof($elty)
                ptrC += size(C, 1) * size(C, 2) * sizeof($elty)
            end
            C
        end
        function gemm(transA::AbstractChar, transB::AbstractChar, alpha::($elty), A::AbstractArray{$elty, 3}, B::AbstractArray{$elty, 3})
            gemm!(transA, transB, alpha, A, B, zero($elty), similar(B, $elty, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1), size(B, 3))))
        end
        function gemm(transA::AbstractChar, transB::AbstractChar, A::AbstractArray{$elty, 3}, B::AbstractArray{$elty, 3})
            gemm(transA, transB, one($elty), A, B)
        end
    end
end

gemm(A::AbstractArray{T, 3}, B::AbstractArray{T, 3}) where T = gemm('N', 'N', A, B)

function tr!(A::AbstractArray{T, 3}, out::AbstractArray{T, 1}) where T
    @inbounds for i in eachindex(out)
        out[i] = LinearAlgebra.tr(view(A, :, :, i))
    end
    out
end

tr(A::AbstractArray{T, 3}) where T = tr!(A, Vector{T}(undef, size(A, 3)))

################################################################################
##########################    Transpose Operations   ###########################

export Transpose

"""
    Transpose{B, T, AT <: AbstractArray{T, 3}} <: AbstractArray{T, 3}

Batched transpose. Transpose a batch of matrix.
"""
struct Transpose{B, T, AT <: AbstractArray{T, 3}} <: AbstractArray{T, 3}
    parent::AT

    Transpose(A::AT) where {T, AT <: AbstractArray{T, 3}} = new{size(A, 3), T, AT}(A)
end

Base.size(m::Transpose) = (size(m.parent, 2), size(m.parent, 1), size(m.parent, 3))
Base.axes(m::Transpose) = (axes(m.parent, 2), axes(m.parent, 1), axes(m.parent, 3))
Base.@propagate_inbounds Base.getindex(m::Transpose, i::Int, j::Int, k::Int) = getindex(m.parent, j, i, k)
Base.@propagate_inbounds Base.setindex!(m::Transpose, v, i::Int, j::Int, k::Int) = setindex!(m.parent, v, j, i, k)
Base.IndexStyle(::Type{<:Transpose}) = IndexCartesian()

# FIXME: implement a different generic function in case this becomes a pirate?
Base.transpose(A::AbstractArray{T, 3}) where T = Transpose(A)

gemm(A::AbstractArray{T, 3}, B::Batched.Transpose) where T = gemm('N', 'T', A, B.parent)
gemm(A::Batched.Transpose, B::AbstractArray{T, 3}) where T = gemm('T', 'N', A.parent, B)
gemm(A::Batched.Transpose, B::Batched.Transpose) = gemm('T', 'T', A.parent, B.parent)

export ScaleMatrix

"""
    ScaleMatrix{B, K, T} <: AbstractArray{T, 3}

A batch of scalar multiplies a batch of identities, where batch size is
`B`, each identity's size is `K`.
"""
struct ScaleMatrix{B, K, T, VT <: AbstractVector{T}} <: AbstractArray{T, 3}
    scalars::VT
    ScaleMatrix{B, K}(scalars::VT) where {B, K, T, VT <: AbstractVector{T}} = new{B, K, T, VT}(scalars)
end

Base.size(x::ScaleMatrix{B, K, T}) where {B, K, T} = (K, K, B)
Base.@propagate_inbounds Base.getindex(m::ScaleMatrix{B, K, T}, i::Int, j::Int, k::Int) where {B, K, T} =
    i == j ? getindex(m.scalars, k) : zero(T)
Base.IndexStyle(::Type{<:ScaleMatrix}) = IndexCartesian()
Base.transpose(A::ScaleMatrix) = A

function gemm!(A::ScaleMatrix{NBatch, K, T}, B::AbstractArray{T, 3}, C::AbstractArray{T, 3}) where {NBatch, K, T}
    @inbounds for k in 1:NBatch, j in 1:size(B, 2), i in 1:size(B, 1)
        C[i, j, k] += A.scalars[k] * B[i, j, k]
    end
    C
end

gemm!(A::AbstractArray{T, 3}, B::ScaleMatrix{NBatch, K, T}, C::AbstractArray{T, 3}) where {NBatch, K, T} = gemm!(B, A, C)
gemm!(A::Transpose{NBatch, T}, B::ScaleMatrix{NBatch, K, T}, C::AbstractArray{T, 3}) where {NBatch, K, T} = gemm!(B, A, C)

end # Batched

Batched.gemm(A::Value, B::Value) = register(Batched.gemm, A, B)

function gradient(::typeof(Batched.gemm), grad, output, A::AbstractArray{T, 3}, B::AbstractArray{T, 3}) where T
    Batched.gemm(grad, transpose(B)), Batched.gemm(transpose(A), grad)
end

Batched.tr(A::Value) = register(Batched.tr, A)

function gradient(::typeof(Batched.tr), grad, output, A::AbstractArray{T, 3}) where T
    (Batched.ScaleMatrix{size(A, 3), size(A, 1)}(grad), )
end


text:
module FEMDriver

using Test
using Gridap
using GridapPardiso

model = CartesianDiscreteModel(
  domain=(0,1,0,1,0,1), partition=(10,10,10))

fespace = FESpace(
  reffe=:Lagrangian, order=1, valuetype=Float64,
  conformity=:H1, model=model, diritags="boundary")

V = TestFESpace(fespace)
U = TrialFESpace(fespace,0.0)

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

t_Ω = AffineFETerm(
  (v,u) -> inner(∇(v),∇(u)),
  (v) -> inner(v, (x) -> x[1]*x[2] ),
  trian, quad)

op = LinearFEOperator(V,U,t_Ω)

ls = PardisoSolver()
solver = LinearFESolver(ls)

uh = solve(solver,op)

x = free_dofs(uh)
A = op.mat
b = op.vec

r = A*x - b

tol = 1e-10
@test maximum(abs.(r)) < tol

end #module