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

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.

-
Can you provide a minimal working example and give your desired output? :) One example and one output are better than 7 lines of text :) –  Öskå May 23 '14 at 16:01
In general pattern matching is going to be slower than more explicit code. A method that scans row by row will likely be the fastest approach. Partition would do most of the work for you. –  mfvonh May 23 '14 at 16:28
Would it work for you to compute a correlation using the highly optimized ListCorrelate/L:istConvolve? Looking at 1D (i.e. a row) now, take a vector (row), and compute the correlation/convolution with {1, 1, ..., 1}. The positions where the correlation is high give you high densities of 1s. Then take this correlation vector, apply UnitStep to it with a threshold, so you get 1 at high densities and 0 elsewhere. FInally subtract it form the original vector to eliminate high densities of 1s and then apply UnitStep again to get rid of potential -1s. –  Szabolcs May 23 '14 at 17:20
kernel=ConstantArray[1., 20]; result = UnitStep[# - 0.5 - N@Unitize@ListCorrelate[kernel, UnitStep[ListCorrelate[kernel, #, {1, 1}, 0.] - 15.], {1, 1}, 0.]] & /@ data; –  Szabolcs May 23 '14 at 17:52
Ah, this is not as interesting as I'd expected (I thought from your prior description you wanted to treat interspersed zeroes as "noise", and ones separated by some number of "noise" or less were to be considered part of the same sequence.) That's an interesting problem. This one submits to Szabolcs correlation suggestion, or look at array filtering (MinFilter/MaxFilter,etc.). You can probably get sneaky, convert array to image, and use image processing for ludicrous speed, then convert back to numeric... –  ciao May 23 '14 at 22:20

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).

-
It does appear to accomplish what I originally intended. What if the "chunk" were to move one increment at a time versus an entire "chunk" size at a time? That way, any given "chunk" in the data would never exceed "lim". –  Kane May 24 '14 at 16:08
@Kane: Take a look at update, give it a whirl... –  ciao May 24 '14 at 22:13
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 *)

-
I knew that using Position was not smart enough.. :) +1 –  Öskå May 23 '14 at 22:32
@Oska, Position-based method was 10x faster than the first few methods I tried:) –  kglr May 23 '14 at 22:45
Only tried latest version, but it blows up if chunksize is not a divisor of row length. –  ciao May 23 '14 at 23:41
@rasher, thanks... I added the divisibility condition in the definitions. –  kglr May 24 '14 at 0:50
@OrdersOfMagnitudeFasterThanMeClub, I'm quite curious how you smarties would implement a moving window (see my edit & comments) –  mfvonh May 24 '14 at 15:41

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.

-

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];
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 :)

-
No, not quite. I think it is my fault for not being clear enough. Let me revise my posting and give a small scale example. –  Kane May 23 '14 at 19:35
Do you want to total the occurrences over fixed positions or with a moving window? That is, supposing a window of 3 instead of 20, should the row {0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0} become {0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0} or {0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0}? –  mfvonh May 23 '14 at 21:26
I like the moving window idea. However, I think you have what I am trying to do backwards. If the density is too high (meaning that if there are too many 1s close together in a given region), I want to make everything zero. I want less ones, not more 1s. –  Kane May 23 '14 at 21:38
Ah yes I did have that backwards. So then my example should map to a row of only 0's? –  mfvonh May 23 '14 at 21:43
If the window is 3 and my max density say 2 or more and I am moving the window 3 indices at a time then the result would be {0,0,0,0,0,0,0,0,0,1,0,0}. –  Kane May 23 '14 at 22:03