text:
#==================== The Macros ====================#

export @cast, @reduce, @matmul

struct CallInfo
    mod::Module
    src::LineNumberNode
    string::String
    flags::Set{Symbol}
    CallInfo(mod::Module, src, str, flags=Set{Symbol}()) = new(mod, src, str, flags)
    CallInfo(syms::Symbol...) = new(Main, LineNumberNode(0), "", Set([syms...]))
end

"""
    @cast Z[i,j,...] := f(A[i,j,...], B[j,k,...])  options

Macro for broadcasting, reshaping, and slicing of arrays in index notation.
Understands the following things:
* `A[i,j,k]` is a three-tensor with these indices.
* `B[(i,j),k]` is the same thing, reshaped to a matrix. Its first axis (the bracket) is indexed
  by `n = i + (j-1) * N` where `i ∈ 1:N`. This may also be written `B[i\\j,k]` or `B[i⊗j,k]`.
* `C[k][i,j]` is a vector of matrices, either created by slicing (if on the left)
  or implying glueing into (if on the right) a 3-tensor `A[i,j,k]`.
* `D[j,k]{i}` is an ordinary matrix of `SVector`s, which may be reinterpreted from `A[i,j,k]`.
* `E[i,_,k]` has two nontrivial dimensions, and `size(E,2)==1`. On the right hand side
  (or when writing to an existing array) you may also write `E[i,3,k]` meaning `view(E, :,3,:)`,
  or `E[i,\$c,j]` to use a variable `c`.
* `f(x)[i,j,k]` is allowed; `f(x)` must return a 3-tensor (and will be evaluated only once).
* `g(H[:,k])[i,j]` is a generalised `mapslices`, with `g` mapping columns of `H`
    to matrices, which are glued into a 3-tensor `A[i,j,k]`.
* `h(I[i], J[j])[k]` expects an `h` which maps two scalars to a vector,
  which gets broadcasted `h.(I,J')`, then glued to make a 3-tensor.
* `K[i,j]'` conjugates each element, equivalent to `K'[j,i]` which is the
  conjugate-transpose of the matrix.
* `M[i,i]` means `diag(M)[i]`, but only for matrices: `N[i,i,k]` is an error.
* `P[i,i']` is normalised to `P[i,i′]` with unicode \\prime.
* `R[i,-j,k]` means roughly `reverse(R, dims=2)`, and `Q[i,~j,k]` similar with `shuffle`.
* `S[i,T[j,k]]` is the 3-tensor `S[:,T]` created by indexing a matrix `S` with `T`,
  where these integers are `all(1 .<= T .<= size(S,2))`.

The left and right hand sides must have all the same indices,
and the only repeated index allowed is `M[i,i]`, which is a diagonal not a trace.
See `@reduce` and `@matmul` for related macros which can sum over things.

If a function of one or more tensors appears on the right hand side,
then this represents a broadcasting operation,
and the necessary re-orientations of axes are automatically inserted.

The following actions are possible:
* `=` writes into an existing array, overwriting its contents,
  while `+=` adds (precisely `Z .= Z .+ ...`) and `*=` multiplies.
* `:=` creates a new array. This need not be named: `Z = @cast [i,j,k] := ...` is allowed.
* `|=` insists that this is a copy, not a view.

Re-ordering of indices `Z[k,j,i]` is done lazily with `PermutedDimsArray(A, ...)`.
Reversing of an axis `F[i,-j,k]` is also done lazily, by `Reverse{2}(F)` which makes a `view`.
Using `|=` (or broadcasting) will produce a simple `Array`.

Options can be specified at the end (if several, separated by `,` i.e. `options::Tuple`)
* `i:3` supplies the range of index `i`. Variables and functions like `j:Nj, k:length(K)`
  are allowed.
* `assert` will turn on explicit dimension checks of the input.
  (Providing any ranges will also turn these on.)
* `cat` will glue slices by things like `hcat(A...)` instead of the default `reduce(hcat, A)`,
  and `lazy` will instead make a `LazyStack.Stacked` container.
* `nolazy` disables `PermutedDimsArray` and `Reverse` in favour of `permutedims` and `reverse`,
  and `Diagonal` in favour of `diagm` for `Z[i,i]` output.
* `strided` will place `@strided` in front of broadcasting operations,
  and use `@strided permutedims(A, ...)` instead of `PermutedDimsArray(A, ...)`.
  For this you need `using Strided` to load that package.
* `avx` will place `@avx` in front of broadcasting operations,
  which needs `using LoopVectorization` to load that package.

To create static slices `D[k]{i,j}` you should give all slice dimensions explicitly.
You may write `D[k]{i:2,j:2}` to specify `Size(2,2)` slices.
They are made most cleanly from the first indices of the input,
i.e. this `D` from `A[i,j,k]`. The notation `A{:,:,k}` will only work in this order,
and writing `A{:2,:2,k}` provides the sizes.
"""
macro cast(exs...)
    call = CallInfo(__module__, __source__, TensorCast.unparse("@cast", exs...))
    _macro(exs...; call=call)
end

"""
    @reduce A[i] := sum(j,k) B[i,j,k]             # A = vec(sum(B, dims=(2,3)))
    @reduce A[i] := prod(j) B[i] + ε * C[i,j]     # A = vec(prod(B .+ ε .* C, dims=2))
    @reduce A[i] = sum(j) exp( C[i,j] / D[j] )    # sum!(A, exp.(C ./ D') )

Tensor reduction macro:
* The reduction function can be anything which works like `sum(B, dims=(1,3))`,
  for instance `prod` and `maximum` and `Statistics.mean`.
* In-place operations `Z[j] = sum(...` will construct the banged version of the given
  function's name, which must work like `sum!(Z, A)`.
* The tensors can be anything that `@cast` understands, including gluing of slices `B[i,k][j]`
  and reshaping `B[i\\j,k]`. See `? @cast` for the complete list.
* If there is a function of one or more tensors on the right,
  then this is a broadcasting operation.
* Index ranges may be given afterwards (as for `@cast`) or inside the reduction `sum(i:3, k:4)`.
* All indices appearing on the right must appear either within `sum(...)` etc, or on the left.

    F = @reduce sum(i,j)  B[i] + γ * D[j]         # sum(B .+ γ .* D')
    @reduce G[] := sum(i,j)  B[i] + γ * D[j]      # F == G[]

Complete reduction to a scalar output `F`, or a zero-dim array `G`.
`G[]` involves `sum(A, dims=(1,2))` rather than `sum(A)`.

    @reduce Z[k] := sum(i,j) A[i] * B[j] * C[k]  lazy, i:N, j:N, k:N

The option `lazy` replaces the broadcast expression with a `BroadcastArray`,
to avoid `materialize`ing the entire array before summing. In the example this is of size `N^3`.

The option `strided` will place `@strided` in front of the broadcasting operation.
You need `using Strided` for this to work.

    @reduce sum(i) A[i] * log(@reduce [i] := sum(j) A[j] * exp(B[i,j]))
    @cast W[i] := A[i] * exp(- @reduce S[i] = sum(j) exp(B[i,j]) lazy)

Recursion like this is allowed, inside either `@cast` or `@reduce`.
The intermediate array need not have a name, like `[i]`,
unless writing into an existing array, like `S[i]` here.
"""
macro reduce(exs...)
    call = CallInfo(__module__, __source__, TensorCast.unparse("@reduce", exs...),
        Set([:reduce]))
    _macro(exs...; call=call)
end

"""
    @matmul M[a,c] := sum(b)  A[a,b] * B[b,c]

Matrix multiplication macro. Uses the same syntax as `@reduce`, but instead of broadcasting
the expression on the right out to a 3-tensor before summing it along one dimension,
this calls `*` or `mul!`, which is usually much faster.
But it only works on expressions of suitable form.

With more than two tensors on the right, it proceeds left to right, and each summed index
must appear on two tensors, probably neighbours.

Note that unlike `@einsum` and `@tensor`, you must explicitly specify
what indices to sum over. Both this macro and `@reduce` could in principal infer this,
and perhaps I will add that later... but then I find myself commenting `# sum_μ` anyway.

    @matmul Z[a⊗b,z] = sum(i,j,k)  D[i,a][j] * E[i⊗j,_,k,b] * F[z,3,k]

Each tensor will be pre-processed exactly as for `@cast` / `@reduce`,
here glueing slices of `D` together, reshaping `E`, and taking a view of `F`.
Once this is done, the right hand side must be of the form `(tensor) * (tensor) * ...`,
which becomes `mul!(ZZ, (DD * EE), FF)`.

    @reduce V[i] := sum(k) W[k] * exp(@matmul [i,k] := sum(j) A[i,j] * B[j,k])

You should be able to use this within the other macros, as shown.
"""
macro matmul(exs...)
    call = CallInfo(__module__, __source__, TensorCast.unparse("@matmul", exs...),
        Set([:matmul, :nolazy]))
    _macro(exs...; call=call)
end

macro mul(exs...)
    call = CallInfo(__module__, __source__, TensorCast.unparse("@matmul", exs...),
        Set([:matmul, :nolazy]))
    @warn "please replace @mul with @matmul, and ensure that it has explicit sum()" call.string maxlog=1 _id=hash(exs)
    _macro(exs...; call=call)
end

#==================== The Main Functions ====================#

function _macro(exone, extwo=nothing, exthree=nothing; call::CallInfo=CallInfo(), dict=Dict())
    store = (dict=dict, assert=[], mustassert=[], seen=[], need=[], top=[], main=[])
    # TODO use OrderedDict() for main? To allow duplicate removal

    if (:reduce in call.flags) || (:matmul in call.flags)
        # First the options:
        optionparse(exthree, store, call)
        # Then the LHS, to get canonical list of indices:
        canon, parsed = reduceparse(exone, extwo, store, call)
    else
        # Much the same for @cast:
        isnothing(exthree) || throw(MacroError("@cast doesn't take three expressions: $exthree", call))
        optionparse(extwo, store, call)
        canon, parsed = castparse(exone, store, call)
    end

    # First pass over RHS just to read sizes, prewalk sees A[i][j] before A[i]
    MacroTools.prewalk(x -> rightsizes(x, store, call), parsed.right)

    # To look for recursion, we need another prewalk:
    right2 = MacroTools.prewalk(x -> recursemacro(x, store, call), parsed.right)

    # Third pass to standardise & then glue, postwalk sees A[i] before A[i][j]
    right3 = MacroTools.postwalk(x -> standardglue(x, canon, store, call), right2)

    if !(:matmul in call.flags)
        # Then finally broadcasting if necc (or just permutedims etc. if not):
        right4 = targetcast(right3, canon, store, call)
    else
        # or, for matmul, can I change just this? outputinplace is also going to need changing.
        right4 = matmultarget(right3, canon, parsed, store, call)
    end

    checkallseen(canon, store, call) # this must be run before inplaceoutput()

    # Return to LHS, build up what it requested:
    if :inplace in call.flags
        rightlist = inplaceoutput(right4, canon, parsed, store, call)
    else
        right5 = newoutput(right4, canon, parsed, store, call)
        rightlist = [:( $(parsed.name) = $right5 )]
    end

    # Sew all these pieces into output:
    outex = quote end
    append!(outex.args, store.top)
    append!(outex.args, findsizes(store, call)) # this must be run after newoutput() etc.
    append!(outex.args, store.main)
    append!(outex.args, rightlist)

    if :recurse in call.flags
        return (name=parsed.name, ind=parsed.outer, scalar=(:scalar in call.flags), steps=outex.args)
    else
        return esc(outex) # could use MacroTools.unblock()?
    end
end

