gaioguys / GAIO.jl

A Julia package for set oriented computations.
MIT License
9 stars 4 forks source link

optional Cartesian indices #71

Closed April-Hannah-Lena closed 1 year ago

April-Hannah-Lena commented 1 year ago

Here is a version of GAIO.jl that offers use of Cartesian indices, which can be chosen by the user. The only real changes are in the partition_regular.jl code:

struct BoxPartition{N,T,I<:Integer,A<:Union{IndexCartesian,IndexLinear}} <: AbstractBoxPartition{Box{N,T}}
    domain::Box{N,T}
    left::SVector{N,T}
    scale::SVector{N,T}
    dims::SVector{N,I}
    dimsprod::SVector{N,I}
    indextype::A
end

From the type parameters one can see that BoxPartition now uses either IndexLinear or IndexCartesian. The standard BoxPartition constructor now offers the choice via a keyword argument that defaults to linear only if the dimension of the system is less that 10:

function BoxPartition(domain::Box{N,T}, dims::NTuple{N,I}; indextype=(N < 10 ? IndexLinear() : IndexCartesian())) where {N,T,I}
    # same code as before within the function ...
    return BoxPartition{N,T,I,typeof(indextype)}(domain, left, scale, dims, dimsprod, indextype)
end

I'm not entirely sure if the extra code complexity is necessarily worth it. On the one hand, linear indexing is more memory efficient for small systems, so linear indexing optimizes storage. The whole indexing system should be entirely abstracted behind the functions point_to_key, key_to_box, so the user should never have to worry about it. However, it does lead to higher potential for bugs since there are two different functionalities working in the background simultaneously.

PS sorry I only opened this PR today, I had the code ready on Saturday but I updated my home laptop and it completely bricked my computer... I've been busy working to save all the files on my hard drive 🙃

gaioguy commented 1 year ago

Nice! But P[:] does not work anymore with IndexCartesian() ...

April-Hannah-Lena commented 1 year ago

Nice! But P[:] does not work anymore with IndexCartesian() ...

That can't be, tests/unstable_manifold.jl uses P[x] and no errors were thrown. I can also run examples/unstable_manifold.jl and examples/attractor.jl without problems. What is the error message?

gaioguy commented 1 year ago
julia> include("examples/attractor.jl")
ERROR: LoadError: iteration is deliberately unsupported for CartesianIndex. Use `I` rather than `I...`, or use `Tuple(I)...`
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] iterate(#unused#::CartesianIndex{2})
    @ Base.IteratorsMD ./multidimensional.jl:167
  [3] _totuple
    @ ./tuple.jl:328 [inlined]
  [4] Tuple{Int64, Int64}(itr::CartesianIndex{2})
    @ Base ./tuple.jl:317
  [5] (::GAIO.var"#12#13"{Int64, 2})(i::CartesianIndex{2})
    @ GAIO ./none:0
  [6] iterate
    @ ./generator.jl:47 [inlined]
  [7] union!
    @ ./abstractset.jl:96 [inlined]
  [8] Set{Union{}}(itr::Base.Generator{CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, GAIO.var"#12#13"{Int64, 2}})
    @ Base ./set.jl:10
  [9] _Set
    @ ./set.jl:30 [inlined]
 [10] Set
    @ ./set.jl:23 [inlined]
 [11] getindex(partition::BoxPartition{2, Float64, Int64, IndexCartesian}, #unused#::Colon)
    @ GAIO ~/Documents/April-Hannah-Lena/GAIO.jl/src/boxset.jl:112
 [12] top-level scope
    @ ~/Documents/April-Hannah-Lena/GAIO.jl/examples/attractor.jl:11
 [13] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [14] top-level scope
    @ REPL[6]:1
 [15] top-level scope
    @ ~/.julia/packages/CUDA/tTK8Y/src/initialization.jl:52
in expression starting at /home/junge/Documents/April-Hannah-Lena/GAIO.jl/examples/attractor.jl:11
April-Hannah-Lena commented 1 year ago

Strange. I suspect the issue lies in Base.keys(partition::BoxPartition{N,T,I,IndexCartesian}) where {N,T,I}. I made it use a more specific construction that I think will lift that error. Which version of julia are you using?

gaioguy commented 1 year ago

1.7.2 - yes, the error is gone now :+1: