text:
abstract type Constraint end
# expected fields (mutable)
# - score

# reduce Base.show() to Base.string()
Base.show(io::IO, cnstr::Constraint) = print(io, string(cnstr))
Base.show(io::IO, ::MIME"application/prs.juno.inline", cnstr::Constraint) = print(io, string(cnstr))

Base.string(cnstr::Constraint) = @sprintf "ToDo: implement Base.string for %s" typeof(cnstr)

integrand_x1(cnstr::Constraint, args...; kwargs...) = integrand_x1(cnstr.score, args...; kwargs...)

update!(cnstr::Constraint, x::TI, n::TI) where {TS<:Score,TI<:Integer} = update!(cnstr.score, x, n)

function update(cnstr::Constraint, prior::Prior)
    cnstr       = deepcopy(cnstr)
    cnstr.score = update(cnstr.score, prior)
    return cnstr
end

function conditional(cnstr::Constraint, bounds::Tuple{Real,Real})
    cnstr = deepcopy(cnstr)
    cnstr.conditional = bounds
    return cnstr
end

mutable struct PowerConstraint <: Constraint
    score::Power
    β::Real
    conditional::Tuple{Real,Real}
end

>=(power::Power, threshold::Real) = PowerConstraint(power, 1 - threshold, (.5, .99))

function Base.string(cnstr::PowerConstraint)

    return @sprintf "%s >= %5.1f%% (given x1: %5.1f%% - %5.1f%%)" string(cnstr.score) 100*(1 - cnstr.β) 100*cnstr.conditional[1] 100*cnstr.conditional[2]
end

mutable struct TypeOneErrorRateConstraint <: Constraint
    score::Power
    α::Real
    conditional::Tuple{Real,Real}
end

<=(power::Power, threshold::Real) = TypeOneErrorRateConstraint(power, threshold, (.001, .99))

function Base.string(cnstr::TypeOneErrorRateConstraint)

    return @sprintf "%s <= %5.1f%% (given x1: %5.1f%% - %5.1f%%)" string(cnstr.score) 100*cnstr.α 100*cnstr.conditional[1] 100*cnstr.conditional[2]
end


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

"""
cli.jl defines the `jltest.cli` module containing functions for the `jltest` CLI.

Notes
-----
* CLI functions are defined in a .jl file so that testing and code quality tools can by
  applied to the CLI source code.
"""
module cli

# --- Exports

export parse_args, run

# --- Imports

# Standard library
using Test

# External packages
using ArgParse: ArgParse

# Local modules
using ..jltest

# --- Functions/Methods

"""
    parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

Parse and return CLI arguments contained in `raw_args`. By default, `raw_args` is set to
`ARGS`, the command-line arguments provided to the executable that called `parse_args()`.
"""
function parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

    # Define the command-line interface
    description = "Run unit tests."
    arg_table = ArgParse.ArgParseSettings(; prog="jltest", description=description)
    ArgParse.@add_arg_table! arg_table begin
        "--fail-fast", "-x"
        help = "stop testing at first failure"
        action = :store_true

        "--no-wrapper", "-W"
        help = """
               run tests without wrapping them in an EnhancedTestSet. Ignored if fail-fast
               is enabled. Note: running tests in this mode manner will disable
               EnhancedTestSet functionality within each test files unless the test file
               (1) explicitly sets the test set when using the `@testset` macro or
               (2) calls `jltest.run_tests()` (without setting the `test_set_type` keyword
               argument to `nothing`)
               """
        action = :store_true

        "--no-recursion", "-R"
        help = "do not recursively search directories for tests"
        action = :store_true

        "--verbose", "-v"
        help = "display more output to the console"
        action = :store_true

        "--version", "-V"
        help = "show version and exit"
        action = :store_true

        "tests"
        help =
            "tests to run. Any directories present in `tests` are searched " *
            "(recursively) for tests to run. If `tests` is omitted, the current " *
            "directory is searched (recursively) for tests to run."
        nargs = '*'
    end

    # Parse command-line arguments
    args = ArgParse.parse_args(raw_args, arg_table)
    args["tests"] = convert(Vector{String}, args["tests"])

    return args
end

"""
    run(tests::Vector; <keyword arguments>)

Run unit tests defined in the list of files or modules provided in `tests`.

# Keyword Arguments

* `fail_fast::Bool`: flag indicating whether or not to stop testing at first failure.
  Default: `false`

* `use_wrapper::Bool`: flag indicating whether or not to run tests without first wrapping
  them in an EnhancedTestSet. Default: `true`

  !!! note
      Ignored if `fail_fast` is set to `true`.

  !!! note
      Setting `use_wrapper` to `false` will disable EnhancedTestSet functionality within
      each test files unless the test file

      * explicitly sets the test set when using the `@testset` macro or

      * calls `jltest.run_tests()` (without setting the `test_set_type` keyword argument
        to `nothing`).

* `recursive::Bool`: flag indicating whether or not to run tests found in subdirectories
  of directories in `tests`. Default: `true`

* `verbose::Bool`: print more output to the console. Default: `false`
"""
function run(
    tests::Vector;
    fail_fast::Bool=false,
    use_wrapper::Bool=true,
    recursive::Bool=true,
    verbose::Bool=false,
)
    # --- Check arguments

    # Ensure that `tests` contains strings
    tests = convert(Vector{String}, tests)

    # --- Preparations

    # Set test set type
    if fail_fast
        test_set_type = EnhancedTestSet{Test.FallbackTestSet}
    else
        if use_wrapper
            test_set_type = EnhancedTestSet{Test.DefaultTestSet}
        else
            test_set_type = nothing
        end
    end

    # Set test options
    # TODO: figure out how to pass verbose option to @testset when the test set type
    #       is explicitly specified
    #test_set_options = ""
    #if verbose
    #    test_set_options *= "verbose=true"
    #end

    # --- Run tests

    # Unit tests
    if length(tests) == 0
        tests = find_tests(pwd())
    end

    run_tests(tests; desc="All tests", test_set_type=test_set_type, recursive=recursive)

    return nothing
end

end  # End of jltest.cli module

ext: jl


text:
module Custom

using MLJ, MLJScikitLearnInterface, DataFrames

export MLJ, MLJScikitLearnInterface, transform

function transform(data, model)
    if "Species" in names(data)
        data = select!(data, Not(:Species))
    elseif "MEDV" in names(data)
        data = select!(data, Not(:MEDV))
    elseif "class" in names(data)
        data = select!(data, Not(:class))
    end
    return data
end

end

max_stars_repo_name: amperie/user-models

max_issues_repo_name: amperie/user-models

max_forks_repo_name: amperie/user-models


text:
@doc raw"""
    moffat([T=Float64], point; x, y, fwhm, alpha=1, amp=1, theta=0, bkg=0)
    moffat([T=Float64], px, py; x, y, fwhm, alpha=1, amp=1, theta=0, bkg=0)

Two dimensional Moffat model. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. If `theta` is given, the PSF will be rotated by `theta` degrees counter-clockwise from the x-axis. If `bkg` is given it will be added as a scalar to the PSF.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that `theta` has no effect for isotropic distributions and is degenerate with the `fwhm` parameters (i.e., theta=90 is the same as reversing the `fwhm` tuple)

# Functional form
```
f(x | x̂, FWHM, α) = A / (1 + ||x - x̂|| / (FWHM / 2)^2)^α
```
where `x̂` and `x` are position vectors (indices) `||⋅||` represents the square-distance, and `FWHM` is the full width at half-maximum. If `FWHM` is a vector or tuple, the weighting is applied along each axis.
"""
moffat(T, px, py; x, y, fwhm, alpha=1, amp=one(T), theta=0, bkg=0) = convert(T, _moffat(px, py, x, y, fwhm, alpha, amp, theta, bkg))

# isotropic
function _moffat(px, py, x, y, fwhm, alpha, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    !iszero(theta) && @warn "isotropic moffat is not affected by non-zero rotation angle $theta"
    # unnormalized moffat
    dist = dx^2 + dy^2
    dist /= (fwhm / 2)^2
    return amp / (1 + dist)^alpha + background
end

# bivariate
function _moffat(px, py, x, y, fwhm::BivariateLike, alpha, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    if !iszero(theta)
        dx, dy = rotate_point(dx, dy, theta)
    end
    # unnormalized moffat
    fwhmx, fwhmy = fwhm
    # unnormalized moffat
    dist = (dx / (fwhmx/2))^2 + (dy / (fwhmy/2))^2
    return amp / (1 + dist)^alpha + background
end


text:
function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

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

function get_converter(::Type{ITRFVersion}, arg0::ITRFVersion, arg1::ITRFVersion)
    return jcall(ITRFVersion, "getConverter", ITRFVersion_Converter, (ITRFVersion, ITRFVersion), arg0, arg1)
end

function get_converter(::Type{ITRFVersion}, arg0::ITRFVersion, arg1::ITRFVersion, arg2::TimeScale)
    return jcall(ITRFVersion, "getConverter", ITRFVersion_Converter, (ITRFVersion, ITRFVersion, TimeScale), arg0, arg1, arg2)
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

function get_itrf_version(::Type{ITRFVersion}, arg0::JString)
    return jcall(ITRFVersion, "getITRFVersion", ITRFVersion, (JString,), arg0)
end

function get_itrf_version(::Type{ITRFVersion}, arg0::jint)
    return jcall(ITRFVersion, "getITRFVersion", ITRFVersion, (jint,), arg0)
end

function get_name(obj::ITRFVersion)
    return jcall(obj, "getName", JString, ())
end

function get_year(obj::ITRFVersion)
    return jcall(obj, "getYear", jint, ())
end

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

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

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

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function value_of(::Type{ITRFVersion}, arg0::JString)
    return jcall(ITRFVersion, "valueOf", ITRFVersion, (JString,), arg0)
end

function values(::Type{ITRFVersion})
    return jcall(ITRFVersion, "values", Vector{ITRFVersion}, ())
end


text:
using NNlib: pad_constant, pad_repeat, pad_zeros, pad_reflect

@testset "padding constant" begin
  x = rand(2, 2, 2)

  p = NNlib.gen_pad((1,2,3,4,5,6), (1,2,3), 4)
  @test p == ((1, 2), (3, 4), (5, 6), (0, 0))

  @test_throws ArgumentError NNlib.gen_pad((1,2,3,4,5,), (1,2,3), 4)

  p = NNlib.gen_pad((1,3), (1,3), 4)
  @test p == ((1, 1), (0, 0), (3, 3), (0, 0))

  p = NNlib.gen_pad(1, (1,2,3), 4)
  @test p == ((1, 1), (1, 1), (1, 1), (0, 0))

  p = NNlib.gen_pad(3, :, 2)
  @test p == ((3, 3), (3, 3))

  y = pad_constant(x, (3, 2, 4))
  @test size(y) == (8, 6, 10)
  @test y[4:5, 3:4, 5:6] ≈ x
  y[4:5, 3:4, 5:6] .= 0
  @test all(y .== 0)

  @test pad_constant(x, (3, 2, 4)) ≈ pad_zeros(x, (3, 2, 4))
  @test pad_zeros(x, 2) ≈ pad_zeros(x, (2,2,2))

  y = pad_constant(x, (3, 2, 4, 5), 1.2, dims = (1,3))
  @test size(y) == (7, 2, 11)
  @test y[4:5, 1:2, 5:6] ≈ x
  y[4:5, 1:2, 5:6] .= 1.2
  @test all(y .== 1.2)

  @test pad_constant(x, (2,2,2,2), 1.2, dims = (1,3)) ≈
          pad_constant(x, 2, 1.2, dims = (1,3))

  @test pad_constant(x, 1, dims = 1:2) ==
          pad_constant(x, 1, dims = (1,2))

  @test size(pad_constant(x, 1, dims = 1)) == (4,2,2)

  @test all(pad_zeros(randn(2), (1, 2))[[1, 4, 5]] .== 0)

  gradtest(x -> pad_constant(x, 2), rand(2,2,2))
  gradtest(x -> pad_constant(x, (2, 1, 1, 2)), rand(2,2))
  gradtest(x -> pad_constant(x, (2, 1,)), rand(2))
end

@testset "padding repeat" begin
  x = rand(2, 2, 2)

  # y = @inferred pad_repeat(x, (3, 2, 4, 5))
  y = pad_repeat(x, (3, 2, 4, 5))
  @test size(y) == (7, 11, 2)
  @test y[4:5, 5:6, :] ≈ x

  # y = @inferred pad_repeat(x, (3, 2, 4, 5), dims=(1,3))
  y = pad_repeat(x, (3, 2, 4, 5), dims=(1,3))
  @test size(y) == (7, 2, 11)
  @test y[4:5, :, 5:6] ≈ x

  @test pad_repeat(reshape(1:9, 3, 3), (1,2)) ==
        [1  4  7
        1  4  7
        2  5  8
        3  6  9
        3  6  9
        3  6  9]

  @test pad_repeat(reshape(1:9, 3, 3), (2,2), dims=2) ==
       [1  1  1  4  7  7  7
        2  2  2  5  8  8  8
        3  3  3  6  9  9  9]

  @test pad_repeat(x, (2, 2, 2, 2), dims=(1,3)) ≈
          pad_repeat(x, 2, dims=(1,3))

  gradtest(x -> pad_repeat(x, (2,2,2,2)), rand(2,2,2))
end

@testset "padding reflect" begin
  y = pad_reflect(reshape(1:9, 3, 3), (2,2), dims=2)
  @test y ==  [ 7  4  1  4  7  4  1
                8  5  2  5  8  5  2
                9  6  3  6  9  6  3]

  y = pad_reflect(reshape(1:9, 3, 3), (2,2,2,2))
  @test y ==   [9  6  3  6  9  6  3
                8  5  2  5  8  5  2
                7  4  1  4  7  4  1
                8  5  2  5  8  5  2
                9  6  3  6  9  6  3
                8  5  2  5  8  5  2
                7  4  1  4  7  4  1]

  x = rand(4, 4, 4)

  @test pad_reflect(x, (2, 2, 2, 2), dims=(1,3)) ≈
          pad_reflect(x, 2, dims=(1,3))

  gradtest(x -> pad_repeat(x, (2,2,2,2)), rand(2,2,2))
end


text:
"""
    rand_graph(n, m; bidirected=true, seed=-1, kws...)

Generate a random (Erdós-Renyi) `GNNGraph` with `n` nodes
and `m` edges.

If `bidirected=true` the reverse edge of each edge will be present.
If `bidirected=false` instead, `m` unrelated edges are generated.
In any case, the output graph will contain no self-loops or multi-edges.

Use a `seed > 0` for reproducibility.

Additional keyword arguments will be passed to the [`GNNGraph`](@ref) constructor.

# Examples

```juliarepl
julia> g = rand_graph(5, 4, bidirected=false)
GNNGraph:
    num_nodes = 5
    num_edges = 4

julia> edge_index(g)
([1, 3, 3, 4], [5, 4, 5, 2])

# In the bidirected case, edge data will be duplicated on the reverse edges if needed.
julia> g = rand_graph(5, 4, edata=rand(16, 2))
GNNGraph:
    num_nodes = 5
    num_edges = 4
    edata:
        e => (16, 4)

# Each edge has a reverse
julia> edge_index(g)
([1, 3, 3, 4], [3, 4, 1, 3])

```
"""
function rand_graph(n::Integer, m::Integer; bidirected=true, seed=-1, kws...)
    if bidirected
        @assert iseven(m) "Need even number of edges for bidirected graphs, given m=$m."
    end
    m2 = bidirected ? m÷2 : m
    return GNNGraph(Graphs.erdos_renyi(n, m2; is_directed=!bidirected, seed); kws...)
end

"""
    knn_graph(points::AbstractMatrix,
              k::Int;
              graph_indicator = nothing,
              self_loops = false,
              dir = :in,
              kws...)

Create a `k`-nearest neighbor graph where each node is linked
to its `k` closest `points`.

# Arguments

- `points`: A num_features × num_nodes matrix storing the Euclidean positions of the nodes.
- `k`: The number of neighbors considered in the kNN algorithm.
- `graph_indicator`: Either nothing or a vector containing the graph assigment of each node,
                     in which case the returned graph will be a batch of graphs.
- `self_loops`: If `true`, consider the node itself among its `k` nearest neighbors, in which
                case the graph will contain self-loops.
- `dir`: The direction of the edges. If `dir=:in` edges go from the `k`
         neighbors to the central node. If `dir=:out` we have the opposite
         direction.
- `kws`: Further keyword arguments will be passed to the [`GNNGraph ](@ref) constructor.

# Examples

```juliarepl
julia> n, k = 10, 3;

julia> x = rand(3, n);

julia> g = knn_graph(x, k)
GNNGraph:
    num_nodes = 10
    num_edges = 30

julia> graph_indicator = [1,1,1,1,1,2,2,2,2,2];

julia> g = knn_graph(x, k; graph_indicator)
GNNGraph:
    num_nodes = 10
    num_edges = 30
    num_graphs = 2

```
"""
function knn_graph(points::AbstractMatrix, k::Int;
        graph_indicator = nothing,
        self_loops = false,
        dir = :in,
        kws...)

    if graph_indicator !== nothing
        d, n = size(points)
        @assert graph_indicator isa AbstractVector{<:Integer}
        @assert length(graph_indicator) == n
        # All graphs in the batch must have at least k nodes.
        cm = StatsBase.countmap(graph_indicator)
        @assert all(values(cm) .>= k)

        # Make sure that the distance between points in different graphs
        # is always larger than any distance within the same graph.
        points = points .- minimum(points)
        points = points ./ maximum(points)
        dummy_feature = 2d .* reshape(graph_indicator, 1, n)
        points = vcat(points, dummy_feature)
    end

    kdtree = NearestNeighbors.KDTree(points)
    if !self_loops
        k += 1
    end
    sortres = false
    idxs, dists = NearestNeighbors.knn(kdtree, points, k, sortres)

    g = GNNGraph(idxs; dir, graph_indicator, kws...)
    if !self_loops
        g = remove_self_loops(g)
    end
    return g
end


text:
using NFFT, DataFrames, LinearAlgebra, Statistics, LaTeXStrings, DelimitedFiles, CuNFFT
using Plots; pgfplotsx()

include("../../Wrappers/NFFT3.jl")
include("../../Wrappers/FINUFFT.jl")

const packagesCtor = [NFFTPlan, NFFTPlan, NFFTPlan, NFFTPlan,  NFFT3Plan, NFFT3Plan, FINUFFTPlan]
const packagesStr = ["NFFT.jl/FULL", "NFFT.jl/LINEAR", "NFFT.jl/TENSOR", "NFFT.jl/POLY", "NFFT3/LINEAR", "NFFT3/TENSOR", "FINUFFT"]
const precomp = [NFFT.FULL, NFFT.LINEAR, NFFT.TENSOR, NFFT.POLYNOMIAL, NFFT.LINEAR, NFFT.TENSOR, NFFT.LINEAR]
const blocking = [false, true, true, true, false, false, false]

#const packagesCtor = [NFFTPlan, CuNFFT.CuNFFTPlan, NFFT3Plan, FINUFFTPlan ]
#const packagesStr = ["NFFT.jl", "CuNFFT.jl", "NFFT3", "FINUFFT", ]
#const precomp = [NFFT.TENSOR,  NFFT.FULL, NFFT.TENSOR, NFFT.LINEAR, ]
#const blocking = [true, true, true, false]

const σs = range(1.25, 4, length=12)
const ms = 3:10
const NBase = [4096, 64, 16]
const Ds = 1:3

function nfft_accuracy_comparison(Ds, σs, ms)
  println("\n\n ##### nfft_accuracy_comparison ##### \n\n")

  df = DataFrame(Package=String[], D=Int[], M=Int[], N=Int[], m = Int[], σ=Float64[],
                   ErrorTrafo=Float64[], ErrorAdjoint=Float64[] )

  for D in Ds
      N = ntuple(d->NBase[D], D)
      M = prod(N)

      for σ in σs
        for m in ms
          @info "m=$m D=$D σ=$σ "
          x = rand(D,M) .- 0.5
          fHat = randn(ComplexF64, M)

          # ground truth (numerical)
          pNDFT = NDFTPlan(x, N)
          f = adjoint(pNDFT) * fHat
          gHat = pNDFT * f

          for pl = 1:length(packagesStr)
            planner = packagesCtor[pl]
            if planner != FINUFFT || σ == 2.0 # FINUFFT is not included in sigma sweep
              p = planner(x, N; m, σ, precompute=precomp[pl], blocking=blocking[pl])

              #if planner == CuNFFT.CuNFFTPlan
              #  fApprox = Array(adjoint(p) * CuNFFT.CuArray(fHat))
              #  gHatApprox = Array(p * CuNFFT.CuArray(f))
              #else
                fApprox = adjoint(p) * fHat
                gHatApprox = p * f
              #end

              eadjoint = norm(f[:] - fApprox[:]) / norm(f[:])

              etrafo = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])

              push!(df, (packagesStr[pl], D, M, N[D], m, σ, etrafo, eadjoint))
            end
        end
      end
    end
  end
  return df
end

function plot_accuracy_m(df, packagesStr, packagesStrShort, filename, D=1)

  σs = range(1.25, 4, length=12)

  df1_ = df[df.σ.==2.0 .&& df.D.==D,:]
  #df2_ = df[df.m.==8 .&& df.D.==D,:]

  Plots.scalefontsizes()
  Plots.scalefontsizes(1.5)

  c#olors = [:black, :orange, :green, :gray, :brown,  :blue, :purple, :yellow ]
  #ls = [:solid, :dashdot, :dash, :dashdotdot, :solid, :dash, :solid, :dash, :solid]
  #shape = [:xcross, :circle, :xcross, :cross, :circle, :xcross, :xcross, :circle]
  colors = [RGB(0.0,0.29,0.57), RGB(0.3,0.5,0.7), RGB(0.95,0.59,0.22), RGB(1.0,0.87,0.0)]
  ls = [:solid, :solid, :solid, :solid]
  shape = [:xcross, :circle, :xcross, :cross]

  p1 = plot(ms, df1_[df1_.Package.==packagesStr[1],:ErrorTrafo],
            yscale = :log10, label=packagesStrShort[1], lw=2, xlabel = "m", ylabel="Relative Error",
            legend = (:topright), title=L"\textrm{NFFT}", shape=:circle, c=:black)

  for p=2:length(packagesStr)
    plot!(p1, ms, df1_[df1_.Package.==packagesStr[p],:ErrorTrafo],
            yscale = :log10, label=packagesStrShort[p], lw=2, shape=shape[p], ls=ls[p],
            c=colors[p], msc=colors[p], mc=colors[p], ms=5, msw=2)
  end

  p2 = plot(ms, df1_[df1_.Package.==packagesStr[1],:ErrorAdjoint],
            yscale = :log10, lw=2, xlabel = "m", #ylabel="Relative Error",
            legend = nothing, title=L"\textrm{NFFT}^H", shape=:circle, c=:black)

  for p=2:length(packagesStr)
    plot!(p2, ms, df1_[df1_.Package.==packagesStr[p],:ErrorAdjoint],
            yscale = :log10, lw=2, shape=shape[p], ls=ls[p],
            c=colors[p], msc=colors[p], mc=colors[p], ms=5, msw=2)
  end

  p = plot(p1, p2, layout=(1,2), size=(800,300), dpi=200)
  #p = plot(p1, layout=(1,2), size=(800,450), dpi=200)

  mkpath("./img/")
  savefig(p, joinpath("./img/",filename))
  return p
end

