text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Gives a reinterpreted view (of element type T) of the underlying array (of element type S).
If the size of `T` differs from the size of `S`, the array will be compressed/expanded in
the first dimension. The variant `reinterpret(reshape, T, a)` instead adds or consumes the first dimension
depending on the ratio of element sizes.
"""
struct ReinterpretArray{T,N,S,A<:AbstractArray{S},IsReshaped} <: AbstractArray{T, N}
    parent::A
    readable::Bool
    writable::Bool

    function throwbits(S::Type, T::Type, U::Type)
        @_noinline_meta
        throw(ArgumentError("cannot reinterpret `$(S)` as `$(T)`, type `$(U)` is not a bits type"))
    end
    function throwsize0(S::Type, T::Type, msg)
        @_noinline_meta
        throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a $msg size"))
    end

    global reinterpret
    function reinterpret(::Type{T}, a::A) where {T,N,S,A<:AbstractArray{S, N}}
        function thrownonint(S::Type, T::Type, dim)
            @_noinline_meta
            throw(ArgumentError("""
                cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
                The resulting array would have non-integral first dimension.
                """))
        end
        function throwaxes1(S::Type, T::Type, ax1)
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first."))
        end
        isbitstype(T) || throwbits(S, T, T)
        isbitstype(S) || throwbits(S, T, S)
        (N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T, "different")
        if N != 0 && sizeof(S) != sizeof(T)
            ax1 = axes(a)[1]
            dim = length(ax1)
            rem(dim*sizeof(S),sizeof(T)) == 0 || thrownonint(S, T, dim)
            first(ax1) == 1 || throwaxes1(S, T, ax1)
        end
        readable = array_subpadding(T, S)
        writable = array_subpadding(S, T)
        new{T, N, S, A, false}(a, readable, writable)
    end
    reinterpret(::Type{T}, a::AbstractArray{T}) where {T} = a

    # With reshaping
    function reinterpret(::typeof(reshape), ::Type{T}, a::A) where {T,S,A<:AbstractArray{S}}
        function throwintmult(S::Type, T::Type)
            @_noinline_meta
            throw(ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got $(sizeof(T))) and `sizeof(eltype(a))` (got $(sizeof(S))) be an integer multiple of the other"))
        end
        function throwsize1(a::AbstractArray, T::Type)
            @_noinline_meta
            throw(ArgumentError("`reinterpret(reshape, $T, a)` where `eltype(a)` is $(eltype(a)) requires that `axes(a, 1)` (got $(axes(a, 1))) be equal to 1:$(sizeof(T) ÷ sizeof(eltype(a))) (from the ratio of element sizes)"))
        end
        isbitstype(T) || throwbits(S, T, T)
        isbitstype(S) || throwbits(S, T, S)
        if sizeof(S) == sizeof(T)
            N = ndims(a)
        elseif sizeof(S) > sizeof(T)
            rem(sizeof(S), sizeof(T)) == 0 || throwintmult(S, T)
            N = ndims(a) + 1
        else
            rem(sizeof(T), sizeof(S)) == 0 || throwintmult(S, T)
            N = ndims(a) - 1
            N > -1 || throwsize0(S, T, "larger")
            axes(a, 1) == Base.OneTo(sizeof(T) ÷ sizeof(S)) || throwsize1(a, T)
        end
        readable = array_subpadding(T, S)
        writable = array_subpadding(S, T)
        new{T, N, S, A, true}(a, readable, writable)
    end
    reinterpret(::typeof(reshape), ::Type{T}, a::AbstractArray{T}) where {T} = a
end

ReshapedReinterpretArray{T,N,S,A<:AbstractArray{S}} = ReinterpretArray{T,N,S,A,true}
NonReshapedReinterpretArray{T,N,S,A<:AbstractArray{S, N}} = ReinterpretArray{T,N,S,A,false}

"""
    reinterpret(reshape, T, A::AbstractArray{S}) -> B

Change the type-interpretation of `A` while consuming or adding a "channel dimension."

If `sizeof(T) = n*sizeof(S)` for `n>1`, `A`'s first dimension must be
of size `n` and `B` lacks `A`'s first dimension. Conversely, if `sizeof(S) = n*sizeof(T)` for `n>1`,
`B` gets a new first dimension of size `n`. The dimensionality is unchanged if `sizeof(T) == sizeof(S)`.

# Examples

