Generate random non-zero "island" in matrix of zeros
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_] :=
Thread[Clip[nElements, {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 *)
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}}
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.