function plot_accuracy_sigma(df, packagesStr, packagesStrShort, filename,  D=1)

  σs = range(1.25, 4, length=12)

  df1_ = df[df.m.==4 .&& df.D.==D,:]

  Plots.scalefontsizes()
  Plots.scalefontsizes(1.5)

  colors = [:black, :orange, :green, :gray, :brown,  :blue, :purple, :yellow ]
  ls = [:solid, :dashdot, :dash, :dashdotdot, :solid, :dash, :solid, :dash, :solid]
  shape = [:xcross, :circle, :xcross, :cross, :circle, :xcross, :xcross, :circle]

  p1 = plot(σs, df1_[df1_.Package.==packagesStr[1],:ErrorTrafo],
            yscale = :log10, label=packagesStrShort[1], lw=2, xlabel = L"\sigma", ylabel="Relative Error",
            legend = (:topright), title=L"\textrm{NFFT}", shape=:circle, c=:black)

  for p=2:length(packagesStr)
    plot!(p1, σs, df1_[df1_.Package.==packagesStr[p],:ErrorTrafo],
            yscale = :log10, label=packagesStrShort[p], lw=2, shape=shape[p], ls=ls[p],
            c=colors[p], msc=colors[p], mc=colors[p]) #ms=5, msw=2
  end

  p2 = plot(σs, df1_[df1_.Package.==packagesStr[1],:ErrorAdjoint],
            yscale = :log10, lw=2, xlabel = L"\sigma", #ylabel="Relative Error",
            legend = nothing, title=L"\textrm{NFFT}^H", shape=:circle, c=:black)

  for p=2:length(packagesStr)
    plot!(p2, σs, df1_[df1_.Package.==packagesStr[p],:ErrorAdjoint],
            yscale = :log10, lw=2, shape=shape[p], ls=ls[p],
            c=colors[p], msc=colors[p], mc=colors[p]) #ms=5, msw=2
  end

  p = plot(p1, p2, layout=(1,2), size=(800,300), dpi=200)

  mkpath("./img/")
  savefig(p, joinpath("./img/",filename))
  return p
end

#dfm = nfft_accuracy_comparison(2, [2.0], ms)
#dfσ = nfft_accuracy_comparison(2, σs, [4])

#writedlm("data/accuracy_m.csv", Iterators.flatten(([names(dfm)], eachrow(dfm))), ',')
#writedlm("data/accuracy_sigma.csv", Iterators.flatten(([names(dfσ)], eachrow(dfσ))), ',')

data, header = readdlm("data/accuracy_m.csv", ',', header=true);
dfm = DataFrame(data, vec(header))
data, header = readdlm("data/accuracy_sigma.csv", ',', header=true);
dfσ = DataFrame(data, vec(header))

plot_accuracy_m(dfm, ["NFFT.jl/TENSOR", "NFFT3/TENSOR", "FINUFFT"],
                     ["NFFT.jl", "NFFT3", "FINUFFT"], "accuracy_m_D2.pdf", 2)
plot_accuracy_m(dfm, ["NFFT.jl/FULL", "NFFT.jl/TENSOR", "NFFT.jl/LINEAR", "NFFT.jl/POLY"],
                     ["FULL", "TENSOR", "LINEAR", "POLYNOMIAL"],
                      "accuracy_m_pre_D2.pdf", 2)
plot_accuracy_sigma(dfσ, ["NFFT.jl/TENSOR", "NFFT3/TENSOR"],
                         ["NFFT.jl", "NFFT3"], "accuracy_sigma_D2.pdf", 2)

@info "Mean error deviation  NFFT.jl / NFFT3"
mean((dfm[dfm.Package.=="NFFT.jl/POLY",:ErrorTrafo] ./ dfm[dfm.Package.=="NFFT3/TENSOR",:ErrorTrafo])[1:5])

@info "Mean error deviation  FINUFFT / NFFT.jl"
mean((dfm[dfm.Package.=="FINUFFT",:ErrorTrafo] ./ dfm[dfm.Package.=="NFFT.jl/POLY",:ErrorTrafo])[1:5])


text:
#=
This code implements the VEM of order 0 for the computation of the resonance
modes of the Helmholtz equation. The method may be found in

"A virtual element method for the acoustic vibration problem"
Lourenço Beirão da Veiga, David Mora, Gonzalo Rivera, Rodolfo Rodríguez
https://arxiv.org/abs/1601.04316
=#

include("../utils.jl")  # some useful functions
include("../mesh.jl")   # defines a very general 2d mesh structure
include("../connectique.jl") # builds connectivity for the lower order method
include("../buildmatrices.jl") # builds the VEM matrices
include("../visu.jl") # visualization tools
include("../melem_rt0.jl") # elementary matrices with Raviart-Thomas of order 0

using LinearAlgebra  # package centralizing all linear algebra tools
using Arpack         # for the computation of the eigen values/vectors

# compute the pressure from u
function pressure(msh::mesh, u::Array{Float64, 1}, conn::Connectivity)
    p = zeros(Float64, length(msh.elt))
    for e in eachindex(msh.elt)
        E = polygonArea(msh.vtx[:,msh.elt[e]])
        for ee in conn.e2dof[e]
            sgnI = sign(ee)
            absI = abs(ee)
            p[e] += sgnI*u[absI]
        end
        p[e] /= -E
    end
    return p
end

function main()
    # read mesh
    println(" -> Reading mesh")
    # msh = readMesh("../meshes/voronoi.txt")
    # msh = readMesh("../meshes/squares.txt")
    # msh = readMesh("../meshes/triangles.txt")
    msh = readMesh("../meshes/non-convex.txt")
    # msh = readMesh("../meshes/smoothed-voronoi.txt")
    msh.vtx[2, :] .*= 1.1

    # connectivity
    println(" -> Building connectivity")
    conn = connectivity_k0(msh)

    # matrices
    println(" -> Assembling the matrices")
    A, B = buildmatrices_k0(msh, conn)
    # comment out below if you want to assemble with RT0; will only give a
    # correct solution for the mesh "triangles.txt"
    # A, B = buildRT0matrices(msh, conn)

    # apply boundary conditions
    mask = findall(msh.freeedg .== false)
    A = A[mask, mask]; # I vectorize because I am too lazy to write a loop
    B = B[mask, mask];

    # EV problem
    println(" -> Solving the generalized EV problem")
    Λ = eigen(A, B) # we get a structure with the eigenvalues and eigenvectors

    # plot eigenvectors
    it  = 0
    # index of the non-zero EV you wish to plot
    ind = 15
    for λ in sort(Λ.values)
        it += 1
        if λ > 1.05 # we only seek the non-zero eigenvalues
            uloc = Λ.vectors[:, it+ind]
            u    = zeros(Float64, conn.ndof)
            u[mask] .= uloc
            p    = pressure(msh, u, conn)
            plotSolutionElt(msh, p,
                myTitle="Pressure field - \$\\lambda=$(Λ.values[it+ind]-1)\$")
            #
            savefig("./EV_nonConvex.png")
            break
        end
    end
    println(Λ.values[it:(it+10)] .- 1)
end

main()


text:
# Covert from cli into KernelFunctions kernels

_default_range(::SqExponentialKernel) = 1:10
_default_range(::LinearKernel)        = -3:0.1:3
_default_range(::CategoricalKernel)   = [1,2,3]

_convert2kernel(k::SExp) = k.lengthscale == 1 ? SqExponentialKernel() : with_lengthscale(SqExponentialKernel(), k.lengthscale)
_convert2kernel(k::OU) = k.lengthscale == 1 ? ExponentialKernel() : with_lengthscale(ExponentialKernel(), k.lengthscale)
_convert2kernel(k::Linear) = LinearKernel(c=k.intercept)
_convert2kernel(::Cat) = CategoricalKernel()

_convert2kernel(::SExp, l) = l == 1 ? SqExponentialKernel() : with_lengthscale(SqExponentialKernel(), l)
_convert2kernel(::OU, l) = l == 1 ? ExponentialKernel() : with_lengthscale(ExponentialKernel(), l)
_convert2kernel(::Linear, c) = LinearKernel(; c)

_walk_kernel(ks::Union{KernelTensorProduct, KernelSum}) = reduce(vcat, [_walk_kernel(k) for k in ks.kernels])
_walk_kernel(ks::TransformedKernel) = _walk_kernel(ks.kernel)
_walk_kernel(k::KernelFunctions.SimpleKernel) = [k]

function _formula_pull_varnames(c::GPCompnent)
    return varname(c)
end

function _formula_pull_varnames(c::Tuple)
    reduce(vcat, _formula_pull_varnames.(c[2:3]))
end

function _convertop(op::Symbol)
    if op == :add
        return +
    elseif op == :multiply
        return *
    else
        throw(ArgumentError("Operation $op not yet supported"))
    end
end

_convert2eq(op::Symbol, c1, c2; hyperparams=Dict()) = _convert2eq(op, _convert2eq(c1; hyperparams), _convert2eq(c2; hyperparams))
_convert2eq(t::Tuple; hyperparams=Dict()) = _convert2eq(t...; hyperparams)
_convert2eq(c::GPCompnent; hyperparams=Dict()) = haskey(hyperparams, varname(c)) ? _convert2kernel(c, hyperparams[varname(c)]) : _convert2kernel(c)

_convert2eq(op::Symbol, c1::Kernel, c2::Kernel; hyperparams=Dict()) = _convertop(op)(c1, c2)

function _apply_select(formula::Tuple;  hyperparams=Dict())
    vars = _formula_pull_varnames(formula)
    varmap = Dict(var=>i for (i, var) in enumerate(unique(vars)))

end

function _apply_vars(formula::Tuple; hyperparams=Dict())
    vars = _formula_pull_varnames(formula)
    ks = _convert2eq(formula; hyperparams)
    retkernel = nothing
    counter = 1
    current_k = 1
    while current_k <= length(vars)
        k = ks.kernels[counter]
        if k isa KernelTensorProduct
            kernel = k ∘ SelectTransform(current_k:(current_k + length(k.kernels) - 1))
            retkernel = isnothing(retkernel) ? kernel : retkernel + kernel
            current_k += length(k.kernels)
            counter += 1
        elseif k isa KernelFunctions.SimpleKernel || k isa KernelFunctions.TransformedKernel
            kernel = k ∘ SelectTransform([current_k])
            retkernel = isnothing(retkernel) ? kernel : retkernel + kernel
            current_k += 1
            counter += 1
        else
            @show typeof(k)
            error()
        end
    end
    return retkernel, vars
end

_apply_vars(formula::GPCompnent; hyperparams=Dict()) = (_convert2eq(formula; hyperparams), [_formula_pull_varnames(formula)])

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
module Transforms
export  dwt, idwt, dwt!, idwt!,
        wpt, iwpt, wpt!, iwpt!,
        modwt, imodwt, cwt, icwt
using ..Util, ..WT
using FFTW

# TODO Use StridedArray instead of AbstractArray where writing to array.
# TODO change integer dependent wavelets to parametric types (see "Value types", https://docs.julialang.org/en/v1/manual/types/index.html#%22Value-types%22-1)
const DWTArray = AbstractArray
const WPTArray = AbstractVector
const ValueType = Union{AbstractFloat, Complex}

# DWT

"""
    dwt(x, wt[, L=maxtransformlevels(x)])

Perform a discrete wavelet transform of the array `x`.
The wavelet type `wt` determines the transform type
(filter or lifting) and the wavelet class, see `wavelet`.

The number of transformation levels `L` can be any non-negative
integer such that the size of `x` is divisible by `L`.
Performs the identity transformation if `L==0`.

# Examples
```julia
dwt(x, wavelet(WT.coif6))
```

**See also:** `idwt`, `dwt!`, `wavelet`
"""
function dwt end

"""
    idwt(x, wt[, L=maxtransformlevels(x)])

Perform an inverse discrete wavelet transform of the array `x`,
the inverse of `dwt(x, wt, L)`.

**See also:** `dwt`, `idwt!`
"""
function idwt end

"""
    dwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)])

    dwt!(y, wt::GLS[, L=maxtransformlevels(x)])

Same as `dwt` but without array allocation.
Perform "out of place" transform with a filter, or
a inplace transform with a lifting scheme. The difference
between the filter and lifting methods is due to the
structure of the transform algorithms.

**See also:** `idwt!`
"""
function dwt! end

"""
    idwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)])

    idwt!(y, wt::GLS[, L=maxtransformlevels(x)])

The inverse of `dwt!`.

**See also:** `dwt!`
"""
function idwt! end

# WPT

"""
    wpt

Perform a discrete wavelet packet transform of the array `x`.
**See also:** `dwt`, `wavelet`
"""
function wpt end

"""
    iwpt

Perform an inverse discrete wavelet packet transform of the array `x`.
**See also:** `idwt`, `wavelet`
"""
function iwpt end

"""
    wpt!

Same as `wpt` but without array allocation.
"""
function wpt! end

"""
    iwpt!

Same as `iwpt` but without array allocation.
"""
function iwpt! end