```jldoctest
julia> A = [1 2; 3 4]
2×2 Matrix{$Int}:
 1  2
 3  4

julia> reinterpret(reshape, Complex{Int}, A)    # the result is a vector
2-element reinterpret(reshape, Complex{$Int}, ::Matrix{$Int}) with eltype Complex{$Int}:
 1 + 3im
 2 + 4im

julia> a = [(1,2,3), (4,5,6)]
2-element Vector{Tuple{$Int, $Int, $Int}}:
 (1, 2, 3)
 (4, 5, 6)

julia> reinterpret(reshape, Int, a)             # the result is a matrix
3×2 reinterpret(reshape, $Int, ::Vector{Tuple{$Int, $Int, $Int}}) with eltype $Int:
 1  4
 2  5
 3  6
```
"""
reinterpret(::typeof(reshape), T::Type, a::AbstractArray)

reinterpret(::Type{T}, a::NonReshapedReinterpretArray) where {T} = reinterpret(T, a.parent)
reinterpret(::typeof(reshape), ::Type{T}, a::ReshapedReinterpretArray) where {T} = reinterpret(reshape, T, a.parent)

# Definition of StridedArray
StridedFastContiguousSubArray{T,N,A<:DenseArray} = FastContiguousSubArray{T,N,A}
StridedReinterpretArray{T,N,A<:Union{DenseArray,StridedFastContiguousSubArray},IsReshaped} = ReinterpretArray{T,N,S,A,IsReshaped} where S
StridedReshapedArray{T,N,A<:Union{DenseArray,StridedFastContiguousSubArray,StridedReinterpretArray}} = ReshapedArray{T,N,A}
StridedSubArray{T,N,A<:Union{DenseArray,StridedReshapedArray,StridedReinterpretArray},
    I<:Tuple{Vararg{Union{RangeIndex, ReshapedUnitRange, AbstractCartesianIndex}}}} = SubArray{T,N,A,I}
StridedArray{T,N} = Union{DenseArray{T,N}, StridedSubArray{T,N}, StridedReshapedArray{T,N}, StridedReinterpretArray{T,N}}
StridedVector{T} = StridedArray{T,1}
StridedMatrix{T} = StridedArray{T,2}
StridedVecOrMat{T} = Union{StridedVector{T}, StridedMatrix{T}}

# the definition of strides for Array{T,N} is tuple() if N = 0, otherwise it is
# a tuple containing 1 and a cumulative product of the first N-1 sizes
# this definition is also used for StridedReshapedArray and StridedReinterpretedArray
# which have the same memory storage as Array
stride(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}, i::Int) = _stride(a, i)

function stride(a::ReinterpretArray, i::Int)
    a.parent isa StridedArray || ArgumentError("Parent must be strided.") |> throw
    return _stride(a, i)
end

function _stride(a, i)
    if i > ndims(a)
        return length(a)
    end
    s = 1
    for n = 1:(i-1)
        s *= size(a, n)
    end
    return s
end

function strides(a::ReinterpretArray)
    a.parent isa StridedArray || ArgumentError("Parent must be strided.") |> throw
    size_to_strides(1, size(a)...)
end
strides(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}) = size_to_strides(1, size(a)...)

function check_readable(a::ReinterpretArray{T, N, S} where N) where {T,S}
    # See comment in check_writable
    if !a.readable && !array_subpadding(T, S)
        throw(PaddingError(T, S))
    end
end

function check_writable(a::ReinterpretArray{T, N, S} where N) where {T,S}
    # `array_subpadding` is relatively expensive (compared to a simple arrayref),
    # so it is cached in the array. However, it is computable at compile time if,
    # inference has the types available. By using this form of the check, we can
    # get the best of both worlds for the success case. If the types were not
    # available to inference, we simply need to check the field (relatively cheap)
    # and if they were we should be able to fold this check away entirely.
    if !a.writable && !array_subpadding(S, T)
        throw(PaddingError(T, S))
    end
end

## IndexStyle specializations

# For `reinterpret(reshape, T, a)` where we're adding a channel dimension and with
# `IndexStyle(a) == IndexLinear()`, it's advantageous to retain pseudo-linear indexing.
struct IndexSCartesian2{K} <: IndexStyle end   # K = sizeof(S) ÷ sizeof(T), a static-sized 2d cartesian iterator

IndexStyle(::Type{ReinterpretArray{T,N,S,A,false}}) where {T,N,S,A<:AbstractArray{S,N}} = IndexStyle(A)
function IndexStyle(::Type{ReinterpretArray{T,N,S,A,true}}) where {T,N,S,A<:AbstractArray{S}}
    if sizeof(T) < sizeof(S)
        IndexStyle(A) === IndexLinear() && return IndexSCartesian2{sizeof(S) ÷ sizeof(T)}()
        return IndexCartesian()
    end
    return IndexStyle(A)
end
IndexStyle(::IndexSCartesian2{K}, ::IndexSCartesian2{K}) where {K} = IndexSCartesian2{K}()

struct SCartesianIndex2{K}   # can't make <:AbstractCartesianIndex without N, and 2 would be a bit misleading
    i::Int
    j::Int
end
to_index(i::SCartesianIndex2) = i

struct SCartesianIndices2{K,R<:AbstractUnitRange{Int}} <: AbstractMatrix{SCartesianIndex2{K}}
    indices2::R
end
SCartesianIndices2{K}(indices2::AbstractUnitRange{Int}) where {K} = (@assert K::Int > 1; SCartesianIndices2{K,typeof(indices2)}(indices2))

eachindex(::IndexSCartesian2{K}, A::AbstractArray) where {K} = SCartesianIndices2{K}(eachindex(IndexLinear(), parent(A)))

size(iter::SCartesianIndices2{K}) where K = (K, length(iter.indices2))
axes(iter::SCartesianIndices2{K}) where K = (Base.OneTo(K), iter.indices2)

first(iter::SCartesianIndices2{K}) where {K} = SCartesianIndex2{K}(1, first(iter.indices2))
last(iter::SCartesianIndices2{K}) where {K}  = SCartesianIndex2{K}(K, last(iter.indices2))

@inline function getindex(iter::SCartesianIndices2{K}, i::Int, j::Int) where {K}
    @boundscheck checkbounds(iter, i, j)
    return SCartesianIndex2{K}(i, iter.indices2[j])
end

function iterate(iter::SCartesianIndices2{K}) where {K}
    ret = iterate(iter.indices2)
    ret === nothing && return nothing
    item2, state2 = ret
    return SCartesianIndex2{K}(1, item2), (1, item2, state2)
end

function iterate(iter::SCartesianIndices2{K}, (state1, item2, state2)) where {K}
    if state1 < K
        item1 = state1 + 1
        return SCartesianIndex2{K}(item1, item2), (item1, item2, state2)
    end
    ret = iterate(iter.indices2, state2)
    ret === nothing && return nothing
    item2, state2 = ret
    return SCartesianIndex2{K}(1, item2), (1, item2, state2)
end

SimdLoop.simd_outer_range(iter::SCartesianIndices2) = iter.indices2
SimdLoop.simd_inner_length(::SCartesianIndices2{K}, ::Any) where K = K
@inline function SimdLoop.simd_index(::SCartesianIndices2{K}, Ilast::Int, I1::Int) where {K}
    SCartesianIndex2{K}(I1+1, Ilast)
end

_maybe_reshape(::IndexSCartesian2, A::ReshapedReinterpretArray, I...) = A

# fallbacks
function _getindex(::IndexSCartesian2, A::AbstractArray{T,N}, I::Vararg{Int, N}) where {T,N}
    @_propagate_inbounds_meta
    getindex(A, I...)
end
function _setindex!(::IndexSCartesian2, A::AbstractArray{T,N}, v, I::Vararg{Int, N}) where {T,N}
    @_propagate_inbounds_meta
    setindex!(A, v, I...)
end
# fallbacks for array types that use "pass-through" indexing (e.g., `IndexStyle(A) = IndexStyle(parent(A))`)
# but which don't handle SCartesianIndex2
function _getindex(::IndexSCartesian2, A::AbstractArray{T,N}, ind::SCartesianIndex2) where {T,N}
    @_propagate_inbounds_meta
    I = _to_subscript_indices(A, ind.i, ind.j)
    getindex(A, I...)
end
function _setindex!(::IndexSCartesian2, A::AbstractArray{T,N}, v, ind::SCartesianIndex2) where {T,N}
    @_propagate_inbounds_meta
    I = _to_subscript_indices(A, ind.i, ind.j)
    setindex!(A, v, I...)
end

## AbstractArray interface

parent(a::ReinterpretArray) = a.parent
dataids(a::ReinterpretArray) = dataids(a.parent)
unaliascopy(a::ReinterpretArray{T}) where {T} = reinterpret(T, unaliascopy(a.parent))

function size(a::NonReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
    psize = size(a.parent)
    size1 = div(psize[1]*sizeof(S), sizeof(T))
    tuple(size1, tail(psize)...)
end
function size(a::ReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
    psize = size(a.parent)
    sizeof(S) > sizeof(T) && return (div(sizeof(S), sizeof(T)), psize...)
    sizeof(S) < sizeof(T) && return Base.tail(psize)
    return psize
end
size(a::NonReshapedReinterpretArray{T,0}) where {T} = ()

function axes(a::NonReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
    paxs = axes(a.parent)
    f, l = first(paxs[1]), length(paxs[1])
    size1 = div(l*sizeof(S), sizeof(T))
    tuple(oftype(paxs[1], f:f+size1-1), tail(paxs)...)
end
function axes(a::ReshapedReinterpretArray{T,N,S} where {N}) where {T,S}
    paxs = axes(a.parent)
    sizeof(S) > sizeof(T) && return (Base.OneTo(div(sizeof(S), sizeof(T))), paxs...)
    sizeof(S) < sizeof(T) && return Base.tail(paxs)
    return paxs
end
axes(a::NonReshapedReinterpretArray{T,0}) where {T} = ()

elsize(::Type{<:ReinterpretArray{T}}) where {T} = sizeof(T)
unsafe_convert(::Type{Ptr{T}}, a::ReinterpretArray{T,N,S} where N) where {T,S} = Ptr{T}(unsafe_convert(Ptr{S},a.parent))

@inline @propagate_inbounds getindex(a::NonReshapedReinterpretArray{T,0}) where {T} = reinterpret(T, a.parent[])
@inline @propagate_inbounds getindex(a::ReinterpretArray) = a[1]

@inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, inds::Vararg{Int, N}) where {T,N,S}
    check_readable(a)
    _getindex_ra(a, inds[1], tail(inds))
end

@inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, i::Int) where {T,N,S}
    check_readable(a)
    if isa(IndexStyle(a), IndexLinear)
        return _getindex_ra(a, i, ())
    end
    # Convert to full indices here, to avoid needing multiple conversions in
    # the loop in _getindex_ra
    inds = _to_subscript_indices(a, i)
    isempty(inds) ? _getindex_ra(a, 1, ()) : _getindex_ra(a, inds[1], tail(inds))
end

@inline @propagate_inbounds function getindex(a::ReshapedReinterpretArray{T,N,S}, ind::SCartesianIndex2) where {T,N,S}
    check_readable(a)
    n = sizeof(S) ÷ sizeof(T)
    t = Ref{NTuple{n,T}}()
    s = Ref{S}(a.parent[ind.j])
    GC.@preserve t s begin
        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
        _memcpy!(tptr, sptr, sizeof(S))
    end
    return t[][ind.i]
end

@inline _memcpy!(dst, src, n) = ccall(:memcpy, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), dst, src, n)

@inline @propagate_inbounds function _getindex_ra(a::NonReshapedReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return reinterpret(T, a.parent[i1, tailinds...])
    else
        @boundscheck checkbounds(a, i1, tailinds...)
        ind_start, sidx = divrem((i1-1)*sizeof(T), sizeof(S))
        t = Ref{T}()
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            # Optimizations that avoid branches
            if sizeof(T) % sizeof(S) == 0
                # T is bigger than S and contains an integer number of them
                n = sizeof(T) ÷ sizeof(S)
                for i = 1:n
                    s[] = a.parent[ind_start + i, tailinds...]
                    _memcpy!(tptr + (i-1)*sizeof(S), sptr, sizeof(S))
                end
            elseif sizeof(S) % sizeof(T) == 0
                # S is bigger than T and contains an integer number of them
                s[] = a.parent[ind_start + 1, tailinds...]
                _memcpy!(tptr, sptr + sidx, sizeof(T))
            else
                i = 1
                nbytes_copied = 0
                # This is a bit complicated to deal with partial elements
                # at both the start and the end. LLVM will fold as appropriate,
                # once it knows the data layout
                while nbytes_copied < sizeof(T)
                    s[] = a.parent[ind_start + i, tailinds...]
                    nb = min(sizeof(S) - sidx, sizeof(T)-nbytes_copied)
                    _memcpy!(tptr + nbytes_copied, sptr + sidx, nb)
                    nbytes_copied += nb
                    sidx = 0
                    i += 1
                end
            end
        end
        return t[]
    end
end

@inline @propagate_inbounds function _getindex_ra(a::ReshapedReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return reinterpret(T, a.parent[i1, tailinds...])
    end
    @boundscheck checkbounds(a, i1, tailinds...)
    if sizeof(T) >= sizeof(S)
        t = Ref{T}()
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            if sizeof(T) > sizeof(S)
                # Extra dimension in the parent array
                n = sizeof(T) ÷ sizeof(S)
                if isempty(tailinds) && IndexStyle(a.parent) === IndexLinear()
                    offset = n * (i1 - firstindex(a))
                    for i = 1:n
                        s[] = a.parent[i + offset]
                        _memcpy!(tptr + (i-1)*sizeof(S), sptr, sizeof(S))
                    end
                else
                    for i = 1:n
                        s[] = a.parent[i, i1, tailinds...]
                        _memcpy!(tptr + (i-1)*sizeof(S), sptr, sizeof(S))
                    end
                end
            else
                # No extra dimension
                s[] = a.parent[i1, tailinds...]
                _memcpy!(tptr, sptr, sizeof(S))
            end
        end
        return t[]
    end
    # S is bigger than T and contains an integer number of them
    n = sizeof(S) ÷ sizeof(T)
    t = Ref{NTuple{n,T}}()
    s = Ref{S}()
    GC.@preserve t s begin
        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
        s[] = a.parent[tailinds...]
        _memcpy!(tptr, sptr, sizeof(S))
    end
    return t[][i1]
end

@inline @propagate_inbounds setindex!(a::NonReshapedReinterpretArray{T,0,S} where T, v) where {S} = (a.parent[] = reinterpret(S, v))
@inline @propagate_inbounds setindex!(a::ReinterpretArray, v) = (a[1] = v)

@inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, inds::Vararg{Int, N}) where {T,N,S}
    check_writable(a)
    _setindex_ra!(a, v, inds[1], tail(inds))
end

@inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, i::Int) where {T,N,S}
    check_writable(a)
    if isa(IndexStyle(a), IndexLinear)
        return _setindex_ra!(a, v, i, ())
    end
    inds = _to_subscript_indices(a, i)
    _setindex_ra!(a, v, inds[1], tail(inds))
end

@inline @propagate_inbounds function setindex!(a::ReshapedReinterpretArray{T,N,S}, v, ind::SCartesianIndex2) where {T,N,S}
    check_writable(a)
    v = convert(T, v)::T
    t = Ref{T}(v)
    s = Ref{S}(a.parent[ind.j])
    GC.@preserve t s begin
        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
        _memcpy!(sptr + (ind.i-1)*sizeof(T), tptr, sizeof(T))
    end
    a.parent[ind.j] = s[]
    return a
end

@inline @propagate_inbounds function _setindex_ra!(a::NonReshapedReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
    v = convert(T, v)::T
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return setindex!(a.parent, reinterpret(S, v), i1, tailinds...)
    else
        @boundscheck checkbounds(a, i1, tailinds...)
        ind_start, sidx = divrem((i1-1)*sizeof(T), sizeof(S))
        t = Ref{T}(v)
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            # Optimizations that avoid branches
            if sizeof(T) % sizeof(S) == 0
                # T is bigger than S and contains an integer number of them
                n = sizeof(T) ÷ sizeof(S)
                for i = 0:n-1
                    _memcpy!(sptr, tptr + i*sizeof(S), sizeof(S))
                    a.parent[ind_start + i + 1, tailinds...] = s[]
                end
            elseif sizeof(S) % sizeof(T) == 0
                # S is bigger than T and contains an integer number of them
                s[] = a.parent[ind_start + 1, tailinds...]
                _memcpy!(sptr + sidx, tptr, sizeof(T))
                a.parent[ind_start + 1, tailinds...] = s[]
            else
                nbytes_copied = 0
                i = 1
                # Deal with any partial elements at the start. We'll have to copy in the
                # element from the original array and overwrite the relevant parts
                if sidx != 0
                    s[] = a.parent[ind_start + i, tailinds...]
                    nb = min((sizeof(S) - sidx) % UInt, sizeof(T) % UInt)
                    _memcpy!(sptr + sidx, tptr, nb)
                    nbytes_copied += nb
                    a.parent[ind_start + i, tailinds...] = s[]
                    i += 1
                    sidx = 0
                end
                # Deal with the main body of elements
                while nbytes_copied < sizeof(T) && (sizeof(T) - nbytes_copied) > sizeof(S)
                    nb = min(sizeof(S), sizeof(T) - nbytes_copied)
                    _memcpy!(sptr, tptr + nbytes_copied, nb)
                    nbytes_copied += nb
                    a.parent[ind_start + i, tailinds...] = s[]
                    i += 1
                end
                # Deal with trailing partial elements
                if nbytes_copied < sizeof(T)
                    s[] = a.parent[ind_start + i, tailinds...]
                    nb = min(sizeof(S), sizeof(T) - nbytes_copied)
                    _memcpy!(sptr, tptr + nbytes_copied, nb)
                    a.parent[ind_start + i, tailinds...] = s[]
                end
            end
        end
    end
    return a
end

@inline @propagate_inbounds function _setindex_ra!(a::ReshapedReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
    v = convert(T, v)::T
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return setindex!(a.parent, reinterpret(S, v), i1, tailinds...)
    end
    @boundscheck checkbounds(a, i1, tailinds...)
    t = Ref{T}(v)
    s = Ref{S}()
    GC.@preserve t s begin
        tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
        sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
        if sizeof(T) >= sizeof(S)
            if sizeof(T) > sizeof(S)
                # Extra dimension in the parent array
                n = sizeof(T) ÷ sizeof(S)
                if isempty(tailinds) && IndexStyle(a.parent) === IndexLinear()
                    offset = n * (i1 - firstindex(a))
                    for i = 1:n
                        _memcpy!(sptr, tptr + (i-1)*sizeof(S), sizeof(S))
                        a.parent[i + offset] = s[]
                    end
                else
                    for i = 1:n
                        _memcpy!(sptr, tptr + (i-1)*sizeof(S), sizeof(S))
                        a.parent[i, i1, tailinds...] = s[]
                    end
                end
            else
                # No extra dimension
                _memcpy!(sptr, tptr, sizeof(S))
                a.parent[i1, tailinds...] = s[]
            end
        else
            # S is bigger than T and contains an integer number of them
            s[] = a.parent[tailinds...]
            _memcpy!(sptr + (i1-1)*sizeof(T), tptr, sizeof(T))
            a.parent[tailinds...] = s[]
        end
    end
    return a
end

# Padding
struct Padding
    offset::Int
    size::Int
end
function intersect(p1::Padding, p2::Padding)
    start = max(p1.offset, p2.offset)
    stop = min(p1.offset + p1.size, p2.offset + p2.size)
    Padding(start, max(0, stop-start))
end

struct PaddingError
    S::Type
    T::Type
end

function showerror(io::IO, p::PaddingError)
    print(io, "Padding of type $(p.S) is not compatible with type $(p.T).")
end

"""
    CyclePadding(padding, total_size)

Cylces an iterator of `Padding` structs, restarting the padding at `total_size`.
E.g. if `padding` is all the padding in a struct and `total_size` is the total
aligned size of that array, `CyclePadding` will correspond to the padding in an
infinite vector of such structs.
"""
struct CyclePadding{P}
    padding::P
    total_size::Int
end
eltype(::Type{<:CyclePadding}) = Padding
IteratorSize(::Type{<:CyclePadding}) = IsInfinite()
isempty(cp::CyclePadding) = isempty(cp.padding)
function iterate(cp::CyclePadding)
    y = iterate(cp.padding)
    y === nothing && return nothing
    y[1], (0, y[2])
end
function iterate(cp::CyclePadding, state::Tuple)
    y = iterate(cp.padding, tail(state)...)
    y === nothing && return iterate(cp, (state[1]+cp.total_size,))
    Padding(y[1].offset+state[1], y[1].size), (state[1], tail(y)...)
end

"""
    Compute the location of padding in a type.