"""
    standardise(f(x)[i,3]) -> view(A,:,3)[i]
    standardise(A[(i,j)]) -> reshape(A, sz_i, sz_j)[i,j]

This mostly aims to re-work the given expression into `some(steps(A))[i,j]`,
but also pushes `A = f(x)` into `store.top`.
"""
function standardise(ex, store::NamedTuple, call::CallInfo; LHS=false)
    # This acts only on single indexing expressions:
    if @capture(ex, A_{ijk__})
        static=true
        push!(call.flags, :staticslice)
    elseif @capture(ex, A_[ijk__])
        static=false
    else
        return ex
    end

    # Ensure that f(x)[i,j] will evaluate once, including in size(A)
    if A isa Symbol || @capture(A, AA_.ff_) # caller has ensured !containsindexing(A)
    else
        Asym = Symbol(A,"_val") # exact same symbol is used by rightsizes()
        push!(store.top,  :( local $Asym = $A ) )
        A = Asym
    end

    # Constant indices A[i,3], A[i,_], A[i,$c], and also A[i,3:5]
    if all(isCorR, ijk)
        needview!(ijk)                                         # replace _ with 1, if any
        Asym = maybepush(:( $A[$(ijk...)] ), store, :allconst) # protect from later processing
        return Asym                                            # and nothing more to do here
    elseif any(isCorR, ijk)
        ijcolon = map(i -> isCorR(i) ? i : (:), ijk)
        if needview!(ijcolon)
            A = :( view($A, $(ijcolon...)) ) # TODO nolazy?
        else
            A = :( TensorCast.rview($A, $(ijcolon...)) )
        end
        ijk = filter(!isconstant, ijk)              # remove actual constants from list,
        ijk = map(i -> isrange(i) ? :(:) : i, ijk)  # but for A[i,3:5] replace with : symbol, for slicing
    end

    # Slicing operations A[i,:,j], and A[i,3:5] after the above treatment
    if any(iscolon, ijk)
        code = Tuple(map(i -> iscolon(i) ? (:) : (*), ijk))
        if !static && (:collect in call.flags)
            A = :( TensorCast.slicecopy($A, $code) )
        elseif !static
            A = :( TensorCast.sliceview($A, $code) )
        else
            sizeorcode = maybestaticsizes(ijk, code, call)
            A = :( TensorCast.static_slice($A, $sizeorcode) )
        end
        ijk = filter(!iscolon, ijk)
    elseif static
        throw(MacroError("shouldn't use curly brackets here", call))
    end

    # Nested indices A[i,j,B[k,l,m],n] or worse A[i,B[j,k],C[i,j]]
    if any(i -> @capture(i, B_[klm__]), ijk)
        newijk, beecolon = [], [] # for simple case
        # listB, listijk = [], []
        for i in ijk
            if @capture(i, B_[klm__])
                append!(newijk, klm)
                push!(beecolon, B)
                # push!(listijk, klm)
                # push!(listB, B)
            else
                push!(newijk, i)
                push!(beecolon, (:))
            end
        end
        if allunique(newijk) # then we are in simple case
            ijk = newijk
            A = :( view($A, $(beecolon...)) )
        else
            error("can't handle this indexing yet, sorry")
            # LHS && error("can't do this indexing on LHS sorry")
            # A = maybepush(A, store, :preget)
            # target = guesstarget(??)
            # Bs = [ targetcast(B, target, store, call) for B in listB ]
            # Aex = getindex.(Ref(A), i)
            # A = maybepush(Aex, store, :getindex)
        end
    end

    # Diagonal extraction A[i,i]
    if length(ijk)==2 && ijk[1]==ijk[2]
        if (:nolazy in call.flags) && !LHS # don't do this for in-place output
            A = :( TensorCast.diag($A) ) # LinearAlgebra really
        else
            A = :( TensorCast.diagview($A) )
        end
        pop!(ijk)
    end

    # Parse the remaining indices; saving of their sizes was done elsewhere:
    flat, _, reversed, shuffled = indexparse(A, ijk)
    append!(store.seen, flat)

    # Combined indices A[i,(j,k)]
    if any(istensor, ijk)
        flatsize = map(szwrap, flat)
        A = :( reshape($A, ($(flatsize...),)) )
        append!(store.need, flat)
        # push!(store.mustassert, :( TensorCast.@assert_ !(A isa OffsetArray) "can't combine indices of an OffsetArray") ) # ??
    end

    # Reversed A[-i,j] and shuffled A[i,~j]
    for i in reversed
        di = findfirst(isequal(i), flat)
        if (:nolazy in call.flags) && !LHS
            A = :( reverse($A, dims=$di) )
        else
            A = :( TensorCast.Reverse{$di}($A) )
        end
    end

    if length(shuffled) > 0
        if (:nolazy in call.flags) && !LHS
            if length(flat) == 1
                A = :( TensorCast.shuffle(A) )
            else # sadly shuffle(A, dims=...) doesn't exist!
                code = Tuple(Any[ i in shuffled ? (*) : (:) for i in flat ])
                A = :( TensorCast.red_glue(TensorCast.shuffle(
                    TensorCast.slicecopy($A, $code)),$code) )
            end
        else
            for i in shuffled
                di = findfirst(isequal(i), flat)
                A = :( TensorCast.Shuffle{$di}($A) )
            end
        end
    end

    # Construct final expression
    return :( $A[$(flat...)] )
end

"""
    standardglue(expr) -> func(old)[i,j,k]

This gets walked over the RHS and:
* standardises any simple expression `f(x)[i⊗j,3,k]` by calling `standardise()`,
* finds anything which needs glueing either directly `A[j,k][i]`,
  or after broadcasting `(A[j]+B[k])[i]` using `targetcast()` first.

target dims not correctly handled yet -- what do I want? TODO
Simple glue / stand. does not permutedims, but broadcasting may have to... avoid twice?
"""
function standardglue(ex, target, store::NamedTuple, call::CallInfo)

    # The sole target here is indexing expressions:
    if @capture(ex, A_[inner__])
        static=false
    elseif @capture(ex, A_{inner__})
        static=true
    else
        return ex
    end

    # Do we have something simeple like M[i,j], no worse than f(x)[i,:,j⊗k,3]?
    if !containsindexing(A)
        return standardise(ex, store, call)
    end

    # Otherwise there are two options, (brodcasting...)[k] or simple B[i,j][k]
    needcast = !@capture(A, B_[outer__])

    if needcast
        outer = unique(reduce(vcat, listindices(A)))
        outer = sort(outer, by = i -> findcheck(i, target, call)) # , " in target list"
        Bex = targetcast(A, outer, store, call)
        B = maybepush(Bex, store, :innercast)
    end

    # Before gluing, deal with any constant inner indices
    Bsym = gensym(:innerfix)
    if all(isconstant, inner)
        B = maybepush(B, store, :preget) # protect any previous operations from @.
        needview!(inner) # replace _ with 1, if any
        push!(store.main, :( $Bsym =  @__dot__ getindex($B, $(inner...)) ) )
        return  :( $Bsym[$(outer...)] ) # actually, no glueing at all to do!

    elseif any(isconstant, inner)
        B = maybepush(B, store, :preconst)
        ijcolon = map(i -> isconstant(i) ? i : (:), inner)
        if needview!(ijcolon)
            ex = :( $Bsym =  @__dot__ view($B, $(ijcolon...)) ) # TODO nolazy here
        else
            ex = :( $Bsym = @__dot__ TensorCast.rview($B, $(ijcolon...)) )
        end
        push!(store.main, ex)
        B = Bsym
        inner = filter(!isconstant, inner)
    end

    # And some quick error checks; standsrdise() already put inner into store.seen
    append!(store.seen, inner)
    checknorepeats(vcat(inner, outer), call,
        " in gluing " * string(:([$(outer...)])) * string(:([$(inner...)])))

    # Now we glue, always in A[k,l][i,j] -> B[i,j,k,l] order to start with:
    ijk = vcat(inner, outer)
    code = Tuple(Any[ i in inner ? (:) : (*) for i in ijk ])

    if static
        AB = :( TensorCast.static_glue($B) )
        pop!(call.flags, :collected, :ok)
    elseif :glue in call.flags
        # allow maximum freedom? TODO
        code, ijk = gluereorder(code, ijk, target)

        AB = :( TensorCast.copy_glue($B, $code) )
    elseif :cat in call.flags
        AB = :( TensorCast.cat_glue($B, $code) )
        push!(call.flags, :collected)
    elseif :lazy in call.flags
        AB = :( TensorCast.lazy_glue($B, $code) )
        pop!(call.flags, :collected, :ok)
    else
        # allow a little freedom?
        if code == (*,:) && ijk == reverse(target)
            code = (:,*)
            ijk = target
        end

        AB = :( TensorCast.red_glue($B, $code) )
        push!(call.flags, :collected)
    end

    return :( $AB[$(ijk...)] )
end

"""
    targetcast(A[i] + B[j], [i,j]) -> @__dot__(A + B')
    targetcast(A[j,i], [i,j]) -> transpose(A)

This beings the expression to have target indices,
by permutedims and if necessary broadcasting, always using `readycast()`.
"""
function targetcast(ex, target, store::NamedTuple, call::CallInfo)

    # If just one naked expression, then we won't broadcast:
    if @capture(ex, A_[ijk__])
        containsindexing(A) && error("that should have been dealt with")
        return readycast(ex, target, store, call)
    end

    # But for anything more complicated, we do:
    ex = MacroTools.prewalk(x -> readycast(x, target, store, call), ex)
    if :lazy in call.flags
        ex = :( TensorCast.lazy($ex) )
        pop!(call.flags, :collected, :ok)
    else
        push!(call.flags, :collected)
    end
    if :strided in call.flags
        ex = Broadcast.__dot__(ex) # @strided does not work on @.
        ex = :( Strided.@strided $ex )
    elseif :avx in call.flags
        ex = :( LoopVectorization.@avx @__dot__($ex) )
    else
        ex = :( @__dot__($ex) )
    end
    return ex

end