# DWT (discrete wavelet transform)

for (Xwt, Xwt!, _Xwt!, fw) in ((:dwt, :dwt!, :_dwt!, true),
                                (:idwt, :idwt!, :_dwt!, false))
@eval begin
    # filter
    function ($Xwt)(x::DWTArray{T}, filter::OrthoFilter,
                    L::Integer=maxtransformlevels(x)) where T<:ValueType
        y = Array{T}(undef, size(x))
        return ($_Xwt!)(y, x, filter, L, $fw)
    end
    function ($Xwt!)(y::DWTArray{T}, x::DWTArray{T}, filter::OrthoFilter,
                    L::Integer=maxtransformlevels(x)) where T<:ValueType
        return ($_Xwt!)(y, x, filter, L, $fw)
    end
    # lifting
    function ($Xwt)(x::DWTArray{T}, scheme::GLS,
                    L::Integer=maxtransformlevels(x)) where T<:ValueType
        y = Array{T}(undef, size(x))
        copyto!(y, x)
        return ($_Xwt!)(y, scheme, L, $fw)
    end
    function ($Xwt!)(y::DWTArray{T}, scheme::GLS,
                    L::Integer=maxtransformlevels(y)) where T<:ValueType
        return ($_Xwt!)(y, scheme, L, $fw)
    end
end # begin
end # for

# CWT (continuous wavelet transform)
# cwt(Y::AbstractVector, ::ContinuousWavelet)

@doc """
     wave = cwt(Y::AbstractArray{T,N}, c::CFW{W, S, WT}, daughters, rfftPlan =
             plan_rfft([1]), fftPlan = plan_fft([1])) where {N, T<:Real,
                                                             S<:Real,
                                                             U<:Number,
                                                             W<:WT.WaveletBoundary,
                                                             WT<:Union{<:WT.Morlet,
                                                                       <:WT.Paul}}

  return the continuous wavelet transform along the first axis with averaging.
  `wave`, is (signalLength)×(nscales+1)×(previous dimensions), of type T of
  Y. averagingLength defines the number of octaves (powers of 2) that are
  replaced by an averaging function. This has form averagingType, which can be
  one of `Mother()` or `Dirac()`- in the `Mother()` case, it uses the same form
  as for the wavelets, while the `Dirac` uses a constant window. J1 is the
  total number of scales; default (when J1=NaN, or is negative) is just under
  the maximum possible number, i.e. the log base 2 of the length of the signal,
  times the number of wavelets per octave. If you have sampling information,
  you will need to scale wave by δt^(1/2).

  """
function cwt(Y::AbstractArray{T,N}, c::CFW{W, S, WaTy}, daughters, rfftPlan::AbstractFFTs.Plan =
             plan_rfft([1]), fftPlan = plan_fft([1])) where {N, T<:Real,
                                                             S<:Real,
                                                             W<:WT.WaveletBoundary,
                                                             WaTy<:Union{<:WT.Morlet,
                                                                         <:WT.Paul}}
    # This is for analytic wavelets, so we need to treat the positive and
    # negative frequencies differently, even for real data

    # TODO: complex input version of this
    @assert typeof(N)<:Integer
    # vectors behave a bit strangely, so we reshape them
    if N==1
        Y= reshape(Y,(length(Y), 1))
    end
    n1 = size(Y, 1);

    _, nScales, _ = getNWavelets(n1, c)
    @debug "" nScales
    #....construct time series to analyze, pad if necessary
    x = reflect(Y, boundaryType(c)()) #this function is defined below

    # check if the plans we were given are dummies or not
    if size(rfftPlan)==(1,)
        rfftPlan = plan_rfft(x, 1)
    end
    if size(fftPlan)==(1,)
        fftPlan = plan_fft(x, 1)
    end
    n = size(x, 1)

    x̂ = rfftPlan * x
    # If the vector isn't long enough to actually have any other scales, just
    # return the averaging
    if nScales <= 0 || size(daughters,2) == 1
        daughters = daughters[:,1:1]
        nScales = 1
    end

    wave = zeros(Complex{T}, size(x)..., nScales);  # result array;
    # faster if we put the example index on the outside loop through all scales
    # and compute transform
    actuallyTransform!(wave, daughters,x̂, fftPlan, c.waveType, c.averagingType)
    wave = permutedims(wave, [1, ndims(wave), (2:(ndims(wave)-1))...])
    ax = axes(wave)
    wave = wave[1:n1, ax[2:end]...]
    if N==1
        wave = dropdims(wave, dims=3)
    end

    return wave
end

function cwt(Y::AbstractArray{T,N}, c::CFW{W, S, WaTy}, daughters, rfftPlan =
             plan_rfft([1])) where {N, T<:Real, S<:Real, U<:Number,
                                    W<:WT.WaveletBoundary, WaTy<:WT.Dog}
    # Dog doesn't need a fft because it is strictly real
    # TODO: complex input version of this
    @assert typeof(N)<:Integer
    # vectors behave a bit strangely, so we reshape them
    if N==1
        Y= reshape(Y,(length(Y), 1))
    end

    n1 = size(Y, 1);

    _, nScales, _ = getNWavelets(n1, c)
    #....construct time series to analyze, pad if necessary
    x = reflect(Y, boundaryType(c)())

    # check if the plans we were given are dummies or not
    if size(rfftPlan)==(1,)
        rfftPlan = plan_rfft(x, 1)
    end
    n = size(x, 1)

    # If the vector isn't long enough to actually have any other scales, just
    # return the averaging. Or if there's only averaging
    if nScales <= 1 || size(daughters,2) == 1
        daughters = daughters[:,1:1]
        nScales = 1
    end

    x̂ = rfftPlan * x

    wave = zeros(Complex{T}, size(x)..., nScales);  # result array;
    # faster if we put the example index on the outside loop through all scales
    # and compute transform

    actuallyTransform!(wave, daughters,x̂, rfftPlan, c.waveType)
    wave = permutedims(wave, [1, ndims(wave), (2:(ndims(wave)-1))...])
    ax = axes(wave)
    wave = wave[1:n1, ax[2:end]...]

    if N==1
        wave = dropdims(wave, dims=3)
    end

    return real.(wave)
end

function reflect(Y, bt)
    n1 = size(Y, 1)
    if bt == WT.padded
        base2 = round(Int,log(n1)/log(2));   # power of 2 nearest to N
        x = cat(Y, zeros(2^(base2+1)-n1, size(Y)[2:end]...), dims=1)
    elseif bt == WT.DEFAULT_BOUNDARY
        x = cat(Y, reverse(Y,dims = 1), dims = 1)
    else
        x = Y
    end
    return x
end

function actuallyTransform!(wave, daughters, x̂, fftPlan, analytic::Union{<:WT.Morlet,
                                                                         <:WT.Paul},
                            averagingType::Union{WT.Father, WT.Dirac})
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    isSourceOdd = mod(size(wave,1)+1,2)
    # the averaging function isn't analytic, so we need to do both positive and
    # negative frequencies
    tmpWave = x̂ .* daughters[:,1]
    wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceOdd,
                                                          outer...]),dims=1)
    wave[1:n1, outer..., 1] = tmpWave
    wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # wavelet transform
    for j in 2:size(daughters,2)
        wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end
function actuallyTransform!(wave, daughters, x̂, fftPlan, analytic::Union{<:WT.Morlet,
                                                                         <:WT.Paul},
                            averagingType::WT.NoAve)
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    # the no averaging version
    for j in 1:size(daughters,2)
        wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function actuallyTransform!(wave, daughters, x̂, rfftPlan, analytic::Union{<:WT.Dog})
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    for j in 1:size(daughters,2)
        wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
        wave[:, outer..., j] = rfftPlan \ (wave[1:n1, outer..., j])  # wavelet transform
    end
end

function cwt(Y::AbstractArray{T}, c::CFW{W}; varArgs...) where {T<:Number, S<:Real, V<: Real,
                                                                                W<:WT.WaveletBoundary}
    daughters,ω = computeWavelets(size(Y, 1), c; varArgs...)
    return cwt(Y, c, daughters)
end

"""
period,scale, coi = caveats(Y::AbstractArray{T}, c::CFW{W}; J1::S=NaN) where {T<:Real, S<:Real, W<:WT.WaveletBoundary}

returns the period, the scales, and the cone of influence for the given wavelet transform. If you have sampling information, you will need to scale the vector scale appropriately by 1/δt, and the actual transform by δt^(1/2).
"""
function caveats(n1, c::CFW{W}; J1::Int64=-1, dt::S=NaN, s0::V=NaN) where {S<:Real, W<:WT.WaveletBoundary, V <: Real}
    # don't alter scaling with sampling information if it doesn't exists
    fλ = (4*π) / (c.σ[1] + sqrt(2 + c.σ[1]^2))
    if isnan(dt) || (dt<0)
        dt = 1
    end

    # smallest resolvable scale
    if isnan(s0) || (s0<0)
        s0 = 2 * dt / fλ
    end
    sj = s0 * 2.0.^(collect(0:J1)./c.scalingFactor)
    # Fourier equivalent frequencies
    freqs = 1 ./ (fλ .* sj)

    # Determines the cone-of-influence. Note that it is returned as a function
    # of time in Fourier periods. Uses triangualr Bartlett window with
    # non-zero end-points.
    coi = (n1 / 2 .- abs.(collect(0:n1-1) .- (n1 - 1) ./ 2))
    coi = (fλ * dt / sqrt(2)).*coi

    n1 = length(Y);
    # J1 is the total number of elements
    if isnan(J1) || (J1<0)
        J1=floor(Int,(log2(n1))*c.scalingFactor);
    end
    #....construct time series to analyze, pad if necessary
    if boundaryType(c) == WT.ZPBoundary
        base2 = round(Int,log(n1)/log(2));   # power of 2 nearest to N
        n = length(Y)+2^(base2+1)-n1
    elseif boundaryType(c) == WT.PerBoundary
        n = length(Y)*2
    end
    ω = [0:floor(Int, n/2); -floor(Int,n/2)+1:-1]*2π
    period = c.fourierFactor*2 .^((0:J1)/c.scalingFactor)
    scale = [1E-5; 1:((n1+1)/2-1); reverse((1:(n1/2-1)),dims=1); 1E-5]
    coi = c.coi*scale  # COI [Sec.3g]
    return sj, freqs, period, scale, coi
end
cwt(Y::AbstractArray{T}, w::WT.ContinuousWaveletClass; J1::Int64=-1, dt::S=NaN, s0::V=NaN) where {T<:Real, S<:Real, V<:Real} = cwt(Y,CFW(w),J1=J1,dt=dt,s0=s0)
caveats(Y::AbstractArray{T}, w::WT.ContinuousWaveletClass; J1::S=NaN) where {T<: Real, S<: Real} = caveats(Y,CFW(w),J1=J1)
cwt(Y::AbstractArray{T}) where T<:Real = cwt(Y,WT.Morlet())
caveats(Y::AbstractArray{T}) where T<:Real = caveats(Y,WT.Morlet())

