text:
using GaussianProcessODEs
using Test
using BenchmarkTools

using KernelFunctions
using QuadGK

rot(phi) = GaussianProcessODEs.Krot(phi)

struct altrotinvKernel{K <: KernelFunctions.Kernel, GKP }  <: GaussianProcessODEs.MatrixKernel
    kernel::K
    gkparams::GKP
end

# constructor
function altrotinvKernel(ker::K; N::Int) where K
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, 0, 2pi))
    gkparams = zip(gkparams.x, gkparams.weights)
    gkparams = Iterators.product(gkparams, gkparams)
    altrotinvKernel{K, typeof(gkparams)}(ker, gkparams)
end

function (rk::altrotinvKernel)(z1, z2)
    gkp = rk.gkparams
    ker = rk.kernel

    mapreduce( g -> g[1][2]*ker(rot(g[1][1]) * z1, rot(g[2][1]) * z2)*g[2][2], +, gkp)
end

function (mker::altrotinvKernel)(w::Array{T, 1}) where {T <: Real}
    typeof(mker)(mker.kernel(w), mker.gkparams)
end

import KernelFunctions: kernelmatrix
function kernelmatrix(ker::altrotinvKernel, X)
    K = zeros(length(X), length(X))
    for i in 1:length(X), j in 1:i
        tmp = ker(X[i], X[j])
        K[i,j] = K[j,i] = tmp
    end
#     for (i, z1) in enumerate(X), (j,z2) in enumerate(X)
# #         K[i,j] = threadriker(ker, z1, z2)
#         K[i,j] = ker(z1,z2)
#     end
    return K
end

### setup

using LinearAlgebra
H(q,p) = norm(p)^2/2 - inv(norm(q))

borders = [[-2,2], [-2,2], [-2, 2], [-2, 2]]
ngp = 4
Z = rectanglegrid(borders, ngp) .+ [0.1*(rand(4) .- 0.5) for i in 1:ngp^length(borders)]
U = H.(getindex.(Z, Ref([1,2])), getindex.(Z, Ref([3,4])) );

### testing implementation above
@testset "invariance of reference" begin
    N = 55
    ker = pskernel(ones(5))
    riker = altrotinvKernel(ker; N);

    let phi = rand()*2pi
        @test (riker(Z[1], Z[2]) ≈ riker(rot(phi)*Z[1], Z[2]) ) atol = 1e-5
        @test (riker(Z[1], Z[2]) ≈ riker(Z[1], rot(phi)*Z[2]) ) atol = 1e-5
    end

end

### testing package implementation
@testset "invariance of package implementation" begin
    N = 55
    ker = pskernel(ones(5))
    tmpZ = Z[1:60];
    criker = crotinvKernel(ker; N)

    let phi = rand()*2pi
        @test (criker(Z[1], Z[2]) ≈ criker(rot(phi)*Z[1], Z[2]) ) atol = 1e-5
        @test (criker(Z[1], Z[2]) ≈ criker(Z[1], rot(phi)*Z[2]) ) atol = 1e-5
    end

end

### comparison
let
    N = 55
    ker = pskernel(ones(5))
    riker = altrotinvKernel(ker; N)
    criker = crotinvKernel(ker; N)
    tmpZ = Z[1:60]

    println("Comparison")
    println("reference")
    @time K1 = kernelmatrix(riker, tmpZ);
    K1 = kernelmatrix(riker, tmpZ);

    println("package")
    println("\r setup")
    @time tmpKx = Kx(criker, tmpZ);
    tmpKx = Kx(criker, tmpZ)
    println("\r kernel matrix")
    @time K2 = tmpKx();
    # @btime K2 = $tmpKx();
    K2 = tmpKx();

    @testset "accuracy" begin
        @test sum(K1 .- K2) < 1e-8
    end

    @testset "kernelmatrix()" begin
        @time K1 = kernelmatrix(criker, tmpZ)
        @time K2 = kernelmatrix(riker, tmpZ)
        @test K1 ≈ K2 atol=1e-6
    end
end

### test package implementation

## complete prediction
@testset "Kx" begin
    N = 55
    ker = pskernel(ones(5))
    riker = altrotinvKernel(ker; N)
    criker = crotinvKernel(ker; N)

    tmpZ = Z[1:60]
    tmpKx = Kx(criker, tmpZ)

    x = rand(4)
    @time Kx1 = tmpKx(x);
    @time Kx2 = riker.([x], tmpZ)';

    @test sum(Kx1[:] .- Kx2[:]) < 1e-6
end

@testset "Complete Prediction" begin
    N = 55
    ker = pskernel(ones(5))
    criker = crotinvKernel(ker; N)

    nZ = 6
    minZ = Z[1:nZ]
    minU = U[1:nZ]

    testker = pskernel([1., 12 .* ones(4)...])
    criker = crotinvKernel(testker; N)
    tmpKx2 = Kx(criker, minZ)
    K = tmpKx2()

    KinvU = K \ minU;
    # gpH(q,p) = (riker.([vcat(q,p)], permutedims(Zcl)) * KinvU)[1];
    gpf(z) = (tmpKx2(z) * KinvU)[1]

    @testset "recover data points" begin
        @test maximum(abs.(gpf.(minZ) .- minU)) < 1e-6
    end

    r1 = criker(Z[1], Z[2])
    phis = criker.gkparams.x
    rz1 = rot.(phis) .* [Z[1]]
    rz2 = rot.(phis) .* [Z[2]]

    r2 = criker(rz1, rz2)

    @testset "kernel itself" begin
        @test r1 ≈ r2 atol=1e-6
    end
end

### Derivative kernel

using Flux
@testset "Base derivative kernel" begin
    ker = pskernel(ones(5))
    dker = dKernel(ker)
    r1 = dker(Z[1], Z[2])
    r2 = gradient(x1 -> ker(x1, Z[2]), Z[1])[1]

    @test r1 ≈ r2 atol = 1e-6
end

### Invariant derivative kernel

@testset "Invariant derivative kernel" begin
    N = 60
    ker = pskernel(ones(5))
    dker = dKernel(ker)
    criker = crotinvKernel(ker; N)
    dcriker = dcrotinvKernel(dker; N)

    @time r1 = dcriker(Z[1], Z[2])
    @time r2 = gradient(x1 -> criker(x1, Z[2]), Z[1])[1]

    @test r1 ≈ r2 atol=1e-6
end

println("done")


text:
module ImplicitAdjoints

    export genlasso, fake_genlasso, L1, TV
    export FISTA
    export run_reconstruction_test
    export Identity
    export proximal, transform_op, projection_op, support, L1, TV, L1Project, TVProject, Gradient

    using ChainRulesCore
    import ChainRulesCore.rrule
    using LinearAlgebra
    import LinearAlgebra: mul!
    using IterativeSolvers: cg, powm
    using LinearMaps
    using LinearMaps: CompositeMap, TransposeMap, BlockMap, AdjointMap
    import LinearMaps._unsafe_mul!
    using Zygote: @adjoint # TODO: don't want this

    include("linear_maps.jl")
    include("regularizers.jl")
    include("fista.jl")
    include("lasso.jl")

end

ext: jl


text:
export BicategoryRelations, FreeBicategoryRelations,
  AbelianBicategoryRelations, FreeAbelianBicategoryRelations,
  mplus, mzero, coplus, cozero

# Bicategory of relations
#########################

""" Doctrine of *bicategory of relations*

TODO: The 2-morphisms are missing. I haven't decided how to handle them yet.

References:

- Carboni & Walters, 1987, "Cartesian bicategories I"
- Walters, 2009, blog post, "Categorical algebras of relations",
  http://rfcwalters.blogspot.com/2009/10/categorical-algebras-of-relations.html
"""
@signature MonoidalCategoryWithDiagonals(Ob,Hom) => BicategoryRelations(Ob,Hom) begin
  # Dagger category.
  dagger(f::Hom(A,B))::Hom(B,A) <= (A::Ob,B::Ob)

  # Self-dual compact closed category.
  dunit(A::Ob)::Hom(munit(), otimes(A,A))
  dcounit(A::Ob)::Hom(otimes(A,A), munit())
end