"""
    readycast(A[j,i], [i,j,k]) -> transpose(A)
    readycast(A[k], [i,j,k]) -> orient(A, (*,*,:))

This is walked over the expression to prepare for `@__dot__` etc, by `targetcast()`.
"""
function readycast(ex, target, store::NamedTuple, call::CallInfo)

    # Scalar functions can be protected entirely from broadcasting:
    # TODO this means A[i,j] + rand()/10 doesn't work, /(...,10) is a function!
    if @capture(ex, fun_(arg__) | fun_(arg__).field_ ) && !containsindexing(fun) &&
            length(arg)>0 && !any(containsindexing, arg) # pull out log(2) but not randn()
        return maybepush(ex, store, :scalarfun)
    end

    # Some things ought to apply elementwise: conjugation,
    @capture(ex, A_[ijk__]') && return :( Base.adjoint($A[$(ijk...)]) )
    # .fields... only one deep for now,
    @capture(ex, A_[ijk__].field_ ) &&
        return :( getproperty($A[$(ijk...)], $(QuoteNode(field))) )
    @capture(ex, fun_(arg__).field_ ) && any(containsindexing, arg) &&
        return :( getproperty($fun($(arg...)), $(QuoteNode(field))) )
    # tuple creation... now including namedtuples
    @capture(ex, (args__,) ) && any(containsindexing, args) &&
        if any(a -> @capture(a, sym_ = val_), args)
            syms, vals = [], []
            map(args) do a
                @capture(a, sym_ = val_ ) || throw(MacroError("invalid named tuple element $a", call))
                push!(syms, QuoteNode(sym))
                push!(vals, val)
            end
            return :( NamedTuple{($(syms...),)}(tuple($(vals...))) )
        else
            return :( tuple($(args...)) )
        end
    # and arrays of functions, using apply:
    @capture(ex, funs_[ijk__](args__) ) &&
        return :( Core._apply($funs[$(ijk...)], $(args...) ) )

    # Apart from those, readycast acts only on lone tensors:
    @capture(ex, A_[ijk__]) || return ex

    dims = Int[ findcheck(i, target, call, " on the left") for i in ijk ]

    # Do we need permutedims, or equvalent?
    perm = sortperm(dims)
    if perm != 1:length(dims)
        tupleperm = Tuple(perm)
        if :strided in call.flags
            A = :( Strided.@strided permutedims($A, $tupleperm) )
        elseif :nolazy in call.flags
            A = :( permutedims($A, $tupleperm) )
            push!(call.flags, :collected)
        elseif tupleperm == (2,1)
            A = :( TensorCast.PermuteDims($A) )
        else
            A = :( PermutedDimsArray($A, $tupleperm) )
            pop!(call.flags, :collected, :ok)
        end
    end

    # Do we need orient()?
    if length(dims)>0 && maximum(dims) != length(dims)
        code = Tuple(Any[ d in dims ? (:) : (*) for d=1:maximum(dims) ])
        A = :( TensorCast.orient($A, $code) )
    end
    # Those two steps are what might be replaced with TransmuteDims

    # Does A need protecting from @__dot__?
    A = maybepush(A, store, :nobroadcast)

    return A
end

"""
    matmultarget(A[i,j] * B[j,k], [i,k]) -> A * B

For inplace, instead it returns tuple `(A,B)`, which `inplaceoutput()` can use.
If there are more than two factors, it recurses, and you get `(A*B) * C`,
or perhaps tuple `(A*B, C)`.
"""
function matmultarget(ex, target, parsed, store::NamedTuple, call::CallInfo)

    @capture(ex, A_ * B_ * C__ | *(A_, B_, C__) ) || throw(MacroError("can't @matmul that!", call))

    # Figure out what to sum over, and make A,B into matrices ready for *
    iA = guesstarget(A)
    iB = guesstarget(B)

    isum = sort(intersect(iA, iB, parsed.reduced),
        by = i -> findfirst(isequal(i), target)) # or target? parsed.reduced
    iAnosum = setdiff(iA, isum)
    iBnosum = setdiff(iB, isum)

    Aex = targetcast(A, vcat(iAnosum, isum), store, call)
    Bex = targetcast(B, vcat(isum, iBnosum), store, call)

    Aex = matrixshape(Aex, iAnosum, isum, store, call)
    Bex = matrixshape(Bex, isum, iBnosum, store, call)

    Aex = maybepush(Aex, store, :mulA)
    Bex = maybepush(Bex, store, :mulB)

    # But don't actually * if you are going to mul!(Z,A,B) instead later:
    if (:inplace in call.flags) && length(C)==0
        return (Aex, Bex, iAnosum, iBnosum)
    end

    ABex = unmatrixshape(:( $Aex * $Bex ), iAnosum, iBnosum, store, call)
    ABsym = gensym(:matmul)
    push!(store.main, :( local $ABsym = $ABex ))
    ABind = :( $ABsym[$(iAnosum...), $(iBnosum...)] )

    if length(C) == 0
        # Last step for Z := A * B case is reshape & permutedims etc:
        subtarget = setdiff(target, isum)
        return targetcast(ABind, subtarget, store, call)

    else
        # But if we have (A*B) * C then continue, whether in- or out-of-place:
        @warn "@matmul has bugs in how it handles more than two factors" call.string maxlog=1 _id=hash(call.string)
        ABCex = :( *($ABind, $(C...)) )
        return matmultarget(ABCex, target, parsed, store, call)
    end
end

"""
    recursemacro(@reduce sum(i) A[i,j]) -> G[j]

Walked over RHS to look for `@reduce ...`, and replace with result,
pushing calculation steps into store.

Also a convenient place to tidy all indices, including e.g. `fun(M[:,j],N[j]).same[i']`.
"""
function recursemacro(ex, store::NamedTuple, call::CallInfo)

    # Actually look for recursion
    if @capture(ex, @reduce(subex__) )
        subcall = CallInfo(call.mod, call.src,
            TensorCast.unparse("innder @reduce", subex...), Set([:reduce, :recurse]))
        name, ind, scalar, steps = _macro(subex...; call=subcall, dict=store.dict)
        append!(store.main, steps)
        ex = scalar ? :($name) :  :($name[$(ind...)])

    elseif @capture(ex, @matmul(subex__) )
        subcall = CallInfo(call.mod, call.src,
            TensorCast.unparse("innder @matmul", subex...), Set([:matmul, :recurse]))
        name, ind, scalar, steps = _macro(subex...; call=subcall, dict=store.dict)
        append!(store.main, steps)
        ex = scalar ? :($name) :  :($name[$(ind...)])

    elseif @capture(ex, @cast(subex__) )
        subcall = CallInfo(call.mod, call.src,
            TensorCast.unparse("inner @cast", subex...), Set([:recurse]))
        name, ind, scalar, steps = _macro(subex...; call=subcall, dict=store.dict)
        append!(store.main, steps)
        ex = scalar ? :($name) :  :($name[$(ind...)])
    end

    # Tidy up indices, A[i,j][k] will be hit on different rounds...
    if @capture(ex, A_[ijk__])
        return :( $A[$(tensorprimetidy(ijk)...)] )
    elseif @capture(ex, A_{ijk__})
        return :( $A{$(tensorprimetidy(ijk)...)} )
    else
        return ex
    end
end

"""
    rightsizes(A[i,j][k])

This saves to `store` the sizes of all input tensors, and their sub-slices if any.
* Deals with `A[i,j][k]` at once, before seeing `A[i,j]`, hence `prewalk` & destruction.
* For `fun(x)[i,j]` it saves `sz_i` under name `Symbol(fun(x),"_val")` used later by standardise
* But for `fun(M[:,j],N[j]).same[i]` it can't save `sz_i` as this isn't calculated yet,
  however it should not destroy this so that `sz_j` can be got later.
"""
function rightsizes(ex, store::NamedTuple, call::CallInfo)
    :recurse in call.flags && return nothing # outer version took care of this

    if @capture(ex, A_[outer__][inner__] | A_[outer__]{inner__} )
        field = nothing
    elseif @capture(ex, A_[outer__].field_[inner__] | A_[outer__].field_{inner__} )
    elseif  @capture(ex, A_[outer__] | A_{outer__} )
        field = nothing
    else
        return ex
    end

    # Special treatment for  fun(x)[i,j], goldilocks A not just symbol, but no indexing
    if A isa Symbol || @capture(A, AA_.ff_)
    elseif !containsindexing(A)
        A = Symbol(A,"_val") # the exact same symbol is used by standardiser
    end

    # When we can save the sizes, then we destroy so as not to save again:
    if A isa Symbol || @capture(A, AA_.ff_) && !containsindexing(A)
        indexparse(A, outer, store, call; save=true)
        if field==nothing
            innerparse(:(first($A)), inner, store, call; save=true)
        else
            innerparse(:(first($A).$field), inner, store, call; save=true)
        end
        return nothing
    end

    return ex
end

#==================== Various Parsing Functions ====================#

"""
    canon, parsed = castparse(Z[out][inner] := right)

For `@cast`, this digests the LHS.
`parsed.outer` is the indices exactly as given, `parsed.flat` cleans up,
and `canon` prepends inner indices too.
"""
function castparse(ex, store::NamedTuple, call::CallInfo; reduce=false)
    Z = gensym(:left)

    # Do we make a new array? With or without collecting:
    if @capture(ex, left_ := right_ )
    elseif @capture(ex, left_ == right_ )
        @warn "using == no longer does anything" call.string maxlog=1 _id=hash(call.string)
    elseif @capture(ex, left_ |= right_ )
        push!(call.flags, :collect)

    # Do we write into an exising array? Possibly updating it:
    elseif @capture(ex, left_ = right_ )
        push!(call.flags, :inplace)
    elseif @capture(ex, left_ += right_ )
        push!(call.flags, :inplace)
        right = :( $left + $right )
        reduce && throw(MacroError("can't use += with @reduce", call))
    elseif @capture(ex, left_ -= right_ )
        push!(call.flags, :inplace)
        right = :( $left - ($right) )
        reduce && throw(MacroError("can't use -= with @reduce", call))
    elseif @capture(ex, left_ *= right_ )
        push!(call.flags, :inplace)
        right = :( $left * ($right) )
        reduce && throw(MacroError("can't use *= with @reduce", call))

    # @reduce might not have a LHS, otherwise some errors:
    elseif @capture(ex, right_ => left_ )
        throw(MacroError("anonymous functions using => currently not supported", call))
        # @capture(right, AA_[ii__] | AA_[ii__].ff_ | AA_{ii__} ) ||
        #     throw(MacroError("anonymous functions => only accept simpler input A[...] => ...", call))
        # push!(call.flags, :anonfunc) # TODO make the output understand this
    elseif reduce
        return (Any[], nothing)
    else
        error("wtf is $ex")
    end

    static = @capture(left, ZZ_{ii__})

    if @capture(left, Z_[outer__][inner__] | [outer__][inner__] | Z_[outer__]{inner__} | [outer__]{inner__} )
        isnothing(Z) && (:inplace in call.flags) && throw(MacroError("can't write into a nameless tensor", call))
        Z = isnothing(Z) ? gensym(:output) : Z
        parsed = indexparse(Z, outer, store, call, save=(:inplace in call.flags))
        innerflat = innerparse(:(first($Z)), inner, store, call, save=(:inplace in call.flags))
        canon = vcat(innerflat, parsed.flat)
        checknorepeats(canon, call, " on the left")

    elseif @capture(left, Z_[outer__] | [outer__] )
        isnothing(Z) && (:inplace in call.flags) && throw(MacroError("can't write into a nameless tensor", call))
        Z = isnothing(Z) ? gensym(:output) : Z
        parsed = indexparse(Z, outer, store, call, save=(:inplace in call.flags))
        canon = parsed.flat
        inner, innerflat = [], []

    elseif left isa Symbol # only for @reduce A := sum(i) ...
        (:reduce in call.flags) || throw(MacroError("@cast really needs an indexing expression on left!", call))
        (:inplace in call.flags) && throw(MacroError("scalar output needs @reduce Z := ...", call))
        Z = left
        push!(call.flags, :scalar)
        parsed = indexparse(Z, [], store, call) # just to get the right namedtuple
        canon, outer, inner, innerflat = [], [], [], []

    else
        throw(MacroError("don't know what to do with left = $left", call))
    end

    return canon, (parsed..., name=Z, outer=outer, inner=inner, innerflat=innerflat,
        static=static, left=left, right=right)
end

"""
    canon, parsed = reduceparse(Z[i,j] := sum(j), right)

Similar to `castparse()`, which it uses before adding on what `@reduce` needs.
Tries to be smart about the order of `canon` in the hope that reduction over `parsed.rdims`
can be done without `permutedims`.
"""
function reduceparse(ex1, ex2, store::NamedTuple, call::CallInfo)

    # Parse the LHS if there is one, else only @reduce sum(i,j) A[i,j]
    leftcanon, parsed = castparse(ex1, store, call; reduce=true)
    if parsed==nothing
        @capture(ex1, redfun_(redlist__))
        push!(call.flags, :scalar)
        parsed = (indexparse(nothing, [])..., name=gensym(:noleft), outer=[], inner=[], static=false)
    else
        @capture(parsed.right, redfun_(redlist__) ) || throw(MacroError("how do I reduce over $(parsed.right) ?", call))
    end

    # Then parse redlist, decoding ranges like sum(i:10,j) which specify sizes
    reduced = []
    for item in tensorprimetidy(redlist) # normalise θ'
        i = @capture(item, j_:s_) ? saveonesize(j, s, store) : item
        push!(reduced, i)
    end
    checknorepeats(reduced, call, " in the reduction") # catches sum(i,i) B[i,j,k]

    # Now work out canonical list. For sum!(), reduced indices must come last.
    if :inplace in call.flags
        canon = vcat(leftcanon, reduced)
    else
        # But for Z = sum(A, dims=...) can try to avoid permutedims, not sure it matters.
        guess = guesstarget(ex2) # TODO make guess smarter, use leftcanon as a target
        # @show guess reduced
        [ deleteat!(guess, findcheck(i, guess, call, " on the right")) for i in reduced ]
        if leftcanon == guess
            canon = guesstarget(ex2)
        else
            canon = vcat(leftcanon, reduced)
        end
    end
    checknorepeats(canon, call, " in the reduction") # catches A[i] := sum(i)

    # Finally, record positions in canon of reductions
    rdims = sort([ findfirst(isequal(i), canon) for i in reduced ])

    return canon, (parsed..., redfun=redfun, reduced=reduced, rdims=rdims, right=ex2) # this replaces parsed.right
end

"""
    p = indexparse(A, [i,j,k])

`p.flat` strips constants, colons, tensor/brackets, minus, tilde, and doubled indices.
`p.outsize` is a list like `[sz_i, 1, star(...)]` for use on LHS.
`p.reversed` has all those with a minus.
Sizes are saved to `store` only with keyword `save=true`, i.e. only when called by `rightsizes()`
"""
function indexparse(A, ijk::Vector, store=nothing, call=nothing; save=false)
    flat, outsize, reversed, shuffled = [], [], [], []

    ijk = tensorprimetidy(ijk) # un-wrap i⊗j to tuples, and normalise i' to i′

    stripminustilde!(ijk, reversed, shuffled)

    for (d,i) in enumerate(ijk)

        if iscolon(i) || isrange(i)
            if i isa QuoteNode
                str = "fixed size in $A[" * join(ijk, ", ") * "]" # DimensionMismatch("fixed size in M[i, \$(QuoteNode(5))]: size(M, 2) == 5") TODO print more nicely
                push!(store.mustassert, :( TensorCast.@assert_ size($A,$d)==$(i.value) $str) )
            end
            continue
        end

        if isconstant(i)
            push!(outsize, 1)
            if i == :_ && save
                str = "underscore in $A[" * join(ijk, ", ") * "]"
                push!(store.mustassert, :( TensorCast.@assert_ size($A,$d)==1 $str) )
            end
            continue
        end

        if @capture(i, (ii__,))
            stripminustilde!(ii, reversed, shuffled)
            append!(flat, ii)
            push!(outsize, szwrap(ii))
            save && saveonesize(ii, :(size($A, $d)), store)

        elseif @capture(i, B_[klm__])
            innerparse(B, klm, store, call) # called just for error on tensor/colon/constant
            sub = indexparse(B, klm, store, call; save=save) # I do want to save size(B,1) etc.
            append!(flat, sub.flat)

        elseif i isa Symbol
            push!(flat, i)
            push!(outsize, szwrap(i))
            save && saveonesize(i, :(size($A, $d)), store)

        else
            throw(MacroError("don't understand index $i", call))
        end
    end

    if save && length(ijk)>0
        N = length(ijk)
        str = "expected a $N-tensor $A[" * join(ijk, ", ") * "]"
        push!(store.assert, :( TensorCast.@assert_ ndims($A)==$N $str) )
    end

    if length(flat)==2 && flat[1]==flat[2] # allow for diag, A[i,i]
        pop!(flat)
    end

    checknorepeats(flat, call)

    reversed = filter(i -> isodd(count(isequal(i), reversed)), reversed)
    shuffled = unique(shuffled)

    return (flat=flat, outsize=outsize, reversed=reversed, shuffled=shuffled)
end

function stripminustilde!(ijk::Vector, reversed, shuffled)
    # @show ijk
    for (d,i) in enumerate(ijk)
        # if @capture(i, -(j, jj__))
        #     append!(reversed, vcat(j,jj))
        #     ijk[d] = :( ($j,$(jj...),) )

        if @capture(i, -j_)
            push!(reversed, j)
            ijk[d] = j
        elseif @capture(i, ~j_)
            push!(shuffled, j)
            ijk[d] = j
        end
    end
end

"""
    innerparse(firstA, [i,j])
Checks that inner indices are kosher, and returns the list.
If `save=true` then it now expects expr. `first(A)` to allow `first(A).field` too,
otherwise it ignores first arg.
"""
function innerparse(firstA, ijk, store::NamedTuple, call::CallInfo; save=false)
    isnothing(ijk) && return []

    ijk = tensorprimetidy(ijk) # only for primes really

    innerflat = []
    for (d,i) in enumerate(ijk)
        iscolon(i) && throw(MacroError("can't have a colon in inner index!", call))
        isrange(i) && @capture(i, alpha_Int:omega_)  && throw(MacroError("can't have a range in inner index!", call)) # this is an imperfect check, must not be triggered by @cast R2[j]{i:3} := M[i,j]  which is another kind of range!
        istensor(i) && throw(MacroError("can't tensor product inner indices", call))
        (@capture(i, -j_) || @capture(i, ~j_)) && throw(MacroError("can't reverse or shuffle along inner indices", call))

        if @capture(i, j_:s_)
            push!(innerflat, j)
            saveonesize(j, s, store) # save=true on LHS only for in-place, save this anyway
        elseif isconstant(i)
            i == :_ && save && push!(store.mustassert, :(TensorCast.@assert_ size($firstA, $d)==1 "inner underscore") )
        else
            push!(innerflat, i)
        end
        save && saveonesize(i, :(size($firstA, $d)), store)
    end

    checknorepeats(innerflat, call)

    return innerflat
end

function optionparse(opt, store::NamedTuple, call::CallInfo)
    if isnothing(opt)
        return
    elseif @capture(opt, (opts__,) )
        [ optionparse(o, store, call) for o in opts ]
        return
    end

    if @capture(opt, i_:s_)
        saveonesize(tensorprimetidy(i), s, store)
        push!(call.flags, :assert)
    elseif opt in (:assert, :lazy, :nolazy, :cat, :strided, :avx)
        push!(call.flags, opt)
    elseif opt == :(!)
        @warn "please replace option ! with assert" call.string maxlog=1 _id=hash(call.string)
        push!(call.flags, :assert)
    else
        throw(MacroError("don't understand option $opt", call))
    end
end

#==================== Smaller Helper Functions ====================#

"""
    saveonesize(:i, size(A,3), store) -> :i

Saves sizes for indices `:i` or products like  `[:i,:j]` to `store`.
It it already has a size for this single index,
then save an assertion that new size is equal to old.
"""
function saveonesize(ind, long, store::NamedTuple)
    if !haskey(store.dict, ind)
        store.dict[ind] = long
    else
        if isa(ind, Symbol)
            str = "range of index $ind must agree"
            push!(store.assert, :(TensorCast.@assert_ $(store.dict[ind]) == $long $str) )
        end
    end
    ind
end

function findsizes(store::NamedTuple, call::CallInfo)
    out = []
    if :assert in call.flags
        append!(out, store.assert)
        empty!(store.assert)
    end
    if length(store.need) > 0
        sizes = sizeinfer(store, call)
        sz_list = map(szwrap, store.need)
        push!(out, :( local ($(sz_list...),) = ($(sizes...),) ) )
    end
    append!(out, store.mustassert) # NB do this after calling sizeinfer()
    out
end

function sizeinfer(store::NamedTuple, call::CallInfo)
    allowedcolons = :strided in call.flags ? 0 : 1

    sort!(unique!(store.need))
    sizes = Any[ (:) for i in store.need ]

    # First pass looks for single indices whose lengths are known directly
    for pair in store.dict
        if isa(pair.first, Symbol)
            d = findfirst(isequal(pair.first), store.need)
            d != nothing && (sizes[d] = pair.second)
        end
    end

    count(isequal(:), sizes) <= allowedcolons && return sizes

    # Second pass looks for tuples where exactly one entry has unknown length
    for pair in store.dict
        if isa(pair.first, Vector)
            known = [ haskey(store.dict, j) for j in pair.first ]

            if sum(.!known) == 1 # bingo! now work out its size:
                num = pair.second
                denfacts = [ store.dict[i] for i in pair.first[known] ]
                if length(denfacts) > 1
                    den = :( *($(denfacts...)) )
                else
                    den = :( $(denfacts[1]) )
                end
                rat = :( $num ÷ $den )

                i = pair.first[.!known][1]

                d = findfirst(isequal(i), store.need)
                d != nothing && (sizes[d] = rat)

                str = "inferring range of $i from range of $(join(pair.first, " ⊗ "))"
                push!(store.mustassert, :( TensorCast.@assert_ rem($num, $den)==0 $str) )
            end
        end
    end

    unknown = store.need[sizes .== (:)]
    str = join(unknown, ", ")
    length(unknown) <= allowedcolons || throw(MacroError("unable to infer ranges for indices $str", call))

    return sizes
end

"""
    maybestaticsizes([:, :, i], (:,:,*)) -> (:,:,*)
    maybestaticsizes([:3, :4, i], (:,:,*)) -> Size(3,4)
Produces the 2nd argument of `static_slice()`, for slicing `A{:3, :4, i}`.
"""
function maybestaticsizes(ijk::Vector, code::Tuple, call::CallInfo)
    iscodesorted(code) || throw(MacroError("static slices need all colons to the left, " *
            "got {$(pretty(ijk))} hence code = $(pretty(code))", call))
    length(ijk) == length(code) || error("wrong length of code!")
    staticsize = Any[ i.value for i in ijk if i isa QuoteNode ]
    if length(staticsize) == count(iscolon, ijk)
        return :( TensorCast.Size($(staticsize...)) ) # really StaticArrays.
    else
        return code
    end
end

"""
    maybestaticsizes([:i,:j,:k], (:,:,*), store) -> Size(3,4)
Produces the 2nd argument of `static_slice()`, using sizes from `store.dict` if available.
"""
function maybestaticsizes(ijk::Vector, code::Tuple, store::NamedTuple, call::CallInfo)
    iscodesorted(code) || throw(MacroError("static slices need all colons to the left, " *
            "got {$(pretty(ijk))} hence code = $(pretty(code))", call))
    length(ijk) == length(code) || error("wrong length of code!")
    staticsize = []
    for d=1:countcolons(code)
        if haskey(store.dict, ijk[d])
            push!(staticsize, store.dict[ijk[d]])
        else
            return code
        end
    end
    return :( TensorCast.Size($(staticsize...)) )
end

"""
    A = maybepush(ex, store, :name)
If `ex` is not just a symbol, then it pushes `:(Asym = ex)` into `store.main`
and returns `Asym`.
"""
maybepush(s::Symbol, any...) = s
function maybepush(ex::Expr, store::NamedTuple, name::Symbol=:A) # TODO make this look for same?
    Asym = gensym(name)
    push!(store.main, :( local $Asym = $ex ) )
    return Asym
end

tensorprimetidy(v::Vector) = Any[ tensorprimetidy(x) for x in v ]
function tensorprimetidy(ex)
    MacroTools.postwalk(ex) do x

        @capture(x, ((ij__,) \ k_) ) && return :( ($(ij...),$k) )
        @capture(x, i_ \ j_ ) && return :( ($i,$j) )

        @capture(x, ((ij__,) ⊗ k_) ) && return :( ($(ij...),$k) )
        @capture(x, i_ ⊗ j_ ) && return :( ($i,$j) )

        @capture(x, ((ij__,), k__) ) && return :( ($(ij...),$(k...)) )

        @capture(x, i_') && return Symbol(i,"′")
        x
    end
end

szwrap(i::Symbol) = Symbol(:sz_,i)
function szwrap(ijk::Vector)
    length(ijk) == 0 && return nothing
    length(ijk) == 1 && return szwrap(first(ijk))
    return :( TensorCast.star($([ Symbol(:sz_,i) for i in ijk ]...)) )
end

isconstant(n::Int) = true
isconstant(s::Symbol) = s == :_
isconstant(ex::Expr) = ex.head == :($)
isconstant(q::QuoteNode) = false

isindexing(s) = false
isindexing(ex::Expr) = @capture(x, A_[ijk__])

isCorI(i) = isconstant(i) || isindexing(ii)

isrange(ex::Expr) = @capture(ex, alpha_:omega_) # this may be a terrible idea
isrange(i) = false

isCorR(i) = isconstant(i) || isrange(i) # ranges are treated a bit like constants!

istensor(n::Int) = false
istensor(s::Symbol) = false
function istensor(ex::Expr)
    @capture(ex, i_' )     && return istensor(i)
    @capture(ex, -(ij__) ) && return length(ij)>1 # TODO maybe reject -(i,j) ... istensor(:( -i⊗j )) is OK without this line
    @capture(ex, i_⊗j_ )   && return true
    @capture(ex, i_\j_ )   && return true
    @capture(ex, (ij__,) ) && return length(ij)>1
end
istensor(::QuoteNode) = false

iscolon(s::Int) = false
iscolon(s::Symbol) = s == :(:)
iscolon(ex::Expr) = false
iscolon(q::QuoteNode) = true

containsindexing(s) = false
function containsindexing(ex::Expr)
    flag = false
    # MacroTools.postwalk(x -> @capture(x, A_[ijk__]) && (flag=true), ex)
    MacroTools.postwalk(ex) do x
        # @capture(x, A_[ijk__]) && !(all(isconstant, ijk)) && (flag=true)
        if @capture(x, A_[ijk__])
            # @show x ijk # TODO this is a bit broken?  @pretty @cast Z[i,j] := W[i] * exp(X[1][i] - X[2][j])
            flag=true
        end
    end
    flag
end

listindices(s::Symbol) = []
function listindices(ex::Expr)
    list = []
    MacroTools.postwalk(ex) do x
        if @capture(x, A_[ijk__])
            flat, _ = indexparse(nothing, ijk)
            push!(list, flat)
        end
        x
    end
    list
end

function guesstarget(ex::Expr)
    list = sort(listindices(ex), by=length, rev=true)
    shortlist = unique(reduce(vcat, list))
end

# function overlapsorted(x,y) # works fine but not in use yet
#     z = intersect(x,y)
#     length(z) ==0 && return true
#     xi = map(i -> findfirst(isequal(i),x), z)
#     yi = map(i -> findfirst(isequal(i),y), z)
#     return sortperm(xi) == sortperm(yi)
# end

"""
    needview!([:, 3, A])   # true, need view(A, :,3,:)
    needview!([:, :_, :])  # false, can use rview(A, :,1,:)

Mutates the given vector, replacing symbol `:_` with `1`.
If the vector contains only colons & underscores, then the result is suitable for use
with `rview`, but if not, we need a real view, so it returns `true`.
"""
function needview!(ij::Vector)
    out = false
    for k = 1:length(ij)
        if ij[k] == :_
            ij[k] = 1
        elseif ij[k] isa Int || ij[k] isa Symbol
            out = true
        elseif ij[k] isa Colon
            nothing
        elseif ij[k] isa Expr && ij[k].head == :($)
            ij[k] = ij[k].args[1]
            out = true
        elseif ij[k] isa Expr && @capture(ij[k], alpha_:omega_) # i.e. isrange(ij[k])
            out = true
        else
            error("this should never happen! needview! got ", ij[k])
        end
    end
    out
end

"""
    matrixshape(A, [i,j], [k,l]) -> reshape(A, sz_i * sz_j, sz_k * sz_l)
    matrixshape(A, [], [k,l]) -> reshape(A, 1, sz_k * sz_l) # rowvector * something
    matrixshape(A, [i,j], []) -> reshape(A, sz_i * sz_j)    # something * vector
"""
function matrixshape(ex, left::Vector, right::Vector, store::NamedTuple, call::CallInfo)
    # Deal with simple matrix, and with empty right of V in M*V
    length(left) == 1 && length(right) <= 1 && return ex
    # and empty right because it's M * vec(T)
    isempty(right) && return :( reshape($ex, :) )

    # Deal with empty left of V in V'*M, or perhaps V=vec(T) first
    if isempty(left)
        if length(right) == 1
            return :( TensorCast.PermuteDims($ex) )
        else
            return :( TensorCast.PermuteDims(reshape($ex, :)) )
        end
    end

    # Otherwise we are going to need to reshape
    left_sz = szwrap(left)
    right_sz = szwrap(right) # this is the product!
    append!(store.need, left)
    append!(store.need, right)
    # push!(call.flags, :reshaped)
    return :( reshape($ex, ($left_sz,$right_sz)) )
end

function unmatrixshape(ex, left::Vector, right::Vector, store::NamedTuple, call::CallInfo)
    # Easy cases M*M and M*V
    length(left) == 1 && length(right) <= 1 && return ex

    # For V' * V, did we want a scalar or not? What we will get is unknown to the macro:
    if length(left) == 0 && length(right) == 0
        if :scalar in call.flags
            return :( first($ex) )
            # If you had arrays of arrays, then PermuteDims would have permutedims-ed, and * would make an array, so this is still OK.
        else
            return :( Base.fill($ex) ) # zero-dim array
        end
    end

    # For V'*M, you may get a Transpose row-vector, for which this is .parent:
    if length(left) == 0
        ex = :( TensorCast.rview($ex, 1,:) )
        length(right) == 1 && return ex # literally V'*M done,
        # but for V'*T we should reshape this .parent
    end

    # For anything more complicated, we will need to reshape:
    sizes = vcat(map(szwrap, left), map(szwrap, right))
    append!(store.need, left)
    append!(store.need, right)
    # push!(call.flags, :reshaped)
    return :( reshape($ex, ($(sizes...),)) )
end

#==================== Nice Errors ====================#

"""
    @assert_ cond str

Like `@assert`, but prints both the given string and the condition.
Throws a `DimensionMismatch` error if `cond` is false.
"""
macro assert_(ex, str)
    msg = str * ": " * string(ex)
    return esc(:($ex || throw(DimensionMismatch($msg))))
end

# m_error(str::String) = @error str
# m_error(str::String, call::CallInfo) =
#     @error str  input=call.str _module=call.mod  _line=call.src.line  _file=string(call.src.file)

struct MacroError <: Exception
    msg::String
    call::Union{Nothing, CallInfo}
    MacroError(msg, call=nothing) = new(msg, call)
end

wherecalled(call::CallInfo) = "@ " * string(call.mod) * " " * string(call.src.file) * ":" * string(call.src.line)

function Base.showerror(io::IO, err::MacroError)
    print(io, err.msg)
    if err.call isa CallInfo
        printstyled(io, "\n    ", err.call.string; color = :blue)
        printstyled(io, "\n    ", wherecalled(err.call); color = :normal)
    end
end

function checknorepeats(flat, call=nothing, msg=nothing)
    msg == nothing && (msg = " in " * string(:( [$(flat...)] )))
    seen = Set{Symbol}()
    for i in flat
        (i in seen) ? throw(MacroError("index $i repeated" * msg, call)) : push!(seen, i)
    end
end

function checkallseen(canon, store, call)
    left = setdiff(canon, unique!(store.seen))
    length(left) > 0 && throw(MacroError("index $(left[1]) appears only on the left", call))
    right = setdiff(store.seen, canon)
    length(right) > 0 && throw(MacroError("index $(right[1]) appears only on the right", call))
end

# this may never be necessary with checkallseen?
function findcheck(i::Symbol, flat::Vector, call=nothing, msg=nothing)
    msg == nothing && (msg = " in " * string(:( [$(flat...)] )))
    res = findfirst(isequal(i), flat)
    res == nothing ? throw(MacroError("can't find index $i" * msg, call)) : return res
end
findcheck(i, flat, call=nothing, msg=nothing) =
    throw(MacroError("expected a single symol not $i", call))

#==================== The Output Functions ====================#

function newoutput(ex, canon, parsed, store::NamedTuple, call::CallInfo)
    isempty(parsed.reversed) || throw(MacroError("can't reverse along indices on LHS right now", call))
    isempty(parsed.shuffled) || throw(MacroError("can't shuffle along on LHS right now", call))
    any(iscolon, parsed.outer) && throw(MacroError("can't have colons on the LHS right now", call))
    any(isrange, parsed.outer) && throw(MacroError("can't have ranges on the LHS right now", call))

    # Is there a reduction?
    if :reduce in call.flags
        if :scalar in call.flags
            ex = :( $(parsed.redfun)($ex) )
        else
            dims = length(parsed.rdims)>1 ? Tuple(parsed.rdims) : parsed.rdims[1]
            ex = :( dropdims($(parsed.redfun)($ex, dims=$dims), dims=$dims) )
            if :strided in call.flags
                pop!(call.flags, :collected, :ok) # dropdims makes reshape(stridedview(...
            end
        end
        canon = deleteat!(copy(canon), sort(parsed.rdims))
    end

    # Were we asked to slice the output?
    if length(parsed.inner) != 0
        code = Tuple(Any[ i in parsed.innerflat ? (:) : (*) for i in canon ])
        if parsed.static
            sizeorcode = maybestaticsizes(canon, code, store, call)
            ex = :( TensorCast.static_slice($ex, $sizeorcode) )
        elseif :collect in call.flags
            ex = :( TensorCast.slicecopy($ex, $code) )
            push!(call.flags, :collected)
        else
            ex = :( TensorCast.sliceview($ex, $code) )
        end
        # Now I allow fixing output indices
        if any(isconstant, parsed.inner)
            any(i -> isconstant(i) && !(i == :_ || i == 1), parsed.inner) && throw(MacroError("can't fix output index to something other than 1", call))
            code = Tuple(map(i -> isconstant(i) ? (*) : (:), parsed.inner))
            Asafe = maybepush(ex, store, :outfix)
            ex = :(TensorCast.orient.($Asafe, Ref($code)) ) # @. would need a dollar
        end
    end

    # Must we collect? Do this now, as reshape(PermutedDimsArray(...)) is awful.
    if :collect in call.flags && !(:collected in call.flags)
        ex = :( identity.($ex) )
    end

    # Do we need to reshape the container? Using orient() avoids needing sz_i
    if any(i -> istensor(i) || isconstant(i), parsed.outer)
        any(i -> isconstant(i) && !(i == :_ || i == 1), parsed.outer) && throw(MacroError("can't fix output index to $i, only to 1", call))
        if any(istensor, parsed.outer)
            ex = :( reshape($ex, ($(parsed.outsize...),)) )
            append!(store.need, parsed.flat)
            # push!(call.flags, :reshaped)
        else
            code = Tuple(map(i -> isconstant(i) ? (*) : (:), parsed.outer))
            ex = :( TensorCast.orient($ex, $code) )
        end
    end

    # Is the result Diagonal or friends? Doesn't allow Z[i,i,1] or Z[i,-i] but that's OK
    if length(parsed.outer)==2 && parsed.outer[1]==parsed.outer[2]
        if :nolazy in call.flags
            ex = :( TensorCast.diagm(0 => $ex) )
        else
            ex = :( TensorCast.Diagonal($ex) )
        end
    end

    return ex
end

function inplaceoutput(ex, canon, parsed, store::NamedTuple, call::CallInfo)
    if !(parsed.static)
        isempty(parsed.inner) || throw(MacroError("can't write in-place into slices right now",call))
        any(iscolon, parsed.outer) && throw(MacroError("can't have colons on the LHS right now",call))
    end

    out = []
    Zsym = gensym(:reverse)

    # We can re-use exactly the same "standardise" function as for RHS terms:
    pop!(call.flags, :nolazy, :ok) # ensure we use diagview(), Reverse{}, etc, not a copy

    if @capture(parsed.left, zed_[]) # special case Z[] = ... else allconst pulls it out
        zed isa Symbol || @capture(zed, ZZ_.field_) || error("wtf")
        newleft = parsed.left
        str = "expected a 0-tensor $zed[]"
        push!(store.mustassert, :( TensorCast.@assert_ ndims($zed)==0 $str) )
    else
        newleft = standardise(parsed.left, store, call)
        @capture(newleft, zed_[ijk__]) || throw(MacroError("failed to parse LHS correctly, $(parsed.left) -> $newleft"))

        if !(zed isa Symbol) # then standardise did something!
            push!(call.flags, :showfinal)
            Zsym = gensym(:reverse)
            push!(out, :( local $Zsym = $zed ) )
            zed = Zsym
        end
    end

    # Now write into that, either sum!(Z,...) or mul!(Z,A,B) or Z .= ..., no copyto!(Z,...)
    if :reduce in call.flags
        redfun! = endswith(string(parsed.redfun),'!') ? parsed.redfun : Symbol(parsed.redfun, '!')
        push!(out, :( $redfun!($zed, $ex) ) )
    elseif :matmul in call.flags
        ex isa Tuple || error("wtf?")

        zmul = targetcast(newleft, vcat(ex[3], ex[4]), store, call)
        zmul = matrixshape(zmul, ex[3], ex[4], store, call)
        zmul isa Symbol || push!(call.flags, :showfinal)
        # zed = maybepush(zed, out, :mul!)
        push!(out, :( TensorCast.mul!($zmul, $(ex[1]), $(ex[2])) ) )
    else
        push!(out, :( $zed .= $ex ) )
    end

    if :showfinal in call.flags
        push!(out, parsed.name)
    end

    return out
end

#==================== The End ====================#


text:
using GitTest
using Test

@testset "testfunctionderivatives.jl" begin
    # Derivatives
    @test dfdx(1.23,-93.2) == 2
end

ext: jl


text:
# imports
include("axelrodnet.jl")

# replication of the original experiment
orig_replication = Axelrodnet.run!(
    agentcount=100,
    n_iter=100,
    nettopology=1,
    networkprops=Dict("grid_height" => 10),
    socialbotfrac=0.0,
    rndseed=151,
    repcount=3,
    export_every_n=1
)

Axelrodnet.export_experiment(orig_replication)

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
const TEST_NAME = "Xoroshiro128Plus"

include("common.jl")

using RNG.Xorshifts

r = Xoroshiro128Plus(123)

test_all(r, 100_000_000)


text:
#
# JSON client with Julia
#
# author: Atsushi Sakai
#

using Requests

function main()

	req = Dict("Name" => "julia json request")
	println("request:")
	println(req)

    url = "http://localhost:8000"
    res = Requests.json(Requests.post(url; json = req))
    println("response:")
    println(res)

    println(PROGRAM_FILE," Done!!")
end

main()

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
"""
    unsigned(T)
    unsigned(x::T)

Obtain the unsigned type or unsigned value that corresponds.
""" unsigned

"""
    signed(T)
    signed(x::T)

Obtain the signed type or signed value that corresponds.
""" signed

"""
    ieeefloat(T)
    ieeefloat(x::T)

Obtain the IEEEFloat type or float value with the same number of bits as T.
""" ieeefloat

Base.unsigned(::Type{Float128}) = UInt128
Base.unsigned(::Type{Float64}) = UInt64
Base.unsigned(::Type{Float32}) = UInt32
Base.unsigned(::Type{Float16}) = UInt16
Base.unsigned(::Type{Float8}) = UInt8

Base.unsigned(x::Float128) = reinterpret(UInt128, x)
Base.unsigned(x::Float64) = reinterpret(UInt64, x)
Base.unsigned(x::Float32) = reinterpret(UInt32, x)
Base.unsigned(x::Float16) = reinterpret(UInt16, x)
Base.unsigned(x::Float8) = reinterpret(UInt8, x)

Base.signed(::Type{Float128}) = Int128
Base.signed(::Type{Float64}) = Int64
Base.signed(::Type{Float32}) = Int32
Base.signed(::Type{Float16}) = Int16
Base.signed(::Type{Float8}) = Int8

Base.signed(x::Float128) = reinterpret(Int128, x)
Base.signed(x::Float64) = reinterpret(Int64, x)
Base.signed(x::Float32) = reinterpret(Int32, x)
Base.signed(x::Float16) = reinterpret(Int16, x)
Base.signed(x::Float8) = reinterpret(Int8, x)

ieeefloat(::Type{UInt128}) = Float128
ieeefloat(::Type{UInt64}) = Float64
ieeefloat(::Type{UInt32}) = Float32
ieeefloat(::Type{UInt16}) = Float16
ieeefloat(::Type{UInt8}) = Float8

ieeefloat(x::UInt128) = reinterpret(Float128, x)
ieeefloat(x::UInt64) = reinterpret(Float64, x)
ieeefloat(x::UInt32) = reinterpret(Float32, x)
ieeefloat(x::UInt16) = reinterpret(Float16, x)
ieeefloat(x::UInt8) = reinterpret(Float8, x)

ieeefloat(::Type{Int128}) = Float128
ieeefloat(::Type{Int64}) = Float64
ieeefloat(::Type{Int32}) = Float32
ieeefloat(::Type{Int16}) = Float16
ieeefloat(::Type{Int8}) = Float8

ieeefloat(x::Int128) = reinterpret(Float128, x)
ieeefloat(x::Int64) = reinterpret(Float64, x)
ieeefloat(x::Int32) = reinterpret(Float32, x)
ieeefloat(x::Int16) = reinterpret(Float16, x)
ieeefloat(x::Int8) = reinterpret(Float8, x)

bybitwidth = Dict{Int, NamedTuple{(:unsigned, :signed, :floating)}}(
    8   => (unsigned=UInt8,   signed=Int8,   floating=isdefined(Main, :Float8) ? Float8 : UInt8),
    16  => (unsigned=UInt16,  signed=Int16,  floating=Float16),
    32  => (unsigned=UInt16,  signed=Int16,  floating=Float32),
    64  => (unsigned=UInt16,  signed=Int16,  floating=Float64),
    128 => (unsigned=UInt128, signed=Int128, floating=isdefined(Main, :Quadmath) ? Float128 : UInt128)
)


text:
using LinearAlgebra, TensorRefinement, Test

function refextdn(L::Int, ℓ::Int, d::Int)
	I = diagm(0 => ones(2^ℓ))
	if ℓ > 0
		S = diagm(-1 => ones(2^ℓ-1))
	else
		S = [0.0]
	end
	η = collect(1:2^(L-ℓ))/2^(L-ℓ)
	ζ = collect(2^(L-ℓ)-1:-1:0)/2^(L-ℓ)
	P = kron(I, η) + kron(S, ζ)
	P .*= 2.0^((ℓ-L)/2)
	PP = P
	for k ∈ 2:d
		PP = kron(P, PP)
	end
	return PP
end

function refdiffdn(ℓ::Int, d::Int, K::Int)
	I = diagm(0 => ones(2^ℓ))
	if ℓ > 0
		S = diagm(-1 => ones(2^ℓ-1))
	else
		S = [0.0]
	end
	M0 = 2.0^(ℓ/2-1) * ( kron(I+S, [1; 0]) + kron(I-S, [0; 1]) )
	M1 = 2.0^(ℓ/2-1+(ℓ+1)) * (I-S)
	M = (K == 1) ? M1 : M0
	for κ ∈ 2:d
		M = kron((κ == K) ? M1 : M0, M)
	end
	return M
end

function refintdn(ℓ::Int, d::Int, K1::Int, K2::Int)
	Λ00 = kron(diagm(0 => ones(2^ℓ))/2^ℓ, diagm(0 => [1,1/3]))
	Λ11 = kron(diagm(0 => ones(2^ℓ))/2^ℓ, [1.0])
	Λ01 = kron(diagm(0 => ones(2^ℓ))/2^ℓ, [1.0; 0.0])
	Λ10 = kron(diagm(0 => ones(2^ℓ))/2^ℓ, [1.0 0.0])
	κ = 1
	if κ == K1 && κ == K2
		Λ = Λ11
	elseif κ == K1 && κ ≠ K2
		Λ = Λ10
	elseif κ ≠ K1 && κ == K2
		Λ = Λ01
	else
		Λ = Λ00
	end
	for κ ∈ 2:d
		if κ == K1 && κ == K2
			Λ = kron(Λ11, Λ)
		elseif κ == K1 && κ ≠ K2
			Λ = kron(Λ10, Λ)
		elseif κ ≠ K1 && κ == K2
			Λ = kron(Λ01, Λ)
		else
			Λ = kron(Λ00, Λ)
		end
	end
	return Λ
end

function refbpxdn(L::Int, d::Int)
	C = zeros(2^(d*L), 2^(d*L))
	for ℓ ∈ 0:L
		P = refextdn(L, ℓ, d)
		C = C + 2.0^(-ℓ) * P*P'
	end
	return C
end

@testset "extdn" begin
	@testset "d = $d, ℓ = $ℓ, L = $L" for (d,L,ℓ) ∈ ((1,7,0),(1,7,4),(1,7,7),(2,5,0),(2,5,3),(2,5,5),(3,4,0),(3,4,3),(3,4,4))
		P = FEM.extdn(L, ℓ, d; major="first")
		P = [ [ reshape(P[k], size(P[k], 1), 2*ones(Int, 2*d)..., size(P[k], 4)) for k ∈ 1:ℓ ]...,
			  [ reshape(P[k], size(P[k], 1), 2*ones(Int, d)..., ones(Int, d)..., size(P[k], 4)) for k ∈ ℓ+1:L ]... ]
		Pf = block(P, 1, 1; major="first")
		Pf = reshape(Pf, 2^(d*L), 2^(d*ℓ))
		P0f = refextdn(L, ℓ, d)
		@test Pf ≈ P0f rtol=1e-14
	end
end

@testset "diffdn" begin
	@testset "d = $d, ℓ = $ℓ" for (d,ℓ) ∈ ((1,7),(2,5),(2,3))
		@testset "K = $K" for K ∈ 0:d
			n = 2*ones(Int, d)
			if K ≠ 0
				n[K] = 1
			end
			M = FEM.diffdn(ℓ, d, K; major="first")
			M = [ [ reshape(M[k], size(M[k], 1), 2*ones(Int, 2*d)..., size(M[k], 4)) for k ∈ 1:ℓ ]...,
				  reshape(M[ℓ+1], size(M[ℓ+1], 1), n..., ones(Int, d)..., size(M[ℓ+1], 4)) ]
			Mf = block(M, 1, 1; major="first")
			Mf = reshape(Mf, 2^(d*ℓ)*prod(n), 2^(d*ℓ))
			M0f = refdiffdn(ℓ, d, K)
			@test Mf ≈ M0f rtol=1e-14
		end
	end
end

@testset "intdn" begin
	@testset "d = $d, ℓ = $ℓ" for (d,ℓ) ∈ ((1,7),(2,4),(3,3))
		@testset "K1 = $K1, K2 = $K2" for K1 ∈ 0:d, K2 ∈ 0:d
			n1 = 2*ones(Int, d)
			n2 = 2*ones(Int, d)
			if K1 ≠ 0
				n1[K1] = 1
			end
			if K2 ≠ 0
				n2[K2] = 1
			end
			Λ = FEM.intdn(ℓ, d, K1, K2; major="first")
			Λ = [ [ reshape(Λ[k], size(Λ[k], 1), 2*ones(Int, 2*d)..., size(Λ[k], 4)) for k ∈ 1:ℓ ]...,
				  reshape(Λ[ℓ+1], size(Λ[ℓ+1], 1), n1..., n2..., size(Λ[ℓ+1], 4)) ]
			Λf = block(Λ, 1, 1; major="first")
			Λf = reshape(Λf, 2^(d*ℓ)*prod(n1), 2^(d*ℓ)*prod(n2))
			Λ0f = refintdn(ℓ, d, K1, K2)
			@test Λf ≈ Λ0f rtol=1e-14
		end
	end
end

@testset "bpxdn" begin
	@testset "d = $d, ℓ = $ℓ" for (d,ℓ) ∈ ((1,7),(2,5),(3,4))
		C = FEM.bpxdn(ℓ, d; major="first")
		C = [ reshape(C[k], size(C[k], 1), 2*ones(Int, 2*d)..., size(C[k], 4)) for k ∈ 1:ℓ ]
		Cf = block(C, 1, 1; major="first")
		Cf = reshape(Cf, 2^(d*ℓ), 2^(d*ℓ))
		C0f = refbpxdn(ℓ, d)
		@test Cf ≈ C0f rtol=1e-14
	end
end

@testset "diffextdn" begin
	@testset "d = $d, ℓ = $ℓ, L = $L" for (d,ℓ,L) ∈ ((1,0,7),(1,1,7),(1,3,7),(1,7,7),(2,3,5),(3,2,4))
		@testset "K = $K" for K ∈ 0:d
			n = 2*ones(Int, d)
			if K ≠ 0
				n[K] = 1
			end
			MP = FEM.diffextdn(L, ℓ, d, K; major="first")
			MP = [ [ reshape(MP[k], size(MP[k], 1), 2*ones(Int, 2*d)..., size(MP[k], 4)) for k ∈ 1:ℓ ]...,
				   [ reshape(MP[k], size(MP[k], 1), 2*ones(Int, d)..., ones(Int, d)..., size(MP[k], 4)) for k ∈ ℓ+1:L ]...,
				  reshape(MP[L+1], size(MP[L+1], 1), n..., ones(Int, d)..., size(MP[L+1], 4)) ]
			MPf = block(MP, 1, 1; major="first")
			MPf = reshape(MPf, 2^(d*L)*prod(n), 2^(d*ℓ))
			MP0f = refdiffdn(L, d, K)*refextdn(L, ℓ, d)
			@test MPf ≈ MP0f rtol=1e-14
		end
	end
end

@testset "diffbpxdn" begin
	@testset "d = $d, ℓ = $ℓ" for (d,ℓ) ∈ ((1,1),(1,7),(2,1),(2,4))
		@testset "K = $K" for K ∈ 0:d
			n = 2*ones(Int, d)
			if K ≠ 0
				n[K] = 1
			end
			Q = FEM.diffbpxdn(ℓ, d, K; major="first")
			Q = [ [ reshape(Q[k], size(Q[k], 1), 2*ones(Int, 2*d)..., size(Q[k], 4)) for k ∈ 1:ℓ ]...,
				  reshape(Q[ℓ+1], size(Q[ℓ+1], 1), n..., ones(Int, d)..., size(Q[ℓ+1], 4)) ]
			Qf = block(Q, 1, 1; major="first")
			Qf = reshape(Qf, 2^(d*ℓ)*prod(n), 2^(d*ℓ))
			M0f = refdiffdn(ℓ, d, K)
			C0f = refbpxdn(ℓ, d)

			Q0f = M0f*C0f
			@test Qf ≈ Q0f rtol=1e-14
		end
	end
end

@testset "stiffness and mass matrices" begin
	@testset "d = 1, ℓ = $ℓ" for ℓ ∈ (1,8,9)
		d = 1
		M0 = FEM.diffdn(ℓ, d, 0; major="first")
		Λ0 = FEM.intdn(ℓ, d, 0, 0; major="first")
		M = decmp(M0, 1, decmp(Λ0, 2, M0, 1), 1)
		decskp!(M, ℓ+1; path="backward")
		Mf = block(M, 1, 1; major="first")

		M1 = FEM.diffdn(ℓ, d, 1; major="first")
		Λ1 = FEM.intdn(ℓ, d, 1, 1; major="first")
		S = decmp(M1, 1, decmp(Λ1, 2, M1, 1), 1)
		decskp!(S, ℓ+1; path="backward")
		Sf = block(S, 1, 1; major="first")

		n = 2^ℓ
		Sf0 = 2^(2*ℓ) * diagm(1 => -ones(Float64, n-1), -1 => -ones(Float64, n-1), 0 => [2*ones(Float64, n-1); 1])
		Mf0 = diagm(1 => 1/6*ones(Float64, n-1), -1 => 1/6*ones(Float64, n-1), 0 => 1/3*[2*ones(Float64, n-1); 1])
		@testset "mass matrix" begin
			@test Mf ≈ Mf0 rtol=1e-14
		end
		@testset "stiffness matrix" begin
			@test Sf ≈ Sf0 rtol=1e-14
		end
	end
	@testset "d = 2, ℓ = $ℓ" for ℓ ∈ (1,4,5)
		d = 2
		M0 = FEM.diffdn(ℓ, d, 0; major="first")
		Λ0 = FEM.intdn(ℓ, d, 0, 0; major="first")
		M = decmp(M0, 1, decmp(Λ0, 2, M0, 1), 1)
		decskp!(M, ℓ+1; path="backward")
		M = [ reshape(M[k], size(M[k], 1), 2*ones(Int, 2*d)..., size(M[k], 4)) for k ∈ 1:ℓ ]
		Mf = block(M, 1, 1; major="first")
		Mf = reshape(Mf, 2^(d*ℓ), 2^(d*ℓ))

		M1 = FEM.diffdn(ℓ, d, 1; major="first")
		Λ1 = FEM.intdn(ℓ, d, 1, 1; major="first")
		M2 = FEM.diffdn(ℓ, d, 2; major="first")
		Λ2 = FEM.intdn(ℓ, d, 2, 2; major="first")
		S1 = decmp(M1, 1, decmp(Λ1, 2, M1, 1), 1)
		decskp!(S1, ℓ+1; path="backward")
		S2 = decmp(M2, 1, decmp(Λ2, 2, M2, 1), 1)
		decskp!(S2, ℓ+1; path="backward")
		S = decadd(S1, S2)

		S1 = [ reshape(S1[k], size(S1[k], 1), 2*ones(Int, 2*d)..., size(S1[k], 4)) for k ∈ 1:ℓ ]
		S2 = [ reshape(S2[k], size(S2[k], 1), 2*ones(Int, 2*d)..., size(S2[k], 4)) for k ∈ 1:ℓ ]
		S  = [ reshape(S[k], size(S[k], 1), 2*ones(Int, 2*d)..., size(S[k], 4)) for k ∈ 1:ℓ ]

		S1f = block(S1, 1, 1; major="first")
		S1f = reshape(S1f, 2^(d*ℓ), 2^(d*ℓ))
		S2f = block(S2, 1, 1; major="first")
		S2f = reshape(S2f, 2^(d*ℓ), 2^(d*ℓ))

		Sf = block(S, 1, 1; major="first")
		Sf = reshape(Sf, 2^(d*ℓ), 2^(d*ℓ))

		n = 2^ℓ
		Sf0 = 2^(2*ℓ) * diagm(1 => -ones(Float64, n-1), -1 => -ones(Float64, n-1), 0 => [2*ones(Float64, n-1); 1])
		Mf0 = diagm(1 => 1/6*ones(Float64, n-1), -1 => 1/6*ones(Float64, n-1), 0 => 1/3*[2*ones(Float64, n-1); 1])
		S1f0 = kron(Mf0, Sf0)
		S2f0 = kron(Sf0, Mf0)
		Mf0 = kron(Mf0, Mf0)
		Sf0 = S1f0 + S2f0
		@testset "mass matrix" begin
			@test Mf ≈ Mf0 rtol=1e-14
		end
		@testset "stiffness matrix" begin
			@test S1f ≈ S1f0 rtol=1e-14
			@test S2f ≈ S2f0 rtol=1e-14
			@test Sf ≈ Sf0 rtol=1e-14
		end
	end
end

@testset "preconditioned stiffness and mass matrices" begin
	@testset "d = 1, ℓ = $ℓ" for ℓ ∈ (1,8,9)
		d = 1
		M0 = FEM.diffbpxdn(ℓ, d, 0; major="first")
		Λ0 = FEM.intdn(ℓ, d, 0, 0; major="first")
		M = decmp(M0, 1, decmp(Λ0, 2, M0, 1), 1)
		decskp!(M, ℓ+1; path="backward")
		Mf = block(M, 1, 1; major="first")

		M1 = FEM.diffbpxdn(ℓ, d, 1; major="first")
		Λ1 = FEM.intdn(ℓ, d, 1, 1; major="first")
		S = decmp(M1, 1, decmp(Λ1, 2, M1, 1), 1)
		decskp!(S, ℓ+1; path="backward")
		Sf = block(S, 1, 1; major="first")

		n = 2^ℓ
		C0f = refbpxdn(ℓ, d)
		Sf0 = 2^(2*ℓ) * diagm(1 => -ones(Float64, n-1), -1 => -ones(Float64, n-1), 0 => [2*ones(Float64, n-1); 1])
		Mf0 = diagm(1 => 1/6*ones(Float64, n-1), -1 => 1/6*ones(Float64, n-1), 0 => 1/3*[2*ones(Float64, n-1); 1])
		Mf0 = C0f*Mf0*C0f
		Sf0 = C0f*Sf0*C0f
		@testset "mass matrix" begin
			@test Mf ≈ Mf0 rtol=1e-14
		end
		@testset "stiffness matrix" begin
			@test Sf ≈ Sf0 rtol=1e-14
		end
	end
	@testset "d = 2, ℓ = $ℓ" for ℓ ∈ (1,3,4)
		d = 2
		M0 = FEM.diffbpxdn(ℓ, d, 0; major="first")
		Λ0 = FEM.intdn(ℓ, d, 0, 0; major="first")
		M = decmp(M0, 1, decmp(Λ0, 2, M0, 1), 1)
		decskp!(M, ℓ+1; path="backward")
		M = [ reshape(M[k], size(M[k], 1), 2*ones(Int, 2*d)..., size(M[k], 4)) for k ∈ 1:ℓ ]
		Mf = block(M, 1, 1; major="first")
		Mf = reshape(Mf, 2^(d*ℓ), 2^(d*ℓ))

		M1 = FEM.diffbpxdn(ℓ, d, 1; major="first")
		Λ1 = FEM.intdn(ℓ, d, 1, 1; major="first")
		M2 = FEM.diffbpxdn(ℓ, d, 2; major="first")
		Λ2 = FEM.intdn(ℓ, d, 2, 2; major="first")
		S1 = decmp(M1, 1, decmp(Λ1, 2, M1, 1), 1)
		decskp!(S1, ℓ+1; path="backward")
		S2 = decmp(M2, 1, decmp(Λ2, 2, M2, 1), 1)
		decskp!(S2, ℓ+1; path="backward")
		S = decadd(S1, S2)

		S1 = [ reshape(S1[k], size(S1[k], 1), 2*ones(Int, 2*d)..., size(S1[k], 4)) for k ∈ 1:ℓ ]
		S2 = [ reshape(S2[k], size(S2[k], 1), 2*ones(Int, 2*d)..., size(S2[k], 4)) for k ∈ 1:ℓ ]
		S  = [ reshape(S[k], size(S[k], 1), 2*ones(Int, 2*d)..., size(S[k], 4)) for k ∈ 1:ℓ ]

		S1f = block(S1, 1, 1; major="first")
		S1f = reshape(S1f, 2^(d*ℓ), 2^(d*ℓ))
		S2f = block(S2, 1, 1; major="first")
		S2f = reshape(S2f, 2^(d*ℓ), 2^(d*ℓ))

		Sf = block(S, 1, 1; major="first")
		Sf = reshape(Sf, 2^(d*ℓ), 2^(d*ℓ))

		n = 2^ℓ
		C0f = refbpxdn(ℓ, d)
		Sf0 = 2^(2*ℓ) * diagm(1 => -ones(Float64, n-1), -1 => -ones(Float64, n-1), 0 => [2*ones(Float64, n-1); 1])
		Mf0 = diagm(1 => 1/6*ones(Float64, n-1), -1 => 1/6*ones(Float64, n-1), 0 => 1/3*[2*ones(Float64, n-1); 1])
		S1f0 = C0f*kron(Mf0, Sf0)*C0f
		S2f0 = C0f*kron(Sf0, Mf0)*C0f
		Mf0 = C0f*kron(Mf0, Mf0)*C0f
		Sf0 = S1f0 + S2f0
		@testset "mass matrix" begin
			@test Mf ≈ Mf0 rtol=1e-14
		end
		@testset "stiffness matrix" begin
			@test S1f ≈ S1f0 rtol=1e-14
			@test S2f ≈ S2f0 rtol=1e-14
			@test Sf ≈ Sf0 rtol=1e-14
		end
	end
end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets
using CSDP

# -----------------------------------------------------------------------------
# Trivial tests
# -----------------------------------------------------------------------------

@testset "Guessing game --------- " begin
	@polyvar x y
	p  = (x - y)^2
	Sx = @set 1 - x^2 >= 0
	Sy = @set 1 - y^2 >= 0
	opt = CSDP.Optimizer

	expected_measure_x = AtomicMeasure(variables(Sx), [
		WeightedDiracMeasure([1], 0.5),
		WeightedDiracMeasure([-1], 0.5)])
	expected_measure_y = AtomicMeasure(variables(Sy), [
		WeightedDiracMeasure([0], 1)])
	expected_value = 1

	value, measure_x, measure_y = solve_game(p, Sx, Sy, opt, iteration=0x1)

	eps = 1e-3
	@test isapprox(expected_value, value, atol=eps)
	@test isapprox(expected_measure_x, measure_x, atol=eps)
	@test isapprox(expected_measure_y, measure_y, atol=eps)
end

# -----------------------------------------------------------------------------
# P. A. Parrilo (2006). Polynomial games and sum of squares optimization.
# In Proceedings of the 45th IEEE Conference on Decision and Control
# (pp. 2855-2860).
# -----------------------------------------------------------------------------

@testset "Parrilo Example 3.1 --- " begin
	@polyvar x y
	p  = 2 * x * y^2 - x^2 - y
	Sx = @set 1 - x^2 >= 0
	Sy = @set 1 - y^2 >= 0
	opt = CSDP.Optimizer

	expected_measure_x = AtomicMeasure(variables(Sx), [
		WeightedDiracMeasure([0.3968], 1)])
	expected_measure_y = AtomicMeasure(variables(Sy), [
		WeightedDiracMeasure([0.6299], 1)])
	expected_value = -0.4724

	value, measure_x, measure_y = solve_game(p, Sx, Sy, opt)

	eps = 1e-3
	@test isapprox(expected_value, value, atol=eps)
	@test isapprox(expected_measure_x, measure_x, atol=eps)
	@test isapprox(expected_measure_y, measure_y, atol=eps)
end

@testset "Parrilo Example 3.2 --- " begin
	@polyvar x y
	p  = 5x * y - 2x^2 - 2x * y^2 - y
	Sx = @set 1 - x^2 >= 0
	Sy = @set 1 - y^2 >= 0
	opt = CSDP.Optimizer

	expected_measure_x = AtomicMeasure(variables(Sx), [
		WeightedDiracMeasure([0.2], 1)])
	expected_measure_y = AtomicMeasure(variables(Sy), [
		WeightedDiracMeasure([ 1], 0.78),
		WeightedDiracMeasure([-1], 0.22)])
	expected_value = -0.48

	value, measure_x, measure_y = solve_game(p, Sx, Sy, opt, iteration=0x1)

	eps = 1e-3
	@test isapprox(expected_value, value, atol=eps)
	@test isapprox(expected_measure_x, measure_x, atol=eps)
	@test isapprox(expected_measure_y, measure_y, atol=eps)
end


text:
import PotentialUQ
using LinearAlgebra: I
using Random
println("Starting test of MAP method")
Random.seed!(1234);

n = 100
β = zeros(10)
A = 1e-1*randn(n, 10)
A[1:10, 1:10] += I(10)
b = randn(n)
β = A\b
σ = 1e-3 .* ones(n)
logσ = log.(σ)
snap = PotentialUQ.Potentials.SNAP(β, 0.9, 6, 1)
dsnap = PotentialUQ.SNAPQDistribution(snap, A, b, logσ)
println(dsnap( (β = β, σ = logσ)))
println( dsnap( (β = A\b, σ = logσ)) )

@time PotentialUQ.MAP(dsnap)
println(dsnap.x)

snap = PotentialUQ.Potentials.SNAP(β, 0.9, 6, 1)
dsnap = PotentialUQ.SNAPQDistribution(snap, A, b, logσ)
println(dsnap.x)
@time PotentialUQ.MAP(dsnap)

println("Solution ", A\b)
println("MAP Solution, Optimization β =  ")
show(stdout, "text/plain", dsnap.x[1])
println("MAP Solution, Optimization σ =  ")
show(stdout, "text/plain", exp.(dsnap.x[2]))
println(" ")

println("End of test.")


text:
# This AIBECS test is derived from the P-model tutorial in AIBECS.jl
# But it uses one of the samll circulations to test derivatives of the cost function
module SubModuleAIBECS

using Test
using F1Method
using AIBECS
import AIBECS: @units, units
import AIBECS: @limits, limits
import AIBECS: @initial_value, initial_value
import AIBECS: @flattenable, flattenable, flatten
import AIBECS: @prior, prior
using Unitful: m, d, s, yr, Myr, mol, mmol, μmol, μM
using Distributions
using WorldOceanAtlasTools
using FiniteDiff

# AIBECS model
grd, T_Circ = Primeau_2x2x2.load()
T_DIP(p) = T_Circ
T_POP(p) = transportoperator(grd, z -> w(z,p))
function w(z,p)
    @unpack w₀, w′ = p
    w₀ + w′ * z
end
z_top = topdepthvec(grd) # uptake only in top layer
function U(x,p)
    @unpack σ, τ_DIP, k = p
    return @. σ * x/τ_DIP * x/(x+k) * (z_top==0) * (x≥0)
end
function R(x,p)
    @unpack τ_POP = p
    return x / τ_POP
end
function G_DIP(DIP, POP, p)
    @unpack DIP_geo, τ_geo = p
    return @. -$U(DIP,p) + $R(POP,p) + (DIP_geo - DIP) / τ_geo
end
function G_POP(DIP, POP, p)
    @unpack τ_geo = p
    return @. $U(DIP,p) - $R(POP,p) - POP / τ_geo
end

∞ = Inf
@initial_value @units @flattenable @limits struct PmodelParameters{U} <: AbstractParameters{U}
    w₀::U       |  0.64 | m/d      | true  | (0,∞)
    w′::U       |  0.13 | m/d/m    | true  | (0,∞)
    τ_DIP::U    | 230.0 | d        | true  | (0,∞)
    k::U        |  6.62 | μmol/m^3 | true  | (0,∞)
    τ_POP::U    |   5.0 | d        | true  | (0,∞)
    τ_geo::U    |   1.0 | Myr      | false | (0,∞)
    DIP_geo::U  |  2.12 | mmol/m^3 | true  | (-∞,∞)
    σ::U        |  0.3  | NoUnits  | true  | (0,1)
end
function prior(::Type{T}, s::Symbol) where {T<:AbstractParameters}
    if flattenable(T, s)
        lb, ub = limits(T, s)
        if (lb, ub) == (0,∞)
            μ = log(initial_value(T, s))
            LogNormal(μ, 1.0)
        elseif (lb, ub) == (-∞,∞)
            μ = initial_value(T, s)
            σ = 10.0 # Assumes that a sensible unit is chosen (i.e., that within 10.0 * U)
            Distributions.Normal(μ, σ)
        else # LogitNormal with median as initial value and bounds
            m = initial_value(T, s)
            f = (m - lb) / (ub - lb)
            LocationScale(lb, ub - lb, LogitNormal(log(f/(1-f)), 1.0))
        end
    else
        nothing
    end
end
prior(::T, s::Symbol) where {T<:AbstractParameters} = prior(T,s)
prior(::Type{T}) where {T<:AbstractParameters} = Tuple(prior(T,s) for s in AIBECS.symbols(T))
prior(::T) where {T<:AbstractParameters} = prior(T)
p = PmodelParameters()
λ = p2λ(p)
nb = sum(iswet(grd))
F = AIBECSFunction((T_DIP, T_POP), (G_DIP, G_POP), nb, PmodelParameters)
@unpack DIP_geo = p
x = DIP_geo * ones(2nb) # initial guess
prob = SteadyStateProblem(F, x, p)
sol = solve(prob, CTKAlg())

# Observations from World Ocean Atlas used to evaluate the
# AIBECS model mismatch with observations
# and generate the objective function
ρSW = 1.035u"kg/L" # approximate mean sea water density
obs = let
    obs = WorldOceanAtlasTools.observations("phosphate")
    obs.value = ustrip.(upreferred.(obs.phosphate * ρSW))
    (obs,)
end
modify(DIP, POP) = (DIP,)
ωs = (1.0,) # the weight for the mismatch (weight of POP = 0)
ωp = 1.0       # the weight for the parameters prior estimates
f, ∇ₓf = f_and_∇ₓf(ωs, ωp, grd, modify, obs, PmodelParameters)

# Now we apply the F1 method
τ = ustrip(u"s", 1e3u"Myr")
mem = F1Method.initialize_mem(F, ∇ₓf, x, λ, CTKAlg(), τstop=τ)
objective(λ) = F1Method.objective(f, F, mem, λ, CTKAlg(), τstop=τ)
gradient(λ) = F1Method.gradient(f, F, ∇ₓf, mem, λ, CTKAlg(), τstop=τ)
hessian(λ) = F1Method.hessian(f, F, ∇ₓf, mem, λ, CTKAlg(), τstop=τ)

# Finally we test the result with the "reliable" FiniteDiff :)
@test FiniteDiff.finite_difference_gradient(objective, 2λ) ≈ gradient(2λ)' rtol=1e-3
@test FiniteDiff.finite_difference_hessian(objective, 2λ) ≈ hessian(2λ) rtol=1e-3
@test FiniteDiff.finite_difference_gradient(objective, λ) ≈ gradient(λ)' rtol=1e-3
@test FiniteDiff.finite_difference_hessian(objective, λ) ≈ hessian(λ) rtol=1e-3

end #submodule


text:
using Colors

# ??? code_warntype not defined?
# @code_warntype colordiff(colorant"red", colorant"darkred")

#=
function pointDifference(targetColor, corpusColor )::Float64
    return colordiff(
        corpusColor,
        targetColor,
        metric )
    # WAS metric=DE_94() i.e. computing the metric deep in the loops
end
=#

function maxPointDifference()
    #=
    According to Wolfram (but not the Julia site)
    the "typical" max for CIE2000 colordiff
    =#
    # return 1.16839  # CIE2000
    return 1.48916  # CIE84 DE_94
end

pointDifference = colordiff

#=
if subtract in loop, only 27 allocations.

=#
function testDiff()
    a = CartesianIndex(1,2)
    b = CartesianIndex(3,4)
    for i in 1:10000
        # This mutates the local a
        # a -= b

        # This does not mutate the local a
        # TODO why not?
        diffCI(a,b)
    end
    println(a)
end

function diffCI(a, b)
    # println itself requires allocations
    # println(typeof(a))
    a -= b
end

#=
Test whether a passed struct containing a CI
allocates ?
=#

struct Foo
    ci::CartesianIndex{N} where N
    b::Int64
end

mutable struct MutableFoo
    ci::CartesianIndex{N} where N
    b::Int64
end

function testDiff2()
    a = Foo(CartesianIndex(1,2), 1)
    b = Foo(CartesianIndex(3,4), 2)
    c = MutableFoo(CartesianIndex(1,2), 1)
    for i in 1:10000
        c.ci = diffCI2(a,b)
    end
    println(a)
end

function diffCI2(a, b)
    # Why isn't this returned without an allocation?
    # !!! Because the struct must be fully annotated with types?
    # return Foo(a.ci - b.ci, 3)
    return a.ci - b.ci
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using IRTools.All
using IRTools: Pipe
using Base: tail

ntail(x, n) = n <= 0 ? x : xcall(:tail, ntail(x, n-1))

function instrument!(pr, v, st)
  ex = st.expr
  if isexpr(ex, :new)
    st = stmt(st, expr = xcall(Zygote, :__new__, ex.args...))
    pr[v] = st
  elseif isexpr(ex, :splatnew)
    st = stmt(st, expr = xcall(Zygote, :__splatnew__, ex.args...))
    pr[v] = st
  end
  return st
end

function dual(ir)
  args = copy(arguments(ir))
  dx = argument!(ir, at = 1)
  Δs = Dict()
  for bl in blocks(ir)[2:end], arg in copy(arguments(bl))
    Δs[arg] = argument!(bl, insert = false)
  end
  pr = Pipe(ir)
  partial(x::Variable) = Δs[x]
  partial(x) = push!(pr, xcall(Forward, :zerolike, x))
  partial(v, x::Variable) = Δs[x]
  partial(v, x) = insert!(pr, v, xcall(Forward, :zerolike, x))
  for (i, x) in enumerate(args)
    if i == length(args) && ir.meta.method.isva
      Δs[x] = push!(pr, ntail(dx, i-1))
    else
      Δs[x] = push!(pr, xcall(:getindex, dx, i))
    end
  end
  branches(pr) do br
    args = arguments(br)
    if isreturn(br)
      args[1] = push!(pr, xcall(:tuple, args[1], partial(args[1])))
    else
      for arg in copy(args)
        push!(args, partial(arg))
      end
    end
    br
  end
  for (v, st) in pr
    st = instrument!(pr, v, st)
    if isexpr(st.expr, :meta, :inbounds, :loopinfo)
      Δs[v] = nothing
    elseif isexpr(st.expr, :boundscheck) ||
           (isexpr(st.expr, :call) && st.expr.args[1] == GlobalRef(Base, :not_int)) ||
           (isexpr(st.expr, :call) && st.expr.args[1] == GlobalRef(Core, :(===))) ||
           (isexpr(st.expr, :call) && st.expr.args[1] == GlobalRef(Main, :(===)))
      Δs[v] = false
    elseif isexpr(st.expr, :call)
      dargs = insert!(pr, v, xcall(:tuple, partial.((v,), st.expr.args)...))
      result = insert!(pr, v, stmt(st, expr = xcall(Forward, :_pushforward, dargs, st.expr.args...)))
      pr[v] = xcall(:getindex, result, 1)
      Δs[v] = push!(pr, xcall(:getindex, result, 2))
    elseif !isexpr(st.expr)
      Δs[v] = push!(pr, xcall(Forward, :zerolike, v))
    else
      error("Unsupported $(st.expr.head) expression")
    end
  end
  ir = finish(pr)
  return ir
end

@dynamo function _pushforward(_, x...)
  ir = IR(x...)
  ir === nothing && return :(error("non-differentiable function $(args[2])"))
  ir = Zygote.instrument(ir)
  ir.meta.code.inlineable = true
  return dual(ir)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export SerialComm, SerialDistributor

"""
    SerialDistributor()

Creates a distributor to work with SerialComm
"""
mutable struct SerialDistributor{GID <: Integer, PID <:Integer, LID <: Integer} <: Distributor{GID, PID, LID}
    post::Union{AbstractArray, Nothing}
    reversePost::Union{AbstractArray, Nothing}

    function SerialDistributor{GID, PID, LID}() where GID <: Integer where PID <: Integer where LID <: Integer
        new(nothing, nothing)
    end
end

function createFromSends(dist::SerialDistributor{GID, PID, LID},
        exportPIDs::AbstractArray{PID})::Integer where GID <: Integer where PID <: Integer where LID <: Integer
    for id in exportPIDs
        if id != 1
            throw(InvalidArgumentError("SerialDistributor can only accept PID of 1"))
        end
    end
    length(exportPIDs)
end

function createFromRecvs(
        dist::SerialDistributor{GID, PID, LID}, remoteGIDs::AbstractArray{GID}, remotePIDs::AbstractArray{PID}
        )::Tuple{AbstractArray{GID}, AbstractArray{PID}} where GID <: Integer where PID <: Integer where LID <: Integer
    for id in remotePIDs
        if id != 1
            throw(InvalidArgumentError("SerialDistributor can only accept PID of 1"))
        end
    end
    remoteGIDs,remotePIDs
end

function resolve(dist::SerialDistributor, exportObjs::AbstractArray{T})::AbstractArray{T} where T
    exportObjs
end

function resolveReverse(dist::SerialDistributor, exportObjs::AbstractArray{T})::AbstractArray{T} where T
    exportObjs
end

function resolvePosts(dist::SerialDistributor, exportObjs::AbstractArray)
    dist.post = exportObjs
end

function resolveWaits(dist::SerialDistributor)::AbstractArray
    if dist.post == nothing
        throw(InvalidStateError("Must post before waiting"))
    end

    result = dist.post
    dist.post = nothing
    result
end

function resolveReversePosts(dist::SerialDistributor, exportObjs::AbstractArray)
    dist.reversePost = exportObjs
end

function resolveReverseWaits(dist::SerialDistributor)::AbstractArray
     if dist.reversePost == nothing
        throw(InvalidStateError("Must reverse post before reverse waiting"))
    end

    result = dist.reversePost
    dist.reversePost = nothing
    result
end

"""
    SerialComm()

Gets an serial communication instance.
Serial communication results in mostly no-ops for the communication operations
"""
struct SerialComm{GID <: Integer, PID <:Integer, LID <: Integer} <: Comm{GID, PID, LID}
end

# most of these functions are no-ops or identify functions since there is only
# one processor

function barrier(comm::SerialComm)
end

function broadcastAll(comm::SerialComm, myVals::AbstractArray{T}, root::Integer)::Array{T} where T
    if root != 1
        throw(InvalidArgumentError("SerialComm can only accept PID of 1"))
    end
    myVals
end

function gatherAll(comm::SerialComm, myVals::AbstractArray{T})::Array{T} where T
    myVals
end

function sumAll(comm::SerialComm, partialsums::AbstractArray{T})::Array{T} where T
    partialsums
end

function maxAll(comm::SerialComm, partialmaxes::AbstractArray{T})::Array{T} where T
    partialmaxes
end

function minAll(comm::SerialComm, partialmins::AbstractArray{T})::Array{T} where T
    partialmins
end

function scanSum(comm::SerialComm, myvals::AbstractArray{T})::Array{T} where T
    myvals
end

function myPid(comm::SerialComm{GID, PID})::PID where GID <: Integer where PID <: Integer
    1
end

function numProc(comm::SerialComm{GID, PID})::PID where GID <: Integer where PID <: Integer
    1
end

function createDistributor(comm::SerialComm{GID, PID, LID})::SerialDistributor{GID, PID, LID} where GID <: Integer where PID <: Integer where LID <: Integer
    SerialDistributor{GID, PID, LID}()
end

getComm(::SerialDistributor{GID, PID, LID}) where {GID, PID, LID} = SerialComm{GID, PID, LID}()


text:
import RowEchelon

#=
Homework 5
Section 2.2
Problem 4
=#

# Initialize the Matrix
A = [2 -4 ; 4 -6]

# Calculate the Inverse
invA = inv(A)

println("The Inverse of A is:")
println(invA)


text:
name         = r"[^\W\d]\w*"                           ^Na
number       = r"(0|[1-9]\d*)(\.\d+)?([eE][+-]?\d+)?"  ^Nu
str          = r"\"([^\"\\]|\\[\"\\])*\""              ^St
operator     = r"[!#%&*\-./:?@\p{S}]+"

literal      = name | number | str

unit()       = literal | group | block

arguments()  = binaries & ( "," & arguments | ϵ )

call         = (( "(" & ")" ) |  "(" & arguments & ")" )
indexing     =                   "[" & arguments & "]"

operand      = unit &̲ ( call ^C | indexing ^I | ϵ )
operand     ^= combine

unary        = operator ^Un &̲ operand

expression() = ( unary | operand ) &̲
               ( nbsp &̲ operator ^Op &̲ sp &̲ expression
               |        operator ^Op &̲    expression | ϵ )

binaries     = expression ^E

separator    = ";" | eol

line()       = binaries &̲ nbspᵋ &̲ (             line  | ϵ )
lines()      = line ^Ln &̲ nbspᵋ &̲ ( separator & lines | ϵ )

group        =                 "(" & line  ^G & ")"
block        = ( "{" & "}" ) | "{" & lines ^B & "}"

program      = bof & lines ^B & eof
program     ^= cleanup


text:
import VoronoiCells: LEFT, RIGHT, LOWER, UPPER
using VoronoiDelaunay

@testset "Clipping lines" begin
	# ------------------------------------------------------------
	# Points inside the bounding box are untouched by clip

	A = Point2D(LEFT+rand(), LOWER+rand())
	B = Point2D(LEFT+rand(), LOWER+rand())
	C, D = clip(A, B)

	@test A == C
	@test B == D

	# ------------------------------------------------------------
	# With one point on the boundary and the other outside (in the correct
	# direction), clip should return two identical points

	A = Point2D(LEFT, LOWER+rand())
	B = Point2D(LEFT-0.5, LOWER+rand())
	C, D = clip(A, B)
	@test C == D

	A = Point2D(LEFT+rand(), UPPER)
	B = Point2D(LEFT+rand(), UPPER+0.5)
	C, D = clip(A, B)
	@test C == D

	# ------------------------------------------------------------
	# When clipping off part of a line segment the slope should not change

	A = Point2D(LEFT+rand(), LOWER+rand())
	B = Point2D(LEFT+rand(), UPPER+rand())
	C, D = clip(A, B)

	slope1 = (gety(A)-gety(B)) / (getx(A)-getx(B))
	slope2 = (gety(C)-gety(D)) / (getx(C)-getx(D))

	@test slope1 ≈ slope2
end


text:
@testset "Newton" begin
    function himmelblau!(∇f, x)
        if !(∇f == nothing)
            ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
                44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
            ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
        end

        fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        return ∇f == nothing ? fx : (fx, ∇f)
    end

    function himmelblau!(∇²f, ∇f, x)
        if !(∇²f == nothing)
            ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 42.0
            ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
            ∇²f[2, 1] = 4.0 * x[1] + 4.0 * x[2]
            ∇²f[2, 2] = 12.0 * x[2]^2 + 4.0 * x[1] - 26.0
        end
        if !(∇f == nothing)
            ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
                44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
            ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
        end

        fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        if ∇f == nothing && ∇²f == nothing
            return fx
        elseif ∇²f == nothing
            return fx, ∇f
        else
            return fx, ∇f, ∇²f
        end
    end

    res = minimize!(himmelblau!, copy([2.0,2.0]), Newton(DirectApprox()))
    @test norm(res[3], Inf) < 1e-8
    @printf("NN! Newton (direct): %2.2e  %d\n", norm(res[3], Inf), res[4])
end


text:
"""
    createtemplate(path=".", template=:book)
Creates the latex files such as `main.tex`,
`preface.tex`, `frontmatter/titlepage.tex` and
`frontmatter/copyright.tex`. The `.tex` files depend on the template used.
At the moment, the available templates are `:book` and
`:mathbook`, which are pretty much the same, but `:mathbook`
imports some more packages specific for mathematics.
"""
function createtemplate(path=".", template=:book)

    if template == :book
        tex = read(String(@__DIR__)*"/../templates/book.tex", String)

    elseif template == :mathbook
        tex = read(String(@__DIR__)*"/../templates/mathbook.tex", String)

    elseif template == :article
        tex = read(String(@__DIR__)*"/../templates/article.tex", String)

    elseif template == :matharticle
        tex = read(String(@__DIR__)*"/../templates/matharticle.tex", String)
    end

    maintex = path * "/main.tex"
    if !isfile(maintex)
        open(maintex, "w") do f
            write(f, tex)
        end
    end

    if template == :book || template == :mathbook
        preface = "\\newpage\n\\chapter*{Preface}\n\\addcontentsline{toc}{chapter}{Preface}\n"
        prefacetex = path * "/preface.tex"

        if !isfile(prefacetex)
            open(prefacetex, "w") do f
                write(f, preface)
            end
        end

        titlepage = read(String(@__DIR__)*"/../templates/frontmatter/titlepage.tex", String)

        titlepagetex = path * "/frontmatter/titlepage.tex"
        if !isfile(titlepagetex)
            open(titlepagetex, "w") do f
                write(f, titlepage)
            end
        end

        copyright = read(String(@__DIR__)*"/../templates/frontmatter/copyright.tex", String)

        copyrighttex = path * "/frontmatter/copyright.tex"
        if !isfile(copyrighttex)
            open(copyrighttex, "w") do f
                write(f, copyright)
            end
        end
    end
end


text:
function update_ticks_guides(d::KW, labs, i, j, n)
    # d[:title]  = (i==1 ? _cycle(labs,j) : "")
    # d[:xticks] = (i==n)
    d[:xguide] = (i==n ? _cycle(labs,j) : "")
    # d[:yticks] = (j==1)
    d[:yguide] = (j==1 ? _cycle(labs,i) : "")
end

function plot!(scene::Scene, ::Type{CorrPlot}, attributes::Attributes, mat)
    mat = cp.args[1]
    n = size(mat, 2)
    C = cor(mat)
    attributes, plot_att = merged_get!(:boxplot, scene, attributes) do
        Theme(
            link = :x,  # need custom linking for y
            layout = (n,n),
            legend = false,
            margin = 1mm,
            titlefont = font(11),
            fillcolor = attributes[:color],
            linecolor = attributes[:],
            markeralpha = 0.4,
            grad = cgrad(get(attributes, :markercolor, cgrad())),
            indices = reshape(1:n^2, n, n)',
            title = "",
        )
    end
    # this is simples way to overwrite plot defaults from the global theme,
    # but still allow user owerwrites
    subscene = Scene(scene, theme = Theme(
        get!(plot_att, k, v)
    ))

    labs = pop!(attributes, :label, [""])
    plotgrid = broadcast(1:n, (1:n)') do i, j
        vi = view(mat,:,i)
        vj = view(mat,:,j)
        if i == j # histograms are on the diagonal
            histogram!(subscene, vi, grid = false)
        elseif i > j
            scatter!(subscene, vj, vi)
        else
            plot!(subscene, vj, vi, seriestype = attributes[:seriestype])
        end
    end
    plot!(scene, rest, plotgrid)
end

function plot!(plot::CorrPlot)
    mat = plot[1]
    n = size(mat, 2)
    C = cor(mat)
    plotgrid = broadcast(1:n, (1:n)') do i, j
        vi = view(mat,:,i)
        vj = view(mat,:,j)
        if i == j # histograms are on the diagonal
            histogram!(plot, vi, grid = false)
        elseif i > j
            scatter!(plot, vj, vi, markercolor = grad[0.5 + 0.5C[i,j]])
        else
            plot!(plot, vj, vi, seriestype = attributes[:seriestype])
        end
    end
    plot!(plotgrid) # matrix of plots -> gridlayout
end