"""
    icwt(W::AbstractArray{T}, c::CFW{W}, sj::AbstractArray; dt::S=NaN,
         dj::V=1/12) where {S<:Real, V<:Real, W<:WT.WaveletBoundary}
    icwt(W::AbstractArray{T}, c::CFW{W}; dt::S=NaN, s0::S
         dj::V=1/12) where {T<:Real, S<:Real, V<:Real,
                            W<:WT.WaveletBoundary}

return the inverse continuous wavelet transform
"""
function icwt(W::AbstractArray, c::CFW, sj::AbstractArray; dt::S=NaN,
              dj::V=1/12) where {S<:Real, V<:Real}
    if isnan(dt) || (dt<0)
        dt = 1
    end

    # Torrence and Compo (1998), eq. (11)
    n = WT.setn(size(W,1), c)
    ω = (0:(n-1))*2π
    ψ = WT.Mother(c, 1, 1, ω)[1:(end-1)]
    println("size(sj) = $(size(sj))")
    println("size(W) = $(size(W))")
    println("size(ψ) = $(size(ψ))")
    iW = (dj * sqrt(dt) / 0.776 * psi(c,0)) .* sum((real.(W) ./ sqrt.(sj')), dims=2)

    return iW
end
function icwt(W::AbstractArray, c::CFW; dt::Real=NaN, dj::Real=1/12, J1::Real=NaN)

    fλ = (4*π) / (c.σ[1] + sqrt(2 + c.σ[1]^2))
    n1 = size(W, 1);
    # J1 is the total number of elements
    if isnan(J1) || (J1<0)
        J1=floor(Int,(log2(n1))*c.scalingFactor);
    end

    sj =  WT.getScales(n1, c)

    if isnan(dt) || (dt<0)
        dt = 1
    end

    return icwt(W, c, sj, dt=dt, dj=1/(c.scalingFactor))
end
icwt(Y::AbstractArray, w::WT.ContinuousWaveletClass; dj::T=1/12, dt::S=NaN,
     s0::V=NaN) where {S<:Real, T<:Real, V<:Real} = icwt(Y,CFW(w))
icwt(Y::AbstractArray) = icwt(Y,WT.Morlet())

function psi(c::CFW{W}, t::Int64) where W<:WT.WaveletBoundary
    return real.(π^(-0.25) * exp.(im*c.σ[1]*t - t^2 / 2))
end

# WPT (wavelet packet transform)

for (Xwt, Xwt!, _Xwt!, fw) in ((:wpt, :wpt!, :_wpt!, true),
                                (:iwpt, :iwpt!, :_wpt!, false))
@eval begin
    function ($Xwt)(x::WPTArray{T}, wt::DiscreteWavelet,
                    L::Integer=maxtransformlevels(x)) where T<:ValueType
        return ($Xwt)(x, wt, maketree(length(x), L, :full))
    end
    # filter
    function ($Xwt)(x::WPTArray{T}, filter::OrthoFilter,
                    tree::BitVector=maketree(x, :full)) where T<:ValueType
        y = Array{T}(undef, size(x))
        return ($_Xwt!)(y, x, filter, tree, $fw)
    end
    function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T},
                    filter::OrthoFilter) where T<:ValueType
        return ($Xwt!)(y, x, filter, maketree(x, :full))
    end
    function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T},
                    filter::OrthoFilter, L::Integer) where T<:ValueType
        return ($Xwt!)(y, x, filter, maketree(length(x), L, :full))
    end
    function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T},
                    filter::OrthoFilter, tree::BitVector) where T<:ValueType
        return ($_Xwt!)(y, x, filter, tree, $fw)
    end
    # lifting
    function ($Xwt)(x::WPTArray{T}, scheme::GLS,
                    tree::BitVector=maketree(x, :full)) where T<:ValueType
        y = Array{T}(undef, size(x))
        copyto!(y, x)
        return ($_Xwt!)(y, scheme, tree, $fw)
    end
    function ($Xwt!)(y::WPTArray{T}, scheme::GLS) where T<:ValueType
        return ($Xwt!)(y, scheme, maketree(y, :full))
    end
    function ($Xwt!)(y::WPTArray{T}, scheme::GLS, L::Integer) where T<:ValueType
        return ($Xwt!)(y, scheme, maketree(length(x), L, :full))
    end
    function ($Xwt!)(y::WPTArray{T}, scheme::GLS, tree::BitVector) where T<:ValueType
        return ($_Xwt!)(y, scheme, tree, $fw)
    end
end # begin
end # for

# DWTC (column-wise discrete wavelet transform)
#dwtc(::AbstractArray, ::DiscreteWavelet)
#idwtc(::AbstractArray, ::DiscreteWavelet)

# SWT (stationary wavelet transform)
#swt(::AbstractVector, ::DiscreteWavelet)
#iswt(::AbstractVector, ::DiscreteWavelet)

# CWT (continuous wavelet transform directly) TODO: direct if sufficiently small

# TODO: continuous inverse, when defined
#icwt(::AbstractVector, ::ContinuousWavelet)

# Int -> Float
for Xwt in (:dwt, :idwt, :dwtc, :idwtc, :wpt, :iwpt)
    @eval $Xwt(x::AbstractArray{<:Integer}, args...) = $Xwt(float(x), args...)
end

# non-exported "out of place" functions
for (Xwt_oop!, Xwt!) in ((:dwt_oop!, :dwt!), (:idwt_oop!, :idwt!))
@eval begin
    # filter
    function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, filter::OrthoFilter,
                        L::Integer=maxtransformlevels(x)) where T<:ValueType
        return ($Xwt!)(y, x, filter, L)
    end
    # lifting
    function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, scheme::GLS,
                        L::Integer=maxtransformlevels(x)) where T<:ValueType
        copyto!(y, x)
        return ($Xwt!)(y, scheme, L)
    end
end # begin
end # for

# Array with shared memory
function unsafe_vectorslice(A::Array{T}, i::Int, n::Int) where T
    return unsafe_wrap(Array, pointer(A, i), n)::Vector{T}
end

# linear indices of start of rows/cols/planes
# 2-D, size(A) = (m,n)
row_idx(i, m) = i
col_idx(i, m) = 1 + (i-1)*m
# 3-D, size(A) = (m,n,d)
row_idx(i, j, m, n=m) = row_idx(i, n) + (j-1)*n*m
col_idx(i, j, m, n=m) = col_idx(i, m) + (j-1)*n*m
plane_idx(i, j, m) = i + (j-1)*m

# filter transforms
include("transforms_filter.jl")

# lifting transforms
include("transforms_lifting.jl")

include("transforms_maximal_overlap.jl")

end # module


text:
using MPI, Test
include("../testhelpers.jl")

@testset "Balance Law Solver" begin
  tests = [(1, "Euler/isentropic_vortex_standalone.jl")
           (1, "Euler/isentropic_vortex_standalone_aux.jl")
           (1, "Euler/isentropic_vortex_standalone_source.jl")
           (1, "Euler/isentropic_vortex_standalone_bc.jl")
           # (1, "Euler/isentropic_vortex_standalone_integral.jl")
           (1, "Euler/isentropic_vortex_standalone_IMEX.jl")
           (1, "util/filter_test.jl")
           (1, "util/grad_test.jl")
           (1, "util/grad_test_sphere.jl")
           (1, "util/integral_test.jl")
           (1, "util/integral_test_sphere.jl")
           (1, "util/nonnegative_test.jl")
           (1, "conservation/sphere.jl")
           (1, "compressible_Navier_Stokes/mms_bc.jl")
           (1, "sphere/advection_sphere_lsrk.jl")
           (1, "sphere/advection_sphere_ssp33.jl")
           (1, "sphere/advection_sphere_ssp34.jl")
          ]

  runmpi(tests, @__FILE__)

  if "linux" != lowercase(get(ENV,"TRAVIS_OS_NAME",""))
    moretests = [(1, "../../examples/DGmethods_old/ex_001_periodic_advection.jl")
                 (1, "../../examples/DGmethods_old/ex_002_solid_body_rotation.jl")
                 (1, "../../examples/DGmethods_old/ex_003_acoustic_wave.jl")
                 (1, "../../examples/DGmethods_old/ex_004_nonnegative.jl")
                ]
    runmpi(moretests, @__FILE__)
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
module DZMisc

export rmk, say, orthonormalize_columns!, linearly_independent_column_indices!

using Base.Threads: lock, unlock, SpinLock
using InteractiveUtils: _dump_function

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

const DZMISC_STDOUT_LOCK = SpinLock()

function rmk(args...; verbose::Bool=true)::Nothing
    if verbose
        lock(DZMISC_STDOUT_LOCK)
        print("\33[2K\r")
        print(args...)
        flush(stdout)
        unlock(DZMISC_STDOUT_LOCK)
    end
end

function say(args...; verbose::Bool=true)::Nothing
    if verbose
        lock(DZMISC_STDOUT_LOCK)
        print("\33[2K\r")
        println(args...)
        flush(stdout)
        unlock(DZMISC_STDOUT_LOCK)
    end
end

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

function next_permutation!(items::Vector{T})::Bool where {T}
    num_items = length(items)
    if num_items == 0; return false; end
    current_item = items[num_items]
    pivot_index = num_items - 1
    while pivot_index != 0
        next_item = items[pivot_index]
        if next_item >= current_item
            pivot_index -= 1
            current_item = next_item
        else; break; end
    end
    if pivot_index == 0; return false; end
    pivot = items[pivot_index]
    successor_index = num_items
    while items[successor_index] <= pivot; successor_index -= 1; end
    items[pivot_index], items[successor_index] =
        items[successor_index], items[pivot_index]
    reverse!(view(items, pivot_index + 1 : num_items))
    return true
end

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

function orthonormalize_columns!(mat::Matrix{T})::Nothing where {T <: Real}
    m, n = size(mat, 1), size(mat, 2)
    for j = 1 : n
        let # normalize j'th column
            acc = zero(T)
            @simd for i = 1 : m
                @inbounds acc += abs2(mat[i, j])
            end
            acc = inv(sqrt(acc))
            @simd ivdep for i = 1 : m
                @inbounds mat[i, j] *= acc
            end
        end
        for k = j + 1 : n # orthogonalize k'th column against j'th column
            acc = zero(T)
            @simd for i = 1 : m
                @inbounds acc += mat[i, j] * mat[i, k]
            end
            @simd ivdep for i = 1 : m
                @inbounds mat[i, k] -= acc * mat[i, j]
            end
        end
    end
end

function linearly_independent_column_indices!(
        mat::Matrix{T}, threshold::T) where {T <: Real}
    indices = Int[]
    lo, hi, m, n = zero(T), T(Inf), size(mat, 1), size(mat, 2)
    for i = 1 : n
        x = zero(T)
        @simd for k = 1 : m
            @inbounds x += abs2(mat[k, i])
        end
        x = sqrt(x)
        if x > threshold
            hi = min(hi, x)
            push!(indices, i)
            y = inv(x)
            @simd ivdep for k = 1 : m
                @inbounds mat[k, i] *= y
            end
            for j = i + 1 : n
                acc = zero(T)
                @simd for k = 1 : m
                    @inbounds acc += mat[k, i] * mat[k, j]
                end
                @simd ivdep for k = 1 : m
                    @inbounds mat[k, j] -= acc * mat[k, i]
                end
            end
        else
            lo = max(lo, x)
        end
    end
    indices, lo, hi
end

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

const MOV_TYPES = ["mov", "movabs", "vmovaps", "vmovups", "vmovapd", "vmovupd"]
const JUMP_TYPES = ["je", "jne", "ja", "jae", "jb", "jbe",
                    "jg", "jge", "jl", "jle"]
const JUMP_PREFIXES = Dict("je" => "", "jne" => "",
                           "ja" => "(unsigned) ", "jae" => "(unsigned) ",
                           "jb" => "(unsigned) ", "jbe" => "(unsigned) ",
                           "jg" => "(signed) ", "jge" => "(signed) ",
                           "jl" => "(signed) ", "jle" => "(signed) ")
const JUMP_SYMBOLS = Dict("je" => " == ", "jne" => " != ",
                          "ja" => " > ", "jae" => " >= ",
                          "jb" => " < ", "jbe" => " <= ",
                          "jg" => " > ", "jge" => " >= ",
                          "jl" => " < ", "jle" => " <= ")
const SHUFFLE_INSTRUCTIONS = ["vunpcklpd", "vunpckhpd", "vperm2f128",
    "vblendpd", "vshufpd", "vpermpd", "vpermilpd"]

