text:
module Test2D

using Test
using Trixi2Vtk

include("test_trixi2vtk.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Create artifacts directory where all files that should be preserved will be stored
artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts")
if !isdir(artifacts_dir)
  mkdir(artifacts_dir)
end

@testset "2D" begin
  @testset "TreeMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("tree_2d_dgsem", "elixir_advection_extended.jl"), maxiters=1)

    @timed_testset "uniform mesh" begin
      if Sys.iswindows()
        # This test fails on Windows due to globbing not working
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"),
                    ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")])
        outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu")

      else
        test_trixi2vtk("solution_00000*.h5", outdir,
            hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"),
                    ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830"),
                    ("solution_00000.pvd", "b9e6742dc2b397b14d8d3964e90dcfadea5d98cb"),
                    ("solution_00000_celldata.pvd", "c12004714a980581450cd4bad16a2541c5ec8f26")])
        outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu",
                    "solution_00000.pvd", "solution_00000_celldata.pvd")
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      testname = "2d-tree-mesh-uniform-mesh"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end

    @timed_testset "uniform mesh with vti output" begin
      if Sys.isapple()
        # This test fails on MacOS due to differing binary VTI files (even though the contents match)
        test_trixi2vtk("restart_000001.h5", outdir,
            format=:vti)
      else
        test_trixi2vtk("restart_000001.h5", outdir,
            hashes=[("restart_000001.vti", "9ade067c71f1f6492242a8aa215bd0d633caf9bc"),
                    ("restart_000001_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")],
            format=:vti)
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("restart_000001.vti", "restart_000001_celldata.vtu")
      testname = "2d-tree-mesh-uniform-mesh-with-vti-output"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "StructuredMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("structured_2d_dgsem", "elixir_advection_waving_flag.jl"), maxiters=1)

    @timed_testset "waving flag" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "a93dbd393647627a861d890568e65598be0062f9")])
      else
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "564701ed0a9a90230f3a67f8bddd0616c818319b")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "2d-curved-mesh-waving-flag"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end

    @timed_testset "waving flag (supersampling)" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir, nvisnodes=6,
            hashes=[("solution_000000.vtu", "c6d74ab831bf4b6de2ba8cf537b6653ad611cfe7")])
      else
        test_trixi2vtk("solution_000000.h5", outdir, nvisnodes=6,
            hashes=[("solution_000000.vtu", "ae8c059c110aaabe2ed7dcfa8516d336c15ba618")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "2d-curved-mesh-waving-flag-supersampling"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "UnstructuredMesh2D" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("unstructured_2d_dgsem", "elixir_euler_basic.jl"), maxiters=1)

    @timed_testset "basic" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "0daedeea99d03d53b925ce5691bd7924abe88861")])
      else
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "acc03f295d2b7daee2cb0b4e29b90035c5e92bd7")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "2d-unstructured-quad-basic"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "P4estMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("p4est_2d_dgsem", "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), initial_refinement_level=0, maxiters=1)

    @timed_testset "nonperiodic" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "476eaaa9c13c3a8dba826b614a68a9d3e7979e6b")])
      else
        test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "a80aadb353ce6ec40baa1b94d278c480f17d0419")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "2d-p4est-mesh-nonperiodic"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end


text:
rawproto(io::IO) = readproto(io, Proto.ModelProto())
rawproto(path::String) = open(rawproto, path)

"""
Retrieve only the useful information from a AttributeProto
object into a Dict format.
"""
function convert_model(x::Proto.AttributeProto)
    if (x._type != 0)
        field = [:f, :i, :s, :t, :g, :floats, :ints, :strings, :tensors, :graphs][x._type]
        return Symbol(x.name) => getfield(x, field)
    end
end

convert_array(as) = Dict(convert_model(a) for a in as)

"""
Convert a ValueInfoProto to  ValueInfo.
"""
function convert_model(model::Proto.ValueInfoProto)
    a = Types.ValueInfo(model.name, model.doc_string)
    return a
end

"""
Convert an OperatorSetIdProto to Dict.
"""
function convert_model(model::KnetOnnx.Proto.OperatorSetIdProto)
    a = Dict{Symbol, Any}()
    fields = [:domain, :version]
    for ele in fields
        a[ele] = getfield(model, ele)
    end
    return a
end

"""
Convert a StringStringEntryProto to Dict.
"""
function convert_model(model::KnetOnnx.Proto.StringStringEntryProto)
    a = Dict{Symbol, Any}()
    fields = [:key, :value]
    for ele in fields
        a[ele] = getfield(model, ele)
    end
    return a
end

"""
Get the array from a TensorProto object.
"""
function get_array(x::Proto.TensorProto)
    if (x.data_type == 1)
        if !isempty(x.float_data)
            x = reshape(reinterpret(Float32, x.float_data), reverse(x.dims)...)
        else
            x = reshape(reinterpret(Float32, x.raw_data), reverse(x.dims)...)
        end
        return x
    end
    if x.data_type == 7
        if !isempty(x.raw_data)
            x = reshape(reinterpret(Int64, x.raw_data), reverse(x.dims)...)
        else
            x = reshape(reinterpret(Int64, x.int64_data), reverse(x.dims)...)
        end
        return x
    end
    if x.data_type == 9
        x = reshape(reinterpret(Int8, x.raw_data), reverse(x.dims)...)
        return x
    end
    if x.data_type == 6
         x = reshape(reinterpret(Int32, x.raw_data), reverse(x.dims)...)
        return x
    end
    if x.data_type == 11
        if !isempty(x.raw_data)
            x = reshape(reinterpret(Float64, x.raw_data), reverse(x.dims)...)
        else
            x = Base.convert(Array{Float32, N} where N, reshape(x.double_data , reverse(x.dims)...))
        end
        return x
    end
    if x.data_type == 10
        x = reshape(reinterpret(Float16, x.raw_data), reverse(x.dims)...)
        return x
    end
end

"""
Convert a ModelProto object to a Model type.
"""
function convert(model::Proto.ModelProto)
    # conversion for opset_import
    arr1 = Array{Any, 1}()
    for ele in model.opset_import
        push!(arr1, convert_model(ele))
    end

    # conversion for stringstringentry proto
    arr2 = Array{Any, 1}()
    for ele in model.metadata_props
        push!(arr2, convert_model(ele))
    end

    m = Types.Model(model.ir_version,
                arr1, model.producer_name,
                model.producer_version,
                model.domain, model.model_version,
                model.doc_string, convert(model.graph),
                arr2)
    return m
end

"""
Convert a GraphProto object to Graph type.
"""
function convert(model::Proto.GraphProto)
    # conversion for vector of nodeproto
    arr1 = Array{Any, 1}()
    for ele in model.node
        push!(arr1, convert(ele))
    end

    # conversion for vector of tensorproto
    arr2 = Dict{Any, Any}()
    for ele in model.initializer
        arr2[ele.name] = get_array(ele)
    end

    #conversion for vector of valueinfoproto
    arr3 = Array{Types.ValueInfo ,1}()
    for ele in model.input
        push!(arr3, convert_model(ele))
    end

    arr4 = Array{Types.ValueInfo ,1}()
    for ele in model.output
        push!(arr4, convert_model(ele))
    end

    arr5 = Array{Types.ValueInfo ,1}()
    for ele in model.value_info
        push!(arr5, convert_model(ele))
    end

    m = Types.Graph(arr1,
            model.name,
            arr2, model.doc_string,
            arr3, arr4, arr5)
    return m
end

"""
Convert a Proto.NodeProto to Node type.
"""
function convert(model::Proto.NodeProto)
    # Conversion of attribute
    arr1 = convert_array(model.attribute)

    m = Types.Node(model.input,
            model.output,
            model.name,
            model.op_type,
            model.domain,
            arr1,
            model.doc_string)
    return m
end

function parent(path)
    temp = split(path, "/")
    res = ""
    for element in temp
        if (element != temp[end])
            res = res * element * "/"
        end
    end
    return res
end

ext: jl


text:
using SampleChainsDynamicHMC
using TransformVariables
using Test

function ℓ(nt)
	z = nt.x/nt.σ
	return -z^2 - nt.σ - log(nt.σ)
end

t = as((x=asℝ, σ=asℝ₊))

@testset "single chain" begin
	chain = newchain(dynamichmc(), ℓ, t)
	@test chain isa SampleChainsDynamicHMC.DynamicHMCChain
	@test length(chain) == 1

	sample!(chain, 9)
	@test length(chain) == 10

	sample!(chain, 90)
	@test length(chain) == 100
end

@testset "multichain" begin
	chains = newchain(4, dynamichmc(), ℓ, t)
	@test chains isa MultiChain
	chains_chains = getfield(chains, :chains)
	@test length(chains_chains) == 4
	@test all(x -> length(x) == 1, chains_chains)

	sample!(chains, 9)
	@test all(x -> length(x) == 10, chains_chains)
	sample!(chains, 90)
	@test all(x -> length(x) == 100, chains_chains)

	samples(chains)
end

using LinearAlgebra
using ReverseDiff
using SampleChainsDynamicHMC.LogDensityProblems

@testset "config options" begin
	config = dynamichmc(
		warmup_stages=default_warmup_stages(
			M=Symmetric, # adapt dense positive definite metric
			stepsize_adaptation=DualAveraging(δ=0.9), # target acceptance rate 0.9
			doubling_stages=7, # 7-stage metric adaptation
		),
		reporter=LogProgressReport(), # log progress using Logging
		ad_backend=Val(:ReverseDiff), # use ReverseDiff AD package
	)

	chain = newchain(config, ℓ, t)
	@test length(chain) == 1

	meta = getfield(chain, :meta)
	@test meta.H.κ.M⁻¹ isa Symmetric
	@test meta.H.ℓ isa LogDensityProblems.ReverseDiffLogDensity

	sample!(chain, 9)
	@test length(chain) == 10
	sample!(chain, 90)
	@test length(chain) == 100