"""
function padding(T)
    padding = Padding[]
    last_end::Int = 0
    for i = 1:fieldcount(T)
        offset = fieldoffset(T, i)
        fT = fieldtype(T, i)
        if offset != last_end
            push!(padding, Padding(offset, offset-last_end))
        end
        last_end = offset + sizeof(fT)
    end
    padding
end

function CyclePadding(T::DataType)
    a, s = datatype_alignment(T), sizeof(T)
    as = s + (a - (s % a)) % a
    pad = padding(T)
    s != as && push!(pad, Padding(s, as - s))
    CyclePadding(pad, as)
end

using .Iterators: Stateful
@pure function array_subpadding(S, T)
    checked_size = 0
    lcm_size = lcm(sizeof(S), sizeof(T))
    s, t = Stateful{<:Any, Any}(CyclePadding(S)),
           Stateful{<:Any, Any}(CyclePadding(T))
    isempty(t) && return true
    isempty(s) && return false
    while checked_size < lcm_size
        # Take padding in T
        pad = popfirst!(t)
        # See if there's corresponding padding in S
        while true
            ps = peek(s)
            ps.offset > pad.offset && return false
            intersect(ps, pad) == pad && break
            popfirst!(s)
        end
        checked_size = pad.offset + pad.size
    end
    return true
end

# Reductions with IndexSCartesian2

function _mapreduce(f::F, op::OP, style::IndexSCartesian2{K}, A::AbstractArrayOrBroadcasted) where {F,OP,K}
    inds = eachindex(style, A)
    n = size(inds)[2]
    if n == 0
        return mapreduce_empty_iter(f, op, A, IteratorEltype(A))
    else
        return mapreduce_impl(f, op, A, first(inds), last(inds))
    end
end

@noinline function mapreduce_impl(f::F, op::OP, A::AbstractArrayOrBroadcasted,
                                  ifirst::SCI, ilast::SCI, blksize::Int) where {F,OP,SCI<:SCartesianIndex2{K}} where K
    if ifirst.j + blksize > ilast.j
        # sequential portion
        @inbounds a1 = A[ifirst]
        @inbounds a2 = A[SCI(2,ifirst.j)]
        v = op(f(a1), f(a2))
        @simd for i = ifirst.i + 2 : K
            @inbounds ai = A[SCI(i,ifirst.j)]
            v = op(v, f(ai))
        end
        # Remaining columns
        for j = ifirst.j+1 : ilast.j
            @simd for i = 1:K
                @inbounds ai = A[SCI(i,j)]
                v = op(v, f(ai))
            end
        end
        return v
    else
        # pairwise portion
        jmid = (ifirst.j + ilast.j) >> 1
        v1 = mapreduce_impl(f, op, A, ifirst, SCI(K,jmid), blksize)
        v2 = mapreduce_impl(f, op, A, SCI(1,jmid+1), ilast, blksize)
        return op(v1, v2)
    end
end

mapreduce_impl(f::F, op::OP, A::AbstractArrayOrBroadcasted, ifirst::SCartesianIndex2, ilast::SCartesianIndex2) where {F,OP} =
    mapreduce_impl(f, op, A, ifirst, ilast, pairwise_blocksize(f, op))

max_stars_repo_name: jishnub/julia

max_issues_repo_name: jishnub/julia

max_forks_repo_name: jishnub/julia


text:
abstract type AbstractTenArray{T<:AbstractTen,L,N} <: AbstractArray{T,N} end

@inline Base.size(u::AbstractTenArray) =
    size(u.xx)

@inline Base.length(::Type{<:AbstractTenArray{T,L,N}}) where {T,L,N} =
    L

@inline Base.length(u::AbstractTenArray) =
    length(typeof(u))

@inline Base.linearindices(::Type{<:AbstractTenArray{T,L,N}}) where {T,L,N} =
    Base.OneTo(L)

@inline Base.linearindices(u::AbstractTenArray) =
    linearindices(typeof(u))

Base.IndexStyle(u::AbstractTenArray) =
    IndexLinear()

struct Ten2DArray{T,L,N} <: AbstractTenArray{Ten2D{T},L,N}
    xx::Array{T,N}
    yx::Array{T,N}
    xy::Array{T,N}
    yy::Array{T,N}

    function Ten2DArray{T,L,N}(xx,yx,xy,yy) where {T,L,N}
        ((length(xx) == L && length(xy) == L && length(yx) == L && length(yy) == L) &&
            (size(xx) == size(yx) == size(xy) == size(yy))) || throw(ArgumentError("Arrays must have the same dimensions"))
        return new{T,L,N}(xx,yx,xy,yy)
    end

end

Ten2DArray(xx,yx,xy,yy) =
    Ten2DArray{eltype(xx),length(xx),ndims(xx)}(xx,yx,xy,yy)

@inline function Base.getindex(u::Ten2DArray{T,L,N},I) where {T,L,N}
    xx = u.xx
    yx = u.yx
    xy = u.xy
    yy = u.yy
    @boundscheck checkbounds(xx,I)
    @inbounds begin
        xxv = xx[I]
        yxv = yx[I]
        xyv = xy[I]
        yyv = yy[I]
    end
    return Ten2D{T}(xxv,yxv,xyv,yyv)
end

@inline function Base.setindex!(u::Ten2DArray{T,L,N},v::Ten2D,I) where {T,L,N}
    xx = u.xx
    yx = u.yx
    xy = u.xy
    yy = u.yy
    @boundscheck checkbounds(xx,I)
    @inbounds begin
        xx[I] = v.xx
        yx[I] = v.yx
        xy[I] = v.xy
        yy[I] = v.yy
    end
    return u
end

ext: jl


text:
using Pathfinder
using Random
using Test
using Transducers

@testset "transducers integration" begin
    @testset "_default_executor" begin
        rng = MersenneTwister(42)
        @test Pathfinder._default_executor(MersenneTwister(42); basesize=2) ===
            SequentialEx()
        if VERSION ≥ v"1.7.0"
            @test Pathfinder._default_executor(Random.GLOBAL_RNG; basesize=1) ===
                PreferParallel(; basesize=1)
            @test Pathfinder._default_executor(Random.default_rng(); basesize=1) ===
                PreferParallel(; basesize=1)
        else
            @test Pathfinder._default_executor(Random.GLOBAL_RNG; basesize=1) ===
                SequentialEx()
            @test Pathfinder._default_executor(Random.default_rng(); basesize=1) ===
                SequentialEx()
        end
    end

    @testset "_findmax" begin
        x = randn(100)
        @test Pathfinder._findmax(x) == findmax(x)
        @test Pathfinder._findmax(x |> Map(sin)) == findmax(sin.(x))
        @test Pathfinder._findmax([NaN, 3.0, 1.0]) === (3.0, 2)
        @test Pathfinder._findmax([NaN, NaN, NaN]) === (NaN, 1)
        @test Pathfinder._findmax([2.0, NaN, 4.0]) === (4.0, 3)
    end
end


text:
const XPR = Xpress
using LinQuadOptInterface
const LQOI = LinQuadOptInterface
const MOI = LQOI.MOI

const SUPPORTED_OBJECTIVES = [
    LQOI.SinVar,
    LQOI.Linear,
    LQOI.Quad
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.Linear, LQOI.IV),
    (LQOI.Quad, LQOI.EQ),
    (LQOI.Quad, LQOI.LE),
    (LQOI.Quad, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, LQOI.SOS1),
    # (LQOI.VecVar, LQOI.SOS2),
    # (LQOI.SinVar, MOI.Semicontinuous{Float64}),
    # (LQOI.SinVar, MOI.Semiinteger{Float64}),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct Optimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase(Model)
    params::Dict{Any,Any}
    l_rows::Vector{Int}
    q_rows::Vector{Int}
    conflict::Union{Nothing, IISData}
    Optimizer(::Nothing) = new()
end

LQOI.LinearQuadraticModel(::Type{Optimizer}, env) = XPR.Model()

function Optimizer(; kwargs...)
    env = nothing
    m = Optimizer(nothing)
    m.params = Dict{Any,Any}()
    MOI.empty!(m)
    for (name,value) in kwargs
        m.params[name] = value
        XPR.setparam!(m.inner, XPR.XPRS_CONTROLS_DICT[name], value)
    end
    m.conflict = nothing
    return m
end

function MOI.empty!(m::Optimizer)
    MOI.empty!(m,nothing)
    # m.constraint_primal_solution = m.qconstraint_primal_solution
    # m.constraint_dual_solution = m.qconstraint_dual_solution
    m.l_rows = Int[]
    m.q_rows = Int[]
    m.conflict = nothing
    for (name,value) in m.params
        XPR.setparam!(m.inner, XPR.XPRS_CONTROLS_DICT[name], value)
    end
end

LQOI.supported_constraints(s::Optimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(s::Optimizer) = SUPPORTED_OBJECTIVES

backend_type(m::Optimizer, ::MOI.GreaterThan{T}) where T = Cchar('G')
backend_type(m::Optimizer, ::MOI.LessThan{T}) where T    = Cchar('L')
backend_type(m::Optimizer, ::MOI.EqualTo{T}) where T     = Cchar('E')
# Implemented separately
# backend_type(m::Optimizer, ::MOI.Interval{T}) where T    = Cchar('R')

backend_type(m::Optimizer, ::MOI.Zeros)        = Cchar('E')
backend_type(m::Optimizer, ::MOI.Nonpositives) = Cchar('L')
backend_type(m::Optimizer, ::MOI.Nonnegatives) = Cchar('G')

#=
    not in LinQuad
=#

setparam!(instance::Optimizer, name, val) = XPR.setparam!(instance.inner, XPR.XPRS_CONTROLS_DICT[name], val)

setlogfile!(instance::Optimizer, path) = XPR.setlogfile(instance.inner, path::String)

cintvec(v::Vector) = convert(Vector{Int32}, v)

#=
    Constraints
=#

LQOI.change_variable_bounds!(instance::Optimizer, colvec, valvec, sensevec) = XPR.chgbounds!(instance.inner, cintvec(colvec), sensevec, valvec)

LQOI.get_variable_lowerbound(instance::Optimizer, col) = XPR.get_lb(instance.inner, col, col)[1]
LQOI.get_variable_upperbound(instance::Optimizer, col) = XPR.get_ub(instance.inner, col, col)[1]

LQOI.get_number_linear_constraints(instance::Optimizer) = XPR.num_linconstrs(instance.inner)

LQOI.get_number_quadratic_constraints(instance::Optimizer) = XPR.num_qconstrs(instance.inner)

# function LQOI.shift_references_after_delete_quadratic!(instance::Optimizer, row)
#     LQOI.shift_references_after_delete_quadratic_base!(instance, row)
#     LQOI.shift_references_after_delete_affine_base!(instance, row)
# end
# function LQOI.shift_references_after_delete_affine!(instance::Optimizer, row)
#     LQOI.shift_references_after_delete_quadratic_base!(instance, row)
#     LQOI.shift_references_after_delete_affine_base!(instance, row)
# end

function LQOI.add_linear_constraints!(instance::Optimizer, A::LQOI.CSRMatrix{Float64}, sensevec, rhsvec)
    newrows = length(rhsvec)
    rows = XPR.num_constrs(instance.inner)
    addedrows = collect((rows+1):(rows+newrows))
    # quadind
    append!(instance.l_rows,addedrows)
    XPR.add_constrs!(instance.inner, A.row_pointers, A.columns, A.coefficients, sensevec, rhsvec)
end

function LQOI.add_ranged_constraints!(instance::Optimizer, A::LQOI.CSRMatrix{Float64}, lowerbound, upperbound)
    newrows = length(lowerbound)
    rows = XPR.num_constrs(instance.inner)
    addedrows = collect((rows+1):(rows+newrows))
    # quadind
    append!(instance.l_rows,addedrows)
    sensevec = fill(Cchar('E'),newrows)
    XPR.add_constrs!(instance.inner, A.row_pointers, A.columns, A.coefficients, sensevec, upperbound)
    XPR.chg_rhsrange!(instance.inner, cintvec(addedrows), upperbound .- lowerbound)
end

function LQOI.modify_ranged_constraints!(instance::Optimizer, rows::Vector{Int}, lowerbound::Vector{Float64}, upperbound::Vector{Float64})
    # quadind
    _rows = instance.l_rows[rows]
    XPR.set_rhs!(instance.inner, _rows, upperbound)
    XPR.chg_rhsrange!(instance.inner, cintvec(_rows), upperbound .- lowerbound)
end

function LQOI.get_rhs(instance::Optimizer, row)
    # quadind
    _row = instance.l_rows[row]
    return XPR.get_rhs(instance.inner, _row, _row)[1]
end
function LQOI.get_quadratic_rhs(instance::Optimizer, row)
    # quadind
    _row = instance.q_rows[row]
    return XPR.get_rhs(instance.inner, _row, _row)[1]
end

function LQOI.get_range(instance::Optimizer, row)
    # quadind
    _row = instance.l_rows[row]
    ub = XPR.get_rhs(instance.inner, _row, _row)[1]
    r = XPR.get_rhsrange(instance.inner, _row, _row)[1]
    return ub .- r, ub
end

# TODO improve
function LQOI.get_linear_constraint(instance::Optimizer, row)
    # quadind
    _row = instance.l_rows[row]
    A = sparse(XPR.get_rows(instance.inner, _row, _row)')
    return rowvals(A), nonzeros(A)
end

function LQOI.get_quadratic_constraint(instance::Optimizer, row)
    # quadind
    _row = instance.q_rows[row]
    A = sparse(XPR.get_rows(instance.inner, _row, _row)')
    Q = XPR.get_qrowmatrix(instance.inner, _row)
    return rowvals(A), nonzeros(A), Q
end

# notin LQOI
# TODO improve
function getcoef(instance::Optimizer, row, col)
    # quadind
    _row = instance.l_rows[row]
    A = sparse(XPR.get_rows(instance.inner, _row, _row)')
    cols = rowvals(A)
    vals = nonzeros(A)

    pos = findfirst(cols, col)
    if pos > 0
        return vals[pos]
    else
        return 0.0
    end
end

function LQOI.change_matrix_coefficient!(instance::Optimizer, row, col, coef)
    # quadind
    _row = instance.l_rows[row]
    XPR.chg_coeffs!(instance.inner, _row, col, coef)
end

function LQOI.change_objective_coefficient!(instance::Optimizer, col, coef)
    XPR.set_objcoeffs!(instance.inner, Int32[col], Float64[coef])
end
function LQOI.change_rhs_coefficient!(instance::Optimizer, row, coef)
    # quadind
    _row = instance.l_rows[row]
    XPR.set_rhs!(instance.inner, Int32[_row], Float64[coef])
end
function LQOI.delete_linear_constraints!(instance::Optimizer, rowbeg, rowend)
    _rows = instance.l_rows[collect(rowbeg:rowend)]
    XPR.del_constrs!(instance.inner, cintvec(_rows))
    for i in rowend:-1:rowbeg
        val = instance.l_rows[i]
        deleteat!(instance.l_rows,i)
        for j in eachindex(instance.l_rows)
            if instance.l_rows[j] > val
                instance.l_rows[j] -= 1
            end
        end
    end
end
function LQOI.delete_quadratic_constraints!(instance::Optimizer, rowbeg, rowend)
    _rows = instance.q_rows[collect(rowbeg:rowend)]
    XPR.del_constrs!(instance.inner, cintvec(_rows))
    for i in rowend:-1:rowbeg
        val = instance.q_rows[i]
        deleteat!(instance.q_rows,i)
        for j in eachindex(instance.q_rows)
            if instance.q_rows[j] > val
                instance.q_rows[j] -= 1
            end
        end
    end
end
LQOI.change_variable_types!(instance::Optimizer, colvec, typevec) = XPR.chgcoltypes!(instance.inner, colvec, typevec)

function LQOI.change_linear_constraint_sense!(instance::Optimizer, rowvec, sensevec)
    _rows = instance.l_rows[rowvec]
    XPR.set_rowtype!(instance.inner, _rows, sensevec)
end
LQOI.add_sos_constraint!(instance::Optimizer, colvec, valvec, typ) = XPR.add_sos!(instance.inner, typ, colvec, valvec)

LQOI.delete_sos!(instance::Optimizer, idx1, idx2) = XPR.del_sos!(instance.inner, cintvec(collect(idx1:idx2)))

# TODO improve getting processes
function LQOI.get_sos_constraint(instance::Optimizer, idx)
    A, types = XPR.get_sos_matrix(instance.inner)
    line = A[idx,:] #sparse vec
    cols = line.nzind
    vals = line.nzval
    typ = types[idx] == Cchar('1') ? :SOS1 : :SOS2
    return cols, vals, typ
end

function LQOI.add_quadratic_constraint!(instance::Optimizer, cols, coefs, rhs, sense, I, J, V)
    @assert length(I) == length(J) == length(V)
    XPR.add_qconstr!(instance.inner, cols, coefs, I, J, V/2, sense, rhs)
    push!(instance.q_rows, num_constrs(instance.inner))
end

#=
    Objective
=#

function LQOI.set_quadratic_objective!(instance::Optimizer, I, J, V)
    @assert length(I) == length(J) == length(V)
    XPR.delq!(instance.inner)
    XPR.add_qpterms!(instance.inner, I, J, V)
    return nothing
end

function LQOI.set_linear_objective!(instance::Optimizer, colvec, coefvec)
    nvars = XPR.num_vars(instance.inner)
    obj = zeros(Float64, nvars)

    for i in eachindex(colvec)
        obj[colvec[i]] += coefvec[i]
    end

    XPR.set_obj!(instance.inner, obj)
    nothing
end

function LQOI.change_objective_sense!(instance::Optimizer, symbol)
    if symbol == :min
        XPR.set_sense!(instance.inner, :minimize)
    else
        XPR.set_sense!(instance.inner, :maximize)
    end
end

function LQOI.get_linear_objective!(instance::Optimizer, x::Vector{Float64})
    obj = XPR.get_obj(instance.inner)
    @assert length(x) == length(obj)
    for i in 1:length(obj)
        x[i] = obj[i]
    end
end

function LQOI.get_quadratic_terms_objective(instance::Optimizer)
    Q = XPR.getq(instance.inner)
    return Q
end

function LQOI.get_objectivesense(instance::Optimizer)
    s = XPR.model_sense(instance.inner)
    if s == :maximize
        return MOI.MaxSense
    else
        return MOI.MinSense
    end
end

#=
    Variables
=#

LQOI.get_number_variables(instance::Optimizer) = XPR.num_vars(instance.inner)

LQOI.add_variables!(instance::Optimizer, int) = XPR.add_cvars!(instance.inner, zeros(int))

LQOI.delete_variables!(instance::Optimizer, col, col2) = XPR.del_vars!(instance.inner, col)

# LQOI.lqs_addmipstarts(m, colvec, valvec)
function LQOI.add_mip_starts!(instance::Optimizer, colvec, valvec)
    x = zeros(XPR.num_vars(instance.inner))
    for i in eachindex(colvec)
        x[colvec[i]] = valvec[i]
    end
    XPR.loadbasis(instance.inner, x)
end

#=
    Solve
=#

LQOI.solve_mip_problem!(instance::Optimizer) = XPR.mipoptimize(instance.inner)

LQOI.solve_quadratic_problem!(instance::Optimizer) = LQOI.solve_linear_problem!(instance)

LQOI.solve_linear_problem!(instance::Optimizer) = XPR.lpoptimize(instance.inner)

function LQOI.get_termination_status(instance::Optimizer)
    stat_lp = XPR.get_lp_status2(instance.inner)
    if XPR.is_mip(instance.inner)
        stat_mip = XPR.get_mip_status2(instance.inner)
        if stat_mip == XPR.MIP_NotLoaded
            return MOI.OTHER_ERROR
        elseif stat_mip == XPR.MIP_LPNotOptimal
            # MIP search incomplete but there is no linear sol
            return MOI.OTHER_ERROR
        elseif stat_mip == XPR.MIP_NoSolFound
            # MIP search incomplete but there is no integer sol
            other = xprsmoi_stopstatus(instance.inner)
            if other == MOI.OTHER_ERROR
                return MOI.SLOW_PROGRESS#OtherLimit
            else
                return other
            end

        elseif stat_mip == XPR.MIP_Solution
            # MIP search incomplete but there is a solution
            other = xprsmoi_stopstatus(instance.inner)
            if other == MOI.OTHER_ERROR
                return MOI.OTHER_LIMIT
            else
                return other
            end

        elseif stat_mip == XPR.MIP_Infeasible
            return MOI.INFEASIBLE
        elseif stat_mip == XPR.MIP_Optimal
            return MOI.OPTIMAL
        elseif stat_mip == XPR.MIP_Unbounded
            return MOI.DUAL_INFEASIBLE
        end
        return MOI.OTHER_ERROR
    else
        if stat_lp == XPR.LP_Unstarted
            return MOI.OTHER_ERROR
        elseif stat_lp == XPR.LP_Optimal
            return MOI.OPTIMAL
        elseif stat_lp == XPR.LP_Infeasible
            return MOI.INFEASIBLE
        elseif stat_lp == XPR.LP_CutOff
            return MOI.OBJECTIVE_LIMIT
        elseif stat_lp == XPR.LP_Unfinished
            return xprsmoi_stopstatus(instance.inner)
        elseif stat_lp == XPR.LP_Unbounded
            return MOI.DUAL_INFEASIBLE
        elseif stat_lp == XPR.LP_CutOffInDual
            return MOI.OBJECTIVE_LIMIT
        elseif stat_lp == XPR.LP_Unsolved
            return MOI.OTHER_ERROR
        elseif stat_lp == XPR.LP_NonConvex
            return MOI.INVALID_MODEL
        end
        return MOI.OTHER_ERROR
    end
end

function xprsmoi_stopstatus(instance::Model)
    ss = XPR.get_stopstatus(instance)
    if ss == XPR.StopTimeLimit
        return MOI.TIME_LIMIT
    elseif ss == XPR.StopControlC
        return MOI.INTERRUPTED
    elseif ss == XPR.StopNodeLimit
        # should not be here
        Compat.@warn("should not be here")
        return MOI.NODE_LIMIT
    elseif ss == XPR.StopIterLimit
        return MOI.ITERATION_LIMIT
    elseif ss == XPR.StopMIPGap
        return MOI.OBJECTIVE_LIMIT
    elseif ss == XPR.StopSolLimit
        return MOI.SOLUTION_LIMIT
    elseif ss == XPR.StopUser
        return MOI.INTERRUPTED
    end
    return MOI.OTHER_ERROR
end

function LQOI.get_primal_status(instance::Optimizer)
    if XPR.is_mip(instance.inner)
        stat_mip = XPR.get_mip_status2(instance.inner)
        if stat_mip in [XPR.MIP_Solution, XPR.MIP_Optimal]
            return MOI.FEASIBLE_POINT
        elseif stat_mip == XPR.MIP_Unbounded && XPR.hasprimalray(instance.inner)
            return MOI.INFEASIBILITY_CERTIFICATE
        elseif stat_mip in [XPR.MIP_LPOptimal, XPR.MIP_NoSolFound]
            return MOI.INFEASIBLE_POINT
        end
        return MOI.NO_SOLUTION
    else
        stat_lp = XPR.get_lp_status2(instance.inner)
        if stat_lp == XPR.LP_Optimal
            return MOI.FEASIBLE_POINT
        elseif stat_lp == XPR.LP_Unbounded && XPR.hasprimalray(instance.inner)
            return MOI.INFEASIBILITY_CERTIFICATE
        # elseif stat_lp == LP_Infeasible
        #     return MOI.InfeasiblePoint - xpress wont return
        # elseif cutoff//cutoffindual ???
        else
            return MOI.NO_SOLUTION
        end
    end
end

function LQOI.get_dual_status(instance::Optimizer)
    if XPR.is_mip(instance.inner)
        return MOI.NO_SOLUTION
    else
        stat_lp = XPR.get_lp_status2(instance.inner)
        if stat_lp == XPR.LP_Optimal
            return MOI.FEASIBLE_POINT
        elseif stat_lp == XPR.LP_Infeasible && XPR.hasdualray(instance.inner)
            return MOI.INFEASIBILITY_CERTIFICATE
        # elseif stat_lp == LP_Unbounded
        #     return MOI.InfeasiblePoint - xpress wont return
        # elseif cutoff//cutoffindual ???
        else
            return MOI.NO_SOLUTION
        end
    end
end

LQOI.get_variable_primal_solution!(instance::Optimizer, place) = XPR.get_solution!(instance.inner, place)

function LQOI.get_linear_primal_solution!(instance::Optimizer, place)
    if num_qconstrs(instance.inner) == 0
        XPR.get_slack_lin!(instance.inner, place)
        rhs = XPR.get_rhs(instance.inner)
        for i in eachindex(place)
            place[i] = -place[i]+rhs[i]
        end
        return nothing
    else
        XPR.get_slack_lin!(instance.inner, place)
        rhs = XPR.get_rhs(instance.inner)
        lrows = XPR.get_lrows(instance.inner)
        for (i,v) in enumerate(lrows)
            place[i] = -place[i]+rhs[v]
        end
    end
    nothing
end

function LQOI.get_quadratic_primal_solution!(instance::Optimizer, place)
    if num_qconstrs(instance.inner) == 0
        return nothing
    else
        qrows = XPR.get_qrows(instance.inner)
        newplace = zeros(num_constrs(instance.inner))
        XPR.get_slack!(instance.inner, newplace)
        rhs = XPR.get_rhs(instance.inner)
        for (i,v) in enumerate(qrows)
            place[i] = -newplace[v]+rhs[v]
        end
    end
    nothing
end

LQOI.get_variable_dual_solution!(instance::Optimizer, place) = XPR.get_reducedcost!(instance.inner, place)

LQOI.get_linear_dual_solution!(instance::Optimizer, place) = XPR.get_dual_lin!(instance.inner, place)

function LQOI.get_quadratic_dual_solution!(instance::Optimizer, place)
    if num_qconstrs(instance.inner) == 0
        return nothing
    else
        qrows = XPR.get_qrows(instance.inner)
        newplace = zeros(num_constrs(instance.inner))
        XPR.get_dual!(instance.inner, newplace)
        for (i,v) in enumerate(qrows)
            place[i] = newplace[v]
        end
    end
end

LQOI.get_objective_value(instance::Optimizer) = XPR.get_objval(instance.inner)

function LQOI.get_objective_bound(instance::Optimizer)
    if XPR.is_mip(instance.inner)
        return XPR.get_bestbound(instance.inner)+LQOI.get_constant_objective(instance)
    else
        return XPR.get_objval(instance.inner)+LQOI.get_constant_objective(instance)
    end
end
function LQOI.get_relative_mip_gap(instance::Optimizer)
    L = XPR.get_mip_objval(instance.inner)
    U = XPR.get_bestbound(instance.inner)
    return abs(U-L)/U
end

LQOI.get_iteration_count(instance::Optimizer)  = XPR.get_simplex_iter_count(instance.inner)

LQOI.get_barrier_iterations(instance::Optimizer) = XPR.get_barrier_iter_count(instance.inner)

LQOI.get_node_count(instance::Optimizer) = XPR.get_node_count(instance.inner)

LQOI.get_farkas_dual!(instance::Optimizer, place) = XPR.getdualray!(instance.inner, place)

LQOI.get_unbounded_ray!(instance::Optimizer, place) = XPR.getprimalray!(instance.inner, place)

finalize(m::Optimizer) = XPR.free_model(m.inner)

function MOI.write_to_file(instance::Optimizer, lp_file_name::String)
    flags = ""
    if endswith(lp_file_name, ".lp")
        flags = "l"
    end
    XPR.write_model(instance.inner, filename, flags)
end

"""
    compute_conflict(model::Optimizer)