function view_asm(@nospecialize(func), @nospecialize(types))

    code_lines = asm_lines(func, types)
    for i = 1 : length(code_lines)
        if (i < length(code_lines) && code_lines[i][1] == "cmp"
                                   && code_lines[i+1][1] in JUMP_TYPES)
            args = split(join(code_lines[i][2:end], ' '), ", ")
            @assert length(args) == 2
            @assert length(code_lines[i+1]) == 2
            prefix = JUMP_PREFIXES[code_lines[i+1][1]]
            symbol = JUMP_SYMBOLS[code_lines[i+1][1]]
            code_lines[i] = [">>>>if (" * prefix * args[1] * symbol * args[2] *
                             ") goto " * code_lines[i+1][2] * ";"]
            code_lines[i+1] = ["nop"]
        end
        if (i < length(code_lines) && code_lines[i][1] == "test"
                                   && code_lines[i+1][1] in JUMP_TYPES)
            args = split(join(code_lines[i][2:end], ' '), ", ")
            @assert length(args) == 2
            @assert length(code_lines[i+1]) == 2
            prefix = JUMP_PREFIXES[code_lines[i+1][1]]
            symbol = JUMP_SYMBOLS[code_lines[i+1][1]]
            if (args[1] == args[2])
                code_lines[i] = [">>>>if (" * prefix * args[1] * symbol *
                                 "0) goto " * code_lines[i+1][2] * ";"]
                code_lines[i+1] = ["nop"]
            end
        end
    end

    unknown_instrs = Dict{String,Int}()
    for line in code_lines

        # Preprocessed lines
        if length(line) == 1 && startswith(line[1], ">>>>")
            say("    ", line[1][5:end])

        # Unprinted lines
        elseif length(line) == 1 && line[1] == ".text"
            # Ignore this line.
        elseif line[1] in ["nop", "push", "pop", "vzeroupper"]
            # Ignore this line.

        # Labels
        elseif length(line) == 1 && endswith(line[1], ':')
            say(line[1])

        # Moves
        elseif line[1] in MOV_TYPES
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            say("    ", line[1], " = ", line[2], ';')
        elseif line[1] == "lea"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            @assert startswith(line[2], '[') && endswith(line[2], ']')
            say("    ", line[1], " = ", line[2][2:end-1], ';')
        elseif line[1] == "seto"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 1
            say("    ", line[1], " = OF;")

        # Increment and decrement
        elseif line[1] == "inc"
            if length(line) == 2
                say("    ++", line[2], ";")
            else
                say("    ++(", join(line[2:end], ' '), ");")
            end
        elseif line[1] == "dec"
            if length(line) == 2
                say("    --", line[2], ";")
            else
                say("    --(", join(line[2:end], ' '), ");")
            end

        # Scalar arithmetic and shifts
        elseif line[1] == "add"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            say("    ", line[1], " += ", line[2], ';')
        elseif line[1] == "imul"
            line = split(join(line[2:end], ' '), ", ")
            if length(line) == 1
                say("    rdx:rax = (signed) ", line[1], " * rax;")
            elseif length(line) == 2
                say("    ", line[1], " *= (signed) ", line[2], ';')
            elseif length(line) == 3
                say("    ", line[1], " = (signed) ", line[2],
                    " * ", line[3], ';')
            else
                @assert false
            end
        elseif line[1] == "sub"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            say("    ", line[1], " -= ", line[2], ';')
        elseif line[1] == "and"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            say("    ", line[1], " &= ", line[2], ';')
        elseif line[1] == "xor"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 2
            if line[1] == line[2]
                say("    ", line[1], " = 0;")
            else
                say("    ", line[1], " ^= ", line[2], ';')
            end
        elseif line[1] == "shl" || line[1] == "sal"
            line = split(join(line[2:end], ' '), ", ")
            if length(line) == 1
                say("    ", line[1], " <<= 1;")
            elseif length(line) == 2
                say("    ", line[1], " <<= ", line[2], ';')
            else
                @assert false
            end
        elseif line[1] == "shr"
            line = split(join(line[2:end], ' '), ", ")
            if length(line) == 1
                say("    ", line[1], " >>>= 1;")
            elseif length(line) == 2
                say("    ", line[1], " >>>= ", line[2], ';')
            else
                @assert false
            end
        elseif line[1] == "sar"
            line = split(join(line[2:end], ' '), ", ")
            if length(line) == 1
                say("    ", line[1], " >>= 1;")
            elseif length(line) == 2
                say("    ", line[1], " >>= ", line[2], ';')
            else
                @assert false
            end
        elseif line[1] == "andn"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ~", line[2], " & ", line[3], ';')

        # Vector arithmetic
        elseif line[1] == "vmovsd"
            line = split(join(line[2:end], ' '), " # ")
            @assert 1 <= length(line) <= 2
            line = split(line[1], ", ")
            @assert length(line) == 2
            say("    ", line[1], " = ", line[2], "; // scalar")
        elseif line[1] == "vaddsd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " + ", line[3], "; // scalar")
        elseif line[1] == "vaddpd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " + ", line[3], ';')
        elseif line[1] == "vsubsd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " - ", line[3], "; // scalar")
        elseif line[1] == "vsubpd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " - ", line[3], ';')
        elseif line[1] == "vmulsd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " * ", line[3], "; // scalar")
        elseif line[1] == "vmulpd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " * ", line[3], ';')
        elseif line[1] == "vdivpd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            say("    ", line[1], " = ", line[2], " / ", line[3], ';')
        elseif line[1] == "vxorpd"
            line = split(join(line[2:end], ' '), ", ")
            @assert length(line) == 3
            if line[2] == line[3]
                say("    ", line[1], " = 0.0;")
            else
                say("    ", line[1], " = ", line[2], " / ", line[3], ';')
            end
        elseif line[1] in SHUFFLE_INSTRUCTIONS
            line = split(join(line[2:end], ' '), " # ")
            @assert length(line) == 2
            say("    ", replace(line[2], "]," => "], "), ';')

        # Control flow
        elseif line[1] == "jmp"
            @assert length(line) == 2
            say("    goto ", line[2], ';')
        elseif line[1] == "call"
            if length(line) == 2
                say("    ", line[2], "();")
            else
                say("    (", join(line[2:end], ' '), ")();")
            end
        elseif line[1] == "ret"
            say("    return;")
        elseif line[1] == "ud2"
            say("    <unreachable code>")

        # Unknown instructions
        else
            say("    {", join(line, ' '), '}')
            unknown_instrs[line[1]] = get(unknown_instrs, line[1], 0) + 1
        end
    end

    if length(unknown_instrs) > 0
        say()
        freqs = sort!(reverse.(collect(unknown_instrs)), rev=true)
        total = sum(f[1] for f in freqs)
        say("$(total) unknown instructions:")
        for i = 1 : min(length(freqs), 10)
            freq, name = freqs[i]
            say("    ", name, " (", freq, ")")
        end
    end
end

end # module DZMisc


text:
function _do_serial_mem(::Type{M}, ::Type{D}, filename::String, count_mode::CountMode, i) where {M<:AbstractMer,D<:ReadDatastore}
    ds = open(D, filename)
    part = i:nworkers():length(ds)
    counts = serial_mem(M, buffer(ds), count_mode, part)
    close(ds)
    return counts
end

function _dist_mem(::Type{M}, ::Type{D}, filename::String, count_mode::CountMode) where {M<:AbstractMer,D<:ReadDatastore}
    nw = nworkers()
    @info "Splitting all reads in $filename, across $nw worker processes and counting kmers"
    if nw === 1
        @warn "There are not any extra worker processes available. This method will not provide any speed benefit over serial_mem."
    end
    local_counts = Vector{Vector{MerCount{M}}}(undef, nw)
    @sync begin
        @inbounds for (i, w) in enumerate(workers())
            @async local_counts[i] = remotecall_fetch(_do_serial_mem, w, M, D, filename, count_mode, i)
        end
        @info "Waiting for workers to finish"
    end
    @info "Merging all counts from $nw worker processes"
    final_count = local_counts[1]
    @inbounds for i in 2:lastindex(local_counts)
        unsafe_merge_into!(final_count, local_counts[i])
    end
    return final_count
end

function _dist_mem(::Type{M}, input::ReadDatastore, count_mode::CountMode) where {M<:AbstractMer}
    return _dist_mem(M, typeof(input), input.filename, count_mode)
end

function _dist_mem(::Type{M}, input::DatastoreBuffer{<:ReadDatastore}, count_mode::CountMode) where {M<:AbstractMer}
    return _dist_mem(M, ReadDatastores.datastore(input), count_mode)
end

"""
The multi-process parallel version of the [`SerialMemCounter`](@ref) k-mer counter.

Like the [`SerialMemCounter`](@ref) this counter does all counting and processing
in memory, so it is not designed to limit memory use, but if using a cluster of
machines the memory used per worker machine would be less. Of course the data is
copied back over to the master process and the results combined and so the master
node has to have enough RAM available to hold all the result.
"""
struct DistMemCounter{M<:AbstractMer,C<:CountMode} <: AbstractKmerCounter{M}
end

function dist_mem(::Type{M}, mode::C) where {M<:AbstractMer,C<:CountMode}
    return DistMemCounter{M,C}()
end

function (dmc::DistMemCounter{M,C})(input) where {M<:AbstractMer,C}
    return _dist_mem(M, input, C())
end


text:
"""
Mod block
"""

export ModBlock

mutable struct ModBlock <: AbstractBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}

    function ModBlock() where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}()) 2 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::ModBlock)
    println(io, "Mod()")
end

"""
to expr
"""

function _toexpr(blk::ModBlock)
    Expr(:call, :%, _toexpr(blk.inport[1]), _toexpr(blk.inport[2]))
end


text:
function get_saved_searches(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/saved_searches/list.json", options))
end

function get_saved_searches_show_id(auth::TwAuth, options::Dict)
  id = options["id"]
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/saved_searches/show/$id.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
using TimeseriesPrediction

using DynamicalSystemsBase

const STE = SpatioTemporalEmbedding

ti = time()

include("localmodeling_tests.jl")
include("st_pred_barkley.jl")
include("st_coupledhenon_tests.jl")
include("reconstruction_tests.jl")
include("symmetry_tests.jl")
ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits=3), " seconds or ", round(ti/60, digits=3), " minutes")


text:
function double_expectation(numsamples=1000;seed=1)

    rg = MersenneTwister(seed)

    N = 5
    F = 2
    idx = ceil.(Int, rand(rg, N)*F)

    μf = rand(rg,N)
    Σf = randn(rg,N,N) ; Σf = Σf*Σf'
    qf = MvNormal(μf, Σf)
    f  = rand(rg, qf)

    μα = randn(rg, F); Σα = Diagonal(rand(rg, F))
    qα = MvNormal(μα, Σα)

    S = Diagonal(rand(rg,N) .+ 0.1);

    # first verify that this holds for diagonal matrices
    A = Diagonal(rand(rg,N))
    f'*A'*inv(S)*A*f, sum(f.^2 .* diag(A).^2 ./ diag(S))  # ✅

    # test inner expectation
    innerexpectation = μf'*A'*inv(S)*A*μf + tr(A*inv(S)*A*Σf)

    # numinnerexpectation = @distributed (+) for i=1:numsamples
    #
    #     local f = rand(qf)
    #
    #     f'*A'*inv(S)*A*f / numsamples
    #
    # end
    #
    # innerexpectation, numinnerexpectation # ✅

    # test double expectation
    μa = μα[idx]
    Σa = Diagonal(Σα.diag[idx])

    doubleexpectation  = sum(μf.^2 .* (μa.^2 .+ Σa.diag) ./ S.diag) + sum((μa.^2 .+ Σa.diag).*diag(Σf) ./S.diag)

    shorterexpectation = sum((μf.^2 .+ diag(Σf)) .* (μa.^2 .+ Σa.diag) ./ S.diag)

    numdoubleexpectation = 0#@distributed (+) for i=1:numsamples
    #
    #     local A = Diagonal(rand(qa))
    #
    #     (μf'*A'*(S\A*μf) + tr(A*inv(S)*A*Σf)) / numsamples
    #
    # end

    altnumdoubleexpectation = @distributed (+) for i=1:numsamples

        local A = Diagonal(rand(qα)[idx])
        local f = rand(qf)

        f'*A'*(S\A*f) / numsamples

    end

    shorterexpectation, doubleexpectation, numdoubleexpectation,altnumdoubleexpectation # ✅

end


text:
import Distributions
import MonteCarloMeasurements
import NamedTupleTools

function Base.rand(sp::SossMLJPredictor{M};
                   response = sp.model.response) where {M}
    pars = rand(sp.post)
    args = merge(sp.args, pars)
    return getproperty(rand(sp.pred(args)), response)
end

function Distributions.logpdf(sp::SossMLJPredictor{M}, y;
        response = sp.model.response) where {M}
    # Get all the distribution mixture components
    dists = Base.Generator(sp.post) do pars
        args = merge(sp.args, pars)
        sp.pred(args)
    end

    y_namedtuple = NamedTupleTools.namedtuple(response)(tuple(y))
    # Evaluate logpdf(d,x) on each component d
    logvals = Base.Generator((d -> Distributions.logpdf(d, y_namedtuple)) ∘ dists.f, dists.iter)
    n = length(sp.post)
    return Distributions.logsumexp(logvals) - log(n)
