Replacing High Densities of 1.s in Very Large Binary Arrays

Question

I am working with very large binary arrays on the order of {100,100000}. It may look something like this.

data = N[Table[RandomInteger[], {i, 100}, {j, 100000}]];

What I am then trying to do with this data is replace high densities of 1.'s with 0.s in each row. For instance, if I look at successive chunks of 20 in each row and find that it has 15 1.s I want to then make them all 0.s. Now, I know that I can index through each row summing up every 20 indices and determine if the sum is 15 or more. If so, then I can make them all 0. However, indexing through large arrays in this fashion can be very time consuming. Here is a small scale example of what I want to do. In this smaller scale of data I am zeroing out densities of 3 or more for every 5 indices and leaving everything else the same.

data = N[Table[RandomInteger[], {i, 2}, {j, 20}]]

{{1., 0., 0., 1., 0., 1., 0., 1., 0., 1., 1., 0., 0., 0., 1., 1., 0., 
  0., 0., 1.}, {0., 1., 0., 0., 1., 1., 0., 1., 0., 1., 0., 0., 1., 
  0., 1., 0., 0., 0., 1., 1.}}

Do[
z = Total[data[[i]][[j ;; j + 4]]];
If[z >= 3, 
data = ReplacePart[data, Table[{i, j + k}, {k, 0, 4}] -> N[0]]], {i,
1, Dimensions[data][[1]]}, {j, 1, Dimensions[data][[2]], 5}]
data

{{1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 1., 1., 0., 
  0., 0., 1.}, {0., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 1., 
  0., 1., 0., 0., 0., 1., 1.}}  

I need to find a way to do this as fast as possible with really big arrays. This way takes a long time for a {20,100000} array looking at every 20 indices zeroing out densities of 15 or more. As you can see...

data = N[Table[RandomInteger[], {i, 20}, {j, 100000}]];
Do[
z = Total[data[[i]][[j ;; j + 19]]];
If[z >= 15, data = ReplacePart[data, Table[{i, j + k}, {k, 0, 19}] -> N[0]]], {i, 1, 
Dimensions[data][[1]]}, {j, 1, Dimensions[data][[2]], 20}] // AbsoluteTiming

{18.767073, Null}

A {100,100000} array would obviously take longer. I've tried some pattern ideas, see the following Replacing Patterns In a List of Varying Length, but that didn't quite do it since the pattern of 1.s had to be exact before replacing.

Answer 1

I'll throw this out, cobbled up on loungebook while cigar smoking. Same caveats as in my comments re: any of the partitioning-based solutions presented so far. It produces same results as your code and the fastest solution posted so far (trnsfrmF3 from Kguler), it is ~8X faster than that on a 100 X 100K test array using the 5 x 3 criteria in the OP, and orders of magnitude faster than your example code:

With[{d = Dimensions@#, w = ArrayReshape[#, {Length@#*Length@#[[1]]/#2, #2}]},
   ArrayReshape[UnitStep[#3 - 1 - Total /@ w] w, d]] &[test, chunk, lim]

Where test,chunk, and lim arguments are the data array, "chunk" size, and cut-off total limit, respectively.

I still wonder as per my comments if this is really what you intended, what with missing boundary cases and chunk size vs row size limitations.

As I noted, just barfed this out while lounging, not heavily tested, no sanity check for allowable chunk size, and the usual caveats re: speed of loungebook.

I'll await clarification on intent before attempting a more sophisticated assault on the problem.

Update: Here's a "moving window" type of solution, which I think gives a better realization of the desired result (though I still think "density" here is a bit ambiguous, and a solution that tracks zero-one spacing would be more "correct", e.g., a row of {0,1,0,1,0,1,0,1,0,1,...} is obviously "dense" in ones, but a sliding-window solution would not treat these ones as "in a row" and you'd have to fine-tune the chunk/limit to catch them):

result = Block[{data = #, a, psize = #2, lim = #3},
    a = Array[Range[#, # + psize - 1] &, {Length@data[[1]] - psize + 1}];
    Block[{tmp = #, m},
       tmp[[DeleteDuplicates@
               Flatten[a[[Pick[Range@Length@(m = MovingAverage[tmp, psize]), 
                          UnitStep[m - lim/psize], 1]]]]]] = 0.;
       tmp] & /@ data] &[test, chunk, lim]

This takes advantage of MovingAverage (performs quite well for machine-precision lists, as is the case here). About twenty times faster on a 100 X 50K array using chunk of 7 and limit of 3 vs mfvonh's interesting solution, bigger problem size makes for even larger delta (did not test beyond that size: mfvonh's solution gets memory-hungry, and I'm on the loungebook). Did a 100 X 100K array in a few seconds on the loungebook, so should be well under a second for same on a "real" machine.

Note, the two produce slightly different results because of the technique used. Mine looks at the sliding window, and any windows exceeding the limit are fully zeroed.

Give it a whirl, let me know what your thoughts are.

Update: Here's another method that might find use here, since it appears the primary goal is to leave "isolated" ones alone while zeroing other ones in "dense" areas.

result1 = Block[{data = #, psize = #2, lim = #3, ma1},
    (ma1 = MovingAverage[ArrayPad[#, psize - 1,0.], psize];
       UnitStep[Subtract[lim, ma1[[;; -psize]] + ma1[[psize ;;]]]]*#) & /@ 
                data] &[test, chunk,threshold]

As above, test is your array. chunk is the size of the "window". threshold is a cut-off threshold between 0 and 2 (2 meaning no filtering - fiddle with values to see how it affects things). Note that by the forward/backward use of the window, the actual window is 2*chunk-1 in size, e.g., a chunk of 3 is a window of 5 (3 forward, 3 back, with one shared (the element).

This works by using the very fast (at least for machine-precision) MovingAverage in a forward and backward scan on appropriately padded rows, summing the results. The resulting values reflect the "density" of ones to both sides of an element.

Quite fast (on the loungebook it's a few hundredths of a second per row with 100K length rows). Note I pad with zeroes on the ends, you might want some other padding (constant, reflection, etc.) - this affects how "end of row" elements behave (since the first element, e.g., has no "left" neighbors, it can't be treated quite the same, only you can decide what's appropriate here).

Answer 2

ClearAll[trnsfrmF, trnsfrmF2];
trnsfrmF[data_, chunksize_, threshold_] := 
    Module[{parts = Partition[Flatten@data, chunksize], dim = Dimensions[data][[2]]}, 
     parts = # (1 - UnitStep[Total[#] - threshold]) & /@ parts;
     Partition[Flatten@parts, dim]];
trnsfrmF2[data_, chunksize_, threshold_]:= 
 Module[{parts = Partition[Flatten@data, chunksize], dim = Dimensions[data][[2]]}, 
  Partition[Flatten[parts  ((1 - UnitStep[Total[#] - threshold]) & /@ parts)], dim]];
(* much faster than the previous ones is *)
trnsfrmF3[data_, chunksize_, threshold_]  := 
   Module[{parts = Partition[Flatten@data, chunksize], dim = Dimensions[data][[2]]}, 
   Partition[Flatten[parts ((1 - UnitStep[Total[parts, {2}] - threshold]))], dim]];

Note: One could add the condition that chunksize divides the column dimension of data using trnsfrmF[data_, chunksize_, threshold_] /; Divisible[Length[data], chunksize] := ....

Small example data

data={{0., 1., 0., 1., 0., 1., 0., 1., 0., 1., 0., 0., 1., 0., 0., 1., 1., 0., 1., 0.}, 
      {1., 1., 1., 0., 1., 0., 1., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 1., 0.,0.}};
trnsfrmF[data, 5, 3]
(* {{0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,  0., 0., 0.},
   {0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}} *)
trnsfrmF[data, 5, 3]== trnsfrmF2[data, 5, 3] ==  trnsfrmF3[data, 5, 3]
(* True *)

Large dataset

randdata = N[Table[RandomInteger[], {100}, {100000}]]

res1 = trnsfrmF[randdata, 20, 15]; // Timing
(* {0.234375, Null} *)
res2 = trnsfrmF2[randdata, 20, 15]; // Timing
(* {0.156250, Null} *)
res3 = trnsfrmF3[randdata, 20, 15]; // Timing
(* {0.078125, Null}  *)

chunks = 20; max = 15;
(parts = Partition[#, chunks] & /@ randdata;
 (parts[[First@#, Last@#]] = parts[[First@#, Last@#]]*0) & /@ 
  Position[Map[Total, parts, {2}] /. {x_ /; x >= max -> 0}, 0];
  res0 = Partition[Flatten@parts, 100000];) // Timing  (* Oska's method)
(* {0.531250, Null} *)

res0 == res1 == res2 == res3
(*  True *)

Answer 3

Trying with OP's data:

Small example:

dataOP = {{1., 0., 0., 1., 0., 1., 0., 1., 0., 1., 1., 0., 0., 0., 1., 1., 0., 0., 0., 1.}, 
  {0., 1., 0., 0., 1., 1., 0., 1., 0., 1., 0., 0., 1., 0., 1., 0., 0., 0., 1., 1.}};
chunks = 5; max = 3;
parts = Partition[#, chunks] & /@ dataOP;
(parts[[First@#, Last@#]] = parts[[First@#, Last@#]]*0) & /@ 
  Position[Map[Total, parts, {2}] /. {x_ /; x >= max -> 0}, 0];
newdata = Partition[Flatten@parts, 20];

outputOP = {{1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 1., 1., 0., 0., 0., 1.}, 
  {0., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 0., 1., 1.}};

outputOP == newdata

True

Bigger one:

SeedRandom@0;
data = N[Table[RandomInteger[], {i, 20}, {j, 100000}]];
Do[z = Total[data[[i]][[j ;; j + 19]]];
  If[z >= 15, 
   data = ReplacePart[data, 
     Table[{i, j + k}, {k, 0, 19}] -> N[0]]], {i, 1, 
   Dimensions[data][[1]]}, {j, 1, Dimensions[data][[2]], 20}] // AbsoluteTiming

{31.608005, Null}

SeedRandom@0;
data1 = N[Table[RandomInteger[], {i, 20}, {j, 100000}]];
chunks = 20; max = 15;
(parts = Partition[#, chunks] & /@ data1;
  (parts[[First@#, Last@#]] = parts[[First@#, Last@#]]*0) & /@ 
   Position[Map[Total, parts, {2}] /. {x_ /; x >= max -> 0}, 0];
newdata = Partition[Flatten@parts, 100000];) // AbsoluteTiming

{0.518846, Null}

newdata == data

True

Then, a 100 x 100000 array takes 1.409655 seconds on my machine.

Answer 4

It sounds like this may be what you want. 5 seconds on my crappy laptop.

(
  data = N[Table[RandomInteger[], {i, 100}, {j, 100000}]];
  parts = Partition[#, n = 20] & /@ data;
  scan = Map[Sequence @@ ConstantArray[If[Mean@# > 15/20, 1, 0], n] &, parts, {2}];
  Tally[Flatten@scan]
  ) // AbsoluteTiming

{4.879750, {{0, 9939360}, {1, 60640}}}

EDIT

Here's a different approach that uses a moving window per the comments.

(
  data = N[Table[RandomInteger[], {i, 100}, {j, 100000}]];
  bandwidth = 21;(*should be odd*)
  threshold = 15;
  kernel = ConstantArray[1, bandwidth];
  corr = ListCorrelate[kernel, #] & /@ data;
  tailwidth = Floor[bandwidth/2];
  padded = PadRight[PadLeft[#, Length@# + tailwidth, #[[1]]], Length@# + 2 tailwidth, #[[-1]]] & /@ corr;
  transformed = ReplacePart[data, Position[Map[UnitStep[# - threshold]&, padded, {2}], 1] -> 0];
  ) // AbsoluteTiming

{2.868025, Null}

Count[data, 1., \[Infinity]]

5001160

Count[transformed, 1., \[Infinity]]

4710877

I'm sure the gurus will show me up in terms of performance in short order :)