@syntax FreeBicategoryRelations(ObExpr,HomExpr) BicategoryRelations begin
  otimes(A::Ob, B::Ob) = associate_unit(Super.otimes(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(Super.otimes(f,g))
  compose(f::Hom, g::Hom) = associate(Super.compose(f,g; strict=true))

  function dagger(f::Hom)
    f = anti_involute(Super.dagger(f), dagger, compose, id)
    distribute_unary(f, dagger, otimes)
  end
end

""" Doctrine of *abelian bicategory of relations*

Unlike Carboni & Walters, we use additive notation and nomenclature.

References:

- Carboni & Walters, 1987, "Cartesian bicategories I", Sec. 5
- Baez & Erbele, 2015, "Categories in control"
"""
@signature BicategoryRelations(Ob,Hom) => AbelianBicategoryRelations(Ob,Hom) begin
  # Second diagonal and codiagonal.
  mplus(A::Ob)::Hom(otimes(A,A),A)
  mzero(A::Ob)::Hom(munit(),A)
  coplus(A::Ob)::Hom(A,otimes(A,A))
  cozero(A::Ob)::Hom(A,munit())
end

@syntax FreeAbelianBicategoryRelations(ObExpr,HomExpr) AbelianBicategoryRelations begin
  otimes(A::Ob, B::Ob) = associate_unit(Super.otimes(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(Super.otimes(f,g))
  compose(f::Hom, g::Hom) = associate(Super.compose(f,g; strict=true))

  function dagger(f::Hom)
    f = anti_involute(Super.dagger(f), dagger, compose, id)
    distribute_unary(f, dagger, otimes)
  end
end

function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Ob{:otimes}; kw...)
  Syntax.show_latex_infix(io, expr, "\\oplus"; kw...)
end
function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Hom{:otimes}; kw...)
  Syntax.show_latex_infix(io, expr, "\\oplus"; kw...)
end
show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Ob{:munit}; kw...) = print(io, "O")

function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Hom{:mplus}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacktriangledown")
end
function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Hom{:coplus}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacktriangle")
end
function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Hom{:mzero}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacksquare")
end
function show_latex(io::IO, expr::FreeAbelianBicategoryRelations.Hom{:cozero}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacklozenge")
end


text:
@testset "Appender Error" begin
    db = DBInterface.connect(DuckDB.DB)
    con = DBInterface.connect(db)

    @test_throws DuckDB.QueryException DuckDB.Appender(db, "nonexistanttable")
    @test_throws DuckDB.QueryException DuckDB.Appender(con, "t")
end

@testset "Appender Usage" begin
    db = DBInterface.connect(DuckDB.DB)

    DBInterface.execute(db, "CREATE TABLE integers(i INTEGER)")

    appender = DuckDB.Appender(db, "integers")
    DuckDB.close(appender)
    DuckDB.close(appender)

    appender = DuckDB.Appender(db, "integers")
    for i in 0:9
        DuckDB.append(appender, i)
        DuckDB.end_row(appender)
    end
    DuckDB.flush(appender)

    results = DBInterface.execute(db, "SELECT * FROM integers")
    df = DataFrame(results)
    @test names(df) == ["i"]
    @test size(df, 1) == 10
    @test df.i == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
end

# @testset "Appender API" begin
#     # Open the database
#     db = DuckDB.open(":memory:")
#     con = DuckDB.connect(db)
#
#     # Create the table the data is appended to
#     DuckDB.execute(con, "CREATE TABLE dtypes(bol BOOLEAN, tint TINYINT, sint SMALLINT, int INTEGER, bint BIGINT, utint UTINYINT, usint USMALLINT, uint UINTEGER, ubint UBIGINT, float FLOAT, double DOUBLE, date DATE, time TIME, vchar VARCHAR, nullval INTEGER)")
#
#     # Create the appender
#     appender = DuckDB.appender_create(con, "dtypes")
#
#     # Append the different data types
#     DuckDB.duckdb_append_bool(appender, true)
#     DuckDB.duckdb_append_int8(appender, 1)
#     DuckDB.duckdb_append_int16(appender, 2)
#     DuckDB.duckdb_append_int32(appender, 3)
#     DuckDB.duckdb_append_int64(appender, 4)
#     DuckDB.duckdb_append_uint8(appender, 1)
#     DuckDB.duckdb_append_uint16(appender, 2)
#     DuckDB.duckdb_append_uint32(appender, 3)
#     DuckDB.duckdb_append_uint64(appender, 4)
#     DuckDB.duckdb_append_float(appender, 1.0)
#     DuckDB.duckdb_append_double(appender, 2.0)
#     DuckDB.duckdb_append_date(appender, 100)
#     DuckDB.duckdb_append_time(appender, 200)
#     DuckDB.duckdb_append_varchar(appender, "Foo")
#     DuckDB.duckdb_append_null(appender)
#     # End the row of the appender
#     DuckDB.duckdb_appender_end_row(appender)
#     # Destroy the appender and flush the data
#     DuckDB.duckdb_appender_destroy(appender)
#
#     # Retrive the data from the table and store it in  a vector
#     df = DuckDB.toDataFrame(con, "select * from dtypes;")
#     data = Matrix(df)
#
#     # Test if the correct types have been appended to the table
#     @test data[1] === true
#     @test data[2] === Int8(1)
#     @test data[3] === Int16(2)
#     @test data[4] === Int32(3)
#     @test data[5] === Int64(4)
#     @test data[6] === UInt8(1)
#     @test data[7] === UInt16(2)
#     @test data[8] === UInt32(3)
#     @test data[9] === UInt64(4)
#     @test data[10] === Float32(1.0)
#     @test data[11] === Float64(2.0)
#     @test data[12] === Dates.Date("1970-04-11")
#     @test data[13] === Dates.Time(0, 0, 0, 0, 200)
#     @test data[14] === "Foo"
#     @test data[15] === missing
#
#     # Disconnect and close the database
#     DuckDB.disconnect(con)
#     DuckDB.close(db)
# end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
function plot_multiclasses(X,y; ms=8,fs=20, disp_legend=1, disp_dim=2, edge_factor=.25)
"""
 plot_multiclasses - display data for classification

   plot_multiclasses(X,y,options);

   Copyright (c) 2017 Gabriel Peyre
"""

Xm = X -> X - repeat(mean(X,1), outer=(size(X,1), 1))
Cov = X -> Xm(X)'*Xm(X)
col = [ [1 0 0]; [0 1 0]; [0 0 1]; [0 0 0]; [0 1 1]; [1 0 1]; [1 1 0];
    [1 .5 .5]; [.5 1 .5]; [.5 .5 1]  ]'

p,n = size(X);
# list of classes
CL = unique(y)
k = length(CL)

# dimensionality reduction
if n>disp_dim
    U,D,V = svd(Xm(X),thin=true)
    Z = Xm(X) * V
else
    Z = X
end

lgd = []
for i=1:min(k,size(col,2))
    I = findall(y.==CL[i])
    if disp_dim==2
        plot(Z[I,1], Z[I,2], "o", c=col[:,i], ms=ms,label=CL[i])
    elseif disp_dim==3
        plot3D(Z[I,1], Z[I,2], Z[I,3], "o", c=col[:,i], ms=ms, label=CL[i])
    else
        error("Works only in 2D and 3D.")
    end
    append!(lgd,CL[i]);
end

axis("tight"); axis("equal"); box("on");gca()[:grid](false)
if disp_dim==3
    gca()[:view_init](50, 70)
end
if disp_legend==1
    legend(loc="best");
end

end


text:
"""
A shape where anything inside could cause a disturbance at the listener position
from a planar wavefront parallel to the y axis starting at the listener. Also
everything inside has a positive `x` coordinate.

More precisely, if the listener is at (l_x,l_y) then the interior of the shape
is defined as
x-l_x+sqrt((x-l_x)^2+(y-l_y)^2)<time and x>0
"""
struct TimeOfFlight{T <: AbstractFloat} <: Shape{T,2}
    listener_position::Vector{T}
    time::T
end

TimeOfFlight(pos::AbstractVector{T}, time::T) where T <:AbstractFloat = TimeOfFlight(Vector{T}(pos), time)

name(shape::TimeOfFlight) = "Time of flight from planar source"

# Let T = shape.time, xr = shape.listener_position[1], yr = shape.listener_position[2]
# and assume all particles are placed in x>0.
# Let (x,y) be a point on the curved part of the shape, then
# x - xr + sqrt((x - xr)^2 + (y - yr)^2) == T => x == T/2 + xr - (y - yr)^2/(2T)
# then the area = 2*Integrate[ T/2 + xr - (y - yr)^2/(2T), {y,yr, yr + sqrt(T^2 + 2xr*T)}]
function volume(shape::TimeOfFlight{T}) where T <: AbstractFloat
    l_x = shape.listener_position[1]
    return 2/(3*shape.time)*(shape.time^2 + 2*l_x*shape.time)^(3//2)
end

import Base.issubset
function issubset(circle::Sphere{T,2}, shape::TimeOfFlight{T}) where T
    l_to_p = origin(circle) - shape.listener_position
    return (origin(circle)[1] > 0) && (l_to_p[1] + norm(l_to_p) <= (shape.time - 2circle.radius))
end

function bounding_box(shape::TimeOfFlight{T}) where T <: AbstractFloat
    t = shape.time
    l = shape.listener_position
    x_max = max(t/2 + l[1], zero(T))
    return Box([SVector(zero(T), l[2] - sqrt(t^2 + 2t*l[1])), SVector(x_max, l[2] + sqrt(t^2 + 2t*l[1]))])
end

function boundary_functions(shape::TimeOfFlight)

    function x(τ)
        check_boundary_coord_range(τ)
        if τ <= 1//2
            return zero(τ)
        else
            t = shape.time
            l = shape.listener_position
            y_sq = (4*(3//4-τ))^2*(t^2 + 2t*l[1])
            return l[1] + (t^2-y_sq)/(2t)
        end
    end

    function y(τ)
        check_boundary_coord_range(τ)
        t = shape.time
        l = shape.listener_position
        if τ <= 1//2
            return l[2] + 4*(τ-1//4)*sqrt(t^2 + 2t*l[1])
        else
            return l[2] + 4*(3//4-τ)*sqrt(t^2 + 2t*l[1])
        end
    end

    return x, y
end


text:
function edit_distance(a, b, submat::AbstractSubstitutionMatrix{T}, insertion_cost::T, deletion_cost::T) where T
    m = length(a)
    n = length(b)
    trace = Matrix{Trace}(undef, (m + 1, n + 1))
    D = Vector{T}(undef, m + 1)
    D[1] = T(0)
    for i in 1:m
        D[i+1] = i * insertion_cost
        trace[i+1,1] = TRACE_INSERT
    end
    for j in 1:n
        d_diag = D[1]
        D[1] = j * deletion_cost
        trace[1,j+1] = TRACE_DELETE
        for i in 1:m
            ins = D[i]   + insertion_cost
            del = D[i+1] + deletion_cost
            mat = d_diag + submat[a[i],b[j]]
            # find the best score and its trace
            best = min(del, ins, mat)
            t = TRACE_NONE
            del == best && (t |= TRACE_DELETE)
            ins == best && (t |= TRACE_INSERT)
            mat == best && (t |= TRACE_MATCH)
            d_diag = D[i+1]
            D[i+1] = best
            trace[i+1,j+1] = t
        end
    end
    return D[m+1], trace, (m, n)
end

function edit_traceback(a, b, trace, endpos)
    anchors = Vector{AlignmentAnchor}()
    i, j = endpos
    @start_traceback
    while i ≥ 1 || j ≥ 1
        t = trace[i+1,j+1]
        if t & TRACE_MATCH > 0
            if a[i] == b[j]
                @anchor OP_SEQ_MATCH
            else
                @anchor OP_SEQ_MISMATCH
            end
        elseif t & TRACE_DELETE > 0
            @anchor OP_DELETE
        elseif t & TRACE_INSERT > 0
            @anchor OP_INSERT
        else
            error("failed to trace back")
        end
    end
    @finish_traceback
    return AlignedSequence(a, anchors)
end


text:
####
# We need different models from the UKF for dispatch

type AugmentedUnscentedModel <: Model
    f::Function
    q::Matrix
end

type AugmentedUnscentedObservationModel <: Model
    h::Function
    r::Matrix
end

type AugmentedUnscentedKalmanFilter <: KalmanFilter
    x::UnscentedState
    f::AugmentedUnscentedModel
    z::AugmentedUnscentedObservationModel
end

function ap(f::AugmentedUnscentedModel,s::UnscentedState)
    n = length(s.x)
    (x,p) = augment(s,f.q)
    (σn,wm,wc) = sigma(x,p,s.α,s.β,s.κ)
    σs = zeros(n,size(σn,2))
    xn = zeros(s.x)
    for i in 1:size(σn,2)
        σs[:,i] = f.f(σn[:,i])
        xn += wm[i] * σs[:,i]
    end
    pn = zeros(s.p)
    for i in 1:size(σn,2)
        pn += wc[i] * (σs[:,i]-xn)*(σs[:,i]-xn)'
    end
    return UnscentedState(xn,pn,s.α,s.β,s.κ)
end

function covs(kf::AugmentedUnscentedKalmanFilter,y::Observation)
    n = size(kf.x.x,1)
    m = size(kf.z.r,1)
    (x,p) = augment(kf.x,kf.z.r)
    (σn,wm,wc) = sigma(x,p,kf.x.α,kf.x.β,kf.x.κ)
    yp = zeros(m,size(σn,2))
    yhat = zeros(m)

    for i in 1:size(σn,2)
        yp[:,i] = kf.z.h(σn[:,i])
        yhat += wm[i] * yp[:,i]
    end

    resx = zeros(σn[1:n,:])
    resy = zeros(yp)

    for i in 1:size(σn,2)
        resx[:,i] = σn[1:n,i]-kf.x.x
        resy[:,i] = yp[:,i]-yhat
    end

    res = y.y-yhat
    ph = zeros(n,m)
    s = zeros(m,m)

    for i in 1:size(σn,2)
        ph += wc[i]*(resx[:,i]*resy[:,i]')
        s += wc[i]*(resy[:,i]*resy[:,i]')
    end

    return res,ph,s
end

# Helper function takes a state and a covariance matrix,
# augments them appropriately, and returns the augmented
# state and covariance
function augment(s::UnscentedState,R::Matrix)
    n = length(s.x)
    m = size(R,1)
    x = [s.x,zeros(m)]
    p = [s.p zeros(n,m);
         zeros(m,n) R]
    return x,p
end


text:
using Base.Test

import AffineTransforms

if VERSION.minor == 2
    macro mytest_throws(extype, arg)
        :(@test_throws $(esc(arg)))
    end
else
    macro mytest_throws(extype, arg)
        :(@test_throws $(esc(extype)) $(esc(arg)))
    end
end

# Pure translations
for nd = 1:4
    t1 = rand(nd)
    t2 = rand(nd)
    ident = eye(nd)

    tfm1 = AffineTransforms.AffineTransform(ident, t1)
    x = zeros(nd)
    @test AffineTransforms.tformfwd(tfm1, x) == t1
    @mytest_throws DimensionMismatch AffineTransforms.tformfwd(tfm1, zeros(nd+1))
    @test AffineTransforms.tforminv(tfm1, x) == -t1
    x = rand(nd)
    @test AffineTransforms.tformfwd(tfm1, x) == x+t1
    @test AffineTransforms.tforminv(tfm1, x) == x-t1
    @test tfm1*x == x+t1
    @test tfm1\x == x-t1
    X = reshape(x, nd, 1)
    @test all(AffineTransforms.tformfwd(tfm1, X) .== x+t1)
    @test all(AffineTransforms.tforminv(tfm1, X) .== x-t1)

    tfm2 = AffineTransforms.tformtranslate(t2)
    tfmP1 = tfm1*tfm2
    @test tfmP1.offset == t1+t2
    @test_approx_eq AffineTransforms.tformfwd(tfmP1, x) x+t2+t1

    if 2 <= nd <= 3
        # Tuple interface
        r = AffineTransforms.tformfwd(tfm1, x...)
        @test [r...] == tfm1*x
        r = AffineTransforms.tforminv(tfm1, x...)
        @test [r...] == tfm1\x
    end
end

# Pure rotations
@test_approx_eq AffineTransforms.tformrotate(0).scalefwd eye(2)
@test_approx_eq AffineTransforms.tformrotate(pi/2)*[1,0] [0,1]
@test_approx_eq AffineTransforms.tformrotate(pi/2).scalefwd [0 -1; 1 0]
@test_approx_eq AffineTransforms.tformrotate(pi)*[1,0] [-1,0]
@test_approx_eq AffineTransforms.tformrotate(3pi/2)*[1,0] [0,-1]
@test_approx_eq AffineTransforms.tformrotate([1,0,0],pi/2).scalefwd [1 0 0; 0 0 -1; 0 1 0]
@test_approx_eq AffineTransforms.tformrotate([1,0,0],pi/2)*[1,0,0] [1,0,0]
@test_approx_eq AffineTransforms.tformrotate([0,1,0],pi/2)*[1,0,0] [0,0,-1]
@test_approx_eq AffineTransforms.tformrotate([0,0,1],pi/2)*[1,0,0] [0,1,0]
s2 = 1/sqrt(2)
@test_approx_eq AffineTransforms.tformrotate([0,1,0],pi/4).scalefwd [s2 0 s2; 0 1 0; -s2 0 s2]
@test_approx_eq AffineTransforms.tformrotate([0,0,1],pi/4).scalefwd [s2 -s2 0; s2 s2 0; 0 0 1]

theta1 = pi/2
# 2d
tfm1 = AffineTransforms.tformrotate(theta1)
tfm2 = AffineTransforms.tformrotate(theta1)
x = rand(2)
@test_approx_eq tfm1*(tfm2*x) -x
@test_approx_eq (tfm1*tfm2)*x -x
@test_approx_eq tfm1\(tfm2*x)  x
@test_approx_eq (tfm1\tfm2)*x  x
@test_approx_eq (tfm1*tfm1*tfm1*tfm1)*x  x
r = AffineTransforms.tformfwd(tfm1, x...)
@test [r...] == tfm1*x
r = AffineTransforms.tforminv(tfm1, x...)
@test [r...] == tfm1\x

# 3d
uaxis = rand(3); uaxis = uaxis/norm(uaxis)
tfm1 = AffineTransforms.tformrotate(uaxis*theta1)
tfm2 = AffineTransforms.tformrotate(uaxis*theta1)
x = rand(3)
xparallel = dot(x,uaxis)*uaxis
xperp = x - xparallel
@test_approx_eq tfm1*(tfm2*x) xparallel-xperp
@test_approx_eq (tfm1*tfm2)*x xparallel-xperp
@test_approx_eq tfm1\(tfm2*x)  x
@test_approx_eq (tfm1\tfm2)*x  x
@test_approx_eq (tfm1*tfm1*tfm1*tfm1)*x  x
r = AffineTransforms.tformfwd(tfm1, x...)
@test [r...] == tfm1*x
r = AffineTransforms.tforminv(tfm1, x...)
@test [r...] == tfm1\x

# random transformations
for nd = 1:4
    A = eye(nd) + 0.1*randn(nd,nd)
    b = randn(nd)
    tfm = AffineTransforms.AffineTransform(A, b)
    x = rand(nd)
    @test_approx_eq tfm*x A*x+b
    @test_approx_eq tfm\x A\(x-b)
    @test_approx_eq (tfm\tfm)*x x
    # Scaling transformations
    s = rand(nd)
    y = tfm*x
    @test_approx_eq scale(s, tfm)*x  s.*y
    @test_approx_eq scale(tfm, s)*x  tfm*(s.*x)
    if 2 <= nd <= 3
        r = AffineTransforms.tformfwd(tfm, x...)
        @test [r...] == tfm*x
        r = AffineTransforms.tforminv(tfm, x...)
        @test [r...] == tfm\x
    end
end

# Pure rotations in 2d and 3d
for i = 1:20
    th = pi*(rand()-0.5)
    R = AffineTransforms.rotation2(th)
    @test_approx_eq th AffineTransforms.rotationparameters(R)
end
for i = 1:100
    ax = randn(3)
    n = norm(ax)
    if n == 0
        continue
    end
    ax = ax / n
    th = pi*rand()
    rv = th*ax
    R = AffineTransforms.rotation3(rv)
    @test_approx_eq rv AffineTransforms.rotationparameters(R)
end

# Transformation of arrays
using Interpolations, Images # Images for padarray
import AffineTransforms: center

padwith0(A) = padarray(A, [2,2], [2,2], "value", 0)

for (IT,GT) in ((BSpline(Constant()), OnCell()),
                (BSpline(Linear()), OnGrid()),
                (BSpline(Quadratic(Flat())), OnCell()))
    A = padwith0([1 1; 1 2])
    itp = interpolate(A, IT, GT)
    tfm = AffineTransforms.tformtranslate([1,0])
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [1 2; 0 0]
    @test_approx_eq tA[3,3] 1
    @test_approx_eq tA[3,4] 2
    tfm = AffineTransforms.tformtranslate([0,1])
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [1 0; 2 0]
    tfm = AffineTransforms.tformtranslate([-1,0])
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [0 0; 1 1]
    tfm = AffineTransforms.tformtranslate([0,-1])
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [0 1; 0 1]

    A = padwith0([2 1; 1 1])
    itp = interpolate(A, IT, GT)
    tfm = AffineTransforms.tformrotate(pi/2)
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [1 2; 1 1]
    tfm = AffineTransforms.tformrotate(-pi/2)
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test_approx_eq AffineTransforms.transform(tA)[3:4,3:4] [1 1; 2 1]

    # Check that getindex and transform yield the same results
    A = rand(6,6)
    itp = interpolate(A, IT, GT)
    tfm = AffineTransforms.tformrotate(pi/7)*AffineTransforms.tformtranslate(rand(2))
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    dest = AffineTransforms.transform(tA)
    tfm_recentered = AffineTransforms.AffineTransform(tfm.scalefwd, tfm.offset + center(A) - tfm.scalefwd*center(dest))
    tA_recentered = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm_recentered)
    nbad = 0
    for j = 1:size(dest,2), i = 1:size(dest,1)
        val = tA_recentered[i,j]
        # try
        #     @test isfinite(dest[i,j]) == isfinite(val)
        # catch err
        #     @show dest tA_recentered
        #     rethrow(err)
        # end
        nbad += isfinite(dest[i,j]) != isfinite(val)
        if isfinite(val) && isfinite(dest[i,j])
            @test abs(val-dest[i,j]) < 1e-12
        end
    end
    @test nbad < 3
    dest = AffineTransforms.transform(tA, origin_src=zeros(2), origin_dest=zeros(2))
    nbad = 0
    for j = 1:size(dest,2), i = 1:size(dest,1)
        val = tA[i,j]
        # try
        #     @test isfinite(dest[i,j]) == isfinite(val)
        # catch err
        #     @show dest tA
        #     rethrow(err)
        # end
        nbad += isfinite(dest[i,j]) != isfinite(val)
        if isfinite(val) && isfinite(dest[i,j])
            @test abs(val-dest[i,j]) < 1e-12
        end
    end
    @test nbad < 3
end

A = Float64[1 2 3 4; 5 6 7 8]
tfm = AffineTransforms.tformeye(2)
dest = zeros(2,2)

itp = interpolate(A, BSpline(Linear()), OnGrid())
tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
@test_approx_eq AffineTransforms.transform!(dest, tA) [2 3; 6 7]
dest3 = zeros(3,3)
@test_approx_eq AffineTransforms.transform!(dest3, tA) [NaN NaN NaN; 3.5 4.5 5.5; NaN NaN NaN]

itp = interpolate(A, BSpline(Constant()), OnCell())
tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
@test AffineTransforms.transform!(dest, tA) == [2 3; 6 7]

# Transforms with non-real numbers
using DualNumbers
a = AffineTransforms.tformtranslate([0,dual(1,0)])
A = reshape(1:9, 3, 3)
At = AffineTransforms.transform(A, a)
@test At[:,1:2] == A[:,2:3]


text:
@testset "Mappers" begin
  @testset "SimpleMapper" begin
    @testset "RegularGrid" begin
      grid1D = RegularGrid{Float64}(100)
      mappings = map(data1D, grid1D, [:value], SimpleMapper())
      @test mappings[:value] == Dict(100=>11,81=>9,11=>2,21=>3,91=>10,51=>6,61=>7,71=>8,31=>4,41=>5,1=>1)

      grid2D = RegularGrid{Float64}(100,100)
      mappings = map(data2D, grid2D, [:value], SimpleMapper())
      @test mappings[:value] == Dict(5076=>3,2526=>1,7551=>2)
    end

    @testset "PointSet" begin
      ps2D = PointSet([25. 50. 75.; 25. 75. 50.])
      mappings = map(data2D, ps2D, [:value], SimpleMapper())
      @test mappings[:value] == Dict(2=>2,3=>3,1=>1)
    end
  end
end


text:
module FESpacesTests

using Test
using Gridap

model = CartesianDiscreteModel(partition=(2,2))

order = 1
tags = [1,2,3,4]
fespace = H1ConformingFESpace(Float64,model,order,tags)

@test string(fespace) == "ConformingFESpace object"
s = "ConformingFESpace object:\n physdim: 2\n refdim: 2\n valuetype: Float64\n nfree: 5\n ndiri: 4"
@test sprint(show,"text/plain",fespace) == s

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

cellbasis = CellBasis(fespace)
a(v,u) = varinner(v,u)
bfun(x) = x[2]
b(v) = varinner(v,CellField(trian,bfun))
mmat = integrate(a(cellbasis,cellbasis),trian,quad)
bvec = integrate(b(cellbasis),trian,quad)

fun(x) = x[1] + x[2]

test_fe_space(fespace,5,4,mmat,bvec,fun)

V = TestFESpace(fespace)
test_fe_space(V,5,4,mmat,bvec,fun)

U = TrialFESpace(fespace,fun)
test_fe_space(U,5,4,mmat,bvec,fun)

fun(x) = x[1] + x[2]
z = 2.0
U = TrialFESpace(fespace,z)
test_fe_space(U,5,4,mmat,bvec,fun)

U = TrialFESpace(fespace,[z,z,z,z])
test_fe_space(U,5,4,mmat,bvec,fun)

end # module


text:
# https://github.com/oxinabox/LoggingExtras.jl

using Logging
using LoggingExtras
using Dates

function learn_ActiveFilteredLogger()
    filtered_logger = ActiveFilteredLogger(
        function(log_args)
            #@show log_args
            # log_args = (
            #     level = Info, message = "Boring message",
            #    _module = Main, group = :learn_LoggingExtras,
            #    id = :Main_1a9832fc,
            #    file = "c:\\Users\\ykyan\\work\\org.allnix.julia\\learn\\learn_LoggingExtras.jl",
            #    line = 13, kwargs = Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}()
            # )
            startswith(log_args.message, "Yo Dawg!")
        end,
        global_logger()
    )

    with_logger(filtered_logger) do
        @info "Boring message"
        @warn "Yo Dawg! it is bad"
        @info "Another boring message"
        @info "Yo Dawg! it is all good"
    end
end

function make_throttled_logger(period)
    history = Dict{Symbol, DateTime}()

    logger = EarlyFilteredLogger(global_logger()) do log
        # https://docs.julialang.org/en/v1/stdlib/Logging/
        # message id: fixed identifier for the source code statement where the logging macro appears
        # @show log
        # log = (
        #     level = Info, _module = Main, group = :learn_LoggingExtras,
        #     id = :Main_cc652dec
        # )
        if !haskey(history, log.id)
            history[log.id] = now()
            return true
        elseif now() - history[log.id] > period # has been long enough since last logging
            history[log.id] = now()
            return true
        else
            return false
        end
    end

    return logger
end

function learn_EarlyFilteredLogger()
    logger = make_throttled_logger(Second(3))
    with_logger(logger) do
        for i in 1:10
            sleep(1)
            @info "It happened" i
        end
    end
end

function learn_MinLevelLogger()
    @debug "Debug"
    @info "Info"

    logger = ConsoleLogger(stdout, Logging.Debug)
    debug_logger = MinLevelLogger(logger, Logging.Debug)
    with_logger(debug_logger) do
        @debug "Debug Logger"
    end
    error_logger = MinLevelLogger(logger, Logging.Error)
    with_logger(error_logger) do
        @debug "Debug no show"
        @error "Error Logger"
    end
end

function learn_TransformerLogger()
    truncating_logger = TransformerLogger(global_logger()) do log
        # @show log
        # log = (
        #     level = Info, message = "The truncating logger only truncates long messages",
        #     _module = Main, group = :learn_LoggingExtras,
        #     id = :Main_31932d27,
        #     file = "C:\\Users\\ykyan\\work\\org.allnix.julia\\learn\\learn_LoggingExtras.jl",
        #     line = 85, kwargs = Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}()
        # )

        if length(log.message) > 50
            short_message = log.message[1:(50-3)] * "..."
            return merge(log, (;message=short_message))
        else
            return log
        end
    end

    with_logger(truncating_logger) do
        @info "The truncating logger only truncates long messages"
        @info "Like this one that is this is a long and rambling message, it just keeps going and going and going,  and it seems like it will never end."
        @info "Not like this one, that is is short"
    end
end

function learn_FileLogger()
    demux_logger = TeeLogger(
        MinLevelLogger(FileLogger("info.log"), Logging.Info),
        MinLevelLogger(FileLogger("warn.log"), Logging.Warn),
    )

    with_logger(demux_logger) do
        @warn "Warn-1"
        @info "Info-1"
        @error "Error-1"
        @debug "Debug-1"
    end
end

function learn_FormatLogger()
    logger = FormatLogger() do io, log
        # @show log
        # log = (
        #     level = Info, message = "Info-1",
        #     _module = Main, group = :learn_LoggingExtras,
        #     id = :Main_1a175f96,
        #     file = "C:\\Users\\ykyan\\work\\org.allnix.julia\\learn\\learn_LoggingExtras.jl",
        #     line = 128, kwargs = Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}()
        # )
        println(io, "$(log._module) | [$(log.level)] $(log.message)")
    end

    with_logger(logger) do
        @info "Info-1"
        @warn "Warn-1"
    end
end

function learn_timestamp()
    date_format = "yyyy-mm-dd HH:MM:SS"

    timestamp_logger(logger) = TransformerLogger(logger) do log
        date_str = Dates.format(now(), date_format)
        merge(log, (; message = "$date_str $(log.message)"))
    end

    logger = ConsoleLogger(stdout, Logging.Debug) |> timestamp_logger

    with_logger(logger) do
        @debug "Time stampped logger"
    end
end

function my_logging()
    # const date_format = "yyyy-mm-dd HH:MM:SS" # if global
    date_format = "yyyy-mm-dd HH:MM:SS.sss"

    timestamp_logger(logger) = TransformerLogger(logger) do log
        date_str = Dates.format(now(), date_format)
        merge(log, (; message = "$date_str $(log.message)"))
    end

    logger = TeeLogger(
        ConsoleLogger(stdout, Logging.Warn) |> timestamp_logger,
        MinLevelLogger(FileLogger("allnix.log", append=false, always_flush=true), Logging.Info) |> timestamp_logger,
    )
    global_logger(logger)

    @info "File only"
    @warn "Console and file"
end

# Simple, clean logging with color coded "Info"
function learn_console_logger()
    # Save a copy
    default_logger = global_logger()
    @info "Default global logger"

    cout_logger = ConsoleLogger(stdout, Logging.Info)
    global_logger(cout_logger)
    @info "cout_logger"

    # Restore
    global_logger(default_logger)
end

# SimpleLogger prints the line number.
function learn_simple_logger()
    # Save a copy
    default_logger = global_logger()
    @info "Default global logger"

    logger = SimpleLogger(stdout, Logging.Info)
    global_logger(logger)
    @info "SimpleLogger"

    # Restore
    global_logger(default_logger)
end

#learn_ActiveFilteredLogger()
#learn_EarlyFilteredLogger()
#learn_MinLevelLogger()
#learn_TransformerLogger()
#learn_FileLogger()
#learn_FormatLogger()
#learn_timestamp()
#my_logging()
#learn_console_logger()
learn_simple_logger()

nothing


text:
module OMETIFF

using AxisArrays
using DataStructures
using DocStringExtensions
using EzXML
using FileIO
using ImageMetadata
using JSON
using ProgressMeter
using TiffImages: AbstractDenseTIFF, TiffFile, Iterable, IFD,
                  IMAGEDESCRIPTION, load!, getcache, getstream, fixcolors
using Unitful
using UUIDs

include("utils.jl")
include("files.jl")
include("parsing.jl")
include("mmap.jl")
include("loader.jl")

end # module


text:
using JuMP, EAGO, GAMS

function baron_factory()
    m = GAMS.Optimizer(GAMS.GAMSWorkspace("C:\\GAMS\\37"))

    m.gams_options["nlp"] = "BARON"
    m.gams_options["optca"] = 1E-4
    m.gams_options["optcr"] = 1E-4

    return m
end
                     m = Model(baron_factory)

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -1.9515287248442508 <= q <= 2.66866690087751)

                     add_NL_constraint(m, :((0.7910255593621245 + -0.9427161680575309*(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3]))/(1 + exp(-(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])))) + 0.08805479609820388*(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))/(1 + exp(-(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3])))))/(1 + exp(-(0.7910255593621245 + -0.9427161680575309*(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3]))/(1 + exp(-(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])))) + 0.08805479609820388*(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))/(1 + exp(-(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))))))) + (-0.7742678604770994 + 0.3023856535046492*(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3]))/(1 + exp(-(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])))) + -0.5348059615310992*(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))/(1 + exp(-(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3])))))/(1 + exp(-(-0.7742678604770994 + 0.3023856535046492*(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3]))/(1 + exp(-(-0.9881665794300214 + 0.6330293110094045*$(x[1]) + 0.8786704770929474*$(x[2]) + 0.11417841288574326*$(x[3])))) + -0.5348059615310992*(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))/(1 + exp(-(0.6512390960566696 + -0.0007920608466527312*$(x[1]) + -0.28084831326063275*$(x[2]) + 0.17707194589795927*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)
                     optimize!(m)
                     @show objective_value(m)
                     @show solve_time(m)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
emit_optimized(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false,
         dump_options=[Brutus.DumpCanonicalized])

emit_lowered(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false, # TODO: change to true when ready
         dump_options=[Brutus.DumpLoweredToStd,
                       Brutus.DumpLoweredToLLVM])

emit_translated(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false,
         dump_options=[Brutus.DumpTranslated])

emit_llvm(f, tt...) =
    emit(typeof(f), tt,
         emit_fptr=false, # TODO: change to true when ready
         dump_options=[Brutus.DumpLoweredToLLVM,
                       Brutus.DumpTranslateToLLVM])

function lit(emit_sym)
    @assert length(ARGS) == 1
    temp = Module(:Main)
    expr = quote
            include(fname) = Core.include(@__MODULE__, fname)
            import Brutus
            emit(args...;kwargs...) = Brutus.$(emit_sym)(args...;kwargs...)
            include($(ARGS[1]))
    end
    Base.eval(temp, Expr(:toplevel, expr))
end

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
#========== the macro! ==========#

"""
    @tullio C[i,k] := A[i,j] * B[j,k]
    @tullio F[i,k] := \$α * D[i].field[j] * E[col=k, row=j] + \$β

This is a replacement for `@einsum` which understands a bit more syntax.
The expression on the right is summed over all possible valued of the free index `k`,
and `:=` makes a new array `C`, while `=` and `+=` would write into an existing one.
Scalar arguments should have a dollar sign, like `\$α` or `A[i,\$γ]`.

    @tullio G[i,j] := M[i+x+1, j+y+1] * K[x,y]
    @tullio H[i,j] := M[2i+x, 2j+y]  (x in -1:1, y in -1:1)

Shifts and scaling of indices are allowed, including shifts by other indices.
Ranges can be provided as shown, for under-constrained indices.
If they are over-constrained, shifted indices run over the intersection allowed by all constraints,
while un-shifted indices demand agreement between them (e.g. `axes(A,2) == axes(B,1)` above).

    @tullio (*) L[i] := A[J[k]+2, i] / B[k]^2

This is a product instead of a sum, which could also enabled by writing `L[i] *= ...` (in-place).
You can use any reduction function such as `@tullio (max) M[i,j] := ...`.
Indexing by `J[k]+2` here demands `issubset(J, axes(A,1) .- 2)`.

    @tullio N[j] := sqrt <| M[i,j]^2

Pipe operators `|>` and `<|` apply a function after the sum, here `N ≈ map(norm, eachcol(M))`.
Underscores create functions, e.g. `|> sqrt(_ / V[i])` where clearly `i` must not have been summed.

See the readme for further options.
"""
macro tullio(exs...)
    _tullio(exs...; mod=__module__)
end

function _tullio(exs...; mod=Main)

    opts, ranges, ex = parse_options(exs...)
    if isnothing(ex) # then we simply updated global settings
        return (verbose=_VERBOSE[], fastmath=_FASTMATH[], threads=_THREADS[], grad=_GRAD[], avx=_AVX[], cuda=_CUDA[], tensor=_TENSOR[])
    end

    if opts.tensor && opts.redfun == :+ && isdefined(mod, :TensorOperations) && opts.grad != :Dual
        res = try_tensor(ex, ranges, DotDict(; mod = mod, opts...,
            newarray = false, scalar = false,
            arrays = Symbol[], indices = [], scalars = Symbol[]))
        if res != nothing # then forward & backward both handled by try_tensor
            return Expr(:block, res...) |> esc
        end
    end

    store = DotDict(; mod = mod, opts...,
    # Reduction
        redind = Symbol[],
        init = nothing,
    # Everything writes into leftarray[leftraw...], sometimes with a generated name
        leftraw = [],
        leftind = Symbol[],    # vcat(leftind, redind) is the complete list of loop indices
        leftarray = nothing,
        leftscalar = nothing, # only defined for scalar reduction
        leftnames = Symbol[],  # for NamedDims
        zero = false,
        scalar = false,
        newarray = false,
        plusequals = false,
    # Whole RHS, without finaliser, plus things extracted:
        right = nothing,
        finaliser = nothing,
        rightind = Symbol[],
        sharedind = Symbol[], # indices appearing on every RHS array, safe for ∇thread
        unsafeleft = Symbol[], # k in A[J[k]] never written to by different threads
        unsaferight = Symbol[], # same for gradient
        arrays = Symbol[],
        scalars = Symbol[],
        cost = 1,
    # Index ranges: first save all known constraints
        constraints = Dict{Symbol,Vector}(), # :k => [:(axes(A,2)), :(axes(B,1))] etc.
        constraintpairs = Tuple[], # (:i, :j, entangled range_i, range_j) from A[i+j] etc.
        notfree = Symbol[], # indices assigned values i = clamp(j, 1,3) within RHS
        shiftedind = Symbol[],
        axisdefs = Expr[],
        padmodclamp = false,
    # Expressions:
        outpre = Expr[],  # preliminary steps
        outex = Expr[],   # the rest!
    )

    parse_input(ex, store)

    parse_ranges(ranges, store)

    index_ranges(store)

    output_array(store)

    ex = action_functions(store)

    opts.verbose > 1 && verboseprint(store)

    ex |> esc
end

#========== options, etc ==========#

OPTS = Dict(
    :verbose => Any[true, false, 2, 3],
    :fastmath => [true, false],
    :threads => Integer,
    :grad => [false, :Base, :Dual],
    :avx => Integer,
    :cuda => Integer,
    :tensor => [true, false],
    )

_VERBOSE = Ref{Any}(false)
_FASTMATH = Ref(true)
_THREADS = Ref{Any}(true)
_GRAD = Ref{Any}(:Base)
_AVX = Ref{Any}(true)
_CUDA = Ref{Any}(256)
_TENSOR = Ref(true)

function parse_options(exs...)
    opts = Dict{Symbol,Any}(
        :redfun => :+,
        :init => TYP, # this means "auto"
        :pad => TYP,
        :verbose => _VERBOSE[],
        :fastmath => _FASTMATH[],
        :threads => _THREADS[],
        :grad => _GRAD[],
        :avx => _AVX[],
        :cuda => _CUDA[],
        :tensor => _TENSOR[],
        )
    expr = nothing
    nograd = Symbol[]
    ranges = Tuple[]
    for ex in exs
        # Actual options:
        if isexpr(ex, :(=)) && haskey(OPTS, ex.args[1])
            checklegal(ex.args[1], ex.args[2])
            opts[ex.args[1]] = ex.args[2]

        # Init & pad keyword
        elseif isexpr(ex, :(=)) && ex.args[1] == :init
            opts[:init] = ex.args[2]
        elseif isexpr(ex, :(=)) && ex.args[1] == :pad
            opts[:pad] = ex.args[2]

        # Nograd keyword
        elseif isexpr(ex, :(=)) && ex.args[1] == :nograd
            if ex.args[2] isa Symbol
                push!(nograd, ex.args[2])
            elseif isexpr(ex.args[2], :tuple)
                append!(nograd, ex.args[2].args)
            else
                throw("this accepts nograd=A or nograd=(A,B,C)")
            end

        # Ranges specified outside:
        elseif isexpr(ex, :call) && ex.args[1] in [:in, :∈]
            push!(ranges, (ex.args[2], ex.args[3]))
        elseif isexpr(ex, :tuple) && isexpr(ex.args[1], :call) && ex.args[1].args[1] in [:in, :∈]
            for el in ex.args
                isexpr(el, :call) && el.args[1] in [:in, :∈] || throw("expected (i ∈ 1:3) but got $el")
                push!(ranges, (el.args[2], el.args[3]))
            end

        # Reduction function
        elseif ex isa Symbol
            opts[:redfun] = ex

        # The main course!
        elseif ex isa Expr
            isnothing(expr) || throw("too many expressions! recognised keywords are $(vcat(:nograd, keys(opts)...))")
            expr = ex
        else
            throw("not sure what to do with input $ex")
        end
    end
    if isnothing(expr) # if run with no expression, it updates global options
        _VERBOSE[] = opts[:verbose]
        _FASTMATH[] = opts[:fastmath]
        _THREADS[] = opts[:threads]
        _GRAD[] = opts[:grad]
        _AVX[] = opts[:avx]
        _CUDA[] = opts[:cuda]
        _TENSOR[] = opts[:tensor]
    end
    (redfun=opts[:redfun],
        initkeyword=opts[:init], # surely there is a tidier way...
        padkeyword=opts[:pad],
        verbose=opts[:verbose],
        fastmath=opts[:fastmath],
        threads=opts[:threads],
        grad=opts[:grad],
        avx=opts[:avx],
        cuda=opts[:cuda],
        tensor=opts[:tensor],
        nograd=nograd,
    ), ranges, expr
end

checklegal(opt, val) =
    if OPTS[opt] isa Vector
        val in OPTS[opt] || throw("keyword $opt accepts values [" * join(OPTS[opt], ", ") * "]")
    elseif val isa Expr || val isa Symbol
        # allows threads=64^3 to work
    elseif OPTS[opt] == Integer
        val isa Integer && val >= 0 || throw("keyword $opt accepts false or a positive integer")
    end

#========== symbols ==========#

# These only need not to clash with symbols in the input:
RHS, AXIS = :𝓇𝒽𝓈, :𝒶𝓍
ZED, TYP, ACC, KEEP, FINAL = :ℛ, :𝒯, :𝒜𝒸𝒸, :♻, :💀
EPS, DEL, EXPR = :𝜀, :𝛥, :ℰ𝓍
MAKE, ACT! = :ℳ𝒶𝓀ℯ, :𝒜𝒸𝓉!

# @gensym RHS MAKE ACT!
# @gensym AXIS ZED TYP ACC KEEP FINAL
# @gensym EPS DEL EXPR

SYMBOLS = [
    RHS, MAKE, ACT!, AXIS, ZED, TYP, ACC, KEEP, EPS, DEL, EXPR,
    Symbol(:∇, MAKE), Symbol(:∇, ACT!), Symbol(DEL, ZED), Symbol(AXIS, :i),
    ] # to test for leaks

#========== input parsing ==========#

function parse_input(expr, store)

    # Equals sign & friends:
    if @capture_(expr, left_ := right_ )
        store.newarray = true
    elseif @capture_(expr, left_ = right_ )
    elseif @capture_(expr, left_ += right_ )
        store.plusequals = true
        store.redfun == :+ || throw("can't use += with reduction $(store.redfun)")
    elseif @capture_(expr, left_ *= right_ )
        store.plusequals = true # slightly abusing the name of the flag!
        if store.redfun == :+ # default, then we change it?
            store.verbose>0 && @info "inferring reduction by *, because of lhs *= rhs"
            store.redfun = :*
        elseif store.redfun == :*
        else
            throw("can't use *= with reduction $(store.redfun)")
        end
    elseif @capture_(expr, left_ ^= right_ )
        store.redfun == :+ && throw("can't use ^= with reduction +, please use +=")
        store.redfun == :* && throw("can't use ^= with reduction *, please use *=")
        store.plusequals = true
    else
        throw("can't understand input, expected A[] := B[] (or with =, or +=, *=, ^=) got $expr")
    end

    # Left hand side:
    if @capture_(left, Z_[leftraw__] ) || @capture_(left, [leftraw__] )
    elseif left isa Symbol # complete reduction
        store.newarray = true
        store.scalar = true
        store.leftscalar = left # because store.leftarray will be the array
        leftraw = Any[1,] # the gradient still indexes a fake 1D array
        expr.head == :(+=) && push!(store.scalars, left)
    else
        throw("can't understand LHS, expected A[i,j,k], got $left")
    end
    leftraw2 = tidyleftraw(leftraw, store)
    store.leftind = filter(i -> i isa Symbol, leftraw2) # this gives correct outer loop order

    isnothing(Z) && !(store.newarray) && throw("can't write into an array whose name isn't given!")
    Zed = isnothing(Z) ? ZED : Z
    store.leftarray = Zed

    store.leftraw = finishleftraw(leftraw2, store)
    if store.newarray && !allunique(store.leftind)
        store.zero = true # making diagonals, etc.
    end
    if !(store.newarray)
        saveconstraints(Zed, leftraw, store, false) # this adds to leftind, e.g. A[2i+1] = ..., is that bad??
        store.plusequals && detectunsafe(left, store.unsafeleft, store) # A[J[k]] += is unsafe, A[J[k]] = is not.
    end

    # Right hand side
    detectunsafe(right, store.unsaferight, store)
    right2 = MacroTools_prewalk(rightwalk(store), right)

    if isexpr(right2, :call) && right2.args[1] in (:|>, :<|)
        if right2.args[1] == :|>
            store.finaliser = makefinaliser(right2.args[3], store)
            store.right = MacroTools_postwalk(dollarwalk(store), right2.args[2])
        elseif right.args[1] == :<|
            store.finaliser = makefinaliser(right2.args[2], store)
            store.right = MacroTools_postwalk(dollarwalk(store), right2.args[3])
        end
        if store.scalar
            throw("can't use a finaliser $(right2.args[1]) with scalar output")
        end
    else
        store.right = MacroTools_postwalk(dollarwalk(store), right2)
        store.finaliser = :identity
    end

    unique!(store.scalars)
    unique!(store.arrays)
    unique!(store.leftind)
    store.sharedind = unique!(setdiff(store.sharedind, store.notfree))
    store.rightind = unique!(setdiff(store.rightind, store.notfree))
    any(==(:_), vcat(store.leftind, store.rightind)) && throw("can't use _ as an index name")

    unique!(store.outpre) # kill mutiple assertions, and evaluate any f(A) only once

    if store.newarray && Zed in store.arrays
        throw("can't create a new array $Zed when this also appears on the right")
    end
end

rightwalk(store) = ex -> begin
        @nospecialize ex
        # First, this will detect any assignment before it is used:
        if isexpr(ex, :(=))
            if ex.args[1] isa Symbol
                push!(store.notfree, ex.args[1])
            elseif isexpr(ex.args[1], :tuple)
                for i in ex.args[1].args
                    i isa Symbol && push!(store.notfree, i)
                end
            end
        end
        isexpr(ex, :return) && throw("can't use return inside body")

        # Second, alter indexing expr. to pull out functions of arrays:
        @capture_(ex, A_[inds__]) || return ex

        if isnothing(arrayonly(A))
            Anew = Symbol(string("≪", A, "≫"))
            push!(store.outpre, :(local $Anew = $(dollarstrip(A))))
            A = Anew
        end

        # Third, save letter A, and what axes(A) says about indices:
        push!(store.arrays, arrayonly(A))
        inds3 = primeindices(inds)
        saveconstraints(A, inds3, store, true)

        # Finally, re-assemble with new A etc:
        return :($A[$(inds3...)])
    end

arrayonly(A::Symbol) = A   # this is for RHS(i,j,k, A,B,C)
arrayonly(A::Expr) =
    if @capture_(A, B_[inds__]) || @capture_(A, B_.field_)
        return arrayonly(B)
    end # returns nothing from :(f(A)), signal to pull function out.

saveconstraints(A, inds, store, right=true) = begin
    A1 = arrayfirst(A, store)
    is = Symbol[]
    foreach(enumerate(inds)) do (d,ex)
        is_const(ex) && return
        containsany(ex, store.notfree) && return
        axis_i = length(inds)==1 ? :($eachindex($A1)) : :($axes($A1,$d))
        ex_i, axis_i = padmodclamp_ind(ex, axis_i, store) # this may pad the axis, or may make it nothing
        range_i, i = range_expr_walk(axis_i, ex_i)
        if isnothing(axis_i) # because mod(i) or clamp(i+j). Do save index, don't save range.
            if i isa Symbol
                push!(is, i)
                ex_i isa Symbol || push!(store.shiftedind, i)
            elseif i isa Tuple
                push!(is, filter(!isnothing, collect(i))...)
                push!(store.shiftedind, filter(!isnothing, collect(i))...)
            end
        elseif i isa Symbol
            push!(is, i)
            ex_i isa Symbol || push!(store.shiftedind, i)
            v = get!(store.constraints, i, [])
            push!(v, dollarstrip(range_i))
        elseif i isa Tuple # from things like A[i+j]
            push!(is, filter(!isnothing, collect(i))...) # collect for Julia ⩽ 1.3
            push!(store.shiftedind, filter(!isnothing, collect(i))...)
            push!(store.constraintpairs, (i..., dollarstrip.(range_i)...))
        elseif isnothing(i) # from A[J[k]], but A[J[k]+i] goes via store.constraintpairs, I said.
            str = "extrema of index $ex must fit within $A1"
            # @show range_i axis_i # @tullio C[i,k] := B[J[i]+1,k] verbose=2 grad=false # comes here, wrong check
            push!(store.outpre, :($issubset($range_i, $axis_i) || $throw($str)))
        end
    end
    if right
        append!(store.rightind, is)
        if A1 in store.nograd # then don't care whether it sharesindices
        elseif isassigned(store.sharedind)
            shared = intersect(is, store.sharedind)
            empty!(store.sharedind)
            append!(store.sharedind, shared)
        else
            append!(store.sharedind, is)
        end
    else
        append!(store.leftind, is) # why can's this be the only path for store.leftind??
    end
    n = length(inds)
    if n==1
        str = "expected a 1-array $A1, or a tuple"
        push!(store.outpre, :( $A1 isa Tuple || $ndims($A1) == 1 || $throw($str) ))
    else
        str = "expected a $n-array $A1" # already arrayfirst(A)
        push!(store.outpre, :( $ndims($A1) == $n || $throw($str) ))
    end
end

arrayfirst(A::Symbol, store) = A  # this is for axes(A,d), axes(first(B),d), etc.
arrayfirst(A::Expr, store) =
    if (@capture_(A, Binds_.field_) && @capture_(Binds, B_[inds__]))
        str = "elements $A must be of uniform size"
        push!(store.outpre, :( $all($ZED -> $axes($ZED.$field) == $axes($first($B).$field), $B) || throw($str) ))
        return :( $first($B).$field )
    elseif @capture_(A, B_[inds__])
        str = "elements $A must be of uniform size"
        push!(store.outpre, :( $all($AXIS -> $axes($AXIS) == $axes($first($B)), $B) || $throw($str) ))
        return :( first($B) )
    elseif @capture_(A, B_.field_)
        return A
    end

containsany(ex, list) = begin
    out = false
    MacroTools_postwalk(ex) do x
        if x in list
            out = true
        end
        x
    end
    out
end

primeindices(inds) = map(inds) do ex
    isexpr(ex, Symbol("'")) &&
        return Symbol(ex.args[1], "′") # normalise i''
    ex
end

# This function is for range inference
padmodclamp_ind(i, ax_i, store) = i, ax_i
padmodclamp_ind(ex::Expr, ax_i, store) =
    if ex.head == :call && ex.args[1] in [:mod, :clamp, :pad] && length(ex.args) == 2
        store.padmodclamp = true
        return ex.args[2], nothing # nothing means that range inference is discarded

    elseif ex.head == :call && ex.args[1] == :pad && length(ex.args) == 3
        store.padmodclamp = true
        _, a, p = ex.args
        return ex.args[2], :($padrange($ax_i, $p, $p)) # padrange() is in shifts.jl
    elseif ex.head == :call && ex.args[1] == :pad && length(ex.args) == 4
        store.padmodclamp = true
        _, a, lo, hi = ex.args
        return ex.args[2], :($padrange($ax_i, $lo, $hi))
    else
        return ex, ax_i
    end

padmodclamp_replace(s, store, inside=false) = s
padmodclamp_replace(ex::Expr, store, inside=false) =
    if ex.head == :(=) && @capture_(ex.args[1], A_[inds__])
        # This tricky case is 𝛥A[pad(i,2)] = 𝛥A[pad(i,2)] + ...
        Aex, fun = padmodclamp_pair(A, inds, store)
        right = if fun != identity
            padmodclamp_replace(ex.args[2], store, true)
        else
            padmodclamp_replace(ex.args[2], store, inside)
        end
        return fun(:($Aex = $right))
    elseif @capture_(ex, A_[inds__])
        Aex, fun = padmodclamp_pair(A, inds, store)
        return inside ? Aex : fun(Aex)
    else
        args = map(x -> padmodclamp_replace(x, store, inside), ex.args)
        Expr(ex.head, args...)
    end

padmodclamp_pair(A, inds, store) = begin
    nopadif = []
    inds4 = map(enumerate(inds)) do (d,ex)
        isexpr(ex, :call) || return ex
        if ex.args[1] == :mod && length(ex.args) == 2
            i = ex.args[2]
            return :($mod($i, $axes($A,$d)))
        elseif ex.args[1] == :clamp && length(ex.args) == 2
            i = ex.args[2]
            return :($clamp($i, $first($axes($A,$d)), $last($axes($A,$d))))
        elseif ex.args[1] == :pad && length(ex.args) >= 2
            i = ex.args[2]
            if !all(==(0), ex.args[3:end]) || length(ex.args) == 2
                push!(nopadif, :($i ∈ $axes($A,$d)))
            end
            return i
        end
        ex
    end
    Aex = :($A[$(inds4...)])
    fun = if isempty(nopadif)
        identity
    else
        cond = first(nopadif)
        for c2 in nopadif[2:end]
            cond = :($cond & $c2)
        end
        if store.padkeyword == TYP # default
            ex -> :($cond ? $ex : $zero($eltype($A)))
        else
            ex -> :($cond ? $ex : $convert($eltype($A), $(store.padkeyword)))
        end
    end
    Aex, fun # fun(Aex), but also fun(Aex = ...)
end

dollarwalk(store) = ex -> begin
        if isexpr(ex, :call)
            callcost(ex.args[1], store) # cost model for threading
        elseif isexpr(ex, :$) # interpolation of $c things:
            ex.args[1] isa Symbol || throw("you can only interpolate single symbols, not $ex")
            push!(store.scalars, ex.args[1])
            return ex.args[1]
        end
        ex
    end

dollarstrip(expr) = MacroTools_postwalk(expr) do ex
        isexpr(ex, :$) && return ex.args[1]
        ex
    end

tidyleftraw(leftraw, store) = begin
    step1 = map(leftraw) do i
        if isexpr(i, :kw) && store.newarray # then NamedDims wrapper is put on later
                push!(store.leftnames, i.args[1])
                return i.args[2]
        elseif i === :_ # underscores on left denote trivial dimensions
            return 1
        end
        i
    end
    primeindices(step1) # normalise i' to i′
end

finishleftraw(leftraw, store) = map(enumerate(leftraw)) do (d,i)
    if isexpr(i, :$)
        store.newarray && throw("can't fix indices on LHS when making a new array")
        i.args[1] isa Symbol || throw("you can only interpolate single symbols, not $ex")
        push!(store.scalars, i.args[1])
        return i.args[1]

    elseif isexpr(i, :call) && i.args[1] == :+ &&
            length(i.args)==3 && i.args[3] == :_ # magic un-shift A[i+_, j] := ...
        i = primeindices(i.args)[2]
        i isa Symbol || throw("index ($i + _) is too complicated, sorry")
        push!(store.leftind, i)
        deli = Symbol(DEL, i)
        push!(store.scalars, deli) # calculating this must be done later
        return :($i + $deli)

    elseif @capture_(i, J_[inds__]) # scatter operation, A[i,J[j,k]] := ...
        push!(store.nograd, J)
        rightwalk(store)(i) # array J viewed as part of RHS, and provides a range for j,k
        inds2 = filter(j->j isa Symbol, tidyleftraw(inds, store))
        append!(store.leftind, inds2) # but j,k aren't to be summed

        ex = :($J[$(tidyleftraw(inds, store)...)])
        if store.newarray
            ax_i = Symbol(AXIS, string("≪", ex, "≫")) # fake index name, to which to attach a size?
            push!(store.axisdefs, :(local $ax_i = $extremerange($J)))
            store.zero = true
        end

        return ex # has primes dealt with
    end
    i
end

detectunsafe(expr, list, store) = MacroTools_postwalk(expr) do ex
        @capture_(ex, A_[inds__]) || return ex
        for i in inds
            MacroTools_postwalk(i) do x
                @capture_(x, B_[inner__]) || return x
                # Now we have found an array which indexes another one, mark its indices unsafe
                append!(list, filter(j -> j isa Symbol, inner))
                unique!(list)
                # and don't compute a gradient for the inner array
                B isa Symbol && push!(store.nograd, B)
                x
            end
        end
        ex
    end

makefinaliser(s::Symbol, store) = s
makefinaliser(expr::Expr, store) = begin
    underscore = false
    out = MacroTools_postwalk(expr) do ex
        if ex == :_
            underscore = true
            return RHS
        elseif @capture_(ex, A_[inds__])
            for i in inds
                i isa Symbol || continue
                i in store.leftind || throw("index $i can't be used in finaliser")
            end
        end
        ex
    end
    if underscore
        return dollarstrip(:($RHS -> $out))
    else
        return dollarstrip(ex)
    end
end

function parse_ranges(ranges, store) # now runs after parse_input
    for (i,r) in ranges
        if isexpr(i, Symbol("'")) # catch primes!
            i = Symbol(i.args[1], "′")
        end
        push!(store.rightind, i)
        v = get!(store.constraints, i, [])
        if isexpr(r, :call) && r.args[1] == :(:) && length(r.args) == 3
            # for a literal range, write OneTo(10) or 0:9 directly into constraints
            if r.args[2] == 1 && r.args[3] isa Integer
                push!(v, :(Base.OneTo($(r.args[3]))))
                continue
            elseif r.args[2] isa Integer && r.args[3] isa Integer
                push!(v, r)
                continue
            end
        end
        # for axes(A,2) where A is already available, just save it
        if isexpr(r, :call) && r.args[1] in (:axes, :eachindex) && r.args[2] in store.arrays
            push!(v, r)
            continue
        end
        # for anything else, treat it as a scalar argument
        if r isa Symbol
            push!(store.scalars, r)
            push!(v, r)
        else
            s = Symbol(string("≪", r, "≫"))
            push!(store.outpre, :(local $s = $r))
            str = "expected a range for ($i in $r), got "
            push!(store.outpre, :($s isa AbstractRange || throw($str * string($r))))
            push!(store.scalars, s)
            push!(v, s)
        end
    end
    unique!(store.rightind)
    unique!(store.scalars)
    store.redind = setdiff(store.rightind, store.leftind)
end

#========== index ranges ==========#

function index_ranges(store)

    todo = Set(vcat(store.leftind, store.redind))
    done = Dict{Symbol,Any}()

    for (i,j,r_i,r_j) in store.constraintpairs

        if isnothing(i) # case of A[j + I[k]]
            v = get!(store.constraints, j, [])
            push!(v, r_j)
        elseif isnothing(j)
            v = get!(store.constraints, i, [])
            push!(v, r_i)

        elseif haskey(store.constraints, i) && i in todo
            resolveintersect(i, store, done) # use existing knowledge to fix i's range
            pop!(todo, i)
            v = get!(store.constraints, j, []) # and then allow j's range to depend on that
            push!(v, r_j)
        elseif haskey(store.constraints, j) && j in todo
            resolveintersect(j, store, done)
            pop!(todo, j)
            v = get!(store.constraints, i, [])
            push!(v, r_i)
        end
    end

    for i in todo
        haskey(store.constraints, i) || throw("unable to infer range of index $i")
        if i in store.shiftedind
            resolveintersect(i, store, done)
        else
            resolvestrict(i, store, done)
        end
        deli = Symbol(DEL,i)
        if deli in store.scalars # magic shift on LHS
            axi, axi_del = Symbol(AXIS,i), Symbol(AXIS,i,DEL)
            push!(store.axisdefs, :($deli = 1 - $first($axi)),
                :(local $axi_del = Base.OneTo($length($axi))))
            # You can't compute deli inside Act! as doesn't always see full range of i.
            # But if you make it a scalar argument, then it's an argument of Make, hence
            push!(store.outpre, :(local $deli = 0)) # ... this awful hack.
        end
    end

    append!(store.outex, store.axisdefs)

    if store.verbose > 0
        if !isempty(store.leftind)
            lex = map(i -> Expr(:(=), i, done[i]), store.leftind)
            push!(store.outex, :(@info "left index ranges" $(lex...)))
        end
        if !isempty(store.redind)
            rex = map(i -> Expr(:(=), i, done[i]), store.redind)
            push!(store.outex, :(@info "reduction index ranges" $(rex...)))
        end
    end
end

resolvestrict(i, store, done) = begin
    res = first(store.constraints[i])
    ax_i = Symbol(AXIS, i)
    push!(store.axisdefs, :( local $ax_i = $res ))
    done[i] = res
    for alt in store.constraints[i][2:end] # in which case it shouldn't be a Set
        str = "range of index $i must agree"
        push!(store.axisdefs, :( $alt == $res || $throw($str) ))
    end
end

resolveintersect(i, store, done) = begin
    res = if isempty(store.constraints[i])
        throw("unable to infer range of index $i")
    elseif length(store.constraints[i])==1
        first(store.constraints[i])  # because intersect(1:3) isa Vector, wtf?
    else
        :( intersect($(store.constraints[i]...)) )
    end
    ax_i = Symbol(AXIS, i)
    push!(store.axisdefs, :( local $ax_i = $res ))
    done[i] = res
end

#========== output array + eltype ==========#

function output_array(store)

    # Initialisation needs to be worked out somewhere...
    if store.initkeyword == TYP # then auto
        store.init = store.redfun == :+ ? :(zero($TYP)) :
                    store.redfun == :* ? :(one($TYP)) :
                    store.redfun == :max ? :(typemin($TYP)) :
                    store.redfun == :min ? :(typemax($TYP)) :
                    store.redfun == :& ? :(true) :
                    store.redfun == :| ? :(false) :
                    begin
                        store.verbose>0 && @warn "guessing init=zero(T) for unknown reduction $(store.redfun)"
                        :(zero($TYP))
                    end
    else
        if store.initkeyword isa Number
            store.init = store.initkeyword
        else
            init_sy = Symbol(string("≪", store.initkeyword, "≫"))
            push!(store.outpre, :(local $init_sy = $(store.initkeyword)))
            push!(store.scalars, init_sy)
            store.init = init_sy
        end
    end

    # And some not-compltely-unrelated errors:
    if isempty(store.redind) && !(store.plusequals)
        store.redfun == :+ || throw("nothing to reduce over using $(store.redfun)")
        store.finaliser == :identity || throw("can't apply finaliser without a reduction")
    end
    if isempty(store.redind)
        store.initkeyword == TYP || throw("nothing to reduce over, so won't use init = $(store.initkeyword)")
    elseif store.plusequals && !(store.scalar)
        store.initkeyword == TYP || throw("in-place update will not use init = $(store.initkeyword)")
    end

    if store.newarray # this includes scalar case!

        ex_right = padmodclamp_replace(:($(store.finaliser)($(store.right))), store)
        push!(store.outex, :( local $RHS($(store.arrays...), $(store.rightind...)) = $ex_right))

        # Try inference first, usually fine, and avoids scalar evaluation on GPU
        allfirst = map(i -> :($first($(Symbol(AXIS, i)))), store.rightind)
        T1 = Symbol(TYP,1)
        T2 = Symbol(TYP,2)
        T3 = Symbol(TYP,3)
        warn = store.verbose>0 ? :(@warn "unable to infer eltype from RHS") : nothing
        push!(store.outex, quote
            local $T1 = Core.Compiler.return_type($RHS, $typeof(($(store.arrays...), $(allfirst...))))
            local $T2 = if Base.isconcretetype($T1)
                $T1
            else
                $warn
                $typeof($RHS($(store.arrays...), $(allfirst...)))
            end
        end)

        # Init. usually depends on type, but sometimes widens type
        if store.initkeyword == TYP
            push!(store.outex, :(local $T3 = $T2))
        else
            push!(store.outex, :(local $T3 = Base.promote_type($T2, $typeof($(store.init)))))
        end

        # Oh, also scalar += might widen type...
        if store.scalar && store.plusequals
            push!(store.outex, :(local $TYP = Base.promote_type($T3, $typeof($(store.leftscalar)))))
        else
            push!(store.outex, :(local $TYP = $T3))
        end

        # This now checks for OffsetArrays, and allows A[i,1] := rhs. Pulls out scatterers.
        outaxes = map(store.leftraw) do i
            i isa Integer && i==1 && return :(Base.OneTo(1))
            i isa Symbol && return Symbol(AXIS, i)
            i isa Expr && @capture_(i, J_[inds__]) && return Symbol(AXIS, string("≪", i, "≫"))
            i isa Expr && i.head == :call && length(i.args)==3 && i.args[1] == :+ &&
                startswith(string(i.args[3]), string(DEL)) && return Symbol(AXIS, i.args[2], DEL)
            throw("can't use index $i on LHS for a new array")
        end

        if !isdefined(store.mod, :OffsetArrays)
            outaxes = map(store.leftraw, outaxes) do i, ax
                ax == :(Base.OneTo(1)) && return ax
                i in store.shiftedind || @capture_(i, J_[inds__]) || return ax
                push!(store.outex, :( $first($ax) == 1 || $throw("to allow indices not starting at 1, OffsetArrays must be visible in the caller's module. Otherwise write `A[i+_] := ...` to remove offset")))
                return :(Base.OneTo($ax)) # This doesn't apply to offsets caused by pad(i+j,3), sadly?
            end
        end

        simex = if store.scalar && store.plusequals
            :( convert($TYP, $(store.leftscalar)) ) # here init is needed only if threading
        elseif store.scalar
            :( convert($TYP, $(store.init)) )
        elseif isempty(store.arrays)
            :( similar(1:0, $TYP, tuple($(outaxes...))) )
        else
            # parent() is a trick to avoid a NamedDims bug
            :( similar(parent($(store.arrays[1])), $TYP, tuple($(outaxes...),)) )
        end
        if store.scalar
            push!(store.outex, :( local $ZED = $simex ))
        elseif isempty(store.leftnames)
            push!(store.outex, :( local $(store.leftarray) = $simex ))
        else
            nex = :(tuple($(QuoteNode.(store.leftnames)...)))
            push!(store.outex, :( local $(store.leftarray) = NamedDims.NamedDimsArray($simex, $nex) ))
        end

        if store.scalar && store.threads != false && store.initkeyword != TYP
            msg = "init=$(store.init) must be compatible with $(store.redfun), for possibly-threaded scalar reduction"
            push!(store.outex, :($(store.redfun)($(store.init), $(store.init)) == $(store.init) || $throw($msg)))
        end
    end

    if store.zero # allow pad=NaN to control this too
        # push!(store.outex, :( $(store.leftarray) .= false )) # zero($TYP) won't work in-place
        if store.padkeyword == TYP # default
            push!(store.outex, :($(store.leftarray) .= zero(eltype($(store.leftarray)))))
        else
            push!(store.outex, :($(store.leftarray) .= $(store.padkeyword)))
        end
    end

    ex_pre = quote $(store.outpre...) end # before act! gets pushed into store.outpre
    store.verbose==2 && @info ">>>>> Preliminary expressions" verbosetidy(ex_pre)
end

#========== action functions ==========#

function action_functions(store)

    axisleft = map(i -> Symbol(AXIS, i), setdiff(store.leftind, store.unsafeleft))
    axisred = map(i -> Symbol(AXIS, i), union(store.redind, store.unsafeleft))
    axislist = vcat(axisleft, axisred)
    # Order of these is convenient for threader(), which divides axisleft up freely,
    # divides axisred up with re-starts.
    # This is independent of the grouping inner/outer for make_many_actors().

    #===== constructing loops =====#

    zed_arg, zed_one = if store.scalar
        :($ZED::$TYP), ZED
    else
        :($ZED::AbstractArray{$TYP}), :($ZED[$(store.leftraw...)])
    end

    ex_init = if store.plusequals && !isempty(axisleft) # then always keep=true
        :( $ACC = $zed_one )
    elseif store.scalar && !(store.plusequals) # then always keep=false
        :( $ACC = $(store.init) )
    else # for non-numbers, similar() may leave undef, so avoid ifelse here
        :( $ACC = $KEEP===nothing ? $(store.init) : $zed_one )
    end

    ex_iter = :( $ACC = $(store.redfun)($ACC, $(store.right) ) )

    ex_write = if store.scalar # then we return the value instead, ZED is immutable
        :( $ACC )
    elseif store.finaliser == :identity
        :( $ZED[$(store.leftraw...)] = $ACC )
    else # this branch is moved outside @avx by finalsplit(expr), below.
        :( $ZED[$(store.leftraw...)] = $FINAL===nothing ? $ACC : $(store.finaliser)($ACC) )
    end

    ex_nored = if store.plusequals # implies keep=true directly, and final=true since no J indices in threader.
        :( $ZED[$(store.leftraw...)] =  $(store.finaliser)($(store.redfun)($ZED[$(store.leftraw...)] ,$(store.right))) )
    else # using finaliser without reduction, and without +=, is now an error.
        :( $ZED[$(store.leftraw...)] = $(store.right) )
    end

    if isempty(store.redind)
        make_many_actors(ACT!,
            vcat(zed_arg, store.arrays, store.scalars, axislist),
            nothing, store.leftind, nothing, Symbol[], ex_nored, nothing, store)
    else
        make_many_actors(ACT!,
            vcat(zed_arg, store.arrays, store.scalars, axislist),
            nothing, store.leftind, ex_init, store.redind, ex_iter, ex_write, store)
    end

    ∇make = if store.newarray
        # make_many_actors and backward_definitions both push into store.outpre
        backward_definitions(store)
    else
        nothing
    end

    #===== action! =====#

    keep = store.plusequals ? :true : :nothing
    block = store.threads==false ? nothing :
        store.threads==true ? cld(BLOCK[], store.cost) :
        store.threads
    if store.scalar
        ST = :($storage_type($(store.arrays...)))
        push!(store.outex, :(
            $thread_scalar($ACT!, $ST, $ZED,
                tuple($(store.arrays...), $(store.scalars...),),
                tuple($(axisred...),), $(store.redfun), $block, $keep)
        ))
    else
        ST = :($storage_type($(store.leftarray), $(store.arrays...)))
        push!(store.outex, :(
            $threader($ACT!, $ST, $(store.leftarray),
                tuple($(store.arrays...), $(store.scalars...),),
                tuple($(axisleft...),), tuple($(axisred...),), $(store.redfun), $block, $keep);
            $(store.leftarray)
        ))
    end
    store.verbose>0 && block != nothing && @info "threading threshold (from cost = $(store.cost))" block

    if store.newarray
        # then slurp up outex to make a function:
        ex_make = quote
            local @inline function $MAKE($(store.arrays...), $(store.scalars...), )
                $(store.outex...)
            end
        end
        store.verbose==2 && @info ">>>>> Maker function" verbosetidy(ex_make)
        ex = quote
            let $ACT! = $ACT!
                $ex_make
                $Eval($MAKE, $∇make)($(store.arrays...), $(store.scalars...), )
            end
        end

        # wrap pre and out in one let block so that ACT! doesn't escape:
        ex = :(let
            $(store.outpre...)
            $ex
        end)

        # and assign the result if necc:
        if store.leftarray != ZED
            push!(store.outex, :($(store.leftarray) = $ex ))
            return :($(store.leftarray) = $ex )
        elseif store.scalar
             push!(store.outex, :($(store.leftscalar) = $ex))
             return :($(store.leftscalar) = $ex)
        else # case of [i,j] := ... with no name given
            # push!(store.outex, ex)
            return ex
        end

    else
        # in-place, no MAKE function, but still keep ACT! from escaping
        ex_body = quote $(store.outex...) end
        store.verbose==2 && @info "In-place body" verbosetidy(ex_body)
        return :(let
            $(store.outpre...)
            $(store.outex...)
        end)
    end
end

"""
    make_many_actors(f!, args, ex1, [:i,], ex3, [:k,], ex5, ex6, store)

This makes several functions of this form,
decorated as necessary with `@inbouds` or `@avx` etc,
and with appropriate `storage_type` as the first argument.
```
f!(::Type, args..., keep=nothing, final=true) where {T}
    ex1
    ex2 = (for i in axis_i
        ex3
        ex4 = (for k in axis_k
            ex5
        end)
        ex6
    end)
end
```
"""
function make_many_actors(act!, args, ex1, outer::Vector, ex3, inner::Vector, ex5, ex6, store, note="")

    if store.padmodclamp
        ex5 = padmodclamp_replace(ex5, store)
    end

    ex4 = recurseloops(ex5, inner)
    ex2 = recurseloops(:($ex3; $ex4; $ex6), outer)

    ex_act = if store.fastmath && isempty(store.notfree)
        quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                @inbounds @fastmath ($ex1; $ex2)
            end
        end
    elseif isempty(store.notfree)
        quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                @inbounds ($ex1; $ex2)
            end
        end
    else
        quote
            local @inline function $act!(::Type, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                ($ex1; $ex2)
            end
        end
    end
    store.verbose==2 && @info "===== Base actor $note" verbosetidy(ex_act)
    push!(store.outpre, ex_act)

    if act! != ACT! && isempty(store.sharedind) && store.threads != false
        store.verbose>0 && @warn "can't parallelise this gradient, no shared indices $note"
    end

    #===== LoopVectorization =====#

    expre, exloop0, expost = if isempty(outer)
        :($ex1; $ex3), ex4, ex6
    else
        ex1, ex2, nothing
    end
    exloop, exloopfinal = finalsplit(exloop0)

    # Disable @avx for scatter, https://github.com/chriselrod/LoopVectorization.jl/issues/145
    safe = if act! == ACT!
        isempty(store.unsafeleft)
    else # working on ∇act!
        isempty(store.unsaferight)
    end

    if safe && store.avx != false && isdefined(store.mod, :LoopVectorization)
        unroll = store.avx == true ? 0 : store.avx # unroll=0 is the default setting
        info1 = store.verbose>0 ? :(@info "running LoopVectorization actor $($note)" maxlog=3 _id=$(hash(store))) : nothing
        check1 = store.verbose>0 ? :(LoopVectorization.check_args($(store.arrays...)) || @error "rejected by LoopVectorization's check_args! $($note)" maxlog=3 _id=$(hash(store))) : nothing
        try
            lex = if isnothing(exloopfinal)
                quote

                    local @inline function $act!(::Type{<:Array{<:Union{Base.HWReal, Bool}}}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                        $expre
                        $info1
                        $check1
                        LoopVectorization.@avx unroll=$unroll $exloop
                        $expost
                    end

                end
            else # "isnothing(final) ? exp(rhs) : rhs" does not prevent execution of finaliser within @avx
                quote

                    local @inline function $act!(::Type{<:Array{<:Union{Base.HWReal, Bool}}}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                        $expre
                        $info1
                        $check1
                        if $FINAL === nothing
                            LoopVectorization.@avx unroll=$unroll $exloop
                        else
                            LoopVectorization.@avx unroll=$unroll $exloopfinal
                        end
                        $expost
                    end

                end
            end
            store.verbose==2 && @info "=====LV===== LoopVectorization actor $note" verbosetidy(lex)
            push!(store.outpre, macroexpand(store.mod, lex))
            store.verbose==2 && @info "success expanding LoopVectorization.@avx"
        catch err
            store.verbose>0 && @warn "LoopVectorization failed $note" err
        end
    end

    #===== KernelAbstractions =====#

    unsafe = if act! == ACT!
        store.unsafeleft
    else # working on ∇act!
        store.unsaferight
    end
    safeouter = setdiff(outer, unsafe)

    if store.cuda > 0 && isdefined(store.mod, :KernelAbstractions)
        kernel = gensym(:🇨🇺)
        axouter = map(i -> Symbol(AXIS, i), safeouter)
        asserts = map(ax -> :( $first($ax)==1 || $throw("KernelAbstractions can't handle OffsetArrays here")), axouter)
        sizes = map(ax -> :(length($ax)), axouter)

        if isempty(safeouter)
            store.verbose>0 && @warn "using KernelAbstractions with no outer indices, this will be slow"
            safeouter = [Symbol(EPS, 1)] # fake index name, only appears in @index
            sizes = [:(one(Int))]    # iterate over 1:1
        end

        kernelbody = recurseloops(:($ex3; $ex4; $ex6), unsafe)
        try
            # const_args = map(args) do a
            #     a isa Symbol || return a  # this skips output ZED::AbstractArray{TYP}
            #     a == store.leftarray && return a  # case A[i] = A[i]^2 / B[i,j]
            #     :(@Const($a))
            # end
            kex1 = quote
                # @Const removed, see https://github.com/mcabbott/Tullio.jl/pull/32
                # KernelAbstractions.@kernel function $kernel($(const_args...), @Const($KEEP), @Const($FINAL)) where {$TYP}
                KernelAbstractions.@kernel function $kernel($(args...), $KEEP, $FINAL) where {$TYP}
                    ($(safeouter...),) = @index(Global, NTuple)
                    $ex1  # This seems dodgy, shouldn't ex1 be outside?
                    $kernelbody
                end

            end
            store.verbose==2 && @info "=====KA===== KernelAbstractions kernel $note" verbosetidy(kex1)
            push!(store.outpre, macroexpand(store.mod, kex1))
            if isdefined(store.mod, :CUDA) && isdefined(store.mod, :CuArray) # new-style, CUDA.jl, with CUDADevice()
                info2 = store.verbose>0 ? :(@info "running KernelAbstractions + CUDA actor $($note)" maxlog=3 _id=$(hash(store))) : nothing
                kex2 = quote

                    local @inline function $act!(::Type{<:CuArray}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                        $info2
                        cu_kern! = $kernel(CUDADevice(), $(store.cuda))
                        $(asserts...)
                        $ACC = cu_kern!($(args...), $KEEP, $FINAL; ndrange=tuple($(sizes...)), dependencies=Event(CUDADevice()))
                        KernelAbstractions.wait(CUDADevice(), $ACC)
                    end

                end
                store.verbose==2 && @info "=====KA===== KernelAbstractions CUDA actor $note" verbosetidy(kex2)
                push!(store.outpre, kex2)
            end
            info3 = store.verbose>0 ? :(@info "running KernelAbstractions CPU actor $($note)" maxlog=3 _id=$(hash(store))) : nothing
            kex3 = quote

                local @inline function $act!(::Type{<:Array}, $(args...), $KEEP=nothing, $FINAL=true) where {$TYP}
                    $info3
                    cpu_kern! = $kernel(CPU(), 4)
                    $(asserts...)
                    $ACC = cpu_kern!($(args...), $KEEP, $FINAL; ndrange=tuple($(sizes...)))
                    KernelAbstractions.wait($ACC)
                end

            end
            if store.threads==false
                # This CPU kernel can't be called by threader, and so threads=false
                # offers a way to control whether it gets used or not. By default, not.
                push!(store.outpre, kex3)
            end
            store.verbose==2 && @info "success expanding KernelAbstractions.@kernel"
        catch err
            store.verbose>0 && @warn "KernelAbstractions failed $note" err
        end
    end
end

recurseloops(ex, list::Vector) =
    if isempty(list)
        return ex
    else
        i = first(list)
        r = Symbol(AXIS, i)
        ex = :(for $i in $r; $ex; end)
        return recurseloops(ex, list[2:end])
    end

finalsplit(expr) = begin
    yes = false
    ex_1 = MacroTools_postwalk(expr) do ex
        yes |= isifelsefinal(ex)
        isifelsefinal(ex) ? ex.args[2] : ex
    end
    ex_2 = MacroTools_postwalk(expr) do ex
        isifelsefinal(ex) ? ex.args[3] : ex
    end
    if yes
        return ex_1, ex_2
    else
        return ex_1, nothing
    end
end

# This matches ex = :(isnothing(💀) ? 𝒜𝒸𝒸 : tanh(𝒜𝒸𝒸))
# and ex = :(💀===nothing ? 𝒜𝒸𝒸 : tanh(𝒜𝒸𝒸))
isifelsefinal(ex) = isexpr(ex, :if, 3) && isexpr(ex.args[1], :call) &&
        ex.args[1].args[1] in (:isnothing, :(===)) && ex.args[1].args[2] == FINAL

#===== define gradient hooks =====#

function backward_definitions(store)
    store.grad == false && return nothing # no gradient wanted

    axisshared = map(i -> Symbol(AXIS, i), setdiff(store.sharedind, store.unsaferight)) # safe to multi-thread
    loopind = vcat(store.leftind, store.redind)
    axisnonshared = map(i -> Symbol(AXIS, i), union(setdiff(loopind, store.sharedind), store.unsaferight))

    axislist = vcat(axisshared, axisnonshared) # this defines the order of arguments of ∇act!

    ok = false
    if store.grad == :Dual && store.redfun == :+
        try
            insert_forward_gradient(axislist, store)
            ok = true
            store.verbose==2 && @info "using ForwardDiff gradient"
        catch err
            store.verbose>0 && @warn "ForwardDiff gradient failed" err
        end
    elseif store.grad == :Base
        try
            insert_symbolic_gradient(axislist, store)
            ok = true
            store.verbose==2 && @info "success wtih symbolic gradient"
        catch err
            store.verbose>0 && @warn "symbolic gradient failed" err
        end
    end

    ok == false && return nothing # failed to make a gradient

    dZ = Symbol(DEL, ZED)
    ∇make = Symbol(:∇, MAKE)
    ∇act! = Symbol(:∇, ACT!)

    gradarrays = map(A -> Symbol(DEL, A), store.arrays)
    # gradscalars = map(A -> Symbol(DEL, A), store.scalars)
    defineempties = map(store.arrays, gradarrays) do A, dA
        if A in store.nograd
            :(local $dA = nothing)
        else
            :( local $dA = fill!(similar($A, Base.promote_type(eltype($A), $TYP)), 0) )
        end
    end
    # append!(defineempties, map((x,dx) -> :($dx = zero(Base.promote_type(typeof($x), $TYP))), store.scalars, gradscalars))
    returns = vcat(gradarrays, map(_->:nothing, store.scalars)) # ?? needs a test!
    # returns = vcat(gradarrays, gradscalars)

    ST = :($storage_type($(gradarrays...), $(store.arrays...)))
    block = store.threads==false ? nothing :
        store.threads==true ? cld(BLOCK[], store.cost) :
        store.threads
    input, acton = if store.scalar
        :($dZ::$TYP), :( $OneBox($dZ) ) # a hack to minimise changes to ∇Act!, for now??
    else
        :($dZ::AbstractArray{$TYP}), dZ
    end
    ex_make = quote

        local function $∇make($input, $ZED, $(store.arrays...), $(store.scalars...), ) where {$TYP}
            $(defineempties...)
            $(store.axisdefs...)
            $∇threader($∇act!, $ST,
                tuple($(gradarrays...), $acton, $ZED, $(store.arrays...), $(store.scalars...),),
                tuple($(axisshared...),), tuple($(axisnonshared...), ), $block)
            return ($(returns...),)
        end

    end
    store.verbose==2 && @info "<<<<< Gradient maker function" verbosetidy(ex_make)
    push!(store.outpre, quote
        local $∇make = let $∇act! = $∇act!
            $ex_make
        end
    end)

    return ∇make
end

fillarrayreplace(rhs, dZ) = MacroTools_postwalk(rhs) do @nospecialize ex
        @capture_(ex, A_[inds__]) && A==dZ || return ex
        return Symbol(dZ, :_value)
    end

#========== the end ==========#


text:
module TestNeighbors
using Test

include("../src/incidence.jl")

@testset "IncidenceIndex: Order          " begin
    idx1 = IncidenceIndex(3, [(3, 1), (3, 2)])
    idx2 = IncidenceIndex(3, [(1, 3), (2, 3)])
    @test idx1 == idx2
end

@testset "IncidenceIndex: Square         " begin
    #       1  2  3  4
    head = [1, 1, 2, 3]
    tail = [2, 3, 4, 4]
    idx = IncidenceIndex(4, head, tail)
    @test length(idx) == 4
    @test Set(idx[1]) == Set([(2, 1), (3, 2)])
    @test Set(idx[2]) == Set([(1, 1), (4, 3)])
    @test Set(idx[3]) == Set([(1, 2), (4, 4)])
    @test Set(idx[4]) == Set([(2, 3), (3, 4)])
end

@testset "IncidenceIndex: Square -1      " begin
    #       1  2  3  4
    head = [1, 1, 2, 0]
    tail = [2, 3, 4, 4]
    idx = IncidenceIndex(4, head, tail)
    @test length(idx) == 4
    @test Set(idx[1]) == Set([(2, 1), (3, 2)])
    @test Set(idx[2]) == Set([(1, 1), (4, 3)])
    @test Set(idx[3]) == Set([(1, 2)])
    @test Set(idx[4]) == Set([(2, 3)])
end

@testset "IncidenceIndex: Test_1         " begin
    #       1  2  3  4
    head = [1, 2, 3, 4]
    tail = [2, 4, 2, 3]
    idx = IncidenceIndex(4, head, tail)
    @test length(idx) == 4
    @test Set(idx[1]) == Set([(2, 1)])
    @test Set(idx[2]) == Set([(1, 1), (4, 2), (3, 3)])
    @test Set(idx[3]) == Set([(2, 3), (4, 4)])
    @test Set(idx[4]) == Set([(2, 2), (3, 4)])
end

@testset "IncidenceIndex: Petgraph::Csr  " begin
    edges = [   # zero indexed!
        (0, 1),
        (0, 2),
        (1, 0),
        (1, 1), # self-reference is not a problem
        (2, 2),
        (2, 4),
    ]
    edges = reshape(reinterpret(Int, edges), 2, length(edges))' .+ 1
    n = maximum(edges)
    @test n == 5
    edges = vec(mapslices(Tuple, edges, dims = [2]))
    idx = IncidenceIndex(n, edges)
    @test Set(map(first, idx[1])) == Set([2, 3])
    @test Set(map(first, idx[2])) == Set([0, 1] .+ 1)
    @test Set(map(first, idx[3])) == Set([0, 2, 4] .+ 1)  # we take both directions!
    @test Set(map(first, idx[4])) == Set([])
    @test Set(map(first, idx[5])) == Set([3])     # both directions!
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
nondatavaluetype(::Type{T}) where {T} = T
nondatavaluetype(::Type{Union{}}) = Union{}
datavaluetype(::Type{T}) where {T} = T
datavaluetype(::Type{Union{}}) = Union{}

Base.@pure function nondatavaluetype(::Type{NT}) where {NT <: NamedTuple{names}} where {names}
    TT = Tuple{Any[ nondatavaluetype(fieldtype(NT, i)) for i = 1:fieldcount(NT) ]...}
    return NamedTuple{names, TT}
end

Base.@pure function datavaluetype(::Tables.Schema{names, types}) where {names, types}
    TT = Tuple{Any[ datavaluetype(fieldtype(types, i)) for i = 1:fieldcount(types) ]...}
    return NamedTuple{names, TT}
end

unwrap(x) = x

struct DataValueUnwrapper{S}
    x::S
end

Tables.istable(::Type{<:DataValueUnwrapper}) = true
Tables.rowaccess(::Type{<:DataValueUnwrapper}) = true
Tables.rows(x::DataValueUnwrapper) = x

function Tables.schema(dv::DataValueUnwrapper)
    eT = eltype(dv.x)
    (!(eT <: NamedTuple) || eT === Union{}) && return nothing
    return Tables.Schema(nondatavaluetype(eT))
end
Base.eltype(rows::DataValueUnwrapper) = DataValueUnwrapRow{eltype(rows.x)}
Base.IteratorSize(::Type{DataValueUnwrapper{S}}) where {S} = Base.IteratorSize(S)
Base.length(rows::DataValueUnwrapper) = length(rows.x)

function Base.iterate(rows::DataValueUnwrapper, st=())
    x = iterate(rows.x, st...)
    x === nothing && return nothing
    row, st = x
    return DataValueUnwrapRow(row), (st,)
end

struct DataValueUnwrapRow{T}
    row::T
end

Base.getproperty(d::DataValueUnwrapRow, ::Type{T}, col::Int, nm::Symbol) where {T} = unwrap(getproperty(getfield(d, 1), T, col, nm))
Base.getproperty(d::DataValueUnwrapRow, nm::Symbol) = unwrap(getproperty(getfield(d, 1), nm))
Base.propertynames(d::DataValueUnwrapRow) = propertynames(getfield(d, 1))


text:
struct LogData

 timestep :: Int32
 dcdfreq :: Int32
 dcdfirststep :: Int32

 time :: Vector{Float64}

 ts :: Vector{Int64}
 bond :: Vector{Float64}
 angle :: Vector{Float64}
 dihed :: Vector{Float64}
 imprp :: Vector{Float64}
 elect :: Vector{Float64}
 vdw :: Vector{Float64}
 boundary :: Vector{Float64}
 misc :: Vector{Float64}
 kinetic :: Vector{Float64}
 total :: Vector{Float64}
 temperature :: Vector{Float64}
 potential :: Vector{Float64}
 total3 :: Vector{Float64}
 tempavg :: Vector{Float64}
 pressure :: Vector{Float64}
 gpressure :: Vector{Float64}
 volume :: Vector{Float64}
 pressavg :: Vector{Float64}
 gpressavg :: Vector{Float64}

end

LogData() = LogData(0, 0, 0, [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.],
                    [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.], [0.])

function Base.show( io :: IO, log :: LogData )
  ndata = length(log.time)
  tempavg = 0.
  pavg = 0.
  eavg = 0.
  for i in 1:ndata
    tempavg = tempavg + log.temperature[i]
    pavg = pavg + log.pressure[i]
    eavg = eavg + log.total[i]
  end
  tempavg = tempavg / ndata
  pavg = pavg / ndata
  eavg = eavg / ndata
  println(" Simulation log data: ")
  println("    Number of steps printed: ", ndata)
  println("    Average temperature: ", tempavg)
  println("    Average pressure: ", pavg)
  println("    Average total energy: ", eavg)
end