end


text:
"""
```
    LinearStretching <: AbstractHistogramAdjustmentAlgorithm
    LinearStretching(; minval = 0, maxval = 1)

    adjust_histogram([T,] img, f::LinearStretching)
    adjust_histogram!([out,] img, f::LinearStretching)
```

Returns an image where the range of the intensities spans the interval [`minval`, `maxval`].

# Details

Linear stretching (also called *normalization*) is a contrast enhancing
transformation that is used to modify the dynamic range of the image. In
particular, suppose that the input image has gray values in the range [A,B] and
one wishes to change the dynamic range to [a,b] using a linear mapping, then the
necessary transformation is given by the relation

```math
f(x) = (x-A) \\frac{b-a}{B-A} + a.
```

# Options

Various options for the parameters of the `adjust_histogram` and
`LinearStretching` type  are described in more detail below.

## Choices for `img`

The function can handle a variety of input types. The returned
image depends on the input type.

For colored images, the input is converted to the
[YIQ](https://en.wikipedia.org/wiki/YIQ)  type and the intensities of the Y
channel are stretched to the specified range. The modified Y channel is then
combined with the I and Q channels and the resulting image converted to the same
type as the input.

## Choices for `minval` and `maxval`

If minval and maxval are specified then intensities are mapped to the range
[`minval`, `maxval`]. The default values are 0 and 1.

# Example

```julia
using ImageContrastAdjustment, ImageView, TestImages

img = testimage("mandril_gray")
imgo = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))

```

# References
1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)

"""
Base.@kwdef struct LinearStretching{T₁ <: Union{Real,AbstractGray},
                                    T₂ <: Union{Real,AbstractGray}} <: AbstractHistogramAdjustmentAlgorithm
    minval::T₁ = 0.0
    maxval::T₂ = 1.0
end

function (f::LinearStretching)(out::GenericGrayImage, img::GenericGrayImage)
    src_minval = minfinite(img)
    src_maxval = maxfinite(img)
    T = eltype(out)
    out .= img
    map!(out,out) do val
        if isnan(val)
            return val
        else
            newval = linear_stretch(val, src_minval, src_maxval, f.minval, f.maxval)
            return  T <: Integer ? round(Int, newval ) : newval
        end
    end
end

function (f::LinearStretching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(out)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)
    adjust_histogram!(view(yiq_view,1,:,:), f)
    out .= convert.(T, yiq)
end

(f::LinearStretching)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
using Test
using ReferenceFrameRotations
using LinearAlgebra
using StaticArrays

# Available rotations.
const rot_seq_array = [:XYX,
                       :XYZ,
                       :XZX,
                       :XZY,
                       :YXY,
                       :YXZ,
                       :YZX,
                       :YZY,
                       :ZXY,
                       :ZXZ,
                       :ZYX,
                       :ZYZ]

# Number of samples.
const samples = 1000

@time @testset "Direction Cosine Matrices" begin
    include("./test_dcm.jl")
end
println("")

@time @testset "Quaternions" begin
    include("./test_quaternions.jl")
end
println("")

@time @testset "DCM <=> Euler Angles conversions" begin
    include("./test_conversion_dcm_euler_angles.jl")
end
println("")

@time @testset "DCM <=> Quaternions conversions" begin
    include("./test_conversion_dcm_quaternions.jl")
end
println("")

@time @testset "Euler Angles <=> Quaternion conversion" begin
    include("./test_conversion_euler_angles_quaternions.jl")
end
println("")

@time @testset "Euler Angle Axis <=> Quaternion conversions" begin
    include("./test_conversion_euler_angle_axis_quaternions.jl")
end
println("")

@time @testset "Kinematics using DCMs" begin
    include("./test_kinematics_dcm.jl")
end
println("")

@time @testset "Kinematics using Quaternions" begin
    include("./test_kinematics_quaternions.jl")
end
println("")

@time @testset "Compose rotations using DCMs" begin
    include("./test_compose_rotations_dcm.jl")
end
println("")

@time @testset "Compose rotations using Quaternions" begin
    include("./test_compose_rotations_quaternions.jl")
end
println("")


text:
module CWorldMCTSVisualize

using Plots, ContinuumWorld, POMDPToolbox, MCTS

export CWorldMCTSVis

struct CWorldMCTSVis
    mdp::CWorld
    s::Vec2
    actions::Vector{Vector{Vec2}}
    next_states::Vector{Vector{Vec2}}
    ns::Vector{Int}
    radii::Vector{Vector{Float64}}
end
CWorldMCTSVis(mdp::CWorld, s::Vec2=Vec2(0,0)) = CWorldMCTSVis(mdp, s, Vec2[], Vec2[], Int[], Float64[])

function circle(x,y,r)
    c = Shape(Plots.partialcircle(0,2π,20,r))
    translate!(c, x, y)
    c
end

@recipe function f(c::CircularRegion)
    @series begin
        plot(circle(c.center[1], c.center[2], c.radius))
    end
end