Compute a minimal subset of the constraints and variables that keep the model
infeasible.

See also `Xpress.ConflictStatus` and `Xpress.ConstraintConflictStatus`.

Note that if `model` is modified after a call to `compute_conflict`, the
conflict is not purged, and any calls to the above attributes will return values
for the original conflict without a warning.
"""
function compute_conflict(model::Optimizer)
    model.conflict = getfirstiis(model.inner)
    return
end

function _ensure_conflict_computed(model::Optimizer)
    if model.conflict === nothing
        error("Cannot access conflict status. Call `Xpress.compute_conflict(model)` first. " *
              "In case the model is modified, the computed conflict will not be purged.")
    end
end

"""
    ConflictStatus()

Return an `MOI.TerminationStatusCode` indicating the status of the last computed conflict.

If a minimal conflict is found, it will return `MOI.OPTIMAL`. If the problem is feasible, it will
return `MOI.INFEASIBLE`. If `compute_conflict` has not been called yet, it will return
`MOI.OPTIMIZE_NOT_CALLED`.
"""
struct ConflictStatus <: MOI.AbstractModelAttribute  end
MOI.is_set_by_optimize(::ConflictStatus) = true

function MOI.get(model::Optimizer, ::ConflictStatus)
    if model.conflict === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    elseif model.conflict.stat == 0
        return MOI.OPTIMAL
    elseif model.conflict.stat == 1
        return MOI.INFEASIBLE
    else
        return MOI.OTHER_LIMIT
    end
end

function MOI.supports(::Optimizer, ::ConflictStatus)
    return true
end

"""
    ConstraintConflictStatus()