end


text:
#!/bin/bash
#=
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#

using Pkg
#Pkg.activate(".")
#include("srg/MockFBA.jl")
import MockFBA

using YAML

function main(ARGS)
    #load the config file
    config=YAML.load_file(ARGS[1])
    tracer=ARGS[2]
    group =ARGS[3]

    #how many partition we want this pass to run with
    npart_zone=parse(Int32,ARGS[4])
    #which partition should I execute now
    mypart_zone=parse(Int32,ARGS[5])

    #process the zones
    MockFBA.PostFBA_JLD2FITS(config,tracer,group,mypart_zone,npart_zone,;columns=["RA","DEC","FITS_index_julia"])
end

@show ARGS
main(ARGS)

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
using Test
using PolynomialBasis
using Roots
using IntervalArithmetic
# using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function allapprox(v1, v2; tol = 1e-14)
    @assert length(v1) == length(v2)
    flags = [isapprox(v1[i], v2[i], atol = tol) for i = 1:length(v1)]
    return all(flags)
end

f(x) = (x - 1.0) * (x - 2.0) * (x - 3.0)
g(x) = (x - 0.1) * (x - 0.2) * (x - 2.5)

quad1d = IDQ.ReferenceQuadratureRule(3)
f2(x) = f(x[1])
g2(x) = g(x[1])
w0 = 2.0
@test_throws AssertionError IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1],
    1,
    0.0,
    3.5,
    [1.0],
    w0,
    quad1d,
)
quad1 = IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1],
    1,
    0.0,
    3.5,
    [1.0],
    w0,
    quad1d,
)
quad2 = IDQ.one_dimensional_quadrature([f, g], [+1, -1], 0.0, 3.5, quad1d)
@test allapprox(quad1.points, IDQ.extend([1.0], 1, quad2.points))
@test allapprox(quad1.weights, w0 * quad2.weights)

edgequad = IDQ.TemporaryQuadrature(reshape([1.0], 1, 1), [w0])
quad1 = IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1],
    1,
    Interval(0.0, 3.5),
    edgequad,
    quad1d,
)
quad2 = IDQ.one_dimensional_quadrature([f, g], [+1, -1], Interval(0.0, 3.5), quad1d)
@test allapprox(quad1.points, IDQ.extend([1.0], 1, quad2.points))
@test allapprox(quad1.weights, w0 * quad2.weights)

x0 = [1.0 2.0 3.0]
w0 = [3.0, 6.0, 9.0, 10.0]
@test_throws AssertionError IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1, +1],
    1,
    0.0,
    3.5,
    x0,
    w0,
    quad1d,
)
@test_throws AssertionError IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1],
    1,
    0.0,
    3.5,
    x0,
    w0,
    quad1d,
)

w0 = [3.0, 6.0, 9.0]
quad_ext = IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1],
    1,
    0.0,
    3.5,
    x0,
    w0,
    quad1d,
)
quad = IDQ.one_dimensional_quadrature([f, g], [+1, -1], 0.0, 3.5, quad1d)
testp = hcat([IDQ.extend([x0[i]], 1, quad.points) for i = 1:3]...)
testw = vcat([quad.weights * w0[i] for i = 1:3]...)
@test allapprox(quad_ext.points, testp)
@test allapprox(quad_ext.weights, testw)

edgequad = IDQ.TemporaryQuadrature(x0, w0)
quad_ext = IDQ.extend_edge_quadrature_to_area_quadrature(
    [f2, g2],
    [+1, -1],
    1,
    Interval(0.0, 3.5),
    edgequad,
    quad1d,
)
@test allapprox(quad_ext.points, testp)
@test allapprox(quad_ext.weights, testw)

