Conseils pour se débarrasser des boucles

Question

J'ai écrit un programme qui fonctionne avec le problème 3n + 1 (alias "nombres merveilleux" et diverses autres choses). Mais il a une double boucle. Comment pourrais-je le vectoriser ?

le code est

count <- vector("numeric", 100000)
L <- length(count)

for (i in 1:L)
{
x <- i
   while (x > 1)
   {
   if (round(x/2) == x/2) 
     {
     x <- x/2
     count[i] <- count[i] + 1 
     } else
     {
     x <- 3*x + 1
     count[i] <- count[i] + 1
     }
   }
}

Merci !

Answer 1

Je l'ai transformé en créant un vecteur x dont le ième élément est la valeur après chaque itération de l'algorithme. Le résultat est relativement intelligible car

f1 <- function(L) {
    x <- seq_len(L)
    count <- integer(L)
    while (any(i <- x > 1)) {
        count[i] <- count[i] + 1L
        x <- ifelse(round(x/2) == x/2, x / 2, 3 * x + 1) * i
    }
    count
}

Ceci peut être optimisé pour (a) suivre uniquement les valeurs encore en jeu (via idx) et (b) éviter les opérations inutiles, par exemple, ifelse évalue les deux arguments pour toutes les valeurs de x, x/2 évalué deux fois.

f2 <- function(L) {
    idx <- x <- seq_len(L)
    count <- integer(L)
    while (length(x)) {
        ix <- x > 1
        x <- x[ix]
        idx <- idx[ix]
        count[idx] <- count[idx] + 1L
        i <- as.logical(x %% 2)
        x[i] <- 3 * x[i] + 1
        i <- !i
        x[i] <- x[i] / 2
    }
    count
}

avec f0 la fonction originale, j'ai

> L <- 10000
> system.time(ans0 <- f0(L))
   user  system elapsed 
  7.785   0.000   7.812 
> system.time(ans1 <- f1(L))
   user  system elapsed 
  1.738   0.000   1.741 
> identical(ans0, ans1)
[1] TRUE
> system.time(ans2 <- f2(L))
   user  system elapsed 
  0.301   0.000   0.301 
> identical(ans1, ans2)
[1] TRUE

Une astuce consiste à mettre à jour les valeurs impaires en 3 * x[i] + 1, puis à effectuer la division par deux sans condition.

x[i] <- 3 * x[i] + 1
count[idx[i]] <- count[idx[i]] + 1L
x <- x / 2
count[idx] <- count[idx] + 1

Avec ceci en f3 (je ne sais pas pourquoi f2 est plus lent ce matin !) j'obtiens

> system.time(ans2 <- f2(L))
   user  system elapsed 
   0.36    0.00    0.36 
> system.time(ans3 <- f3(L))
   user  system elapsed 
  0.201   0.003   0.206 
> identical(ans2, ans3)
[1] TRUE

Il semble que de plus grandes étapes peuvent être prises à l'étape de division par deux, par exemple, 8 est 2^3 donc nous pourrions prendre 3 étapes (ajouter 3 pour compter) et avoir fini, 20 est 2^2 * 5 donc nous pourrions prendre deux étapes et entrer dans l'itération suivante à 5. Implémentations ?

Answer 2

Parce que vous avez besoin d'itérer sur les valeurs de x on ne peut pas vraiment vectoriser ça. À un moment donné, R doit travailler sur chaque valeur de x séparément et à tour de rôle. Il est possible d'exécuter les calculs sur des cœurs de processeur distincts pour accélérer les choses, en utilisant par exemple foreach dans le paquet du même nom.

Sinon, (et ceci ne fait que vous cacher la boucle), enveloppez le corps principal de votre boucle comme une fonction, par ex :

wonderous <- function(n) {
    count <- 0
    while(n > 1) {
        if(isTRUE(all.equal(n %% 2, 0))) {
            n <- n / 2
        } else {
            n <- (3*n) + 1
        }
        count <- count + 1
    }
    return(count)
}

et ensuite vous pouvez utiliser sapply() pour exécuter la fonction sur un ensemble de chiffres :

> sapply(1:50, wonderous)
 [1]   0   1   7   2   5   8  16   3  19   6  14   9   9  17  17
[16]   4  12  20  20   7   7  15  15  10  23  10 111  18  18  18
[31] 106   5  26  13  13  21  21  21  34   8 109   8  29  16  16
[46]  16 104  11  24  24

Ou vous pouvez utiliser Vectorize pour retourner une version vectorisée de wonderous qui est elle-même une fonction qui vous cache encore plus de choses :

> wonderousV <- Vectorize(wonderous)
> wonderousV(1:50)
 [1]   0   1   7   2   5   8  16   3  19   6  14   9   9  17  17
[16]   4  12  20  20   7   7  15  15  10  23  10 111  18  18  18
[31] 106   5  26  13  13  21  21  21  34   8 109   8  29  16  16
[46]  16 104  11  24  24

Je pense que c'est à peu près tout ce que vous pouvez faire avec les outils R standard pour le moment. @Martin Morgan montre que l'on peut faire beaucoup mieux que cela en résolvant de manière ingénieuse le problème que pose le fait que fait a utilisé les capacités vectorielles de R.

Answer 3

Une approche différente reconnaît que l'on revient souvent sur des chiffres faibles, alors pourquoi ne pas les retenir et économiser le coût du nouveau calcul ?

memo_f <- function() {
    e <- new.env(parent=emptyenv())
    e[["1"]] <- 0L
    f <- function(x) {
        k <- as.character(x)
        if (!exists(k, envir=e))
            e[[k]] <- 1L + if (x %% 2) f(3L * x + 1L) else f(x / 2L)
        e[[k]]
    }
    f
}

ce qui donne

> L <- 100
> vals <- seq_len(L)
> system.time({ f <- memo_f(); memo1 <- sapply(vals, f) })
   user  system elapsed 
  0.018   0.000   0.019 
> system.time(won <- sapply(vals, wonderous))
   user  system elapsed 
  0.921   0.005   0.930 
> all.equal(memo1, won) ## integer vs. numeric
[1] TRUE

Cela pourrait ne pas être bien parallélisé, mais peut-être que ce n'est pas nécessaire avec une vitesse de 50 fois supérieure ? De plus, la récursion pourrait devenir trop profonde, mais elle pourrait être écrite comme une boucle (ce qui est probablement plus rapide, de toute façon).