How to speed up this SparseArray construction

Question

This code is (obviously) getting slower with the growing number of trials, currently running over 5 hours. This stylised example takes 2.5 minutes. It is not entirel illustrative as the real situation is much sparser and less spread out but I would not know how to condense this into the below example.

How can I speed it up? The first section defines the input, a set of two dimensional slabs in a 4 dimensional space (5 racks and 20 shelves are allocated via a random discrete number, where slabs are defined by horStart, horEnd, verStart and verEnd - also randomly - which define the coordinates of the boundaries via 2 dimensional array. The idea is to stack the individual cells of these arrays and measure how high the stack becomes.

rack = RandomInteger[5, trials];
shelf = RandomInteger[20, trials];
horStart = RandomInteger[100, trials];
verStart = RandomInteger[200, trials];
horEnd = 
 RandomChoice[{.8, .15, .04, .01} -> {0, 10, 100, 200}, trials] + 
  horStart;
verEnd = 
 RandomChoice[{.8, .15, .04, .01} -> {0, 10, 100, 200}, trials] + 
  verStart;

The whole set is subsequently assigned to input, which is a List of Integer with length 6, identifying the rack and shelf space, and how much what individual cells are taken up by the slabs.

input = {rack, shelf, horStart, verStart, horEnd, verEnd}\[Transpose];

The sparseArray needs to be initialised with a value for the dimensions.

reach = Max /@ {rack, shelf, horEnd, verEnd};
sa2 = SparseArray[{}, reach, 0];

The second section allocates this to a SparseArray where the number of overlaps are counted to identify how high the stacks will become.

The area of interest is:

sa2=SparseArray[{{ra_, sh_, hor_, ver_} /; 
     Apply[Or, 
      Or[And[ra == #1, sh == #2, Between[hor, {#3, #5}], 
          Between[ver, {#4, #6}]] & @@@ input]] :> (sa2[[ra, sh, 
       hor, ver]] + 1)}, reach]]

which is at the bottom of the below integrated code.

AbsoluteTiming[
sa = Module[{rack, shelf, horStart, verStart, horEnd, verEnd, 
     trials = 1000, input, reach, sa2},

(*first section: define rack, shelf and slab sizes to stack*)

rack = RandomInteger[5, trials];
shelf = RandomInteger[20, trials];
horStart = RandomInteger[100, trials];
verStart = RandomInteger[200, trials];
horEnd = 
 RandomChoice[{.8, .15, .04, .01} -> {0, 10, 100, 200}, trials] + 
  horStart;
verEnd = 
 RandomChoice[{.8, .15, .04, .01} -> {0, 10, 100, 200}, trials] + 
  verStart;
reach = Max /@ {rack, shelf, horEnd, verEnd};
input = {rack, shelf, horStart, verStart, horEnd, verEnd}\[Transpose];

(* second section, allocate to SparseArray, measure size of stacks *)

sa2 = SparseArray[{}, reach, 0];
sa2 = SparseArray[{{ra_, sh_, hor_, ver_} /; 
     Apply[Or, 
      Or[And[ra == #1, sh == #2, Between[hor, {#3, #5}], 
          Between[ver, {#4, #6}]] & @@@ input]] :> (sa2[[ra, sh, 
       hor, ver]] + 1)}, reach]];]


Out (* 140 *)

Answer 1

The code is super inefficient as it uses pattern matching and as to perform many more checks than there are nonzero entries in the final matrix. It is usually better to compute the sparse array pattern, i.e., the list of all nonzero positions. This is done by the following compiled function:

cf = Compile[{{r, _Integer}, {s, _Integer}, {hS, _Integer}, {vS, _Integer}, {hE, _Integer}, {vE, _Integer}},
   Flatten[Table[{r, s, i, j}, {i, hS, hE}, {j, vS, vE}], 1],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Now the matrix is assembled with

A = SparseArray[
   Join @@ cf[rack, shelf, horStart, verStart, horEnd, verEnd] -> 1, 
   reach];

For trials = 1000, this takes less than a 10000th part of the time.

Edit: Additive assembly can be done with the following helper function:

MySparseArray[R_Rule, dims_?VectorQ, fun_: Total, background_: 0] :=
 With[{spopt = SystemOptions["SparseArrayOptions"]},
  Internal`WithLocalSettings[
   SetSystemOptions[
    "SparseArrayOptions" -> {"TreatRepeatedEntries" -> fun}],
   SparseArray[R, dims, background],
   SetSystemOptions[spopt]]
  ]

and then

A = MySparseArray[
   Join @@ cf[rack, shelf, horStart, verStart, horEnd, verEnd] -> 1, 
   reach, Total, 0];