How to store this structure (list of lists of integers) in Matlab?

Question

I need to store a list of lists of integers. For example, X[1] should be able to contain [1 3 5] while X[2] could contain [1 2]. What's the best solution? A cell array?

---

The back story:

For a project, I pre-calculate the intersections between N lines and M cubes. These are retrieved in two ways: given a line index, I want a list of cubes that it passes through, and given a cube index, I want a list of lines that pass through it.

Typical values are N=2^24 and M=2^18, which means an intersection matrix (NxM) is out of the question. Fortunately, the average line passes through only M^(1/3)=2^6 cubes. Currently, I'm storing the structure as an NxM^(1/3) matrix, so that X(n,:) is a vector of cubes that the nth line passes through (padded with zeros).

This works fine for retrieving cubes given a list index, but it turns out that the bottleneck of my code is the retrieval of lines given a cube index. (I do it with find(X==m) where m is the cube index.) I can't create the opposite matrix, as the number of lines passing through a single cube can be very high, even though it is low on average.

Answer 1

Generally a cell array is the right answer for this. This is the simplest case. Some example uses:

%Writes
X = {[1], [1 2 3], [1 2]};
X{4} = [1 2 3 4];

%Reads
a = X{1}
b = cat(2,X{:});
c = X([2 4]);

However, it is not the only answer.

You could use an array of structures, each with a field called .indexes (or an appropriate name based on your problem). This allows a little more flexibility if there is additional information you would like attached to your list of lists, for example cube position could be added as a .position field. Example uses:

%Writes
X(1).indexes = 1;
X(2).indexes = [1 2 3];
X(3).indexes = [1 2];

%Reads
a = X(1).indexes
b = cat(2,X.indexes)
c = X([2 4]);

You could also use a containers.Map object. This has the same advantages as an array of structures, but with greater flexibility in how you reference your objects. Whereas when using an array of structures is structure is references by an index, using a containers.Map object can reference each structure using an arbitrary number (not integers near 1), or a name (not practical for 2^24 cases). This is probably not the best answer for you, but for reference examples uses are below:

%Writes
X = containers.Map('keyType','uint32','valueType','Any');
X(1) = [1];
X(2) = [1 2 3];
X(3) = [1 2];
X(4) = [1 2 3 4];

%Reads
a = X(1);
b = cat(2,X.values);

Finally, it may be worth defining a pair of custom classes for this. This is a bit more work to set up, but is probably the easiest way to get constant time lookups into your precalculated values. Some code to get you started down this path is below.

%A start at cube.m.  Most of the code handles smartly reallocating the list of lines.
classdef cube < handle
    properties (SetAccess = private, GetAccess = public)
        numLines = 0
        intersectingLines = [];
    end
    methods (Access = public)
        function addLine(self, lineToAdd)
            if self.numLines == 0
                self.intersectingLines = lineToAdd;
                self.numLines = 1;
            elseif self.numLines>=length(self.intersectingLines)
                self.intersectingLines(length(self.intersectingLines)*2) = line();
                self.intersectingLines(self.numLines+1) = lineToAdd;
                self.numLines = self.numLines+1;
            else
                self.intersectingLines(self.numLines+1) = lineToAdd;
                self.numLines = self.numLines+1;
            end
        end
    end
end

%A start at line.m.  A near copy/paste of cube.m
    classdef line < handle
    properties (SetAccess = private, GetAccess = public)
        numCubes = 0
        intersectingCubes = [];
    end
    methods (Access = public)
        function addCube(self, cubeToAdd)
            if self.numCubes == 0
                self.intersectingCubes = cubeToAdd;
                self.numCubes = 1;
            elseif self.numCubes>=length(self.intersectingCubes)
                self.intersectingCubes(length(self.intersectingCubes)*2) = cube();
                self.intersectingCubes(self.numCubes+1) = cubeToAdd;
                self.numCubes = self.numCubes+1;
            else
                self.intersectingCubes(self.numCubes+1) = cubeToAdd;
                self.numCubes = self.numCubes+1;
            end
        end
    end 
end

To use these classes as written, you need to call the add methods in pairs (an obvious upgrade for later is to properly cross add. In the meantime (because I'm lazy), we'll define a helper function.

function crossAdd(cube, line)
cube.addLine(line);
line.addCube(cube);

Now example use is:

%Create two class arrays of cubes and lines
allCubes(1) = cube;
allCubes(2) = cube;
allCubes(3) = cube;
allCubes(4) = cube;

allLines(1) = line;
allLines(2) = line;
allLines(3) = line;
allLines(4) = line;

%Define links (matching above "writes" examples)
crossAdd(allCubes(1), allLines(1));
crossAdd(allCubes(2), allLines(1));
crossAdd(allCubes(2), allLines(2));
crossAdd(allCubes(2), allLines(3));
crossAdd(allCubes(3), allLines(1));
crossAdd(allCubes(3), allLines(2));
crossAdd(allCubes(4), allLines(1));
crossAdd(allCubes(4), allLines(2));
crossAdd(allCubes(4), allLines(3));
crossAdd(allCubes(4), allLines(4));

%Use linked values
aLines = allCubes(1).getLines   %Only one intersecting line
bLines = allCubes(2).getLines   %Three intersecting lines
cubesFromSecondLine = bLines(2).getCubes %Three cubes here (2, 3, 4)

BTW, we're really just taking advanatge of the fact that < handle classes behave as pass-by-reference, so we can use complex, cross-linked data structures.