# Generate random non-zero “island” in matrix of zeros

I'm trying to write a code which will generate a 3 dimensional array of 0s, then pick some random location in the array, and subsequently expand contiguously to random neighboring sites until a set of N contiguous locations in the matrix have a value of 1 (rather than the default value of 0).

I've written a code that does this (see below), and it works rather quickly, but I wonder if there is a functional approach which would be faster/more "Mathematica"-esque. Specifically, in my implementation I had to resort to using a for loop which continually sampled from a growing list of "neighboring sites" until the desired number of 1-sites had been populated. However, as I am ultimately generating several thousands of these, I feel this may be acting as a bottleneck, and there could be a more clever way of accomplishing this goal.

Any help/guidance would be very much appreciated.

My working code:

neighbors = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
neighbors = neighbors~Join~(-neighbors);
AvailableNeighbors[curr_, size_] :=
Select[# + curr & /@
neighbors, ! AnyTrue[#, Or[# == 0, # > size] &] &]

GrowArray[percent_, size_] :=
Module[{phase, seed, ii, all, next, island, new},
seed = RandomInteger[{1, size}, 3];
island = {seed};
all = AvailableNeighbors[seed, size];
For[ii = 1, ii < percent*size^3, ii += 1,
all = DeleteCases[all, x_ /; MemberQ[island, x]];
next = RandomChoice[all];
all = DeleteCases[all, next];
island = Append[island, next];
new = AvailableNeighbors[next, size];
all = all~Join~new;
];
phase =
SparseArray[
island -> ConstantArray[1, Length[island]], {size, size, size}]
]


• Would you please show an example of using GrowArray? – Mark R Jun 11 '19 at 22:43
• Edit: I realize from looking more that you need to do something like this: GrowArray[.1, 10] – Mark R Jun 11 '19 at 22:58
• That is correct - GrowArray[0.1,10] generates a 10x10x10 matrix where 1% (so 10 total) contiguous sites are 1 rather than 0 – KHAAAAAAAAN Jun 11 '19 at 23:04
• I think 0.1=10% and when I run it, the total elements is 100. – Mark R Jun 11 '19 at 23:33
• Have you considered using Nest in place of your For loop? – Mark R Jun 12 '19 at 0:18

ClearAll[growArray]
growArray[percent_, size_] := Module[ {island = {RandomInteger[{1, size}, 3]},
length = Ceiling[percent*size^3]},
island = Take[NestWhile[
Function[x, Union[x, Clip[ RandomChoice[x] + # & /@ neighbors, {1, size}]]],
island,
Length[#] < length &], UpTo[length]];
SparseArray[island -> 1, {size, size, size}] ]


This is slightly faster than Mark R's GrowArrayFaster:

First @ RepeatedTiming[growArray[.1, 10]]


0.0014

First @ RepeatedTiming[GrowArrayFaster[.1, 10]]


0.0019

First @ RepeatedTiming[GrowArray[.1, 10]]


0.0579

To get the positions you can use the property "NonzeroPositions":

SeedRandom[1]
growArray[.1, 10]["NonzeroPositions"]


{{1, 3, 3}, {1, 3, 4}, {1, 4, 1}, {1, 4, 2}, {1, 4, 3}, {1, 4, 4}, {1,    4, 5}, {1, 5, 1}, {1, 5, 2}, {1, 5, 3}, {1, 5, 4}, {1, 5, 5}, {1,   6, 3}, {1, 6, 4}, {1, 7, 1}, {1, 7, 2}, {1, 7, 3}, {1, 7, 4}, {1, 8,    1}, {1, 8, 2}, {1, 8, 3}, {1, 9, 1}, {1, 9, 2}, {1, 9, 3}, {1, 9,   4}, {1, 10, 3}, {2, 2, 4}, {2, 3, 2}, {2, 3, 3}, {2, 3, 4}, {2, 4,   1}, {2, 4, 2}, {2, 4, 3}, {2, 4, 4}, {2, 5, 1}, {2, 5, 2}, {2, 5,   3}, {2, 5, 4}, {2, 6, 1}, {2, 6, 2}, {2, 6, 3}, {2, 7, 1}, {2, 7,   2}, {2, 7, 3}, {2, 7, 4}, {2, 8, 1}, {2, 8, 2}, {2, 8, 3}, {2, 9,   1}, {2, 9, 2}, {3, 2, 1}, {3, 3, 1}, {3, 3, 2}, {3, 3, 3}, {3, 4,   1}, {3, 4, 2}, {3, 4, 3}, {3, 4, 4}, {3, 5, 1}, {3, 5, 2}, {3, 5,   3}, {3, 6, 1}, {3, 6, 2}, {3, 7, 1}, {3, 7, 2}, {3, 7, 3}, {3, 8,   1}, {3, 8, 2}, {4, 2, 1}, {4, 2, 2}, {4, 3, 1}, {4, 3, 2}, {4, 4,   1}, {4, 4, 2}, {4, 4, 3}, {4, 4, 4}, {4, 5, 1}, {4, 5, 2}, {4, 5,   3}, {4, 5, 4}, {4, 5, 5}, {4, 6, 2}, {4, 7, 1}, {4, 7, 2}, {4, 8,   2}, {5, 3, 1}, {5, 4, 1}, {5, 4, 2}, {5, 4, 3}, {5, 4, 4}, {5, 5,   1}, {5, 5, 2}, {5, 6, 1}, {5, 6, 2}, {5, 6, 3}, {5, 7, 1}, {5, 7,   2}, {5, 8, 1}, {5, 8, 2}, {6, 4, 1}}

• Nice feature with "NonzeroPositions"! I also love what you did with Clip. I hadn't figured out how to apply it across a list and you did. By using that (instead of my created ClipAll), my version got faster ;-) – Mark R Jun 12 '19 at 17:23

I don't know if this is what you had in mind, but playing around with what you are doing, I came up with a couple of new routines.

ClipAll[nElements_List, min_, max_] :=
GrowArrayFaster[percent_, size_] :=
Block[ {seed, island, new, possibleNextPoints, next,
maxElements = Ceiling[percent*size^3]},
seed = RandomInteger[{1, size}, 3];
island = {seed};
next = seed;
While[Length[island] < maxElements,
possibleNextPoints =
Transpose@ClipAll[next + # & /@ neighbors, 1, size];
island = DeleteDuplicates[island~Join~possibleNextPoints];
next = RandomChoice[island];
];
island = Take[island, UpTo[maxElements]];
SparseArray[
island -> ConstantArray[1, maxElements], {size, size, size}]
];


ClipAll will clip any value that is out of range and it does it on a list of values. GrowArrayFaster starts from a seed and propagates out to all elements as best as it can (yes, all, not just a single one).

Using this

RepeatedTiming[GrowArray[.1, 10]]


And

RepeatedTiming[GrowArrayFaster[.1, 10]]


The difference is a factor of 36.

And to find the indices, right after running one of these, use this:

Position[Normal@%,1]


Inspired by what @kglr wrote, I eliminated ClipAll and just used Clip correctly for the list and now the timing is around his version. Faster sometimes and slower.

GrowArrayFaster2[percent_, size_] :=
Block[{seed, island, new, possibleNextPoints, next,
maxElements = Ceiling[percent*size^3]},
seed = RandomInteger[{1, size}, 3];
island = {seed};
next = seed;
While[Length[island] < maxElements,
possibleNextPoints = Clip[next + # & /@ neighbors, {1, size}];
island = DeleteDuplicates[island~Join~possibleNextPoints];
next = RandomChoice[island];
];
island = Take[island, UpTo[maxElements]];
SparseArray[
island -> ConstantArray[1, maxElements], {size, size, size}]
];

First@RepeatedTiming[GrowArrayFaster2[0.1,10]]
(* 0.00078 *)


Here's a partial solution, one which 'grows' a 3D island of 1s in a volume of 0s. I don't have time right now to refine it. The intention is that you would use this method to grow an island within a sub-volume of your ocean of 0s, a sub-volume which is just large enough to contain the island you want, and at a subsequent step place the island at a (possibly) random location (and orientation) within your ocean.

First, a function to generate a 2D island. This creates an island of m^2 1s in an area of n^2.

make2DIsland[m_, n_] :=
Partition[RandomSample[Join[Table[1, m^2], Table[0, n^2 - m^2]]], n]


By keeping n not much bigger than m (I tested with n==m+1 mostly) you increase the chances of creating a single island at this stage, and one which is kind-of compact. Which may or may not be desirable; one modification you might make is to use two arguments for the size of the sub-volume (ie factors of n) to change the overall shape of islands.

Now, simply generate p such islands and stack them up ...

make3DIsland[m_, n_, p_] := Table[make2DIsland[m, n], p]


Whoaa, you cry, what's the guarantee that this produces a single island !? None whatsoever :-), so let's check that this is a single island ... first create an island

i3 = make3DIsland[5, 7, 3];


then check how many components it has

Max[MorphologicalComponents[i3]]


and throw it away if this produces 2 or greater.

I don't know if:

• this meets your criteria for island-icity; one change I can think you might make is to use the CornerNeighbors -> False option when checking the morphological components; and you might want to apply the morphology test to 2D islands before stacking them up;

• this is faster than your existing method; as I said I'm in a bit of hurry, and I haven't checked. I wouldn't be surprised to learn that this approach, taking into account the number of rejections you might have to make, is slower.

Obviously this first draft only creates islands with p*m^2 1s, it shouldn't be too difficult to modify to work with any 3 factors of the size of the island required, but might get tricky if you want islands with prime size.

But it's a lot less code.

Not saying anything about speed, but below is a graph-based approach, finding random vertices adjacent to those already selected. This performs essentially randomized "flood fill" over a graph (it wouldn't really need to be a GridGraph) until the number of locations is satisfied.

ClearAll@GrowArray;

GrowArray[ratio_, size_] :=
With[{g = GridGraph[{size, size, size}]},
{RandomChoice@VertexList@g}, Ceiling[Min[ratio, 1] size^3 - 1]] //
SparseArray[(1 + IntegerDigits[# - 1, size, 3]) -> 1, {size, size, size}] &]


GridGraph numbers them sequentially in a specific order starting from 1 - the IntegerDigits construct converts these to coordinates. Replacing GridGraph with own construct which would name the vertices differently could... well, move unnecessary hurdles from one point of this code to another.

Adjacency lists could be accumulated for extra speed, but I would assume a generation of adjacency list of a list of vertices is a relatively fast operation, so I didn't start optimizing a relatively clear piece of code. Another optimization would be to incrementally remove those vertices which have no free adjacent vertices left...

• I tried hard to figure out a graph based approach! I'm not ending up using this method, but think its quite cool, and may actually be ideal for some other applications I'm thinking of. – KHAAAAAAAAN Jun 13 '19 at 7:19
• @KHAAAAAAAAN It seems to have a pretty bad time complexity when you grow the amount of locations, indeed. It does demonstrate how otherwise somewhat convoluted tasks can be turned into graph problems, though. – kirma Jun 13 '19 at 7:22

This implementation utilizes the adjacency lists of the adjacency matrix of the underlying graph. I use Associations to store index list because I expected the modification operations (deleting elements and appending new ones) to be faster than with Lists. However, I could not make out a timing difference. =/

Overall, this seems to be way slower as kglr's solution.

ClearAll[GrowArray2]
GrowArray2[percent_, size_, nseeds_] := Module[{A, island, shore},
island = AssociationThread[RandomChoice[1 ;; size^3, nseeds] -> 1];
shore = AssociationThread[Complement[Join @@ A[[Keys[island]]], Keys[island]] -> 1];
Do[
With[{i = RandomChoice[Keys[shore]]},
AssociateTo[island, i -> 1];
KeyDropFrom[shore, i];
];
,
{Floor[percent size^3]}];
ArrayReshape[SparseArray[Partition[Keys[island], 1] -> 1, size^3], {size, size size}]
]


Example:

a = GrowArray2[0.1, 30, 3]; // AbsoluteTiming // First
Image3D[a]


0.538474

• Points at the very least for introducing me to Image3D which works splendidly for visualization for this purpose. – KHAAAAAAAAN Jun 15 '19 at 7:35