@recipe function f(vis::CWorldMCTSVis, h::SimHistory)
    mdp = vis.mdp
    xlim := mdp.xlim
    ylim := mdp.ylim
    aspect_ratio := :equal
    @series begin
        label := "reward region"
        [circle(c.center[1], c.center[2], c.radius) for c in mdp.reward_regions]
    end
    @series begin
        label := "path"
        color := :blue
        x = [s[1] for s in state_hist(h)[1:end-1]]
        y = [s[2] for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        label := "current position"
        marker := :circle
        seriestype := :scatter
        color := :red
        x,y = state_hist(h)[end-1]
        [x], [y]
    end
end

@recipe function f(vis::CWorldMCTSVis, h::SimHistory, index::Int)
    mdp = vis.mdp
    xlim := mdp.xlim
    ylim := mdp.ylim
    aspect_ratio := :equal
    @series begin
        label := "actions"
        alpha := 0.25
        [circle(x,y,r) for ((x,y),r) in zip(vis.next_states[index], vis.radii[index])]
    end
    @series begin
        label := "reward region"
        [circle(c.center[1], c.center[2], c.radius) for c in mdp.reward_regions]
    end
    @series begin
        label := "path"
        color := :blue
        x = [s[1] for s in state_hist(h)[1:end-1]]
        y = [s[2] for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        label := "current position"
        marker := :circle
        seriestype := :scatter
        color := :red
        x,y = state_hist(h)[end-1]
        [x], [y]
    end
end

function MCTS.notify_listener(vis::CWorldMCTSVis,dsb::DSBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(dsb.tree)
        sol = dsb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            push!(vis.next_states, [s+action for action in actions])
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [sol.r0_action/tree.total_n[snode]^sol.lambda_action for action in actions])
        end
    end
end
function MCTS.notify_listener(vis::CWorldMCTSVis,asb::ASBPlanner,s,a,sp,r,snode,sanode,spnode)
    if s == vis.s
        tree = get(asb.tree)
        sol = asb.solver
        actions = [tree.a_labels[c] for c in tree.children[snode]]
        if isempty(vis.actions) || (actions != vis.actions[end])  #only push the deltas
            push!(vis.actions, actions)
            push!(vis.next_states, [s+action for action in actions])
            push!(vis.ns, tree.total_n[snode])
            push!(vis.radii, [tree.a_radius[x] for x in tree.children[snode]])
        end
    end
end

end # module


text:
using Counters

# stdlib
using Test
using Dates

# internal dependencies
using SystemSnoop

# external dependenceis
using Tables

include("pcm.jl")
include("database.jl")

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
# generate a randomized Beehive

export generateCanonicalFG_Beehive!

"""
    $SIGNATURES

Pretend a bee is walking in a hive where each step (pose) follows one edge of an imaginary honeycomb lattice,
and at after each step a new direction left or right is stochastically chosen and the process repeats.

Notes
- The keyword `locality=1` is a positive `::Real` ∈ [0,∞) value, where higher numbers imply direction decisions are more sticky for multiple steps.
- Use keyword callback function `postpose_cb = (fg, lastpose) -> ...` to hook in your own features right after each new pose step.

DevNotes
- TODO rewrite as a recursive generator function instead.

See also: [`generateCanonicalFG_Honeycomb!`](@ref), [`generateCanonicalFG_Hexagonal`](@ref), [`generateCanonicalFG_ZeroPose`](@ref)
"""
function generateCanonicalFG_Beehive!(poseCountTarget::Int=10;
                                      graphinit::Bool = true,
                                      dfg::AbstractDFG = LightDFG{SolverParams}(solverParams=SolverParams(graphinit=graphinit)),
                                      useMsgLikelihoods::Bool=getSolverParams(dfg).useMsgLikelihoods,
                                      solvable::Int = 1,
                                      refKey::Symbol = :simulated,
                                      addLandmarks::Bool = true,
                                      landmarkSolvable::Int=0,
                                      poseRegex::Regex=r"x\d+",
                                      pose0::Symbol=Symbol(match(r"[a-zA-Z_]+", poseRegex.pattern).match, 0),
                                      yaw0::Real = ([0.0;-2pi/3;2pi/3])[rand(1:3)],
                                      μ0::AbstractVector{<:Real} = [0;0;yaw0],
                                      postpose_cb::Function=(fg_,latestpose)->(),
                                      locality::Real=1,
                                      atol::Real=1     )
  #

  # does anything exist in the graph yet
  posecount = if pose0 in ls(dfg)
    # what is the last pose
    lastPose = (ls(dfg, poseRegex) |> sortDFG)[end]
    # get latest posecount number
    match(r"\d+", string(lastPose)).match |> x->parse(Int,x)
  else
    # pick a random staring direction
    # initial zero pose
    generateCanonicalFG_ZeroPose(dfg=dfg, varType=Pose2, μ0=μ0, graphinit=graphinit, postpose_cb=postpose_cb)
    getSolverParams(dfg).useMsgLikelihoods = useMsgLikelihoods

    # add a new landmark (if not yet present)
    !addLandmarks ? nothing : _addLandmarkBeehive!(dfg, pose0, refKey=refKey, solvable=landmarkSolvable, atol=atol, graphinit=false, _doHoneycomb=false)

    # staring posecount (i.e. :x0)
    0
  end

  # keep adding poses until the target number is reached
  direction = rand(1:2) === 1 ? :left : :right
  _locality = 1/(1.0+locality)
  while posecount < poseCountTarget
    #change or keep direction according to locality keyword
    direction = rand() < _locality ? (direction == :left ? :right : :left) : direction
    posecount = _offsetHexLeg(dfg, posecount, direction=direction, graphinit=graphinit,
                              landmarkSolvable=landmarkSolvable, poseCountTarget=poseCountTarget,
                              atol=atol, postpose_cb=postpose_cb, _doHoneycomb=false)
  end

  # NOTE solvable forced for everything at this time
  setSolvable!.(dfg, ls(dfg),  solvable)
  setSolvable!.(dfg, lsf(dfg), solvable)

  return dfg
end

#


text:
#=
    Copyright (C) 2008, Martin Albrecht
    Copyright (C) 2008, 2009 William Hart.
    Copyright (C) 2010, Fredrik Johansson

    This file is part of FLINT.

    FLINT is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 2.1 of the License, or
    (at your option) any later version.  See <http://www.gnu.org/licenses/>.

    to julia:
    Copyright (C) 2020, Claus Fieker
=#

module Strassen
using Hecke
import AbstractAlgebra, Nemo

const cutoff = 64

function Nemo.mul!(C::AbstractArray, A::AbstractArray, B::AbstractArray, add::Bool = false)
  @assert size(C) == (2,2) && size(A) == (2,2) && size(B) == (2,2)
  C[1,1] = A[1,1] * B[1,1] + A[1,2] * B[2,1]
  C[1,2] = A[1,1] * B[1,2] + A[1,2] * B[2,2]
  C[2,1] = A[2,1] * B[1,1] + A[2,2] * B[2,1]
  C[2,2] = A[2,1] * B[1,2] + A[2,2] * B[2,2]
  return C
  if add
    C .+= A*B
  else
    C .= A*B
  end
end

function mul_strassen!(C::AbstractArray, A::AbstractArray, B::AbstractArray)
  sA = size(A)
  sB = size(B)
  sC = size(C)
  a = sA[1]
  b = sA[2]
  c = sB[2]

  @assert a == sC[1] && b == sB[1] && c == sC[2]

  if (a <= cutoff || b <= cutoff || c <= cutoff)
      mul!(C, A, B)
      return
  end

  anr = div(a, 2)
  anc = div(b, 2)
  bnr = anc
  bnc = div(c, 2)

  #nmod_mat_window_init(A11, A, 0, 0, anr, anc);
  #nmod_mat_window_init(A12, A, 0, anc, anr, 2*anc);
  #nmod_mat_window_init(A21, A, anr, 0, 2*anr, anc);
  #nmod_mat_window_init(A22, A, anr, anc, 2*anr, 2*anc);
  A11 = view(A, 1:anr, 1:anc)
  A12 = view(A, 1:anr, anc+1:2*anc)
  A21 = view(A, anr+1:2*anr, 1:anc)
  A22 = view(A, anr+1:2*anr, anc+1:2*anc)

  #nmod_mat_window_init(B11, B, 0, 0, bnr, bnc);
  #nmod_mat_window_init(B12, B, 0, bnc, bnr, 2*bnc);
  #nmod_mat_window_init(B21, B, bnr, 0, 2*bnr, bnc);
  #nmod_mat_window_init(B22, B, bnr, bnc, 2*bnr, 2*bnc);
  B11 = view(B, 1:bnr, 1:bnc)
  B12 = view(B, 1:bnr, bnc+1:2*bnc)
  B21 = view(B, bnr+1:2*bnr, 1:bnc)
  B22 = view(B, bnr+1:2*bnr, bnc+1:2*bnc)

  #nmod_mat_window_init(C11, C, 0, 0, anr, bnc);
  #nmod_mat_window_init(C12, C, 0, bnc, anr, 2*bnc);
  #nmod_mat_window_init(C21, C, anr, 0, 2*anr, bnc);
  #nmod_mat_window_init(C22, C, anr, bnc, 2*anr, 2*bnc);
  C11 = view(C, 1:anr, 1:bnc)
  C12 = view(C, 1:anr, bnc+1:2*bnc)
  C21 = view(C, anr+1:2*anr, 1:bnc)
  C22 = view(C, anr+1:2*anr, bnc+1:2*bnc)

  #nmod_mat_init(X1, anr, FLINT_MAX(bnc, anc), A->mod.n);
  #nmod_mat_init(X2, anc, bnc, A->mod.n);

  #X1->c = anc;

  #=
      See Jean-Guillaume Dumas, Clement Pernet, Wei Zhou; "Memory
      efficient scheduling of Strassen-Winograd's matrix multiplication
      algorithm"; http://arxiv.org/pdf/0707.2347v3 for reference on the
      used operation scheduling.
  =#

  #nmod_mat_sub(X1, A11, A21);
  X1 = A11 .- A21
  #nmod_mat_sub(X2, B22, B12);
  X2 = B22 .- B12
  #nmod_mat_mul(C21, X1, X2);
  mul_strassen!(C21, X1, X2)

  #nmod_mat_add(X1, A21, A22);
  X1 .= A21 .+ A22
  #nmod_mat_sub(X2, B12, B11);
  X2 .= B12 .- B11
  #nmod_mat_mul(C22, X1, X2);
  mul_strassen!(C22, X1, X2)

  #nmod_mat_sub(X1, X1, A11);
  X1 .-= A11
  #nmod_mat_sub(X2, B22, X2);
  X2 .= B22 .- X2
  #nmod_mat_mul(C12, X1, X2);
  mul_strassen!(C12, X1, X2)

  #nmod_mat_sub(X1, A12, X1);
  X1 .= A12 .- X1
  #nmod_mat_mul(C11, X1, B22);
  mul_strassen!(C11, X1, B22)

  #X1->c = bnc;
  #nmod_mat_mul(X1, A11, B11);
  X1 = similar(C11)
  mul_strassen!(X1, A11, B11)

  #nmod_mat_add(C12, X1, C12);
  C12 .+= X1
  #nmod_mat_add(C21, C12, C21);
  C21 .+= C12
  #nmod_mat_add(C12, C12, C22);
  C12 .+= C22
  #nmod_mat_add(C22, C21, C22);
  C22 .+= C21
  #nmod_mat_add(C12, C12, C11);
  C12 .+= C11
  #nmod_mat_sub(X2, X2, B21);
  X2 .-= B21
  #nmod_mat_mul(C11, A22, X2);
  mul_strassen!(C11, A22, X2)

  #nmod_mat_sub(C21, C21, C11);
  C21 .-= C11

  #nmod_mat_mul(C11, A12, B21);
  mul_strassen!(C11, A12, B21)

  #nmod_mat_add(C11, X1, C11);
  C11 .+= X1

  if c > 2*bnc #A by last col of B -> last col of C
    error("missing")
    #=
  {
      nmod_mat_t Bc, Cc;
      nmod_mat_window_init(Bc, B, 0, 2*bnc, b, c);
      nmod_mat_window_init(Cc, C, 0, 2*bnc, a, c);
      nmod_mat_mul(Cc, A, Bc);
      nmod_mat_window_clear(Bc);
      nmod_mat_window_clear(Cc);
  }
    =#
  end

  if a > 2*anr #last row of A by B -> last row of C
    error("missing")
    #=
  {
      nmod_mat_t Ar, Cr;
      nmod_mat_window_init(Ar, A, 2*anr, 0, a, b);
      nmod_mat_window_init(Cr, C, 2*anr, 0, a, c);
      nmod_mat_mul(Cr, Ar, B);
      nmod_mat_window_clear(Ar);
      nmod_mat_window_clear(Cr);
  }
    =#
  end

  if b > 2*anc # last col of A by last row of B -> C
    error("missing")
    #=
  {
      nmod_mat_t Ac, Br, Cb;
      nmod_mat_window_init(Ac, A, 0, 2*anc, 2*anr, b);
      nmod_mat_window_init(Br, B, 2*bnr, 0, b, 2*bnc);
      nmod_mat_window_init(Cb, C, 0, 0, 2*anr, 2*bnc);
      nmod_mat_addmul(Cb, Cb, Ac, Br);
      nmod_mat_window_clear(Ac);
      nmod_mat_window_clear(Br);
      nmod_mat_window_clear(Cb);
  }
    =#
  end
end

function mul_strassen!(A::Generic.Mat, B::Generic.Mat, C::Generic.Mat)
  mul_strassen!(C.entries, A.entries, B.entries)
end

function mul_strassen!(C::MatElem{T}, A::MatElem{T}, B::MatElem{T}) where {T}
  sA = size(A)
  sB = size(B)
  sC = size(C)
  a = sA[1]
  b = sA[2]
  c = sB[2]

  @assert a == sC[1] && b == sB[1] && c == sC[2]

  if (a <= cutoff || b <= cutoff || c <= cutoff)
      mul!(C, A, B)
      return
  end

  anr = div(a, 2)
  anc = div(b, 2)
  bnr = anc
  bnc = div(c, 2)

  #nmod_mat_window_init(A11, A, 0, 0, anr, anc);
  #nmod_mat_window_init(A12, A, 0, anc, anr, 2*anc);
  #nmod_mat_window_init(A21, A, anr, 0, 2*anr, anc);
  #nmod_mat_window_init(A22, A, anr, anc, 2*anr, 2*anc);
  A11 = view(A, 1:anr, 1:anc)
  A12 = view(A, 1:anr, anc+1:2*anc)
  A21 = view(A, anr+1:2*anr, 1:anc)
  A22 = view(A, anr+1:2*anr, anc+1:2*anc)

  #nmod_mat_window_init(B11, B, 0, 0, bnr, bnc);
  #nmod_mat_window_init(B12, B, 0, bnc, bnr, 2*bnc);
  #nmod_mat_window_init(B21, B, bnr, 0, 2*bnr, bnc);
  #nmod_mat_window_init(B22, B, bnr, bnc, 2*bnr, 2*bnc);
  B11 = view(B, 1:bnr, 1:bnc)
  B12 = view(B, 1:bnr, bnc+1:2*bnc)
  B21 = view(B, bnr+1:2*bnr, 1:bnc)
  B22 = view(B, bnr+1:2*bnr, bnc+1:2*bnc)

  #nmod_mat_window_init(C11, C, 0, 0, anr, bnc);
  #nmod_mat_window_init(C12, C, 0, bnc, anr, 2*bnc);
  #nmod_mat_window_init(C21, C, anr, 0, 2*anr, bnc);
  #nmod_mat_window_init(C22, C, anr, bnc, 2*anr, 2*bnc);
  C11 = view(C, 1:anr, 1:bnc)
  C12 = view(C, 1:anr, bnc+1:2*bnc)
  C21 = view(C, anr+1:2*anr, 1:bnc)
  C22 = view(C, anr+1:2*anr, bnc+1:2*bnc)

  #nmod_mat_init(X1, anr, FLINT_MAX(bnc, anc), A->mod.n);
  #nmod_mat_init(X2, anc, bnc, A->mod.n);

  #X1->c = anc;

  #=
      See Jean-Guillaume Dumas, Clement Pernet, Wei Zhou; "Memory
      efficient scheduling of Strassen-Winograd's matrix multiplication
      algorithm"; http://arxiv.org/pdf/0707.2347v3 for reference on the
      used operation scheduling.
  =#

  X1 = A11 - A21
  X2 = B22 - B12
  #nmod_mat_mul(C21, X1, X2);
  mul_strassen!(C21, X1, X2)

  add!(X1, A21, A22);
  sub!(X2, B12, B11);
  #nmod_mat_mul(C22, X1, X2);
  mul_strassen!(C22, X1, X2)

  sub!(X1, X1, A11);
  sub!(X2, B22, X2);
  #nmod_mat_mul(C12, X1, X2);
  mul_strassen!(C12, X1, X2)

  sub!(X1, A12, X1);
  #nmod_mat_mul(C11, X1, B22);
  mul_strassen!(C11, X1, B22)

  #X1->c = bnc;
  #nmod_mat_mul(X1, A11, B11);
  X1 = similar(C11)
  mul_strassen!(X1, A11, B11)

  add!(C12, X1, C12);
  add!(C21, C12, C21);
  add!(C12, C12, C22);
  add!(C22, C21, C22);
  add!(C12, C12, C11);
  sub!(X2, X2, B21);
  #nmod_mat_mul(C11, A22, X2);
  mul_strassen!(C11, A22, X2)

  sub!(C21, C21, C11);

  #nmod_mat_mul(C11, A12, B21);
  mul_strassen!(C11, A12, B21)

  add!(C11, X1, C11);

  if c > 2*bnc #A by last col of B -> last col of C
      #nmod_mat_window_init(Bc, B, 0, 2*bnc, b, c);
      Bc = view(B, 1:b, 2*bnc+1:c)
      #nmod_mat_window_init(Cc, C, 0, 2*bnc, a, c);
      Cc = view(C, 1:a, 2*bnc+1:c)
      #nmod_mat_mul(Cc, A, Bc);
      mul!(Cc, A, Bc)
  end

  if a > 2*anr #last row of A by B -> last row of C
      #nmod_mat_window_init(Ar, A, 2*anr, 0, a, b);
      Ar = view(A, 2*anr+1:a, 1:b)
      #nmod_mat_window_init(Cr, C, 2*anr, 0, a, c);
      Cr = view(C, 2*anr+1:a, 1:c)
      #nmod_mat_mul(Cr, Ar, B);
      mul!(Cr, Ar, B)
  end

  if b > 2*anc # last col of A by last row of B -> C
      #nmod_mat_window_init(Ac, A, 0, 2*anc, 2*anr, b);
      Ac = view(A, 1:2*anr, 2*anc:b)
      #nmod_mat_window_init(Br, B, 2*bnr, 0, b, 2*bnc);
      Br = view(B, 2*bnr+1:b, 1:2*bnc)
      #nmod_mat_window_init(Cb, C, 0, 0, 2*anr, 2*bnc);
      Cb = view(C, 1:2*anr, 1:2*bnc)
      #nmod_mat_addmul(Cb, Cb, Ac, Br);
      mul!(Cb, Ac, Br, true)
  end
end

end # module