A Boolean constraint attribute indicating whether the constraint participates in the last computed conflict.
"""
struct ConstraintConflictStatus <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::ConstraintConflictStatus) = true

function MOI.get(model::Optimizer, ::ConstraintConflictStatus, index::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction, <:Union{LQOI.LE, LQOI.GE, LQOI.EQ}})
    _ensure_conflict_computed(model)
    return (model[index] - 1) in model.conflict.miisrow
end

function MOI.supports(::Optimizer, ::ConstraintConflictStatus, ::Type{MOI.ConstraintIndex{<:MOI.ScalarAffineFunction, <:Union{LQOI.LE, LQOI.GE, LQOI.EQ}}})
    return true
end


text:
@testset "Macros" begin
  @testset "@estimsolver" begin
    s = ESolver(:z => (A=2,B=3), C=false)
    @test s.C == false
    @test variables(s) == [:z]
    names, params = covariables(:z, s)
    @test names == (:z,)
    @test params[(:z,)].A == 2
    @test params[(:z,)].B == 3

    @test sprint(show, s) == "ESolver"
    @test sprint(show, MIME"text/plain"(), s) == "ESolver\n  └─z\n    └─A ⇨ 2\n    └─B ⇨ 3\n"
  end

  @testset "@simsolver" begin
    s = SSolver(:z => (A=2,B=3), C=false)
    @test s.C == false
    @test variables(s) == [:z]
    names, params = covariables(:z, s)
    @test names == (:z,)
    @test params[(:z,)].A == 2
    @test params[(:z,)].B == 3

    @test sprint(show, s) == "SSolver"
    @test sprint(show, MIME"text/plain"(), s) == "SSolver\n  └─z\n    └─A ⇨ 2\n    └─B ⇨ 3\n"

    s = SSolver()
    @test s.C == true
    @test isempty(variables(s))
    names, params = covariables(:z, s)
    @test names == (:z,)
    @test params[(:z,)].A == 1.0
    @test params[(:z,)].B == 2

    s = SSolver(:z=>(A=1,B=2), :w=>(A=2,B=3), (:z,:w)=>(J="bar",))
    @test s.C == true
    @test Set(variables(s)) == Set([:z,:w])
    names, params = covariables(:z, s)
    @test Set(names) == Set((:z,:w))
    @test params[(:z,)].A == 1
    @test params[(:z,)].B == 2
    @test params[(:w,)].A == 2
    @test params[(:w,)].B == 3
    @test params[(:z,:w)].J == "bar"
  end
end


text:
import Dolo
import Dolo: Point,ListOfPoints
using StaticArrays

model = Dolo.yaml_import(joinpath(Dolo.pkg_path, "examples","models","rbc.yaml"));

@time sol_iti = Dolo.improved_time_iteration(model; verbose=true);
@time sol_iti = Dolo.improved_time_iteration(model; verbose=true);

@time sol_ti = Dolo.time_iteration(model, verbose=true; tol_η=1e-10, maxit=2000)
exit()
# @profile Dolo.improved_time_iteration(model);
#
# #
# f = open("out.jll", "r")
# out = deserialize(f)
# close(f)
#
#
# res, dres, jres, fut_S = out

using JLD
ss = load("myfile.jld")
res = ss["res"]
dres = ss["dres"]
jres = ss["jres"]
fut_S = ss["fut_S"]

module NewModule

    import Dolo
    import Dolo: Point,ListOfPoints
    using StaticArrays

end

import NewModule

R_i, D_i, J_ij, S_ij, M_ij = NewModule.reorder_data(dprocess, res, dres, jres, fut_S)

ddr_filt = Dolo.CachedDecisionRule(dprocess, grid, R_i)

@time tot0, n0, lam0 = Dolo.invert_jac(res,dres,jres,fut_S,ddr_filt)

n0

@time tot, n, lam, errors= NewModule.invert_jac(res,dres,jres,fut_S,ddr_filt)

print("Invert me many times")

function invert_me_1(N, res,dres,jres,fut_S,ddr_filt;maxit=1000)
    tt = 0
    for n=1:N
        tot0, n0, lam0 = Dolo.invert_jac(res,dres,jres,fut_S,ddr_filt;maxit=maxit)
        tt += n0
    end
    return tt
end

function invert_me_2(N, res,dres,jres,fut_S,ddr_filt;maxit=1000)
    tt = 0
    for n=1:N
        tot0, n0, lam0 = NewModule.invert_jac(res,dres,jres,fut_S,ddr_filt;maxit=maxit)
        tt += n0
    end
    return tt
end

@time invert_me_1(1,res,dres,jres,fut_S,ddr_filt;maxit=2)
@time invert_me_2(1,res,dres,jres,fut_S,ddr_filt;maxit=2)

@time invert_me_1(1,res,dres,jres,fut_S,ddr_filt;maxit=1000)
@time invert_me_2(1,res,dres,jres,fut_S,ddr_filt;maxit=1000)


text:
module LLLplus

using Base.BLAS
using ArrayViews

export
    lll,
    lllnative,
    seysen,
    vblast,
    hard_sphere

include("lll.jl")
include("lllnative.jl")
include("seysen.jl")
include("vblast.jl")
include("hard_sphere.jl")

end # LLLplus module


text:
module SimSearchManifoldLearning

using Arpack
using Distances
using LinearAlgebra
using LsqFit: curve_fit
using SimilaritySearch
using SparseArrays

import StatsAPI: fit, predict
export fit, predict

include("utils.jl")
include("precomputedknns.jl")
include("layouts.jl")
include("embeddings.jl")
include("umap_.jl")
include("mlapi.jl")

end # module


text:
import FiniteDiff, NLsolve

mutable struct BVPJacobianWrapper{LossType} <: Function
    loss::LossType
end
(p::BVPJacobianWrapper)(resid,u) = p.loss(resid,u)
(p::BVPJacobianWrapper)(u) = (resid = similar(u); p.loss(resid,u); resid)

function ConstructJacobian(f!::BVPJacobianWrapper, S::BVPSystem, y)
    jac_cache = FiniteDiff.JacobianCache(
                                  similar(y),similar(y),similar(y))
    function fj!(F, J, x)
        FiniteDiff.finite_difference_jacobian!(J, f!, x, jac_cache)
        f!(F,x)
    end
    j!(J, x)  = (F = jac_cache.fx; fj!(F, J, x))
    J0 = BandedMatrix(Zeros{eltype(S.y[1])}(S.M*S.N, S.M*S.N), (S.M-1, S.M-1))
    NLsolve.OnceDifferentiable(f!.loss, j!, fj!, jac_cache.x1, jac_cache.fx, sparse(J0))
end

function ConstructJacobian(f!::BVPJacobianWrapper, y)
    jac_cache = FiniteDiff.JacobianCache(
                                    similar(y),similar(y),similar(y))
    function fj!(F, J, x)
        FiniteDiff.finite_difference_jacobian!(J, f!, x, jac_cache)
        f!(F,x)
    end
    j!(J, x)  = (F = jac_cache.fx; fj!(F, J, x))
    NLsolve.OnceDifferentiable(f!.loss, j!, fj!, jac_cache.x1, jac_cache.fx)
end


text:
using Dates, AxisArrays

@testset "Functions" begin
# findmax
d = collect(DateTime(2003,1,1):Day(1):DateTime(2005,12,31))

data = Array{Float64,3}(undef, 2, 2,1096)
data[1,1,:] = collect(1.0:1096.0); data[1,2,:] = collect(1.0:1096.0); data[2,1,:]=collect(1.0:1096.0); data[2,2,:] = collect(1.0:1096.0);
data[1,1,1] = NaN

# ClimGrid based tests
x = 1:2
y = 1:2
axisdata = AxisArray(data, Axis{:lon}(x), Axis{:lat}(y), Axis{:time}(d))
X, Y = ndgrid(x, y)
C = ClimateTools.ClimGrid(axisdata, variable = "pr", longrid=X, latgrid=Y)
@test isnan(findmin(C)[1])
@test findmin(C)[2] == CartesianIndex(1,1,1)
@test isnan(findmax(C)[1])
@test findmax(C)[2] == CartesianIndex(1,1,1)
@test findmin(C, skipnan=true) == (1.0, CartesianIndex(2,1,1))
@test findmax(C, skipnan=true) == (1096.0, CartesianIndex(1,1,1096))

# Monthly means
D = monthmean(C)
@test size(D[1]) == (2, 2, 36)
@test isnan(D[1][1, 1, 1])
@test D[1][2,1,1] == 16.0
@test D[1][2,2,1] == 16.0
@test D[1][2,2,end] == 1081.0
timevec = get_timevec(D)
@test timevec[1] == DateTime(2003,01,01,00,00,00)
@test timevec[end] == DateTime(2005,12,01,00,00,00)

D = monthsum(C)
@test size(D[1]) == (2, 2, 36)
@test isnan(D[1][1, 1, 1])
@test D[1][2,1,1] == 496.0
@test D[1][2,2,1] == 496.0
@test D[1][2,2,end] == 33511.0
timevec = get_timevec(D)
@test timevec[1] == DateTime(2003,01,01,00,00,00)
@test timevec[end] == DateTime(2005,12,01,00,00,00)

# get dims
x, y, timevec = ClimateTools.getdims(C)
@test length(timevec) == length(d)
for itime = 1:length(timevec)
    @test timevec[itime] == d[itime]
end
@test x == 1:2 && y == 1:2

# get grids
longrid, latgrid = ClimateTools.getgrids(C)
@test longrid == X
@test latgrid == Y

# Extension file
@test ClimateTools.extension("test") == ""
@test ClimateTools.extension("test.nc") == ".nc"
@test ClimateTools.extension("test.") == ""

dt = DateTimeNoLeap(2000, 01, 01, 23)
@test yearmonthdayhour(dt) == (2000, 1, 1, 23)

end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ f1d0ed8f-bfd3-4210-84fd-f189bba6ebf9
begin
	using PlutoUI
	PlutoUI.TableOfContents(title = "Contents")
end

# ╔═╡ 08eadc12-f8e9-11eb-1094-b7adc6936aa1
md"""
# Permutations
"""

# ╔═╡ c4d2ff29-2579-4079-a6cb-29015ef5ecd2
md"""
## Objectives
- _**Create a custom type**_ to represent a permutation
- _**Utilise multiple dispatch to create new methods for inbuilt functions**_ with `Permutation` objects as arguments
- Generate the elements of groups from generators in `Permutation` form
"""

# ╔═╡ db919e43-6dc8-4e53-8081-d9312446117a
md"""
## What is multiple dispatch?
Multiple dispatch is a feature central to the design of Julia with incredibly powerful consequences. In a most basic sense, it allows for the same function name to be used for multiple different operations, such as `-` being used for unary minus (``x \to -x``) as well as for subtraction. More generally, if the types of the inputs are specified differently for different declarations of the function, the function can act differently depending on the inputs that it gets. An example of this inbuilt into Julia is `a * b`, which multiplies `a` and `b` if they are numbers, and concatenates `a` and `b` if they are strings:
"""

# ╔═╡ 524dee4b-6687-44de-af79-b99d3724f9ac
3 * 4

# ╔═╡ a16202d3-406a-4b43-807e-b60c673d4132
"3" * "4"

# ╔═╡ 41fc4769-5b67-4503-b1d1-2854a0d9c3a2
md"""
These, or similar examples, are duplicated notation common to most programming languages, but what Julia allows that others don't is for your own custom functions to utilise multiple dispatch as well. For example, say that I want to check if an input is four or not, but I don't mind if the input is a number or a string. To do this, I follow the argument with double colons and the type to determine when this method will run (here `Number` is an abstract type, which cannot itself be the type of an object, but encompasses all numeric types into one keyword). Note that such notation of arguments can still be used for functions with only a single methods, where it can be helpful for debugging purposes.
"""

# ╔═╡ acf690f4-b65b-4576-9145-1e894d5342f2
isfour(x::String) = (lowercase(x) ∈ ["4", "four"])

# ╔═╡ 9f193c51-b6ea-477f-ab0e-1bc5a4be4b86
md"""
`isfour` is now defined for a `Number` and a `String` input, but not currently for anything else. If I decide that nothing else can be considered as four, then I can create a third method using the keyword `Any` (the abstract type representing all possible types) which will deal with everything else:
"""

# ╔═╡ 29b729bb-cd1a-4ee1-ae2e-26b55ce744a4
isfour(x::Any) = false

# ╔═╡ 740f8275-ee0c-4aa2-b0c1-76d46fad3ecd
md"""
Of course, an input like `4` is also of type `Any`, but Julia's type system recognises that `Number` is more restrictive and so is of higher precedence, hence any `Number` input to `isfour` triggers the specialised method for that type.