f2(x) = x[2]
numqp = 5
quad1d = IDQ.ReferenceQuadratureRule(numqp)
xL, xR = [-1.0, -1.0], [1.0, 1.0]
P = InterpolatingPolynomial(1, 2, 2)
interpgrad = InterpolatingPolynomial(2,2,2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
update_interpolating_gradient!(interpgrad,P)
quad = IDQ.area_quadrature(P, interpgrad, +1, xL, xR, numqp)
p, w = IDQ.transform(quad1d, 0.0, 1.0)
p2 = hcat([IDQ.extend([quad1d.points[i]], 2, p) for i = 1:5]...)
w2 = vcat([quad1d.weights[i] * w for i = 1:5]...)
@test allapprox(quad.points, p2)
@test allapprox(quad.weights, w2)

f2(x) = x[2] + 1.5
numqp = 5
quad1d = IDQ.ReferenceQuadratureRule(numqp)
xL, xR = [-1.0, -1.0], [1.0, 1.0]
P = InterpolatingPolynomial(1, 2, 2)
dP = InterpolatingPolynomial(2, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
update_interpolating_gradient!(dP,P)
quad = IDQ.area_quadrature(P,dP, -1, xL, xR, numqp)
@test size(quad.points) == (2, 0)
@test length(quad.weights) == 0

f2(x) = x[2] + 1.5
numqp = 5
quad1d = IDQ.ReferenceQuadratureRule(numqp)
xL, xR = [-1.0, -1.0], [1.0, 1.0]
P = InterpolatingPolynomial(1, 2, 2)
dP = InterpolatingPolynomial(2, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
update_interpolating_gradient!(dP,P)
quad = IDQ.area_quadrature(P,dP, +1, xL, xR, numqp)
p = hcat([IDQ.extend([quad1d.points[i]], 2, quad1d.points) for i = 1:5]...)
w = vcat([quad1d.weights[i] * quad1d.weights for i = 1:5]...)
@test allapprox(quad.points, p)
@test allapprox(quad.weights, w)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
# p5.jl - repetition of p4.jl via FFT
#        For complex v, delete "real" commands.

# Differentiation of a hat function:
N = 24; h = 2*pi/N; x = h*(1:N);
v = @. max(0,1-abs(x-pi)/2); v_hat = fft(v);
w_hat = 1im*[0:N/2-1;0;-N/2+1:-1] .* v_hat;
w = real(ifft(w_hat)); clf();
subplot(2,2,1);  plot(x,v,".-",markersize=6);
axis([0,2pi,-.5,1.5]); grid(true); title("function");
subplot(2,2,2), plot(x,D*v,".-",markersize=6);
axis([0,2pi,-1,1]); grid(true); title("spectral derivative");

# Differentiation of exp(sin(x)):
v = @. exp(sin(x)); vprime = @. cos(x)*v;
v_hat = fft(v);
w_hat = 1im*[0:N/2-1;0;-N/2+1:-1] .* v_hat;
w = real(ifft(w_hat));
subplot(2,2,3), plot(x,v,".-",markersize=6);
axis([0,2pi,0,3]), grid(true);
subplot(2,2,4), plot(x,D*v,".-",markersize=6);
axis([0,2pi,-2,2]), grid(true);
error = signif(norm(D*v-vprime,Inf),4);
text(2.2,1.4,"max error = $error",fontsize=8);


text:
# ---
# title: 1037. Valid Boomerang
# id: problem1037
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/valid-boomerang/description/>
# hidden: true
# ---
#
# A _boomerang_ is a set of 3 points that are all distinct and **not** in a
# straight line.
#
# Given a list of three points in the plane, return whether these points are a
# boomerang.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,1],[2,3],[3,2]]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: [[1,1],[2,2],[3,3]]
#     Output: false
#
#
#
# **Note:**
#
#   1. `points.length == 3`
#   2. `points[i].length == 2`
#   3. `0 <= points[i][j] <= 100`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
# TODO: This all needs a fair bit of rethinking

"""
    accumulate(Δ, ∂)

Return `Δ + ∂` evaluated in a manner that supports ChainRulesCore's
various `AbstractDifferential` types.

See also: [`accumulate!`](@ref), [`store!`](@ref), [`AbstractRule`](@ref)
"""
accumulate(Δ, ∂) = Δ .+ ∂

"""
    accumulate!(Δ, ∂)

Similar to [`accumulate`](@ref), but attempts to compute `Δ + rule(args...)` in-place,
storing the result in `Δ`.

Note: this function may not actually store the result in `Δ` if `Δ` is immutable,
so it is best to always call this as `Δ = accumulate!(Δ, ∂)` just in-case.

This function is overloadable by using a [`InplaceThunk`](@ref).
See also: [`accumulate`](@ref), [`store!`](@ref).
"""
accumulate!(Δ, ∂) = store!(Δ, accumulate(Δ, ∂))

accumulate!(Δ::Number, ∂) = accumulate(Δ, ∂)

"""
    store!(Δ, ∂)

Stores `∂`, in `Δ`, overwriting what ever was in `Δ` before.
potentially avoiding intermediate temporary allocations that might be
necessary for alternative approaches  (e.g. `copyto!(Δ, extern(∂))`)

Like [`accumulate`](@ref) and [`accumulate!`](@ref), this function is intended
to be customizable for specific rules/input types.

See also: [`accumulate`](@ref), [`accumulate!`](@ref), [`AbstractRule`](@ref)
"""
store!(Δ, ∂) = materialize!(Δ, broadcastable(∂))


text:
using Test
using WaveProp.Nystrom
using StaticArrays

@testset "Kernels" begin
    pde  = Helmholtz(;dim=3,k=1)
    G    = SingleLayerKernel(pde)
    dG    = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == ComplexF64
    @test Nystrom.return_type(dG) == ComplexF64
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Laplace(;dim=3)
    G    = SingleLayerKernel(pde)
    dG    = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == Float64
    @test Nystrom.return_type(dG) == Float64
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Elastostatic(;μ=1,λ=2.0,dim=3)
    G    = SingleLayerKernel(pde)
    dG   = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == SMatrix{3,3,Float64,9}
    @test Nystrom.return_type(dG) == SMatrix{3,3,Float64,9}
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Maxwell(;k=1.0,dim=3)
    G    = SingleLayerKernel(pde)
    dG   = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == SMatrix{3,3,ComplexF64,9}
    @test Nystrom.return_type(dG) == SMatrix{3,3,ComplexF64,9}
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# To run this script, `cd` to the `./test/fixtures` directory and then, from the Julia terminal, `include("./runner.jl")`.

import JSON

# Negative small values:
x = linspace( -2.0^-29, -2.0^-54, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./small_negative.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive small values:
x = linspace( 2.0^-54, 2.0^-29, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./small_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Negative tiny values:
x = linspace( -2.0^-54, 0, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./tiny_negative.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive tiny values:
x = linspace( 0, 2.0^-54, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./tiny_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Negative medium values:
x = linspace( -0.41422, -0.2929, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./medium_negative.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive medium values:
x = linspace( 0.2929, 0.41422, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./medium_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Negative large values:
x = linspace( -0.9999, -0.41422, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./large_negative.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive large values:
x = linspace( 0.41422, 100, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./large_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive big values:
x = linspace( 2.0^30, 2.0^52, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./big_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )

# Positive huge values:
x = linspace( 1e200, 1e300, 500 )

y = log1p( x )
println( y )

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open( "./huge_positive.json", "w" )
JSON.json( data )

write( outfile, JSON.json(data) )
close( outfile )


text:
using Test, Pkg
using VizCLIMA

ENV["JULIA_LOG_LEVEL"] = "WARN"

function include_test(_module)
  println("Starting tests for $_module")
  t = @elapsed include(joinpath(_module,"runtests.jl"))
  println("Completed tests for $_module, $(round(Int, t)) seconds elapsed")
  return nothing
end

@testset "VizCLIMA" begin
    @test VizCLIMA.greet() == "Hello World!"
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false

    for submodule in [
                     ]

      if all_tests || "$submodule" in ARGS || "VizCLIMA" in ARGS
        include_test(submodule)
      end
    end

end


text:
import AbstractAlgebra.Generic: FreeModuleElem

@doc Markdown.doc"""
    convert_type(Q::MatrixElem, K::MatrixElem, d::RingElement) -> ZLat, Vector{RingElement}, RingElement

  Converts a quadratic triple QT = [Q, K, d] to the input values required for closest vector problem (CVP).
"""
function convert_type(G::MatrixElem, K::MatrixElem, d::RingElement)
    @assert typeof(G) == typeof(K) "Input matrices should be of the same type"
    if G[1,1] > 0
        Q = G
    else
        Q = -G
    end
    v = -inv(Q) * K
    upperbound = (-transpose(K) * inv(Q) * -K)[1] - d
    Lattice = Zlattice(gram = Q)
    vector = Vector{RingElement}(undef, size(v, 1))
    for i in 1:size(v,1)
        vector[i] = v[i]
    end
    return Lattice, vector, upperbound
end

@doc Markdown.doc"""
    convert_type(L::ZLat, v::Vector{RingElement}, c::RingElement)-> fmpq_mat, Vector{fmpq}, fmpq

  Converts the input values from closest vector enumeration (CVE) to the corresponding quadratic triple QT = [Q, K, d].
"""
function convert_type(L::ZLat, v::Vector{RingElement}, c::RingElement)
    V = ambient_space(L)
    G = gram_matrix(V)
    Q = G; K = -G*matrix(QQ,length(v),1,v); d = inner_product(V,v,v)-c;
    return Q, K, d
end

@doc Markdown.doc"""
    pojective_quadratic_triple(Q::MatrixElem, L::MatrixElem, c::RingElement, k=nothing) -> Tuple

  Returns a tuple, where the first item in the tuple is a kxk (positive definite)
  symmetric matrix, the next entry is a column vector of length k and the last
  entry is a rational number.

  THEORY:
  R, (x1,x2,x3) = PolynomialRing(QQ, ["x1", "x2", "x3"])
  Given a Polynomial ring of n variables we proceed as follows:
  Generate a list of projective quadratic triples of n-1, ..., 1 variables,
  given an n-variabled quadratic triple QT=[Q,L,c], where Q is a nxn symmetric
  matrix, L is a column vector of length n and c is a rational number.
  NOTE: if k is an int < n then it returns the k-variabled
  projective quadratic triple.
"""
function pojective_quadratic_triple(Q::MatrixElem, L::MatrixElem, c::RingElement, k=nothing)
    n = nrows(Q)
    Q1 = []; p1 = []; L1 = []; r1 = []; m1 = [];
    Lqqt = []; pp1 = []; LL1 = []; pL = [];
    #------------------------------------------------
    c1 = c
    pQT = []; ppQT = []; pLL = []; ppLL = []; cc1 = [];
    Lqt = [];
    #------------------------------------------------
    for i in 1:n-1
        Q1 = Q[1:n-i, 1:n-i] #n-i symmetric matrix extracted from the matrix Q
        p1 = Q[1:n-i, n+1-i:n+1-i] #Q[ri:rj, ci:cj] returns enrties from a matrix in the respective rows r_ and columns c_
        L1 = L[1:n-i,1:1] #column vector of length n-i
        #--------------------------------------------
        push!(Lqqt, Q1)
        push!(pp1, p1)
        push!(LL1, L1)
        push!(r1, Q[n+1-i, n+1-i])
        push!(m1, L[n+1-i,1:1])
        pL =[Lqqt, pp1, LL1, r1, m1]
    end
    #------------------------------------------------
    for j in 1:size(pL[1],1)
        pt1 = transpose(pL[2][j])
        pQT = pL[1][j] - pL[4][j]^-1 * (pL[2][j] * pt1)
        pLL = pL[3][j] - ((pL[5][j] * pL[4][j]^-1)[1] * pL[2][j])
        c1 = c1 - (pL[5][j]^2 * pL[4][j]^-1)[1] #rational number
        push!(ppQT, pQT)
        push!(ppLL, pLL)
        push!(cc1, c1)
    end
    Lqt = ppQT, ppLL, cc1 #set of n-1, ..., 1 variabled (projection) quadratic triple
    #------------------------------------------------
    if k === nothing
        return Lqt
    else
        if k < n
            while k >= 1
                return Lqt[1][n-k], Lqt[2][n-k], Lqt[3][n-k]
                break
            end
        else
            return Q, L, c #since k = n the result will be the input quadratic triple
        end
    end
end

@doc Markdown.doc"""
    range_ellipsoid_dim1(Q::MatrixElem, L::MatrixElem, c::RingElement) -> StepRange{Int, Int}

  Returns a finite set of integers for which the inhomogeneous quadratic function,
  given a one variabled quadratic triple, is less than or equal to zero.
"""
function range_ellipsoid_dim1(Q::MatrixElem, L::MatrixElem, c::RingElement)
    if nrows(Q) > 1
        throw(ArgumentError("The input quadratic triple is not one variabled."))
    end
    if L[1]^2 - Q[1]*c < 0
        throw(DomainError("The inhomogeneous quadratic function has no real solutions"))
    end
    #x1, x2 are roots of q s.t x1 <= x2
    #x1 = (-L[1] - sqrt(L[1]^2 - Q[1]*c)) / Q[1]
    #x2 = (-L[1] + sqrt(L[1]^2 - Q[1]*c)) / Q[1]
    # return the integers in the interval [x1,x2]
    #The following was achieved with Simon Brandhorst's and Tommy Hofmann's help
    sqrt_floor(a::fmpz) = isqrt(a)
    a = -L[1] // Q[1]; b = (L[1]^2 - Q[1]*c) // Q[1]^2;
    cc = lcm(denominator(a), denominator(b))
    aa = FlintZZ(a * cc); bb = FlintZZ(b * cc^2);
    l = round(fmpz, (aa - sqrt_floor(bb))//cc, RoundUp)
    r = round(fmpz, (aa + sqrt_floor(bb))//cc, RoundDown)
    return l:r
end

@doc Markdown.doc"""
    positive_quadratic_triple(a::fmpz, Q::MatrixElem, L::MatrixElem, c::RingElement) -> Tuple

  This function computes a n-1 variabled quadraric triple i*(a, QT), given an n>1
  variabled projective quadratic triple QT and an integer a.
  Return value: an n-1 variabled positive quadratic triple i*(a, QT)
"""
function positive_quadratic_triple(a::fmpz, Q::MatrixElem, L::MatrixElem, c::RingElement)
    n = nrows(Q)
    Q1 = []; p1 = []; L1 = []; r1 = []; m1 = [];
    Lqqt = []; pp1 = []; LL1 = []; pL = [];
    #------------------------------------------------
    pQT = []; ppQT = []; pLL = []; ppLL = []; cc1 = [];
    Lqt = [];
    #------------------------------------------------
    if n == 1
        throw(ArgumentError("The input should be a positive quadratic triple of n variables, for n > 1."))
    else
        for i in 2:n
            Q1 =  Q[i:n, i:n] #n-i symmetric matrix extracted from the matrix Q
            p1 = Q[i:n, i-1:i-1] #Q[ri:rj, ci:cj] returns enrties from a matrix in the respective rows r_ and columns c_
            L1 = L[i:n,1:1] #column vector of length n-i
            #--------------------------------------------
            push!(Lqqt, Q1)
            push!(pp1, p1)
            push!(LL1, L1)
            push!(r1, Q[i-1, i-1])
            push!(m1, L[i-1])
            pL =[Lqqt, pp1, LL1, r1, m1]
        end
        #------------------------------------------------
        pQT = pL[1][1]
        pLL = a * pL[2][1] + pL[3][1]
        c1 = a^2 * pL[4][1] + 2 * a * pL[5][1] + c #rational number
        #------------------------------------------------
        return Lqt = pQT, pLL, c1 #n-1 variabled quadratic triple
    end
end

@doc Markdown.doc"""
    positive_quadratic_triple2(aa::Array{fmpz,1}, Q::MatrixElem,
                               L::MatrixElem, c::RingElement) -> Tuple

  NOTE: aa is an array of length k.

  This function computes an n-k variabled quadratic triple i*(aa, QT), given an n
  variabled quadratic triple QT and a tuple aa = {a1, ..., ak}, where k < n.
  Return value: an n-k variabled positive quadratic triple i*(aa, QT).
"""
function positive_quadratic_triple2(aa::Array{fmpz,1}, Q::MatrixElem, L::MatrixElem, c::RingElement)
    QT = Q, L, c
    if typeof(aa) != String
        k = length(aa)
        if k == 1
            R = positive_quadratic_triple(fmpz(aa[1]), QT[1], QT[2], QT[3])
            QT = R[1], R[2], R[3]
        else
            for i in 1:k
                R = positive_quadratic_triple(fmpz(aa[i]), QT[1], QT[2], QT[3])
                QT = R[1], R[2], R[3]
            end
        end
        if typeof(QT[1]) != Char
            return QT
        else
            error("The input QT = {Q, L, c} should be at least 2 variabled.")
        end
    else
        error("The first input value is not an integer array.")
    end
end

@doc Markdown.doc"""
    list_positive_quadratic_triple(aa::Array{fmpz,1}, Q::MatrixElem,
                                   L::MatrixElem, c::RingElement) -> Array

  Computes a list (for v = 1) QT_m^1 == i*(aa, QT_m^0), for aa in E(QT_1^0)
"""
function list_positive_quadratic_triple(aa::Array{fmpz,1}, Q::MatrixElem, L::MatrixElem, c::RingElement)
    n = nrows(Q)
    v = length(aa)+1
    ListI = []
    for m in v:n
        P = pojective_quadratic_triple(Q, L, c, m) #QT_m^0
        Li = positive_quadratic_triple2(aa, P[1], P[2], P[3])
        push!(ListI, Li)
    end
    return ListI
end

@doc Markdown.doc"""
    list_positive_quadratic_triple2(b::fmpz, ListIv::Array{Any,1}) -> Array

  For a fixed v, this function computes a list of QT_m^v := i*(b, QT_m^{v-1}),
  where b is in E(QT_v^{v-1}).
"""
function list_positive_quadratic_triple2(b::fmpz, ListIv::Array{Any,1}) #ListIv is the list of quadratic triples QT_m^{v-1}
    LL = [];
    for i in 1:size(ListIv,1)-1
        L = positive_quadratic_triple(b, ListIv[i+1][1], ListIv[i+1][2], ListIv[i+1][3])
        push!(LL, L)
    end
    return LL
end

@doc Markdown.doc"""
    posQuadTriple_intVectors(QT::Array{Any,1}, E::Array{Any,1}) -> Tuple

  INPUT:
    A list of quadratic triples QT a list of tuples E of the form {a_1, ..., a_k}
  OUTPUT:
    A new list of quadratic triples in the one dimension higher, and an updated
  list E of the form {a_1, ..., a_k, a_{k+1}}
"""
function posQuadTriple_intVectors(QT::Array{Any,1}, E::Array{Any,1})
    QTT1 = [];
    EE = E;
    QTT = Array{Array}(undef, size(QT, 1))
    bbb = [];
    for j in 1:size(QT,1) #for a fixed v
        b = range_ellipsoid_dim1(QT[j][1][1], QT[j][1][2], QT[j][1][3]) #E(QT_v^{v-1}) = {b_1, ..., b_N}
        QT1 = Array{Array}(undef, BigInt(size(b, 1)))
        for k in 1:BigInt(size(b, 1))
            QT1[k] = list_positive_quadratic_triple2(b[k], QT[j]) #QT_m^v for m = v+1, ..., n
            for kk in EE
                push!(bbb,vcat(kk,b[k])) #bb is now a list of tuples [x_1, ..., x_{v-1},b[k]] that satisfiy a inhomogeneous quad function for a v variabled quad triple
            end
        end
        QTT[j] = QT1
    end
    for j in 1:size(QT, 1)
        append!(QTT1, QTT[j])
    end
    EEE = unique!(bbb)
    return QTT1, EEE
end

@doc Markdown.doc"""
    closest_vectors(Q::MatrixElem, L::MatrixElem, c::RingElement)
                                                    -> Array{Array{fmpz, 1}, 1}

  Returns all the integer vectors of length n such that the inhomogeneous
  quadratic function corresponding to an n variabled quadratic triple is
  less than or equal to zero.
"""
function closest_vectors(G::MatrixElem, L::MatrixElem, c::RingElement)
    #1 < v <= n+1, a = [a_1, ..., a_{v-1}] int tuple &
    if G[1,1] > 0
        Q = G
    else
        Q = -G
    end
    # @assert det(Q) != 0 "the symmetric matrix is not positive definite"
    if det(Q) == 0
        throw("the symmetric matrix is not indefinite")
    else
        n = nrows(Q)
        QT = []; QTT = []; bbb = []; EE = Array{Array{fmpz,1},1}();
        P = pojective_quadratic_triple(Q, L, c, 1)
        aa = range_ellipsoid_dim1(P[1], P[2], P[3])#E(QT_1^0)
        for i in 1:size(aa,1)
            qt = list_positive_quadratic_triple(fmpz[aa[i]], Q, L, c) #List  QT_m^1 i*(aa[i], QT_m^0) for m = 2,...,n. Here v = 1
            push!(QT, qt)
        end
        for j in 1:size(QT,1) #v=2
            b = range_ellipsoid_dim1(QT[j][1][1], QT[j][1][2], QT[j][1][3])#E(QT_2^1) = {b_1, ..., b_N}
            QT1 = Array{Array}(undef, BigInt(size(b, 1)))
            for kk in 1:size(aa, 1)
                for k in 1:BigInt(size(b, 1))
                    QT1[k] = list_positive_quadratic_triple2(b[k], QT[j]) #QT_m^2 for m = 3, ..., n
                    bb = [aa[kk], b[k]]
                    push!(bbb,bb) #bbb is a list of tuples [a,b] that satisfiy a inhomogeneous quad function for a 2 variabled quad triple
                end
            end
            QT[j] = QT1
        end
        for j in 1:size(QT, 1)
            append!(QTT, QT[j])
        end
        QTT1 = QTT;
        E = unique!(bbb)
        for v in 3:n
            R1 = posQuadTriple_intVectors(QTT1, E)
            QTT1 = R1[1]
            E = R1[2]
        end
        #------------------------------------------------
        for k in E
            if (transpose(matrix(QQ,size(k,1),1,k))*Q*matrix(QQ,size(k,1),1,k))[1] + (2*transpose(matrix(QQ,size(k,1),1,k))*L)[1] + c <= 0
                push!(EE, k)
            end
        end
        return EE
    end
end

@doc Markdown.doc"""
    closest_vectors(L:ZLat, v:Vector, c::fmpq) -> Vector{Tuple{Vector{Int}, fmpq}} -> Array{Array{fmpz, 1}, 1}

  Returns all vectors `x` in `L` such that `b(v-x,v-x) <= c`, where `b` is the bilinear form on `L`.
  written by: Simon Brandhorst
"""
function closest_vectors(L::ZLat, v::Vector{RingElement} , upperbound::RingElement)
    epsilon = QQ(1//10)   # some number > 0, not sure how it influences performance
    d = size(v)[1]
    if isdefinite(L) == false
        throw("Zlattice is indefinite.")
    else
        if rank(L) != d
            throw("Zlattice must have the same rank as the length of the vector in the second argument.")
        else
            g1 = gram_matrix(L)
            if g1[1,1] > 0
                G1 = g1
            else
                G1 = -g1
            end
            e = matrix(QQ, 1, 1, [upperbound//3+epsilon])
            G = diagonal_matrix(G1, e)
            B = diagonal_matrix(basis_matrix(L),matrix(QQ,1,1,[1]))
            for i in 1:d[1]
                B[end,i] = -v[i]
            end
            N = Zlattice(B,gram=G)

            delta = QQ(4//3)*upperbound + epsilon
            sv = Hecke.short_vectors(N, delta)
            cv = Array{Array{fmpz,1},1}()
            for a in sv
                a = a[1]
                if a[end] == 0
                continue
                end
                if a[end] == -1
                    a = -a
                end
                x = a[1:end-1]
                push!(cv, x)
            end
            V = ambient_space(L)
            # debug remove later when performance matters
            for x in cv
                t = x - v
                dist = inner_product(V,t,t)
                @assert dist<= upperbound
            end
            return sort!(cv)
        end
    end
end


text:
# todo : document and add some other nasty functions
using SparseArrays

"""
    Mesh1d

A type describing a mesh for a one dimensional domain.
"""
struct Mesh1d
    nodes::Array{Float64, 1}
    elements::Array{Int64, 2}
end

"""
    assemble_1d_FE_matrix(elem::Array{Float64, 2}, nbNodes::Int64;
      intNodes = 0, dof1 = 1, dof2 = 1)

Assemble a finite elements matrix corresponding to a 1 dimensional uniform mesh.

# Arguments
  * `elem` : elementary finite elements matrix
  * `nbNodes` : number of nodes in the mesh
  * `intNodes1` : number of interior nodes in the interior of each element for lhs
  * `intNodes2` : number of interior nodes in the interior of each element for rhs
  * `dof1`     : number of degrees of freedom for each node for lhs
  * `dof2`     : number of degrees of freedom for each node for rhs
"""
function assemble_1d_FE_matrix(elem::Array{Float64, 2}, nbNodes::Int64;
                               intNodes1 = 0, intNodes2 = 0, dof1 = 1, dof2 = 1)

    nT = Threads.nthreads()

    nbNodesTotal1 = (nbNodes + (nbNodes - 1) * intNodes1)
    nbNodesTotal2 = (nbNodes + (nbNodes - 1) * intNodes2)

    MM = Array{SparseMatrixCSC{Float64,Int64}}(undef, nT)
    for i = 1:nT
        MM[i] = spzeros(Float64, nbNodesTotal1 * dof1, nbNodesTotal2 * dof2)
    end

    l1 = zeros(UInt64, nT)
    r1 = zeros(UInt64, nT)
    l2 = zeros(UInt64, nT)
    r2 = zeros(UInt64, nT)

    Threads.@threads for i = 1:(nbNodes - 1)
        k = Threads.threadid()

        l1[k] = (i - 1 + (i - 1) * intNodes1) * dof1 + 1
        r1[k] = (i + 1 + i * intNodes1) * dof1
        l2[k] = (i - 1 + (i - 1) * intNodes2) * dof2 + 1
        r2[k] = (i + 1 + i * intNodes2) * dof2
        MM[k][l1[k]:r1[k], l2[k]:r2[k]] += elem
    end
    M = MM[1]
    for i=2:nT
        M += MM[i]
    end
    M
end

"""
    assemble_1d_nu_FE_matrix(elem::Matrix{SymPy.Sym}, nodes::Array{Float64, 1};
      intNodes = 0, dof1 = 1, dof2 = 1)

Assemble a finite elements matrix corresponding to a 1 dimensional non-uniform mesh.

# Arguments
  * `elem` : elementary finite elements matrix (with elements of type SymPy.Sym)
  * `nodes` : vector of nodes composing the non uniform-mesh
  * `intNodes1` : number of interior nodes in the interior of each element for lhs
  * `intNodes2` : number of interior nodes in the interior of each element for rhs
  * `dof1`     : number of degrees of freedom for each node for lhs
  * `dof2`     : number of degrees of freedom for each node for rhs
"""
function assemble_1d_nu_FE_matrix(elem::Matrix{SymPy.Sym}, nodes::Array{Float64, 1};
                                  intNodes1 = 0, intNodes2 = 0, dof1 = 1, dof2 = 1)

    nT = Threads.nthreads()
    global h
    nbNodes = length(nodes)
    nbNodesTotal1 = (nbNodes + (nbNodes - 1) * intNodes1)
    nbNodesTotal2 = (nbNodes + (nbNodes - 1) * intNodes2)

    MM = Array{SparseMatrixCSC{Float64,Int64}}(undef, nT)
    for i = 1:nT
        MM[i] = spzeros(Float64, nbNodesTotal1 * dof1, nbNodesTotal2 * dof2)
    end

    l1 = zeros(UInt64, nT)
    r1 = zeros(UInt64, nT)
    l2 = zeros(UInt64, nT)
    r2 = zeros(UInt64, nT)
    Threads.@threads for i = 1:nbNodes - 1
        k = Threads.threadid()
        l1[k] = (i - 1 + (i - 1) * intNodes1) * dof1 + 1
        r1[k] = (i + 1 + i * intNodes1) * dof1
        l2[k] = (i - 1 + (i - 1) * intNodes2) * dof2 + 1
        r2[k] = (i + 1 + i * intNodes2) * dof2
        elem_loc = elem.subs(h, nodes[i + 1] - nodes[i])
        MM[k][l1[k]:r1[k], l2[k]:r2[k]] += elem_loc
    end
    M = MM[1]
    for i=2:nT
        M += MM[i]
    end
    M
end


text:
using HetaSimulator, Plots

p = load_platform(".", rm_out=false)

scenarios = read_scenarios("data-mumenthaler-2000/scenarios.csv")
add_scenarios!(p, scenarios)

data = read_measurements("data-mumenthaler-2000/data.csv")
data_scn = Dict()
data_scn[:scn1] = filter(:scenario => ==(:scn1), data)
data_scn[:scn2] = filter(:scenario => ==(:scn2), data)
data_scn[:scn3] = filter(:scenario => ==(:scn3), data)

loss_add = Dict()
loss_add[:scn1] = 2*sum(log.(data_scn[:scn1][:,"prob.sigma"])) + length(data_scn[:scn1][:,"prob.sigma"])*log(2π)
loss_add[:scn2] = 2*sum(log.(data_scn[:scn2][:,"prob.sigma"])) + length(data_scn[:scn1][:,"prob.sigma"])*log(2π)
loss_add[:scn3] = 2*sum(log.(data_scn[:scn3][:,"prob.sigma"])) + length(data_scn[:scn1][:,"prob.sigma"])*log(2π)

add_measurements!(p, data)

# initial plot
sim(p) |> plot
sim(p, parameters_upd = [:Vmax=>0.3, :ke=>0., :k_a=>5.]) |> plot

# fitting 1 best -63.45
params_df_1 = read_parameters("./julia/parameters-1.csv")
res_fit_1 = fit(p, params_df_1; scenarios=[:scn1], ftol_abs=1e-6, ftol_rel=0.)
res_fit_1 = fit(p, optim(res_fit_1); scenarios=[:scn1], ftol_abs=1e-6, ftol_rel=0.)
fig = sim(p, scenarios=[:scn1], parameters_upd=optim(res_fit_1)) |> plot
# savefig(fig, "diagnostics/scn1_best.png")

# fitting 2 best -63.64
params_df_2 = read_parameters("./julia/parameters-2.csv")
res_fit_2 = fit(p, params_df_2; scenarios=[:scn2], ftol_abs=1e-6, ftol_rel=0.)
res_fit_2 = fit(p, optim(res_fit_2); scenarios=[:scn2], ftol_abs=1e-6, ftol_rel=0.)
fig = sim(p, scenarios=[:scn2], parameters_upd=optim(res_fit_2)) |> plot
# savefig(fig, "diagnostics/scn2_best.png")

# fitting 3 best -63.92
params_df_3 = read_parameters("./julia/parameters-3.csv")
res_fit_3 = fit(p, params_df_3; scenarios=[:scn3], ftol_abs=1e-6, ftol_rel=0.)
res_fit_3 = fit(p, optim(res_fit_3); scenarios=[:scn3], ftol_abs=1e-6, ftol_rel=0.)
fig = sim(p, scenarios=[:scn3], parameters_upd=optim(res_fit_3)) |> plot
# savefig(fig, "diagnostics/scn3_best.png")

############################## Identification ##########################

using LikelihoodProfiler, CSV

chi_level = 3.84

p_optim_1 = optim(res_fit_1)
p_optim_2 = optim(res_fit_2)
p_optim_3 = optim(res_fit_3)

sim_scn1 = sim(p.scenarios[:scn1], parameters_upd=p_optim_1)
sim_scn2 = sim(p.scenarios[:scn2], parameters_upd=p_optim_2)
sim_scn3 = sim(p.scenarios[:scn3], parameters_upd=p_optim_3)

p_names = Dict(
  :scn1 => first.(p_optim_1),
  :scn2 => first.(p_optim_2),
  :scn3 => first.(p_optim_3),
)

function loss_func(params::Vector{P}, scen) where P<: Pair
  sim_res = sim(p.scenarios[scen]; parameters_upd=params, reltol=1e-6, abstol=1e-8)
  #sim_res = last(sim_vec[1])
  return loss(sim_res, sim_res.scenario.measurements) - loss_add[scen]
end

function loss_func(pvals::Vector{N}, scen) where N <: Number
 # @show pvals
  @assert length(p_names[scen]) == length(pvals) "Number of params values doesn't match params names"
  params = [pn => pv for (pn,pv) in zip(p_names[scen],pvals)]
  loss_func(params,scen)
end

loss_scn1(params) = loss_func(params, :scn1)
loss_scn2(params) = loss_func(params, :scn2)
loss_scn3(params) = loss_func(params, :scn3)

function scan_func(params::Vector{P}, timepoint, scen) where P<: Pair
  sim_vec = sim(p; parameters_upd=params, scenarios=[scen])
  return last(sim_vec[1])(timepoint)[:BrAC]
end

function scan_func(pvals::Vector{N}, timepoint, scen) where N <: Number
  @assert length(p_names[scen]) == length(pvals) "Number of params values doesn't match params names"
  params = [pn => pv for (pn,pv) in zip(p_names[scen],pvals)]
  scan_func(params, timepoint, scen)
end

scan_scn1(params, timepoint) = scan_func(params, timepoint, :scn1)
scan_scn2(params, timepoint) = scan_func(params, timepoint, :scn2)
scan_scn3(params, timepoint) = scan_func(params, timepoint, :scn3)

saveat_1 = saveat(p.scenarios[:scn1])
saveat_2 = saveat(p.scenarios[:scn2])
saveat_3 = saveat(p.scenarios[:scn3])

p_ident_1 = [get_interval(
  last.(p_optim_1),
  i,
  loss_scn1,
  :CICO_ONE_PASS,
  theta_bounds = fill((1e-10,1e10), length(p_names[:scn1])),
  scan_bounds=((last.(p_optim_1)[i])/1e4, (last.(p_optim_1)[i])*1e4),
  scan_tol=1e-5,
  #scale =  fill(:log, length(p_names[:scn1])),
  loss_crit = loss_scn1(p_optim_1) + chi_level
  ) for i in eachindex(p_names[:scn1])]

plb_1 = [iv.result[1].value for iv in p_ident_1]
pub_1 = [iv.result[2].value for iv in p_ident_1]

pliter_1 = [iv.result[1].counter for iv in p_ident_1]
puiter_1 = [iv.result[2].counter for iv in p_ident_1]

df = DataFrame(params = p_names[:scn1], optim = last.(p_optim_1), lower = plb_1, upper = pub_1, liter = pliter_1, uiter = puiter_1)
df.lower = replace(df.lower, nothing => missing)
df.upper = replace(df.upper, nothing => missing)
CSV.write("./julia/scn1_intervals.csv", df)

p_ident_2 = [get_interval(
  last.(p_optim_2),
  i,
  loss_scn2,
  :CICO_ONE_PASS,
  #theta_bounds = fill((-10.,10.), length(p_names[:scn2])),
  #scale = fill(:log, length(p_names[:scn2])),
  loss_crit = loss_scn2(p_optim_2) + chi_level
) for i in eachindex(p_names[:scn2])]

p_ident_3 = [get_interval(
  last.(p_optim_3),
  i,
  loss_scn3,
  :CICO_ONE_PASS,
  scale = fill(:log, length(p_names[:scn3])),
  loss_crit = loss_scn3(p_optim_3) + chi_level
) for i in eachindex(p_names[:scn3])]

BrAC_ident_1 = [get_interval(
  last.(p_optim_1),
  params->scan_scn1(params,t),
  loss_scn1,
  :CICO_ONE_PASS,
  scale = fill(:log, length(p_names[:scn1])),
  loss_crit = loss_scn1(p_optim_1) + chi_level
) for t in saveat_1]

BrAC_ident_2 = [get_interval(
  last.(p_optim_2),
  params->scan_scn2(params,t),
  loss_scn2,
  :CICO_ONE_PASS,
  scale = fill(:log, length(p_names[:scn2])),
  loss_crit = loss_scn2(p_optim_2) + chi_level
) for t in saveat_2]

BrAC_ident_3 = [get_interval(
  last.(p_optim_3),
  params->scan_scn3(params,t),
  loss_scn3,
  :CICO_ONE_PASS,
  scale = fill(:log, length(p_names[:scn3])),
  loss_crit = loss_scn3(p_optim_3) + chi_level
) for t in saveat_3]

lb_1 = [iv.result[1].value for iv in BrAC_ident_1]
ub_1 = [iv.result[2].value for iv in BrAC_ident_1]

liter_1 = [iv.result[1].counter for iv in BrAC_ident_1]
uiter_1 = [iv.result[2].counter for iv in BrAC_ident_1]

df = DataFrame(times = saveat_1, lower = lb_1, upper = ub_1, liter = liter_1, uiter = uiter_1)
df.lower = replace(df.lower, nothing => missing)
df.upper = replace(df.upper, nothing => missing)
CSV.write("./julia/scn1_conf_band.csv", df)

BrAC_1 = sim_scn1.(saveat_1, :BrAC)
plot(sim_scn1, show_measurements=false, vars=[:BrAC])
scatter!(data_scn[:scn1].t, data_scn[:scn1].measurement, yerror=data_scn[:scn1][!,"prob.sigma"], label = "Measurements")
plot!(saveat_1, lb_1, fillrange =  ub_1, fillalpha = 0.35, c = 1, label = "Confidence band")
savefig("./julia/conf_band.png")

lb_2 = [iv.result[1].value for iv in BrAC_ident_2]
ub_2 = [iv.result[2].value for iv in BrAC_ident_2]
BrAC_2 = sim_scn2.(saveat_2, :BrAC)
plot(sim_scn2, show_measurements=false, ribbon = (BrAC_2-lb_1,ub_2-BrAC_2), fc=:orange, fa=0.7)

lb_3 = [iv.result[1].value for iv in BrAC_ident_3]
ub_3 = [iv.result[2].value for iv in BrAC_ident_3]
BrAC_3 = sim_scn3.(saveat_3, :BrAC)
plot(sim_scn3, show_measurements=false, ribbon = (BrAC_3-lb_3,ub_3-BrAC_3), fc=:orange, fa=0.7)

### Validation band

t_scn1 = data_scn[:scn1].t
sigma_scn1 = data_scn[:scn1][!,"prob.sigma"]
sim_scn1 = sim(p.scenarios[:scn1]; parameters_upd=p_optim_1, reltol=1e-6, abstol=1e-8)
BrAC_scn1 = sim_scn1.(t_scn1, :BrAC)

function valid_obj1(params, i)
  d1 = last(params)
  _params = [pn => pv for (pn,pv) in zip(p_names[:scn1],params[1:end-1])]

  sim_res = sim(p.scenarios[:scn1]; parameters_upd=_params, reltol=1e-6, abstol=1e-8)
  d1_sim = sim_res(t_scn1[i])[:BrAC]
  return loss(sim_res, sim_res.scenario.measurements) + (d1 - d1_sim)^2/(sigma_scn1[i])^2 - loss_add[:scn1]
end

valid_1 = []
for i in eachindex(t_scn1)
println(" Calculating CI for $(t_scn1[i]) timepoint")
push!(valid_1,
  get_interval(
    [last.(p_optim_1); BrAC_scn1[i]],
    5,
    p->valid_obj1(p,i),
    :CICO_ONE_PASS,
    theta_bounds = fill((1e-8,1e8), length(p_names[:scn1])+1),
    scan_bounds=(1e-7,1e7),
    scan_tol=1e-5,
    scale =  fill(:log, length(p_names[:scn1])+1),
    loss_crit = loss_scn1(p_optim_1) + chi_level)
  )
end

lb_1 = [iv.result[1].value for iv in valid_1]
ub_1 = [iv.result[2].value for iv in valid_1]

liter_1 = [iv.result[1].counter for iv in valid_1]
uiter_1 = [iv.result[2].counter for iv in valid_1]

df = DataFrame(times = t_scn1, lower = lb_1, upper = ub_1, liter = liter_1, uiter = uiter_1)
df.lower = replace(df.lower, nothing => missing)
df.upper = replace(df.upper, nothing => missing)
CSV.write("./julia/scn1_valid_bans.csv", df)

lb_1[1] = 0.0
lb_1[end-4:end] .= 0.0
plot(sim_scn1, show_measurements=false, vars=[:BrAC])
scatter!(data_scn[:scn1].t, data_scn[:scn1].measurement, yerror=data_scn[:scn1][!,"prob.sigma"], label = "Measurements")
plot!(t_scn1, lb_1, fillrange =  ub_1, fillalpha = 0.35, c = 1, label = "Validation band")
savefig("./julia/valid_band.png")

max_stars_repo_name: ivborissov/alco

max_issues_repo_name: ivborissov/alco

max_forks_repo_name: ivborissov/alco


text:
using UUIDs

@testset "Niancat game state" begin
    @testset "State has puzzle PUSSGUKRA; Get the puzzle; Puzzle is PUSSGUKRA" begin
        # Arrange
        round = uuid4()
        state = "PUSSGUKRA#sv-14#$(round)"
        game = NiancatGame(state, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, GetPuzzle())

        # Assert
        @test response == PuzzleIs("PUSSGUKRA")
    end

    @testset "State has a given round; Get the game round; The round is same as the state" begin
        # Arrange
        round = uuid4()
        state = "PUSSGUKRA#sv-14#$(round)"
        game = NiancatGame(state, TestingGameService())

        # Act
        actualround = gameround(game)

        # Assert
        @test string(round) == actualround
    end

    @testset "State has dictionary sv-14; Requests dictionary sv-14" begin
        # Arrange
        round = uuid4()
        state = "PUSSGUKRA#sv-14#$(round)"
        gameservice = TestingGameService()

        # Act
        game = NiancatGame(state, gameservice)

        # Assert
        @test gameservice.dictionaryrequests == ["sv-14"]
    end

    @testset "State has dictionary en-00; Requests dictionary en-00" begin
        # Arrange
        round = uuid4()
        state = "PUSSGUKRA#en-00#$(round)"
        gameservice = TestingGameService()

        # Act
        game = NiancatGame(state, gameservice)

        # Assert
        @test gameservice.dictionaryrequests == ["en-00"]
    end

    @testset "State has no puzzle or round; No puzzle is set" begin
        # Arrange
        state = "#sv-14#"
        gameservice = TestingGameService()
        game = NiancatGame(state, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, GetPuzzle())

        # Assert
        @test response == NoPuzzleSet()
    end
end


text:
# Things to test
# [x] Construction limited to ImageOperation
# [x] Construction also works for ops like ElasticDistortion

@test (Either <: Augmentor.AffineOperation) == false
@test Either <: Augmentor.ImageOperation

@testset "constructor" begin
    @test_throws ArgumentError Either()
    @test_throws ArgumentError Either(())
    @test_throws ArgumentError Either((NoOp(),),(0,))
    @test_throws ArgumentError Either((NoOp(),),(-1,))
    @test_throws MethodError Either((NoOp(),),(1,1))
    @test_throws MethodError Either(SplitChannels(), NoOp())
    @test_throws MethodError SplitChannels() * NoOp()
    @test_throws MethodError SplitChannels() * SplitChannels()
    let op = @inferred Either(Rotate90(), 0.3)
        @test op.operations === (Rotate90(), NoOp())
        @test op.chances === @SVector([0.3,0.7])
        @test op.cum_chances === @SVector([0.3,1.0])
    end
    let op = @inferred Either(Rotate90())
        @test op.operations === (Rotate90(), NoOp())
        @test op.chances === @SVector([0.5,0.5])
        @test op.cum_chances === @SVector([0.5,1.0])
    end
    let op = @inferred Either(1 => Rotate90(), 2 => Rotate180(), 1 => Crop(1:10,1:10))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances === @SVector([0.25,0.5,0.25])
        @test op.cum_chances === @SVector([0.25,0.75,1.])
    end
    let op = @inferred Either(1 => Rotate90(), 2 => Rotate180(), 1 => Crop(1:10,1:10))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances === @SVector([0.25,0.5,0.25])
        @test op.cum_chances === @SVector([0.25,0.75,1.])
    end
    let op = @inferred Either((Rotate90(),Rotate180(),Crop(1:10,1:10)), (0.1,0.7,0.2))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances === @SVector([0.1,0.7,0.2])
        @test op.cum_chances ≈ @SVector([0.1,0.8,1.])
    end
    let op = @inferred Either((Rotate90(),Rotate180(),Crop(1:10,1:10)))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances ≈ @SVector([1/3,1/3,1/3])
        @test op.cum_chances ≈ @SVector([1/3,2/3,3/3])
    end
    let op = Either(Rotate90(), Rotate180(), Crop(1:10,1:10), chances = [1,7,2])
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances === @SVector([0.1,0.7,0.2])
        @test op.cum_chances ≈ @SVector([0.1,0.8,1.])
    end
    let op = @inferred Either(Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances ≈ @SVector([1/3,1/3,1/3])
        @test op.cum_chances ≈ @SVector([1/3,2/3,3/3])
    end
    let op = @inferred(Rotate90()*Rotate180()*Crop(1:10,1:10))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances ≈ @SVector([1/3,1/3,1/3])
        @test op.cum_chances ≈ @SVector([1/3,2/3,3/3])
    end
    let op = @inferred(Rotate90()*Rotate180()*ElasticDistortion(5))
        @test op.operations === (Rotate90(),Rotate180(),ElasticDistortion(5))
        @test op.chances ≈ @SVector([1/3,1/3,1/3])
        @test op.cum_chances ≈ @SVector([1/3,2/3,3/3])
    end
    let op = @inferred((1=>Rotate90())*(2=>Rotate180())*(1=>Crop(1:10,1:10)))
        @test op.operations === (Rotate90(),Rotate180(),Crop(1:10,1:10))
        @test op.chances ≈ @SVector([1/4,2/4,1/4])
        @test op.cum_chances ≈ @SVector([1/4,3/4,4/4])
    end
end

@testset "show" begin
    @test str_show(Either((Rotate90(),Rotate270(),NoOp()), (0.2,0.3,0.5))) == """
    Either (1 out of 3 operation(s)):
      - 20% chance to: Rotate 90 degree
      - 30% chance to: Rotate 270 degree
      - 50% chance to: No operation"""
    @test str_showcompact(Either((Rotate90(),Rotate270(),NoOp()), (0.2,0.3,0.5))) ==
        "Either: (20%) Rotate 90 degree. (30%) Rotate 270 degree. (50%) No operation."
    @test str_show(Either((Rotate90(),Rotate270(),NoOp()), (0.15,0.8,0.05))) == """
    Either (1 out of 3 operation(s)):
      - 15% chance to: Rotate 90 degree
      - 80% chance to: Rotate 270 degree
      -  5% chance to: No operation"""
    @test str_showcompact(Either((Rotate90(),Rotate270(),NoOp()), (0.15,0.8,0.05))) ==
        "Either: (15%) Rotate 90 degree. (80%) Rotate 270 degree. (5%) No operation."
    @test str_show(Either((Rotate90(),Rotate270(),NoOp()), (0.155,0.8,0.045))) == """
    Either (1 out of 3 operation(s)):
      - 15.5% chance to: Rotate 90 degree
      - 80.0% chance to: Rotate 270 degree
      -  4.5% chance to: No operation"""
    @test str_showcompact(Either((Rotate90(),Rotate270(),NoOp()), (0.155,0.8,0.045))) ==
        "Either: (16%) Rotate 90 degree. (80%) Rotate 270 degree. (4%) No operation."
    @test str_showconst(Either(Rotate90(), Rotate270(), NoOp())) == "Rotate90() * Rotate270() * NoOp()"
    @test str_showconst(Either((Rotate90(), Rotate270(), NoOp()),(1,1,2))) == "(0.25=>Rotate90()) * (0.25=>Rotate270()) * (0.5=>NoOp())"
    @test str_showconst(Either((Rotate90(), NoOp()),(1,2))) == "(0.333=>Rotate90()) * (0.667=>NoOp())"
end

@testset "eager" begin
    @test_throws MethodError Augmentor.applyeager(Either(Rotate90(),NoOp()), nothing)
    imgs = [
        (rect),
        (view(rect, :, :)),
        (Augmentor.prepareaffine(rect)),
        (OffsetArray(rect, -2, -1)),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
    ]
    for img in imgs
        let op = @inferred Either((Rotate90(),ElasticDistortion(5)), (1,0))
            @test_throws MethodError Augmentor.applylazy(op, img)
            @test Augmentor.supports_eager(op) === true
            @test Augmentor.supports_affine(op) === false
            @test Augmentor.supports_lazy(op) === false
            @test Augmentor.supports_affineview(op) === false
            @test Augmentor.supports_view(op) === false
            @test Augmentor.supports_stepview(op) === false
            @test Augmentor.supports_permute(op) === false
            @test @inferred(Augmentor.applyeager(op, img)) == rotl90(rect)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
            res1, res2 = @inferred(Augmentor.applyeager(op, (square2, img)))
            @test res1 == rotl90(square2)
            @test res2 == rotl90(rect)
        end
        let op = @inferred Either((Rotate90(),Rotate270()), (1,0))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == rotl90(rect)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),Rotate270(),Crop(1:2,2:3)), (0,1,0))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == rotr90(rect)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),Rotate270(),NoOp()), (0,0,1))
            @test Augmentor.supports_eager(op) === true
            if img isa Union{Array,OffsetArray}
                @test parent(@inferred(Augmentor.applyeager(op, img))) === rect
            else
                @test @inferred(Augmentor.applyeager(op, img)) == rect
                @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
            end
        end
        let op = @inferred Either((Rotate90(),Rotate270(),Crop(1:2,2:3)), (0,0,1))
            @test Augmentor.supports_eager(op) === true
            @test collect(@inferred(Augmentor.applyeager(op, img))) == rect[1:2,2:3]
            @test_throws MethodError Augmentor.applyaffine(op, rect)
            @test_throws MethodError Augmentor.applyview(op, rect)
            @test_throws MethodError Augmentor.applystepview(op, rect)
            @test_throws MethodError Augmentor.applypermute(op, rect)
        end
        let op = @inferred Either((Rotate90(),Zoom(.8)), (1,0))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == rotl90(rect)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),FlipX()), (1,0))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == rotl90(rect)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),FlipX()), (0,1))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == reverse(rect; dims=2)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),FlipY()), (0,1))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == reverse(rect; dims=1)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),Resize(5,5)), (0,1))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == imresize(rect,5,5)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Crop(1:2,1:2),Resize(5,5)), (0,1))
            @test Augmentor.supports_eager(op) === true
            @test @inferred(Augmentor.applyeager(op, img)) == imresize(rect,5,5)
            @test typeof(Augmentor.applyeager(op, img)) <: OffsetArray
        end
        let op = @inferred Either((Rotate90(),Rotate270(),Crop(1:2,2:3)))
            @test Augmentor.supports_eager(op) === true
            res1, res2 = @inferred(Augmentor.applyeager(op, (N0f8.(img), img)))
            @test res1 == res2
        end
    end
