Nested Loop Running for 12+ hours Now

[Brian Green]

Just checking to see if this is normal - I'm trying to create a grid of values over which I'm going to be searching for a condition, but just setting up the dataset has already taken roughly 12 hours.

clear all
set obs 10240000

gen sigma = 0
gen gamma = 0
gen alpha = 0
gen f = 0
gen f_h = 0
gen tau = 0
gen t = 0


local n = 1
forvalues d = 1(.25)10 {
    forvalues h = 1(.25)10 {
        forvalues tb = 2.5(.25)4.5 {
            forvalues s = 1.25(1)20.25 {
                forvalues tx = .15(.5).35 {
                    forvalues a = .5(.1).8 {
                        replace alpha = `a' if _n==`n'
                        replace tau = `tx' if _n==`n'
                        replace sigma = `s' if _n==`n'
                        replace t = `tb' if _n==`n'
                        replace f_h = `h' if _n==`n'
                        replace f = `d' if _n==`n'
                        local n = `n'+1
                    }
                }
                
            }
            
        }
        
    }
}

replace gamma = 2*(sigma-1)+.1

[Carlo Lazzaro]

Brian:
your code works for me (I would add -quietly-) if I

set obs 10

and gives back the following results

 list

     +---------------------------------------------+
     | sigma   gamma   alpha   f   f_h   tau     t |
     |---------------------------------------------|
  1. |  1.25      .6      .5   1     1   .15   2.5 |
  2. |  1.25      .6      .6   1     1   .15   2.5 |
  3. |  1.25      .6      .7   1     1   .15   2.5 |
  4. |  1.25      .6      .8   1     1   .15   2.5 |
  5. |  2.25     2.6      .5   1     1   .15   2.5 |
     |---------------------------------------------|
  6. |  2.25     2.6      .6   1     1   .15   2.5 |
  7. |  2.25     2.6      .7   1     1   .15   2.5 |
  8. |  2.25     2.6      .8   1     1   .15   2.5 |
  9. |  3.25     4.6      .5   1     1   .15   2.5 |
 10. |  3.25     4.6      .6   1     1   .15   2.5 |
     +---------------------------------------------+

The procedure lasted 59 seconds (as per -rmsg-).

[Leonardo Guizzetti]

Your code is indeed quite slow because Stata has to read the entire dataset with every iteration of each loop and also edit it each time for just one change at a time. Here's how I would think about approaching your problem using -fillin- which is seldom used but is quite handy when trying to make these kinds of factorial/crossed datasets. Alternative solutions using -merge- or -cross- are also possible but I think this probably saves the most typing. Refer to the documentation for -fillin-, but it is doing the heavy lifting of creating every combination among all listed variables. You'll need to drop the results in the final dataset that get crossed with missing values as those aren't of interest to you.

As an aside, I think you have a typo in your -tx- increment which I assume should be 0.05 instead 0.5. The below code runs in <1 minute for me.

clear *
cls

set obs 100

gen double d = 1 + 0.25*(_n-1)
qui replace d = . if !inrange(d, 1, 10)

gen double h = 1 + 0.25*(_n-1)
qui replace h = . if !inrange(h, 1, 10)

gen double tb = 2.5 + 0.25*(_n-1)
qui replace tb = . if !inrange(tb, 2.5, 4.5)

gen double s = 1.25 + 1*(_n-1)
qui replace s = . if !inrange(s, 1.25, 20.25)

gen double tx = 0.15 + 0.05*(_n-1)
qui replace tx = . if !inrange(tx, 0.15, 0.35)

gen double a = 0.5 + 0.1*(_n-1)
qui replace a = . if !inrange(a, 0.5, 0.8)

egen nmiss = rowmiss(d h tb s tx a)
drop if nmiss==6
drop nmiss

list

fillin d h tb s tx a

egen nmiss = rowmiss(d h tb s tx a)
drop if nmiss
drop nmiss

gen double gamma = 2*(s-1)+.1

[Clyde Schechter]

You can speed this up astronomically. Replace -if _n == `n'- everywhere by -in `n'-, and wrap all those -replace- statements in a -quietly { ... } - block. In my setup the whole thing completes in less than 20 seconds:

. clear all

. timer on 1

. set obs 10240000
Number of observations (_N) was 0, now 10,240,000.

.
. gen sigma = 0

. gen gamma = 0

. gen alpha = 0

. gen f = 0

. gen f_h = 0

. gen tau = 0

. gen t = 0

.
.
. local n = 1

. forvalues d = 1(.25)10 {
  2.     forvalues h = 1(.25)10 {
  3.         forvalues tb = 2.5(.25)4.5 {
  4.             forvalues s = 1.25(1)20.25 {
  5.                 forvalues tx = .15(.5).35 {
  6.                     forvalues a = .5(.1).8 {
  7.                                                 quietly {
  8.                                                         replace alpha = `a' in `n'
  9.                                                         replace tau = `tx' in `n'
 10.                                                         replace sigma = `s' in `n'
 11.                                                         replace t = `tb' in `n'
 12.                                                         replace f_h = `h' in `n'
 13.                                                         replace f = `d' in `n'
 14.                                                         local n = `n'+1
 15.                                                 }
 16.                     }
 17.                 }
 18.                
.             }
 19.            
.         }
 20.        
.     }
 21. }

.
. replace gamma = 2*(sigma-1)+.1
(10,240,000 real changes made)

. timer off 1

. timer list 1
   1:     18.43 /        1 =      18.4290

You cannot overestimate how slow -if- is. On every iteration of the innermost loop, which gets iterated over 5,000,000 times, you have 6 commands using -if _n == `n'-. The problem is that each execution of -if- requires that Stata scan every observation in the data set to determine if it meets the condition, and there are over 10 million observations. So your code requires Stata to check 5 * 10<sup>13</sup> observations in order to complete. By contrast, -in `n'- requires no condition checking at all: Stata can internally pick out the `n'th observation without looking at any of them.

Added: Crossed with #2, which offers another solution.

Also Added: If Leonardo Guizzetti is correct that the increment for tx should be .05 rather than .5, the execution time for this version on my setup is 85 seconds.

[Daniel Schaefer]

In addition to if, writing to the console is very slow, so I agree with the advice in #2. I'll go ahead and add my optimization attempt to the pile.

clear all
set obs 10240000
* alpha
egen alpha = fill(.5(.1).8 .5(.1).8)
* tau
egen tau = seq(), from(1) to(5) block(4)
replace tau = (tau / 20) + 0.1
* sigma
egen sigma = seq(), from(1) to(20) block(20)
replace sigma = sigma + 0.25
* t
egen t = seq(), from(1) to(9) block(400)
replace t = (t / 4) + 2.25
* f_h
egen f_h = seq(), from(1) to(37) block(3600)
replace f_h = (f_h / 4) + 0.75
* f
egen f = seq(), from(1) to(37) block(133200)
replace f = (f / 4) + 0.75
* gamma
gen gamma = 2*(sigma-1)+.1

[Daniel Schaefer]

By the way, this loop doesn't behave the way I would expect for some reason on my end. It only iterates once...

forvalues tx = .15(.5).35 {
    display `tx'
}

. forvalues tx = .15(.5).35 {
  2.         display `tx'
  3. }
.15

Edit: Leonardo gives the reason in #3

[Clyde Schechter]

I want to raise another point about the code in #1 that doesn't affect the timing. It is not a good idea to write loops like

forvalues tx = .15(.05).35 {
     forvalues a = .5(.1).8 {

The difficulty is that neither .05 nor .1 can be represented exactly in a finite number of binary digits (just as 1/3 has no finite decimal representation). So the actual increment Stata uses is some number that is as close to .05 (resp. .1) as it can get in the number of bits available to it. This "rounding error" will accumulate at each iteration of the loop. In addition to the fact that the iterations will not actually have tx = .20, .25, .30, and .35 (which is unavoidable as none of these numbers other than .25 has a finite binary representation), it is possible that the .35 iteration will never take place if the rounding errors have been "rounded upwards." This is because the final version of tx may exceed 0.35.

Now, as it happens, with these particular numbers, it works out OK. You get reasonable approximations to the intended values, and the upper limits are not bypassed. But, in general, a safer way to do this is:

forvalues txi = 15(5)35 {
    local tx = `txi'/100
    forvalues ai = 5(1)8 {
        local a = `ai'/10

This doesn't avoid rounding errors on the values of `a' and `tx' at most of the steps of the loop, but it does guarantee that the final iteration will not be bypassed.

The other loops are not problematic because the increments of .25 and 1 do have finite binary representations. While it is not hard to identify which increments have finite binary representations and which ones do not, it is frankly a better habit to just avoid using floating point numbers to bound or increment loops and always use integers for the purpose. That way if you overlook something or miscalculate whether it is finitely representable in binary, everything will still work out correctly.

Added: Crossed with #5 and #6. #5 does not explicitly raise the issue of the impossibility of finitely representing .05 or .1 in binary, but the solution offered there does overcome the problem because the -egen, seq()- function bounds the iteration with integer values in its -from()- and -to()- options.

[Leonardo Guizzetti]

I believe that my solution in #3 will also avoid the issue of rounding precision raised by Clyde in general, only if the fractional values (e.g., 0.25) are replaced by their fractional representation (e.g., 1/4).

[daniel klein]

I believe that my solution in #3 will also avoid the issue of rounding precision raised by Clyde in general, only if the fractional values (e.g., 0.25) are replaced by their fractional representation (e.g., 1/4).

I don't follow. The code in #3 stores all variables in double precision, matching (most of) Stata’s computations. The example here is confusing because 0.25 has an exact binary representation so there’s no rounding issue. In general, whether you type .05 as a decimal or as 5/100 should rarely (if ever) affect results. I know it's getting a it off-topic, but could you elaborate on your thoughts?

[Daniel Schaefer]

(At the risk of getting even more off topic, daniel klein and others following this thread might enjoy this article: https://chadnauseam.com/coding/random/calculator-app)

[Leonardo Guizzetti]

I don't follow. The code in #3 stores all variables in double precision, matching (most of) Stata’s computations. The example here is confusing because 0.25 has an exact binary representation so there’s no rounding issue. In general, whether you type .05 as a decimal or as 5/100 should rarely (if ever) affect results. I know it's getting a it off-topic, but could you elaborate on your thoughts?

Yeah, 0.25 was a poor example. The mitigation is only in regards to the precision that Clyde warned about. Let's simplify the example for argument's sake and use 1/3 or it's approximation 0.3333, which has no exact binary equivalent.

set obs 28

gen double x = 1 + (0.3333)*(_n-1)
qui replace x = . if !inrange(x, 1, 10)

gen double y = 1 + (1/3)*(_n-1)
qui replace y = . if !inrange(y, 1, 10)

list

The final values are 9.9991 vs 10, respectively. The division of _n by 3, instead of multiplication by 0.3333, results in the correct spacing of values. That's all I intended to say.

[daniel klein]

Let's simplify the example for argument's sake and use 1/3 or it's approximation 0.3333, which has no exact binary equivalent.

Fair enough. However,1/3 doesn't have a finite decimal representation either. The issue only arises with numbers like this. Fractions that do have a finite decimal representation, such as 0.1, don't pose a problem, regardless of whether their binary representation is finite:

set obs 28

gen double x = 1 + (0.1)*(_n-1)
qui replace x = . if !inrange(x, 1, 10)

gen double y = 1 + (1/10)*(_n-1)
qui replace y = . if !inrange(y, 1, 10)

list

yields the same results, at least in Stata (other software might differ).


Added: Here is an example of the issue that Clyde Schechter has in mind:

. forvalues i = 4(.1)4.2 {
  2.     display `i'
  3.     display %21x `i'
  4. }
4
+1.0000000000000X+002
4.1
+1.0666666666666X+002
4.2
+1.0ccccccccccccX+002

. forvalues i = 40(1)42 {
  2.     local x = `i'/10
  3.     display `x'
  4.     display %21x `x'
  5. }
4
+1.0000000000000X+002
4.1
+1.0666666666666X+002
4.2
+1.0cccccccccccdX+002

Note that 4.2 does not have a finite binary representation. Incrementing the (first) loop by .1, which also has no finite binary representation, will end at \[i = +1.0ccccccccccccX+002 = (1+ 0/16 + 12/16^2 + 12/16^3 + ... + 12/16^{13}) * 2^2 \approx 4.199...\], whereas using integers in the second loop ends at i=42 and \[x = +1.0cccccccccccdX+002 = (1+ 0/16 + 12/16^2 + 12/16^3 + ... + 13/16^{13}) * 2^2 \approx 4.200...\].

[Brian Green]

Thank you all for your assistance, I greatly appreciate it.