More of the power of multiple dispatch is illustrated later in this case study.
"""

# ╔═╡ 73d93917-b72e-468b-970c-f34ec6d51de4
md"""
## Creating the `Permutation` type

### Declaring `Permutation` as a new type
I begin by declaring a new type with the keyword `struct`. The most basic syntax to do this is:
```julia
struct Permutation
    image::Vector{Int64}
end
```

although I will want more functionality than this. The new type has a single field `image` which is a vector of integers, representing the images of ``1, 2, \dots n`` under the permutation (where ``n`` is the maximum integer permuted). For example, the permutation ``(123)(45)`` has image ``(2, 3, 1, 5, 4)^T``, and can be created as an object by the syntax:
```julia
σ = Permutation([2,3,1,5,4])
```

Its field will also be accessible through dot notation by:
```julia
σ.image
```
"""

# ╔═╡ c18ae801-b186-44cf-9e42-83188c4781aa
md"""
There are currently some issues with the new type, primarily that at the moment any vector of integers could be the image of a `Permutation`, so I need to be able to check that an image is valid, i.e.
- The image contains all of the integers ``1, 2, ..., n`` for some ``n`` exactly once
- The image contains nothing else

I will also impose an additional condition, which is that the largest integer `n` must not be fixed, since there is no point in storing it. However, this can be amended within the code, rather than rejecting such an image vector.

The way to impose these conditions is to use an inner constructor, which is a method within the `struct` block.
```julia
struct Permutation
    image::Vector{Int64}

	function Permutation(image::Vector{Int64})
		# Inner constructor code
	end
end
```

To check that the vector `image` contains the right elements, I sort it and compare it to the vector `1:n`, where `n` is the length of `image`.
```julia
sort(image) == 1:length(image)
```

Then, I find the integers which aren't fixed by using Boolean indexing, which is indexing an array with an array of `true`s and `false`s of the same size, and taking the elements which correspond to `true`. The maximum of them is then found, with `1` added to the list since `maximum` cannot find the maximum element of an empty vector.
```julia
m = maximum([image[image .!= 1:length(image)]..., 1])
```

Now, I can make use of the inner constructor specific keyword `new` to construct the permutation, by ordering the arguments as they would be if I were calling the `Permutation` function. With an error message if the conditions are not met by `image`, this completes my inner constructor
```julia
function Permutation(image::Vector{Int64})
    sortedimage = sort(image)
 	m = maximum([image[image .!= 1:length(image)]..., 1])
    sortedimage == 1:last(sortedimage) && return new(image[1:m])
    error("not a valid permutation")
end
```
"""

# ╔═╡ 99fe42de-2afa-409f-9416-5c7671e4d0f4
md"""
One last thing I want to adjust before I am done with the `Permutation` type is to be able to construct a permutation by just giving it a list of integers and not having to wrap them up in vector form myself. This I do by utilising multiple dispatch, with an outer constructor, which is simply a new method with different inputs to construct a `Permutation` object.
```julia
Permutation(imagevals::Int64...) = Permutation([x for x ∈ imagevals])
```

Putting this all together, I get (including some code from later sections, due to the reactivity of the Pluto notebook in which this is written):
"""

# ╔═╡ b9c535e5-3317-4291-9925-696bc6cb0252
begin
	struct Permutation
    	image::Vector{Int64}

	    function Permutation(image::Vector{Int64})
    	    sortedimage = sort(image)
        	m = maximum([image[image .!= 1:length(image)]..., 1])
	        sortedimage == 1:last(sortedimage) && return new(image[1:m])
    	    error("not a valid permutation")
	    end
	end

	Permutation(imagevals::Int64...) = Permutation([x for x ∈ imagevals])

	# See section Evaluating a Permutation as a function below
	maxarg(σ::Permutation) = length(σ.image)
	(σ::Permutation)(x::Int64) = x ∈ 1:maxarg(σ) ? σ.image[x] : x

	# See section Generating groups with Permutations
	import Base.:(==)
	==(σ::Permutation, τ::Permutation) = (σ.image == τ.image)

	nothing
end

# ╔═╡ 1d56249e-d86a-4d1e-b7d2-6d63a621bfa8
isfour(x::Number) = (x == 4)

# ╔═╡ 9ab255ff-b40a-4f16-bed8-fe14e4495502
isfour(4.0)

# ╔═╡ 98750bd0-84bd-4cb1-900b-9ee26142e4ec
isfour("Twelve")

# ╔═╡ 4e4508db-47c0-462f-9cec-93e52bbffd7c
isfour([4])

# ╔═╡ de16f61d-853a-43ed-adf0-fb77ea3c8f4a
σ = Permutation(4,5,1,3,2); σ.image

# ╔═╡ 322f2dbc-f8f4-4f62-a9bb-ea6b00a333ae
md"""
### Evaluating a `Permutation` as a function
I now have a `Permutation` object, but it doesn't do anything yet. The next step is to be able to use a Permutation object as the function that it represents, i.e. evaluate it at an integer. I can do this by defining a function with its name as a `Permutation` parameter:
```julia
(σ::Permutation)(x::Int64) = σ.image[x]
```

However this will give an error for any integer not between `1` and the maximum unfixed integer (`m` in the inner constructor), since no such index exists. Instead, I want it to fix all other integers, which will make composing permutations later much easier, so I amend the code, defining a new function `maxarg` to go with it:
```julia
maxarg(σ::Permutation) = length(σ.image)
(σ::Permutation)(x::Int64) = x ∈ 1:maxarg(σ) ? σ.image[x] : x
```

This is included in the block above with the constructors. I can test the result of this on the sample permutation `σ` defined above
"""

# ╔═╡ ff8ab6b4-34de-44eb-b98a-2f5e0f1c48fc
σ(4)

# ╔═╡ 5ce89616-326b-4404-9f2f-707c91835fa3
maxarg(σ)

# ╔═╡ 9617f0f2-1b45-418f-8818-26df1abde6cb
md"""
### Creating a custom display format for a `Permutation` object
At the moment, whenever a `Permutation` object is returned, it displays in the default form for a new `struct`, which is the name of the type followed by a list of its fields (in this case `σ` would look like `Permutation([4, 5, 1, 3, 2])`). This isn’t very useful, so be more understandable (and more in keeping with standard mathematical notation), I want to customise this displayed form, which I can do with multiple dispatch by adding a method to the function `show` specifically for the `Permutation` type.

The format that I want to display the permutation in is as its disjoint cycle decomposition. This will first require a function to calculate the orbit generated by acting repeatedly on a given element:
"""

# ╔═╡ ee9d45f3-0a25-47f6-8445-8ba733c3a8f6
function orbit(σ::Permutation, x::Int64)
    orb = [x]
    y = σ(x)
    while y != x
        push!(orb,y)
        y = σ(y)
    end
    return orb
end

# ╔═╡ b119e9f0-b8ba-42a3-a76c-82ed2dc170a4
md"""
To calculate the disjoint cycle decomposition, I build up the output as follows:
- I start out with an empty `Vector{Vector{Int64}}` (that is, a vector whose elements are vectors of integers)
```julia
decomp = Vector{Int64}[]
```

- A vector `unaccounted` tracks which values in the range `1:maxarg(σ)` are yet to be added to the decomposition. I will iterate until this vector is entirely false
```julia
unaccounted = trues(maxarg(σ))

while any(unaccounted)
    # code to iterate
end
```

- Inside the loop, I will look for the first value which is unaccounted for in the decomposition so far, calculate its orbit, and then update `unaccounted` accordingly
```julia
x = findfirst(unaccounted)
xorbit = orbit(σ,x)
unaccounted[xorbit] .= false
```