end

@testset "affine" begin
    let op = @inferred Either((Rotate90(),Rotate(-90)), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 1)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        res1, res2 = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine.((square2, square)))
        @test res1 == rotl90(square2)
        @test res2 == rotl90(square)
    end
    let op = @inferred Either((Rotate90(),Rotate270()), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === true
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 1)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    let op = @inferred Either((Rotate180(),FlipX(),FlipY()), (0,1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 2)) ≈ AffineMap([1. 0.; 0. -1.], [0,4])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == reverse(square; dims=2)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    let op = @inferred Either((Rotate90(),FlipY()), (0,1))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 2)) ≈ AffineMap([-1. 0.; 0. 1.], [3,0])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == reverse(square; dims=1)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    let op = @inferred Either((Rotate90(),Rotate270()), (0,1))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === true
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 2)) ≈ AffineMap([6.12323e-17 1.0; -1.0 6.12323e-17], [-0.5,3.5])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotr90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    for op = (@inferred(Rotate90(1)), @inferred(Either((Rotate90(),Scale(0.8)),(1,0))))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applyaffine(op, nothing)
        @test @inferred(Augmentor.toaffinemap(op, rect, 1)) ≈ AffineMap([6.12323e-17 -1.0; 1.0 6.12323e-17], [3.5,0.5])
        wv = @inferred Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test wv == rotl90(square)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test parent(wv2).itp.coefs === square
        @test wv2 == rotl90(square)
        @test typeof(wv2) == typeof(wv)
    end
    let op = @inferred Either((Rotate180(),FlipX(),FlipY()))
        res1, res2 = @inferred(Augmentor.applyaffine(op, (N0f8.(square), square)))
        @test res1 == res2
        @test typeof(res1) <: InvWarpedView{N0f8,2}
        @test typeof(res2) <: InvWarpedView{eltype(square),2}
    end
