# Looping through 3 million elements of a list to 'Set' values is slow

I am looping through the 3 million elements of a time-series (v1, o1) and corresponding values (v2, o2).

• v1 contains a uniformly increasing timestamp (e.g. 1 to 3,000,000 in increments of 1)
• o1 contains a non-uniformly increasing timestamp (includes gaps, meaning o1 is a shorter list)
• v2 contains zeroes
• o2 contains real-type values

The purpose is to assign the o2 values to v2 for each v1 timestamp. In case of a gap, a previous o2 value will be used.

Could this code be optimized?

Sample input

 v1 = {1, 2, 3, 4, 5, 6}
o1 = {2, 3, 5}
v2 = {0, 0, 0, 0, 0, 0}
o2 = {0.3, 0.5, 0.1}


Sample output

 v2 = {0, 0.3, 0.5, 0.5, 0.1, 0.1}


Code

v1 = Range[1, 3000000, 1];
v2 = ConstantArray[0., 3000000];
o1 = Delete[v1, {#} & /@ RandomInteger[3000000, 500000]];
o2 = RandomReal[{0, 1}, Length[o1]];

Length[o1]
Length[v1]

y = 1;
i = 1;                              (* Revision: Set i to 1 *)
Monitor[
While[i<Length[v1],             (* Revision: substitute Do for While *)
If[SameQ[v1[[i]], o1[[y]]],
v2[[i]] = o2[[y]];
i++;                     (* Revision: i++ *)
y++,
If[v1[[i]] < o1[[y]],
v2[[i]] = o2[[y - 1]];
i++,                   (* Revision: i++ *)
y++]]
], {i, y}] // Timing


Big thank-you to all the participants emphasizing the importance of clear problem formulation, as well as providing with inputs, and expected outputs.

• Can you please provide some example data or a procedure to generate simulated data. – s0rce Sep 2 '15 at 2:23
• Don't do it that way - mutating a step at a time is... slow. Make a zero list of correct length. Populate the indices at o1 with o1. Fold Max over that. Use result to pull values. All are one-step operations, should handle millions of entries in a fraction of a second... – ciao Sep 2 '15 at 4:26
• @ciao Are you assuming o2 is increasing? – Dr. belisarius Sep 2 '15 at 5:26
• @belisarius: you mean o1? Then yes - it's time steps, so unless we're back to the future... if you meant o2 - not sure what you mean, it's just data values from my read. – ciao Sep 2 '15 at 5:28
• @ciao I've read o2. Silly me – Dr. belisarius Sep 2 '15 at 5:32

Now that OP has accepted an answer, implicitly making the results from it correct, here's a much faster method than the fastest there (as in 3x to nearly order of magnitude faster, depending on amount of missing points):

v2 = c2 = ConstantArray[0, Length@v1];
v2[[o1]] = o2;
c2[[o1]] = o1;
c2[[Range[1, First@o1]]] = Range[1, First@o1];
v2 = v2[[With[{xx = c2}, Compile[{}, FoldList[Max, xx]]][]]];


Had this flash when nodding off, should be much faster still:

c2 = ConstantArray[0, Length@v1];
c2[[o1]] = 1;
If[c2[] == 0, c2[] = 1; o2 = Prepend[o2, 0]];
c2 = Accumulate@c2;
v2 = o2[[c2]];


Take a look at SquareOne's neat post of various methods and timings...

• Faster indeed ! ;) Very instructive answer (for me at least) ! (As requested, see timings in my post). – SquareOne Sep 5 '15 at 17:00
• @SquareOne: Glad you like it, give update a whirl - simpler, and much faster again in my quick tests... – ciao Sep 5 '15 at 22:31
• @Matvei Kruglyak :appreciate the accept, glad you find answer useful. I'd humbly ask however that you return accept to SquareOne - their answer is most educational, shows effort, and they were kind enough to mention mine so if a reader has performance as top priority, they'll see mine referenced. I mark my answers as wiki so I get no points, I'd rather a deserving new poster like SquareOne get them... – ciao Sep 6 '15 at 3:23
• Added the new faster solution to my timings. Nice ! Also I had to correct your previous timings because I forgot to include the beginning part of your code, sorry ;( However, it is only a little bit slower. – SquareOne Sep 6 '15 at 10:14
• I did actually ask OP to accept your answer as performance is clearly the goal. You can mention my post for the detailed timings and comparison with other approaches. – SquareOne Sep 6 '15 at 10:18

Here are 3 other approaches.

I'll first test their output with very basic input:

n = 10;
o1 = {3, 7, 9}; o2 = {2.2, 3.3, 4.4}


It is not clear what the OP exactly needs at the extremities of the new list: for example here do the values for time steps 1 and 2 have to be 0 or 2.2 ? Same question for time step 10 ? The different approaches below show some differences regarding to this, but they could be easily adapted.

TimeSeries

Normal@Quiet@
TimeSeriesResample[TimeSeries[o2, {o1}], {1, n, 1},
ResamplingMethod -> {"Interpolation", InterpolationOrder -> 0}]


{{1, 2.2}, {2, 2.2}, {3, 2.2}, {4, 2.2}, {5, 2.2}, {6, 2.2}, {7, 3.3}, {8, 3.3}, {9, 4.4}, {10, 4.4}}

FoldList & SparseArray

(*v2 = FoldList[If[#2 == 0, #1, #2] &, Normal@SparseArray[o1 -> o2, n]]*)
v2 = FoldList[Max, Normal@SparseArray[o1 -> o2, n]]


{0, 0, 2.2, 2.2, 2.2, 2.2, 3.3, 3.3, 4.4, 4.4}

To compare to the previous result with the time series form:

Thread@{Range@n, v2}


{{1, 0}, {2, 0}, {3, 2.2}, {4, 2.2}, {5, 2.2}, {6, 2.2}, {7, 3.3}, {8, 3.3}, {9, 4.4}, {10, 4.4}}

ConstantArray & Differences

v2 = Join @@ MapThread[ConstantArray, Join @@@ {{{0}, o2}, {{First@o1 - 1},
Differences@o1, {n - Last@o1 + 1}}}]


{0, 0, 2.2, 2.2, 2.2, 2.2, 3.3, 3.3, 4.4, 4.4}

Thread@{Range@n, v2}


{{1, 0}, {2, 0}, {3, 2.2}, {4, 2.2}, {5, 2.2}, {6, 2.2}, {7, 3.3}, {8,3.3}, {9, 4.4}, {10, 4.4}}

## Performance

(Comparison of all the user's answers)

Here are some timings on my macbook for the following parameters:

SeedRandom;
n = 3000000;
n1 = 2500000; (*and n1=500000*)
v1 = Range[n];
o1 = Union@RandomInteger[{1, n}, n1];
o2 = Reverse[o1/100 // N];


(Timings are in seconds): Ciao has provided the fastest answer ! It's a very instructive answer as it shows how fast can be working with Part and list of indices ([[{i1,i2,...}]]) and to always deal as much as possible with integers (FoldList[Max,...]).

Also, as commented below, the built-in TimeSeriesResample function is 40.4/0.3 =about 130x slower than Ciao's ...

• I don't think any of these is are the fastest possible way, but seeing as OP code is broken (or question is underspecified), moot at this point... – ciao Sep 4 '15 at 22:15
• Innovative use of Join, MapThread, ConstantArray and Differences. I get it to be about an order of magnitude faster than the original loop method. – Jack LaVigne Sep 4 '15 at 23:13
• If you have time, please add my answer to timings - curious of results on your hardware, it is much faster on mine... – ciao Sep 5 '15 at 5:43
• BTW, if you do re-test, might want to do it reflective of OP example - you have the opposite: OP keeps most of the data, you are dropping most of it - should change n1 to about 2500000 to mimic OP.... I'd venture the order of performance of the different solutions might get juggled if the OP example is followed – ciao Sep 5 '15 at 6:22
• The interesting question is why TimeSeries is so slow. Is this a bug? Does WRI care? – Rolf Mertig Sep 6 '15 at 7:30

Expanding @ciao's comment above:

n = 3 10^6;
n1 = 10000;
v1 = Range[n];
o1 = Union@RandomInteger[{1, n}, n1];
o2 = Reverse[ o1/1000 // N]; (*Anything*)

getV2[v1_, o1_, o2_] := Module[{f, v2},
v2 = 0 v1;
v2[[o1]] = o1;
(f[#1] = #2) & @@@ Transpose[{o1, o2}];
f = o2[];
v2 = f /@ FoldList[Max, v2]
]

getV2[v1, o1, o2]; // Timing

• The beauty of direct references... +1 o/c! – ciao Sep 2 '15 at 5:47
• I don't think this produces what the OP is asking for. I think he is asking for a zero-order interpolation. – m_goldberg Sep 2 '15 at 14:34
• @m_goldberg And I don´t think what you do think :) Who knows! – Dr. belisarius Sep 2 '15 at 19:14
• @m_goldberg This answer substitues the previous valid non-zero value for zero values in v2 (a direct substitution of o2 values is used). I applied Matvei Kruglyak code to a short version of example data and it matched belisarius's result. I believe it answers the question correctly. – Jack LaVigne Sep 2 '15 at 19:51

My interpretation of this question is that the OP is asking for a zero-order interpolation of his data. There for I propose the following solution.

First a proof of concept using about 50 points based on belisarius' data model.

SeedRandom;
Module[{interp, n, n1, o1, o2},
n = 50;
n1 = 15;
o1 = {1, Sequence @@ Union@RandomInteger[{2, n - 1}, n - n1], n};
o2 = Reverse[o1/1000 // N];
interp = Interpolation[Transpose[{o1, o2}], InterpolationOrder -> 0];
v2 = Table[interp[i], {i, 1, n}]];

ListPlot[v2, Filling -> Bottom] Here is timing with 3*^6 points with about 1000 gaps punched into it.

SeedRandom;
Module[{interp, n, n1, o1, o2},
n = 3*^6;
n1 = 1*^3;
o1 = {1, Sequence @@ Union@RandomInteger[{2, n - 1}, n - n1], n};
o2 = Reverse[o1/1000 // N];
Timing[
interp = Interpolation[Transpose[{o1, o2}], InterpolationOrder -> 0];
Table[interp[i], {i, 1, n}];]]

{11.7651, Null}


This generate the required list about 8500 times faster than the OP's estimate of 30 elements per second.

• When I make a table of a short example o1 with gaps and associated o2, it becomes apparent that rather than getting the previous value it uses the next value. Here there is a gap between 3 and 8 Table[{i, Interpolation[{{1, 1}, {2, 5}, {3, 7}, {8, 2}, {9, 3}, {10, 1}}, InterpolationOrder -> 0][i]}, {i, 1, 10}] and the result {{1, 5}, {2, 5}, {3, 7}, {4, 2}, {5, 2}, {6, 2}, {7, 2}, {8, 2}, {9, 3}, {10, 1}} shows that 4 thru 8 are filled with 2's. Kruglyak's question asked for a previous value (check his code on a short sample). Can your answer be modified to provide that? – Jack LaVigne Sep 2 '15 at 19:48