- Then, if the orbit is non-trivial, I add it to the decomposition
```julia
length(xorbit) > 1 && push!(decomp,xorbit)
```

This results in `decomp` being a list of cycles exactly determining the disjoint cycle decomposition.
"""

# ╔═╡ e4ce9583-5c18-4df7-a5fb-a46670fe40af
function dcd(σ::Permutation)
    decomp = Vector{Int64}[]
    unaccounted = trues(maxarg(σ))
    while any(unaccounted)
        x = findfirst(unaccounted)
        xorbit = orbit(σ,x)
        unaccounted[xorbit] .= false
        length(xorbit) > 1 && push!(decomp,xorbit)
    end
    return decomp
end

# ╔═╡ b2ac585a-d84a-4d19-afe5-ad6d288e847f
md"""
I then need to consider how to build up the string from this decomposition
- Each cycle from the decomposition will be expressed as an open parenthesis, followed by the list of integers in order, separated by spaces, and then closed with another parenthesis. To do this, I list out strings that I can then concatenate with `*`.
```julia
*( "( ", ["$y " for y ∈ x]... , ")" )
```

- This needs to be repeated for all of the cycles in the decomposition, with an additional empty string added in case the disjoint cycle decomposition is empty.
```julia
toprint = *([ *( "( ", ["$y " for y ∈ x]... , ")" ) for x ∈ dcd(σ)]..., "")
```

- Finally, if (and only if) the permutation is the identity permutation, then the string will be empty at this point. Instead, I want to use the symbol `ι`.
```julia
toprint == "" && (toprint = "ι")
```

