9
$\begingroup$

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]

enter image description here

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
     },
    Transpose@PadLeft[raggedOnesArray, {cols, totalRows}]
    ]
   ];

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

For my use case this takes about 0.50s per array

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

How can I make this faster?

Improve this question
$\endgroup$
2
  • 5
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    May 25 at 13:58
  • $\begingroup$ 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]. $\endgroup$
    – IntroductionToProbability
    May 25 at 15:17

7 Answers 7

Reset to default
8
$\begingroup$

Plain C-style loop with Compile:

foo1 = Compile[{{startRows, _Integer, 1}, {totalRows, _Integer}},
   
   Table[
    Boole[j >= Compile`GetElement[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[1];
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.

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

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[0]);
  mint startRowsLen = libData->MTensor_getFlattenedLength(startRows);
  mint* startRowsPtr = libData->MTensor_getIntegerData(startRows);
  
  mint totalRows = MArgument_getInteger(Args[1]);
  if (totalRows < 1) return LIBRARY_DIMENSION_ERROR;
  
  MTensor x;
  mint dims[2] = {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 = 
  LibraryFunctionLoad[lib, 
   "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×

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    May 25 at 14:00
  • $\begingroup$ @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. $\endgroup$
    – Roman
    May 25 at 14:14
  • $\begingroup$ @Roman I raise you an OpenCL! (though mine has some reliability problems) $\endgroup$
    – flinty
    May 25 at 20:24
  • $\begingroup$ @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! $\endgroup$
    – Roman
    May 26 at 7:18
10
$\begingroup$

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}
Improve this answer
$\endgroup$
10
$\begingroup$

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]
Improve this answer
$\endgroup$
9
$\begingroup$ 
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×

Improve this answer
$\endgroup$
8
$\begingroup$ 
SeedRandom[1];
startRowsBig = RandomInteger[{1, 11000}, 1000];
totalRowsBig = 10000;

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

{0.127292, Null}

Improve this answer
$\endgroup$
5
$\begingroup$

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;
  }";
fun = OpenCLFunctionLoad[src, 
   "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}} *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.