# How to quickly make a rectangular array of continous 1 and 0 given the start position of 1

I would like to make a rectangular array of 1 and 0 with continuous columns of 1.

We are given the total rows of the array and the starting points of the 1.

For example:

totalRows = 5;
startRows = {1, 3, 10};

desiredArray = {{1, 0, 0}, {1, 0, 0}, {1, 1, 0}, {1, 1, 0}, {1, 1, 0}};
TableForm[desiredArray] I need to make about 10,000 of such arrays of 10,000x1000 size. How can I do this quickly?

My current implementation is

foo[startRows_List, totalRows_Integer] :=
With[
{
cols = Length[startRows],
lengthOfOnes = Min[Max[#, 0], totalRows] & /@ (totalRows - startRows + 1)
},
Module[
{
raggedOnesArray = ConstantArray[1, #] & /@ lengthOfOnes
},
]
];

foo[startRows, totalRows] == desiredArray (* True *)


For my use case this takes about 0.50s per array

SeedRandom;
startRowsBig = RandomInteger[{1, 11000}, 1000];
totalRowsBig = 10000;
foo[startRowsBig, totalRowsBig]; // RepeatedTiming (* {0.496737, Null} *)


How can I make this faster?

• Although it is fun to play with a problem like this, I have got the feeling that this is an XY-problem. Not sure what you are going to do with the matrix afterwards, but the encoding of the actual data is quite inefficient. So I guess, you would be better off by just coding whatever you want to do downstream without building the matrix in the first place. May 25 at 13:58
• I have another matrix with raw data, D, that I multiply with this 1and0s matrix, B. I then work on the clean M matrix. I guess I could encode the process of checking startRow in the creation of D so it comes out clean (as M) but I thought this way could be quicker than having some If[row>stratRow, d, 0]. May 25 at 15:17

Plain C-style loop with Compile:

foo1 = Compile[{{startRows, _Integer, 1}, {totalRows, _Integer}},

Table[
Boole[j >= CompileGetElement[startRows, i]]
, {j, 1, totalRows}, {i, 1, Length[startRows]}],

CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


A solution without Compile and with equal performance:

foo2[startRows_List, totalRows_Integer] := Module[{m, n, B, start},
m = totalRows;
n = Length[startRows];
B = ConstantArray[0, {n, m}];
Do[
start = startRows[[i]];
If[start <= m, B[[i, start ;;]] = 1;], {i, 1, n}];
Transpose[B]
];


And here the obligatory solution using a SparseArray:

foo3[startRows_List, totalRows_Integer] := Accumulate[
SparseArray[
Transpose[{startRows, Range[Length[startRows]]}] -> 1, {totalRows,
Length[startRows]}]
];


Timing comparisons:

SeedRandom;
totalRows = 10000;
startRows = RandomInteger[{1, 10000}, 1000];

A = foo[startRows, totalRows]; // RepeatedTiming // First
A1 = foo1[startRows, totalRows]; // RepeatedTiming // First
A2 = foo2[startRows, totalRows]; // RepeatedTiming // First
A3 = foo3[startRows, totalRows]; // RepeatedTiming // First
A == A1 == A2 == A3


0.484038

0.0162226

0.0161834

0.0183129

True

The SparseArray version was meant as a joke. But it performs surprisingly well.

• Impressive! I can only get another factor of 2 in speed by going to plain C. May 25 at 12:51
• @Roman Well, the problem is embarassingly parallel. So I think it is easy to squeeze out more than that by using OpenMP from C. May 25 at 13:44
• @Roman I have not been able to find a good way to parallelize this with Compile without a non-parallel Transpose or Join operation. May 25 at 13:54
• No need for OpenMP; a bunch of pthreads should be able to handle it! Maybe I’ll give it a try later. May 25 at 14:10
• Impressive performance with foo2 and with SparseArray and but foo3 doesn't handle the case where a startRow element is larger than the totalRows May 25 at 16:55

Obligatory bare-metal C answer for ultimate speed:

Needs["CCompilerDriver"]

code = "
#include \"WolframLibrary.h\"

DLLEXPORT int foo(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
MTensor startRows = MArgument_getMTensor(Args);
mint startRowsLen = libData->MTensor_getFlattenedLength(startRows);
mint* startRowsPtr = libData->MTensor_getIntegerData(startRows);

mint totalRows = MArgument_getInteger(Args);
if (totalRows < 1) return LIBRARY_DIMENSION_ERROR;

MTensor x;
mint dims = {totalRows, startRowsLen};
int err = libData->MTensor_new(MType_Integer, 2, dims, &x);
if (err) return err;
mint* xPtr = libData->MTensor_getIntegerData(x);

for (int i=0; i<totalRows; i++)
for (int j=0; j<startRowsLen; j++)
xPtr[i*startRowsLen+j] = (i >= startRowsPtr[j]-1);

MArgument_setMTensor(Res, x);
return LIBRARY_NO_ERROR;
}
";

lib =
CreateLibrary[code, "testDLL", "ShellOutputFunction" -> Print,
"ShellCommandFunction" -> Print, "CompileOptions" -> "-O3"];

Cfoo =
"foo", {{Integer, 1}, Integer}, {Integer, 2}];


Let's try it out:

Cfoo[{1, 3, 10}, 5] // MatrixForm


$$\left( \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ \end{array} \right)$$

Speed test:

Cfoo[startRowsBig, totalRowsBig]; // RepeatedTiming // First
(*    0.00708623    *)


A speedup by 68×

• Just out of curiosity: Why do you check libData->MTensor_getType(startRows) != MType_Integer and libData->MTensor_getRank(startRows) != 1? The call to the generated LibraryFunction object should already take care of that. May 25 at 14:00
• @HenrikSchumacher it's good old-fashioned paranoia. You're right that it's superfluous, I've removed these checks now to keep the code leaner. May 25 at 14:14
• @Roman I raise you an OpenCL! (though mine has some reliability problems) May 25 at 20:24
• @flinty Awesome! I've tried to multi-thread the above C code, with various block sizes, but did not get any speed improvements (on an Apple M1 Pro with 8+2 cores). Likely the whole thing is RAM-speed constrained. So OpenCL will do much better because of a better memory bus! May 26 at 7:18

Instead of constructing every column separately, note that each column is contained in a vector of n 0's and n 1's, we only need to pick the corresponding part. For ease, we create the matrix row wise and transpose in the end. For a 10^4x10^3 matrix we have:

row = 10^4;
col = 10^3;
start = RandomInteger[{1, row}, row]; start -= 1;
v = Join[ConstantArray[0, row], ConstantArray[1, row]];
org = 1 + row;


With this we may make a timing:

mat = Table[v[[(t = org - start[[i]]) ;; (t + row - 1)]], {i, col}] //
Transpose; // Timing

{0.046875, Null}


This is pretty fast:

generate[offsets_, rows_] :=
Transpose[IntegerDigits[2^(rows - Min[#, rows + 1] + 1) - 1, 2, rows] & /@ offsets]

TableForm[generate[{1, 3, 10}, 5]]

(**
1 0 0
1 0 0
1 1 0
1 1 0
1 1 0
**)


I get repeated timing of 0.126228 on my machine, compared to 0.624356 for your foo

If we compile this, then my repeated timing goes down a little more to 0.101534

generate = Compile[{{offsets, _Integer, 1}, {rows, _Integer}},
Transpose[
IntegerDigits[2^(rows - Min[#, rows + 1] + 1) - 1, 2, rows] & /@
offsets], CompilationTarget -> "C", RuntimeOptions -> "Speed"];


.. and this is a slight improvement over the original: Since IntegerDigits is Listable there's no need for the Map (/@), though I cannot seem to compile this one:

generate[rows_, offsets_] :=
Transpose@
IntegerDigits[(2^(rows - Clip[offsets, {1, rows + 1}] + 1) - 1), 2,
rows]

Rfoo[startRows_List, totalRows_Integer] :=
Transpose[
Join[ConstantArray[0, #],
ConstantArray[1, totalRows - #]] & /@
Clip[startRows - 1, {0, totalRows}]]

Rfoo[startRowsBig, totalRowsBig]; // RepeatedTiming // First
(*    0.0352359    *)


A speedup of 14×

SeedRandom;
startRowsBig = RandomInteger[{1, 11000}, 1000];
totalRowsBig = 10000;

ref = Transpose[
1] & /@ (ConstantArray[0, #] & /@ (startRowsBig - 1))
]; // RepeatedTiming


{0.127292, Null}

If you want to go really, really, absurdly fast, there's always OpenCL. This is a very naïve implementation, and the blockdim is just {1,1}

Needs["OpenCLLink"];
src = "__kernel void kern( __global mint * offsets, __global unsigned char * result, mint width, mint height) {
int row = get_global_id(0);
int col = get_global_id(1);
result[row*width + col] = row >= offsets[col]-1 ? 1u : 0u;
}";
"kern", {{_Integer}, {"UnsignedByte"}, _Integer, _Integer}, {1, 1}];
offsets = RandomInteger[{1, 1500}, 1000];
result = ConstantArray[0, {1000, 1000}];
result = Last[fun[offsets, result, 1000, 1000]];
Image[result]


... but I sometimes encounter command queue and memory limit problems with OpenCL / CUDA, so this might not be the most reliable method if you want to repeatedly call it, and it's certainly overkill.

fun = OpenCLFunctionLoad[src,
"kern", {{_Integer}, {"UnsignedByte"}, _Integer, _Integer}, {1, 1}];
offsets = {1, 3, 10};
result = ConstantArray[0, {5, 3}];
result = Last[fun[offsets, result, 3, 5]]
TableForm[result]
(* {{1, 0, 0}, {1, 0, 0}, {1, 1, 0}, {1, 1, 0}, {1, 1, 0}} *)
`