In order to write my new method for `show`, I first need to `import` it, ensuring that I do not overwrite the inbuilt methods that it has. Then, I construct my function, mirroring the syntax used for `show` in the source code as it is good practice.
"""

# ╔═╡ 7a96d2c8-aa9c-405f-b34f-1f3c7c53b6c6
begin
	import Base.show

	function show(io::IO,σ::Permutation)
	    toprint = *([ *( "( ", ["$y " for y ∈ x]... , ")" ) for x ∈ dcd(σ)]..., "")
    	toprint == "" && (toprint = "ι")
	    print(io, toprint)
	end
end

# ╔═╡ 12346644-107b-4c43-aaca-3e404fecc8b3
md"""
Permutations now have a much more aesthetic and useful display style:
"""

# ╔═╡ 0afa030a-58ad-44f6-a10f-0ae3baabe6f8
σ

# ╔═╡ 1e8659dc-27d5-429d-8431-d5cfa9006411
md"""
### Adding `Permutation` arithmetic
I now want to add some arithmetic for combining `Permutation` objects, starting with a shortcut for the identity. This will be represented by the constant `ι` (mirroring the way it is displayed by `show`), and is given by:
"""

# ╔═╡ 068a9d15-ce46-4aa5-9c13-3f79e70b0eda
const ι = Permutation(1)

# ╔═╡ b3151f52-6dd7-4f63-a158-becb5d1aec2e
md"""
Then, I add two methods to the inbuilt function one allowing it to be obtained either by passing a permutation, or the `Permutation` type. For these functions, the input variables are not given names, only their types specified, since the type is the only relevant property about the input, as needed for multiple dispatch.
"""

# ╔═╡ d614739d-b22f-4991-9722-e8371732bf4b
begin
	import Base.one
	one(::Permutation) = ι
	one(::Type{Permutation}) = ι
end

# ╔═╡ 363481da-d185-4846-82ab-c227e9f668f8
md"""
Next, I want to be able to compose permutations, for which I will write new `Permutation`-specific methods for the inbuilt operator `∘`. In order for an arbitrary number of permutations to be composable, I use an inductive definition which mirrors the definition of composition of functions [from the Julia source code](https://github.com/JuliaLang/julia/blob/master/base/operators.jl#L1096), that is:
```julia
∘(f) = f
∘(f, g) = ComposedFunction(f, g)
∘(f, g, h...) = ∘(f ∘ g, h...)
```

Hence, my methods for `∘` are:
"""

# ╔═╡ 8f925ec6-7cd8-49fe-a467-807711b3425b
begin
	import Base.∘
	∘(σ::Permutation) = σ
	∘(σ::Permutation, τ::Permutation) =
        Permutation([σ(τ(x)) for x ∈ 1:max(maxarg(σ),maxarg(τ))])
	∘(σ::Permutation, τ::Permutation, υ::Permutation...) = ∘(σ ∘ τ, υ...)
end

# ╔═╡ 4552e49c-af66-4fd3-a827-28284215f4f9
md"""
Thirdly, I want an inverse function to be able to find the inverse of a permutation. To do this, I need to find the index of each of `1:maxarg(σ)` in the image vector and let that be the image of the new permutation, which is done by:
"""

# ╔═╡ 0e710b07-b5db-4ea7-986b-8eb7fc0946b1
begin
	import Base.inv
	inv(σ::Permutation) = Permutation([findfirst(σ.image .== x) for x ∈ 1:maxarg(σ)])
end

# ╔═╡ 7c082457-980c-4961-bb10-d91b7248d0f6
md"""
Finally, I want to be able to exponentiate by any integer (including zero and negative integers), which can be done using some logic combined with the three operations above:
"""

# ╔═╡ 59075aa4-dd35-47af-b4d1-40e89b2e599f
begin
	import Base.^

	function ^(σ::Permutation,n::Int64)
	    n == 0 && return one(Permutation)
	    n < 0 && ((n,σ) = (-n,inv(σ)))
	    return ∘(fill(σ,n)...)
	end
end

# ╔═╡ 764fb693-f4ac-4327-93a3-2171621ad039
md"""
I can now test these out:
"""

# ╔═╡ 75e5491a-6536-44b0-a5b9-10cf7f9e5c66
ι

# ╔═╡ 9d43f21c-c60a-4bf8-b1a0-bf3b56b3a97b
σ

# ╔═╡ b421f174-d554-47a8-a72d-4d5f3ec6ba2e
τ = Permutation(4,5,2,3,1)

# ╔═╡ ff82e088-aa60-49fe-bb12-fe39b0430d51
σ ∘ τ

# ╔═╡ 7c0c4e99-ebff-409b-87f3-919a6870cccc
σ^-1

# ╔═╡ 029a0c9f-5d32-409d-b6ad-9e0d4e19e854
τ^5

# ╔═╡ c75f7100-778e-44c6-9d43-675d11175673
md"""
## Generating groups with `Permutation` objects

### Constructing symmetric and alternating groups
The first groups that I would like to be able to list the elements of are symmetric or alternating groups. Before this however, I will define two more functions which will be useful to me in their construction.

First, I would like a way to quickly construct a transposition:
"""

# ╔═╡ 0598e94b-a043-4106-bb7b-d032eb67be0b
function transposition(m::Int64, n::Int64)
    image = collect(1:max(m,n))
    image[m], image[n] = n, m
    return Permutation(image)
end

# ╔═╡ 00a61c28-230d-41b0-81ac-e2f0869d6640
transposition(3,4)

# ╔═╡ 32b52a58-6086-4c2f-a833-e1a7efb69d29
transposition(3,3)

# ╔═╡ 84bd5d92-00eb-4e66-a1c5-c86626b8f7d2
md"""
Note that `transposition(m,m)` returns the identity, which turns out to be exactly what I want.

To find the alternating group, I will need to consider the parity of elements. The easiest way to find the parity of a permutation is to consider it as a product of cycles and consider their parities. I already have such a decomposition, given by `dcd`, so using the `iseven` function, and `sum` to count the number of `true` elements of an array, I can write such a function:
"""

# ╔═╡ 1b10a694-a310-4a4e-9831-654e3e9691d1
parity(σ::Permutation) = (-1)^sum([iseven(length(x)) for x ∈ dcd(σ)])

# ╔═╡ 9812e6fb-9308-4c27-9d57-9fbf71474343
parity(ι)

# ╔═╡ d69ce761-f489-4728-a19d-783b84b436a5
parity(σ)

# ╔═╡ a6be1358-652d-4084-9c8e-f3089fb99169
parity(τ)

# ╔═╡ 41cf3bc2-c70c-4053-bed0-74a11f8a5ad5
md"""
Now, I have all the tools to create the `n`th symmetric group as a vector of `Permutation` objects, which I will do recursively using the inductive formula:
```math
S_1 = \{ \iota \}, \qquad S_n = \{ \sigma \circ (m n) : \sigma \in S_{n-1}, m \in \{ 1, 2, \dots, n \} \}
```
where as in the `transposition` function above, ``(n n) = \iota``.

- First, I will check that `n` is positive to avoid non-terminating loops (and also because the group doesn’t make sense otherwise)
```julia
n ≥ 1 || error("symmetric group must have a positive parameter")
```

- Then, I implement the base case of the trivial group when `n == 1`:
```julia
n == 1 && return [ι]
```

- For the recursive step, I copy `symmetricgroup(n-1)` into a vector `n` times, and then set up the corresponding vector of transpositions that I will multiply each by, using array filling and concatenation. Then, I compose the two elementwise:
```julia
permutations = vcat(fill(symmetricgroup(n-1),n)...)
transpositions = vcat([fill(transposition(i,n),factorial(n-1)) for i ∈ 1:n]...)
return permutations .∘ transpositions
```

The entire function is:
"""

# ╔═╡ c70a297d-adb7-4f70-b6bb-949a3d49d3c0
function symmetricgroup(n::Int64)
    n ≥ 1 || error("symmetric group must have a positive parameter")
    n == 1 && return [ι]

    permutations = vcat(fill(symmetricgroup(n-1),n)...)
    transpositions = vcat([fill(transposition(i,n),factorial(n-1)) for i ∈ 1:n]...)
    return permutations .∘ transpositions
end

# ╔═╡ eb5640fa-0814-497c-b157-da718df1c129
md"""
Now, finding the `n`th alternating group is as simple as finding the `n`th symmetric group and taking the even permutations:
"""

# ╔═╡ d546b70a-eeaa-45db-8b65-d36cb13b67f6
function alternatinggroup(n::Int64)
    Sₙ = symmetricgroup(n)
    return Sₙ[parity.(Sₙ) .== 1]
end

# ╔═╡ 5dea6f25-65fe-4304-a7e7-2d4b830f4724
symmetricgroup(5)

# ╔═╡ 13cacc56-7fcc-4a71-9cc5-be298295939b
alternatinggroup(5)

# ╔═╡ 46eeeed1-a8e4-49cc-8d97-6ddfb1fa2b1e
md"""
### Generating arbitrary finite groups
Instead of just symmetric and alternating groups, I would like to be able to generate a group given a list of `Permutation` objects as generators. By Cayley's theorem, this suffices to construct any finite group, assuming that you can find an appropriate list of generators. This will require one final use of multiple dispatch, which will be for me to define equality of permutations:
```julia
import Base.:(==)
==(σ::Permutation, τ::Permutation) = (σ.image == τ.image)
```

This snippet of code is included in the block with the definition of the `Permutation` type at the end of the section [Declaring `Permutation` as a new type](#73d93917-b72e-468b-970c-f34ec6d51de4), since otherwise the import of `==` here conflicts with its use in the definition of `Permutation` in the Pluto notebook format.

I can now write a function to generate a group from a list of generators. This I will do by repeatedly checking closure under multiplication, which is not a particularly clever or efficient algorithm, but it will work which is all I am looking for in this case.
"""

# ╔═╡ 3e8b7d82-1d63-4aec-8ed6-0b1f11c5cb33
function generate(σs::Permutation...)
    G = [σs...]
    while true
        newelements = Permutation[]

        for σ ∈ G, τ ∈ G
            υ = σ ∘ τ
            υ ∉ G && υ ∉ newelements && push!(newelements, υ)
        end

        if isempty(newelements)
            return G
        else
            G = vcat(G,newelements)
        end
    end
end

# ╔═╡ 79584487-dd19-4014-97e7-7fcc0582efa8
md"""
An interesting point of note here is that nowhere in the function do I use the syntax `==`, and yet the function would not work without it. This demonstrates further the power of multiple dispatch, since the operator `∉` uses `==` to determine whether or not `υ` is an element of `G` (or `newelements`), and I have defined what `==` means for a `Permutation`, so `∉` will automatically work for permutations without me explicitly telling it how to.

Some demonstrations of the `generate` function are shown below, with `C₅`, `D₁₀`, and `G` the group generated by `σ` and `τ`, which happens to be `S₅`:
"""

# ╔═╡ 4725c3bb-a1c1-49cc-9bf7-a7e443db085c
C₅ = generate(Permutation(2,3,4,5,1))

# ╔═╡ 9e84f004-fd88-4a25-bfba-0797d2df0942
D₁₀ = generate(Permutation(2,3,4,5,1),Permutation(5,4,3,2,1))

# ╔═╡ bd677d6a-3b42-497d-bcfb-58543dee7ebd
G = generate(σ,τ)

# ╔═╡ 076885f1-61b4-4913-8530-fd65c531f419
S₅ = symmetricgroup(5)

# ╔═╡ 738df784-eeb7-4d3b-b19f-6317f2459c6d
all([g ∈ G for g ∈ S₅]) && all([g ∈ S₅ for g ∈ G])

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "477bf42b4d1496b454c10cce46645bb5b8a0cf2c"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.2"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─08eadc12-f8e9-11eb-1094-b7adc6936aa1
# ╟─c4d2ff29-2579-4079-a6cb-29015ef5ecd2
# ╟─db919e43-6dc8-4e53-8081-d9312446117a
# ╠═524dee4b-6687-44de-af79-b99d3724f9ac
# ╠═a16202d3-406a-4b43-807e-b60c673d4132
# ╟─41fc4769-5b67-4503-b1d1-2854a0d9c3a2
# ╠═1d56249e-d86a-4d1e-b7d2-6d63a621bfa8
# ╠═acf690f4-b65b-4576-9145-1e894d5342f2
# ╠═9ab255ff-b40a-4f16-bed8-fe14e4495502
# ╠═98750bd0-84bd-4cb1-900b-9ee26142e4ec
# ╟─9f193c51-b6ea-477f-ab0e-1bc5a4be4b86
# ╠═29b729bb-cd1a-4ee1-ae2e-26b55ce744a4
# ╠═4e4508db-47c0-462f-9cec-93e52bbffd7c
# ╟─740f8275-ee0c-4aa2-b0c1-76d46fad3ecd
# ╟─73d93917-b72e-468b-970c-f34ec6d51de4
# ╟─c18ae801-b186-44cf-9e42-83188c4781aa
# ╟─99fe42de-2afa-409f-9416-5c7671e4d0f4
# ╠═b9c535e5-3317-4291-9925-696bc6cb0252
# ╠═de16f61d-853a-43ed-adf0-fb77ea3c8f4a
# ╟─322f2dbc-f8f4-4f62-a9bb-ea6b00a333ae
# ╠═ff8ab6b4-34de-44eb-b98a-2f5e0f1c48fc
# ╠═5ce89616-326b-4404-9f2f-707c91835fa3
# ╟─9617f0f2-1b45-418f-8818-26df1abde6cb
# ╠═ee9d45f3-0a25-47f6-8445-8ba733c3a8f6
# ╟─b119e9f0-b8ba-42a3-a76c-82ed2dc170a4
# ╠═e4ce9583-5c18-4df7-a5fb-a46670fe40af
# ╟─b2ac585a-d84a-4d19-afe5-ad6d288e847f
# ╠═7a96d2c8-aa9c-405f-b34f-1f3c7c53b6c6
# ╟─12346644-107b-4c43-aaca-3e404fecc8b3
# ╠═0afa030a-58ad-44f6-a10f-0ae3baabe6f8
# ╟─1e8659dc-27d5-429d-8431-d5cfa9006411
# ╠═068a9d15-ce46-4aa5-9c13-3f79e70b0eda
# ╟─b3151f52-6dd7-4f63-a158-becb5d1aec2e
# ╠═d614739d-b22f-4991-9722-e8371732bf4b
# ╟─363481da-d185-4846-82ab-c227e9f668f8
# ╠═8f925ec6-7cd8-49fe-a467-807711b3425b
# ╟─4552e49c-af66-4fd3-a827-28284215f4f9
# ╠═0e710b07-b5db-4ea7-986b-8eb7fc0946b1
# ╟─7c082457-980c-4961-bb10-d91b7248d0f6
# ╠═59075aa4-dd35-47af-b4d1-40e89b2e599f
# ╟─764fb693-f4ac-4327-93a3-2171621ad039
# ╠═75e5491a-6536-44b0-a5b9-10cf7f9e5c66
# ╠═9d43f21c-c60a-4bf8-b1a0-bf3b56b3a97b
# ╠═b421f174-d554-47a8-a72d-4d5f3ec6ba2e
# ╠═ff82e088-aa60-49fe-bb12-fe39b0430d51
# ╠═7c0c4e99-ebff-409b-87f3-919a6870cccc
# ╠═029a0c9f-5d32-409d-b6ad-9e0d4e19e854
# ╟─c75f7100-778e-44c6-9d43-675d11175673
# ╠═0598e94b-a043-4106-bb7b-d032eb67be0b
# ╠═00a61c28-230d-41b0-81ac-e2f0869d6640
# ╠═32b52a58-6086-4c2f-a833-e1a7efb69d29
# ╟─84bd5d92-00eb-4e66-a1c5-c86626b8f7d2
# ╠═1b10a694-a310-4a4e-9831-654e3e9691d1
# ╠═9812e6fb-9308-4c27-9d57-9fbf71474343
# ╠═d69ce761-f489-4728-a19d-783b84b436a5
# ╠═a6be1358-652d-4084-9c8e-f3089fb99169
# ╟─41cf3bc2-c70c-4053-bed0-74a11f8a5ad5
# ╠═c70a297d-adb7-4f70-b6bb-949a3d49d3c0
# ╟─eb5640fa-0814-497c-b157-da718df1c129
# ╠═d546b70a-eeaa-45db-8b65-d36cb13b67f6
# ╠═5dea6f25-65fe-4304-a7e7-2d4b830f4724
# ╠═13cacc56-7fcc-4a71-9cc5-be298295939b
# ╟─46eeeed1-a8e4-49cc-8d97-6ddfb1fa2b1e
# ╠═3e8b7d82-1d63-4aec-8ed6-0b1f11c5cb33
# ╟─79584487-dd19-4014-97e7-7fcc0582efa8
# ╠═4725c3bb-a1c1-49cc-9bf7-a7e443db085c
# ╠═9e84f004-fd88-4a25-bfba-0797d2df0942
# ╠═bd677d6a-3b42-497d-bcfb-58543dee7ebd
# ╠═076885f1-61b4-4913-8530-fd65c531f419
# ╠═738df784-eeb7-4d3b-b19f-6317f2459c6d
# ╟─f1d0ed8f-bfd3-4210-84fd-f189bba6ebf9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#TODO: Compare performance to a naive implementation
#TODO: Skip over huge number of invalid boards at beginning (or better yet do a real search)

#Default board size if no argument given
#This also makes it convenient to run inside an IDE
N = 5

if length(ARGS) >= 1
  N = parse(Int64,ARGS[1])
end

function displaysquarebits(square, size)
  for i = 1:size
    for j = 1:size
      print(square & 1)
      square = square >> 1
    end
    println()
  end
end

# Create all the square sizes
function create_bitmasks(n)
  if n < 1
    error("The board must have at least one square")
  end
  if n == 1
    return Array{UInt64,1}()
  end
  if n > 8
    error("Only board sizes of up to 8 are supported due to bitfield size limitations")
  end

  bitmasks = Array{UInt64,1}()
  # Generate "flat" square bitmasks
  for square_size::UInt64=1:n-1
    offset::UInt64 = square_size * n
    for y=0:n-square_size-1
      for x=0:n-square_size-1
        initial_position::UInt64 = x + y*n
        push!(bitmasks,UInt64(2)^initial_position+
                       UInt64(2)^(initial_position+square_size)+
                       UInt64(2)^(initial_position+offset)+
                       UInt64(2)^(initial_position+offset+square_size))
      end
    end
  end

  #Generate "diamond" square bitmasks
  for square_size::UInt64=1:n÷2
    for y=0:n-(square_size*2+1)
      for x=square_size:n-square_size-1
        initial_position::UInt64 = x + y*n
        middle_offset::UInt64 = initial_position + square_size * n
        push!(bitmasks,UInt64(2)^initial_position+
                       UInt64(2)^(UInt64(2)*square_size*n+initial_position)+
                       UInt64(2)^(middle_offset+square_size)+
                       UInt64(2)^(middle_offset-square_size))
      end
    end
  end

  return bitmasks
end

function create_single_bit_boards(n)
  if n < 1
    error("The board must have at least one square")
  end
  if n > 8
    error("Only board sizes of up to 8 are supported due to bitfield size limitations")
  end

  bit_boards = Array{UInt64,1}(undef,n^2)
  for i::UInt64=1:n^2
    bit_boards[i] = UInt64(2)^(i-UInt64(1))
  end
  return bit_boards
end

function has_square(board, bitmasks, inverter)
  board_inverse = board ⊻ inverter
  for bitmask in bitmasks
    if board & bitmask == bitmask
      return true
    end
    if board_inverse & bitmask == bitmask
      return true
    end
  end
  return false
end

function find_last_one(board, start, single_bit_boards)
  index = start-1
  while index > 0
    if board & single_bit_boards[index] == single_bit_boards[index]
      return index
    end
    index -= 1
  end
  return -1
end

function iterate_valid_boards(n)
  n_squared::UInt64 = UInt64(n)^UInt64(2)
  max_board::UInt64 = UInt64(2)^n_squared-UInt64(1)
  min_half_board::UInt64 = UInt64(2)^ceil(UInt64,n_squared/2)-UInt64(1)

  bitmasks = create_bitmasks(n)
  inverter::UInt64 = max_board
  single_bit_boards = create_single_bit_boards(n)

  board::UInt64 = min_half_board
  while true
    if !has_square(board, bitmasks, inverter)
      displaysquarebits(board,N)
      println()
    end

    if board & single_bit_boards[n_squared] == single_bit_boards[n_squared]
      board_inverse = board ⊻ inverter
      scan_from = find_last_one(board_inverse, n_squared, single_bit_boards)
      num_ones = n_squared - scan_from
      for i=1:num_ones
        board = board - single_bit_boards[scan_from+i]
      end

      index = find_last_one(board, scan_from, single_bit_boards)
      if index == -1
        return
      end
      board = board - single_bit_boards[index]
      for i=1:num_ones+1
        board = board + single_bit_boards[index+i]
      end
    else
      index = find_last_one(board, n_squared, single_bit_boards)
      board = board - single_bit_boards[index]
      board = board + single_bit_boards[index+1]
    end
  end
end

println("START")
println()

iterate_valid_boards(N)

function has_square_20_million_times(board, bitmasks, inverter)
  x = false
  for i=1:20000000
    x = has_square(board, bitmasks, inverter)
  end
  return x
end

function time_has_square(board, n)
  inverter::UInt64 = UInt64(2)^(UInt64(n)^UInt64(2))-UInt64(1)
  bitmasks = create_bitmasks(n)
  @timev has_square(board, bitmasks, inverter)
  @timev has_square(board, bitmasks, inverter)
  @timev has_square(board, bitmasks, inverter)
  @timev has_square(board, bitmasks, inverter)
  @timev has_square(board, bitmasks, inverter)
  @timev has_square_20_million_times(board, bitmasks, inverter)
  #@code_native has_square(board, bitmasks, inverter)
end

#Checking one board reportedly takes around 45-90ns, but this is below the
#minimum time resolution so it's actually taking less time
#On one core of my Intel Core i7-5820K CPU @ 3.30GHz
#it can process approximately 40 million 5x5 boards per second (with no matches)
#or 12.5 million 7x7 boards per second
#This implies that it's actually taking 25ns for 5x5 and 80ns for 7x7

#time_has_square(UInt64(0b0000001011011101100110111), 5)
#time_has_square(UInt64(0b1111011010000000011011100110001110010101010110011), 7)


text:
abstract AbstractSparseArray{Tv,Ti,N} <: AbstractArray{Tv,N}

typealias AbstractSparseVector{Tv,Ti} AbstractSparseArray{Tv,Ti,1}
typealias AbstractSparseMatrix{Tv,Ti} AbstractSparseArray{Tv,Ti,2}

issparse(A::AbstractArray) = false
issparse(S::AbstractSparseArray) = true

indtype{Tv,Ti}(S::AbstractSparseArray{Tv,Ti}) = Ti

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# Testing the function D1matrix. The first derivative of a constant is zero.
@testset "D1matrix building" begin
    @test all(D1matrix(5)*ones(5).==0)
end

@testset "D2matrix building" begin
    @test all(D2matrix(4).==D1matrix(3)*D1matrix(4))
end

@testset "D3matrix building" begin
    @test all(D3matrix(5).==D2matrix(4)*D1matrix(5))
end