end

@testset "affineview" begin
    let op = @inferred Either((Rotate90(),Rotate270(),Crop(1:2,1:2)),(1,0,0))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.toaffinemap(op, nothing)
        @test_throws MethodError Augmentor.toaffinemap(op, rect, 1)
        @test_throws MethodError Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        wv = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == rotl90(square)
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test typeof(wv) == typeof(wv2)
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2 == wv
        res1, res2 = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine.((square2,square)))
        @test res1 == rotl90(square2)
        @test res2 == rotl90(square)
    end
    let op = @inferred Either((Rotate90(),Rotate270(),CropSize(2,2)),(0,0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.toaffinemap(op, nothing)
        @test_throws MethodError Augmentor.toaffinemap(op, rect, 1)
        @test_throws MethodError Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        wv = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == view(Augmentor.prepareaffine(square), IdentityRange(1:2), IdentityRange(1:2))
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test typeof(wv) == typeof(wv2)
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2 == wv
    end
    let op = @inferred Either((Rotate90(),Zoom(.8)), (0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.toaffinemap(op, nothing)
        @test_throws MethodError Augmentor.toaffinemap(op, rect, 1)
        @test_throws MethodError Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        wv = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test typeof(wv) == typeof(wv2)
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2 == wv
    end
    let op = @inferred Either((Rotate90(),Resize(2,3)), (0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.toaffinemap(op, nothing)
        @test_throws MethodError Augmentor.toaffinemap(op, rect, 1)
        @test_throws MethodError Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        wv = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == imresize(square, (2,3))
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test typeof(wv) == typeof(wv2)
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2 == wv
    end
    let op = @inferred Either((Crop(1:2,1:2),Resize(2,3)), (1,0))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.toaffinemap(op, nothing)
        @test_throws MethodError Augmentor.toaffinemap(op, rect, 1)
        @test_throws MethodError Augmentor.applyaffine(op, Augmentor.prepareaffine(square))
        wv = @inferred Augmentor.applyaffineview(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == square[1:2,1:2]
        wv2 = @inferred Augmentor.applylazy(op, square)
        @test typeof(wv) == typeof(wv2)
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2 == wv
    end
    let op = @inferred Either((Rotate90(),Rotate(45),Crop(1:2,1:2)))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === false
        res1, res2 = @inferred(Augmentor.applyaffineview(op, Augmentor.prepareaffine.((N0f8.(square), square))))
        @test res1 == res2

        @test typeof(res1) <: SubArray{N0f8,2,<:InvWarpedView}
        @test typeof(res2) <: SubArray{eltype(square),2,<:InvWarpedView}
    end
end

@testset "view" begin
    let op = @inferred Either((NoOp(),Crop(1:2,2:3)), (1,0))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === true
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test @inferred(Augmentor.applyview(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
        @test @inferred(Augmentor.applylazy(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
        @test @inferred(Augmentor.applyaffineview(op, Augmentor.prepareaffine(rect))) == view(Augmentor.prepareaffine(rect), IdentityRange(1:2), IdentityRange(1:3))
        res1, res2 = @inferred(Augmentor.applylazy(op, (square2, rect)))
        @test res1 === view(square2, IdentityRange(1:4), IdentityRange(1:4))
        @test res2 === view(rect, IdentityRange(1:2), IdentityRange(1:3))
    end
    let op = @inferred Either((NoOp(),Crop(1:2,2:3)), (0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === true
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test @inferred(Augmentor.applyview(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(2:3))
        @test @inferred(Augmentor.applylazy(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(2:3))
        @test @inferred(Augmentor.applyaffineview(op, Augmentor.prepareaffine(rect))) == view(Augmentor.prepareaffine(rect), IdentityRange(1:2), IdentityRange(2:3))
    end
    let op = @inferred Either((Crop(1:2,2:4),CropSize(2,3)), (0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === true
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test @inferred(Augmentor.applyview(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
        @test @inferred(Augmentor.applylazy(op, rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
        @test @inferred(Augmentor.applyaffineview(op, Augmentor.prepareaffine(rect))) == view(Augmentor.prepareaffine(rect), IdentityRange(1:2), IdentityRange(1:3))
    end
    let op = @inferred Either((Crop(1:2,2:4),CropSize(2,3),CropRatio(1)))
        tres1, tres2 = @inferred(Augmentor.applylazy(op, (N0f8.(square2), square2)))
        res1, res2 = @inferred(Augmentor.applyview(op, (N0f8.(square2), square2)))
        @test typeof((tres1,tres2)) == typeof((res1,res2))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:Array}
        @test typeof(res2) <: SubArray{eltype(square),2,<:Array}
        res1, res2 = @inferred(Augmentor.applyview(op, (N0f8.(square2), Augmentor.prepareaffine(square2))))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:Array}
        @test typeof(res2) <: SubArray{eltype(square),2,<:InvWarpedView}
    end
end

@testset "stepview" begin
    let op = @inferred Either((Rotate180(),FlipX()), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,2:-1:1,3:-1:1)
        @test v == rot180(rect)
        @test typeof(v) <: SubArray
        res1, res2 = @inferred Augmentor.applystepview(op, (square2, rect))
        @test res1 === view(square2,4:-1:1,4:-1:1)
        @test res2 === view(rect,2:-1:1,3:-1:1)
    end
    let op = @inferred Either((Rotate180(),FlipX()), (0,1))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,1:1:2,3:-1:1)
        @test v == reverse(rect; dims=2)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((FlipY(),FlipX()), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,2:-1:1,1:1:3)
        @test v == reverse(rect; dims=1)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((Rotate180(),NoOp()), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,2:-1:1,3:-1:1)
        @test v == rot180(rect)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((Rotate180(),NoOp(),Crop(1:2,2:3)),(0,1,0))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,1:1:2,1:1:3)
        @test v == rect
        @test typeof(v) <: SubArray
        wv = @inferred Augmentor.applylazy(op, Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == square
    end
    let op = @inferred Either((Rotate180(),NoOp(),Crop(1:2,2:3)),(0,0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applystepview(op, rect))
        @test v === view(rect,1:1:2,2:1:3)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((Rotate180(),CropSize(2,3)),(0,1))
        @test Augmentor.supports_affine(op) === false
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === true
        @test Augmentor.supports_permute(op) === false
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applystepview(op, nothing)
        v = @inferred Augmentor.applylazy(op, square)
        @test v === @inferred(Augmentor.applystepview(op, square))
        @test v === view(square,1:1:2,1:1:3)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((Rotate180(),NoOp(),Crop(1:2,2:3)))
        tres1, tres2 = @inferred(Augmentor.applylazy(op, (N0f8.(square2), square2)))
        res1, res2 = @inferred(Augmentor.applystepview(op, (N0f8.(square2), square2)))
        @test typeof((tres1,tres2)) == typeof((res1,res2))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:Array}
        @test typeof(res2) <: SubArray{eltype(square),2,<:Array}
        res1, res2 = @inferred(Augmentor.applystepview(op, (N0f8.(square2), Augmentor.prepareaffine(square2))))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:Array}
        @test typeof(res2) <: SubArray{eltype(square),2,<:InvWarpedView}
    end
end

@testset "permute" begin
    let op = @inferred Either((Rotate90(),Rotate270()), (1,0))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === true
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applypermute(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applypermute(op, rect))
        @test v === view(PermutedDimsArray(rect,(2,1)),3:-1:1,1:1:2)
        @test v == rotl90(rect)
        @test typeof(v) <: SubArray
        res1, res2 = @inferred Augmentor.applylazy(op, (square2, rect))
        @test res1 == rotl90(square2)
        @test res2 == rotl90(rect)
    end
    let op = @inferred Either((Rotate90(),Rotate270()), (0,1))
        @test Augmentor.supports_affine(op) === true
        @test Augmentor.supports_lazy(op) === true
        @test Augmentor.supports_affineview(op) === true
        @test Augmentor.supports_view(op) === false
        @test Augmentor.supports_stepview(op) === false
        @test Augmentor.supports_permute(op) === true
        @test_throws MethodError Augmentor.applylazy(op, nothing)
        @test_throws MethodError Augmentor.applypermute(op, nothing)
        v = @inferred Augmentor.applylazy(op, rect)
        @test v === @inferred(Augmentor.applypermute(op, rect))
        @test v === view(PermutedDimsArray(rect,(2,1)),1:1:3,2:-1:1)
        @test v == rotr90(rect)
        @test typeof(v) <: SubArray
    end
    let op = @inferred Either((Rotate90(),Rotate270()))
        tres1, tres2 = @inferred(Augmentor.applylazy(op, (N0f8.(square2), square2)))
        res1, res2 = @inferred(Augmentor.applypermute(op, (N0f8.(square2), square2)))
        @test typeof((tres1,tres2)) == typeof((res1,res2))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:PermutedDimsArray}
        @test typeof(res2) <: SubArray{eltype(square),2,<:PermutedDimsArray}
        res1, res2 = @inferred(Augmentor.applypermute(op, (N0f8.(square2), Augmentor.prepareaffine(square2))))
        @test res1 == res2
        @test typeof(res1) <: SubArray{N0f8,2,<:PermutedDimsArray}
        @test typeof(res2) <: SubArray{eltype(square),2,<:PermutedDimsArray}
        @test typeof(parent(parent(res2))) <: InvWarpedView
    end
end