UP | HOME

Séance 6 : Algo & Prog avec R

Table des matières

1 Recopie et clonage

Le problème de la copie et du clonage est important dans les langages de programmation.

  1. Exécutons la session ci-dessous au toplevel R :
    x <- 1:5
    y <- x
    x[1] <- -1
    print(x)
    print(y)
    
    [1] -1  2  3  4  5
    [1] 1 2 3 4 5
    
  2. Voici le résultat de la même session en Python :
    x = range(1,6)
    y= x
    x[0]=-1
    print(x)
    print(y)
    
    [-1, 2, 3, 4, 5]
    [-1, 2, 3, 4, 5]
    

2 Le tri par fusion (mergesort)

  1. Programmez la fonction Fusion(x, y) prenant deux vecteurs croissants d’entiers x et y, et retournant un vecteur croissant mélange des deux.
    Fusion <- function(x,y) {
      ## x et y triées croissantes, resultat trié croissant
      res <- numeric(length(x)+length(y)) ; 
      i <- i1 <- i2 <- 1
      while (i1 <= length(x) && i2 <= length(y) ) {
        if (x[i1] <= y[i2]) {
          res[i] <- x[i1] ; 
          i1 = i1 + 1
        } else {
          res[i] <- y[i2] 
          i2 = i2 + 1
        }
        i = i + 1
      }
      ## je dois verifier pourquoi je suis sorti de la boucle while !
      if(i1 > length(x)) res[i:length(res)] <- y[i2:length(y)]
      else res[i:length(res)] <- x[i1:length(x)]
      return(res)
    }
    
    x <- c(2, 5, 6, 6, 8, 12, 20)
    y <- c(0, 2, 7, 8, 15, 19, 23)
    print(paste('x = [',paste(x,collapse=','),'] et y =',paste(y,collapse=',')))
    print(paste('Fusion : [',paste(Fusion(x,y),collapse=','),']'))
    
    [1] "x = [ 2,5,6,6,8,12,20 ] et y = 0,2,7,8,15,19,23"
    [1] "Fusion : [ 0,2,2,5,6,6,7,8,8,12,15,19,20,23 ]"
    
  2. Le principe du tri par fusion suit le diagramme dichotomique ci-dessous, dans lequel la scission coupe en deux sous-listes de même longueur ou presque [à 1 élément près] :
        L = [8, 5, 2, 3, 1, 5, 4, 9, 2]
                   /             \
                  /               \
    LT1 = [1, 2, 4, 5, 8]   LT2 = [2, 3, 5, 9]
                  \               /
                   \             /
         LT = [1, 2, 2, 3, 4, 5, 5, 8, 9]
    
  3. Utilisez la fonction Fusion pour programmer par récurrence une fonction TriFusion(x) prenant un vecteur de nombres réels x et retournant une copie triée en ordre croissante de ce vecteur.
    TriFusion <- function(x) {
      n <- length(x)
      m =n%/%2
      ## Attention au cas de base, n == 0 ne suffit pas !
      if(n <= 1) return(x)         
      else return(Fusion(TriFusion(x[1:(m)]),TriFusion(x[(m + 1):n])))
    }
    
    x <- sample(1:100, 10, replace = TRUE)
    print(paste('Test du tri fusion : L = [',paste(x,collapse=','),']'))
    print(paste('Test du tri fusion : LT = [',paste(TriFusion(x),collapse=','),']'))
    
    [1] "Test du tri fusion : L = [ 10,22,19,93,76,78,25,53,14,97 ]"
    [1] "Test du tri fusion : LT = [ 10,14,19,22,25,53,76,78,93,97 ]"
    
  4. Comparez les temps d'exécution des fonctions TriFusion et sort grâce à system.time.
    x <- sample(1:10000, 10000, replace = TRUE)
    t <- 5
    system.time(replicate(t, TriFusion(x)))
    system.time(replicate(t, sort(x)))
    
    utilisateur     système      écoulé 
          1.774       0.000       1.774
    utilisateur     système      écoulé 
          0.005       0.000       0.005
    
  5. Calcul de la complexité
    1. Définissez par récurrence la suite C(N) donnant le nombre d'opérations élémentaires en fonction du nombre d'éléments à trier N.
    2. Donnez le terme général de la suite en supposant que N=2p.

3 Crible d’Eratosthène

Il s’agit d’un algorithme célèbre permettant de construire les entiers premiers jusqu’à n :

  1. On commence par construire le vecteur x des entiers de 2 jusqu’à n inclus.

Nous allons maintenant rayer les nombres non premiers de ce vecteur !

  1. Le premier nombre non rayé est premier (pourquoi ?). Je raye tous ses multiples. Et je répète cette étape jusqu’à avoir examiné tous les nombres du vecteur.

Exemple, n = 26 :

2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26
2  3  .  5  .  7  .  9   .  11   .  13   .  15   .  17   .  19   .  21   .  23   .  25   .
2  3  .  5  .  7  .  .   .  11   .  13   .   .   .  17   .  19   .   .   .  23   .  25   .
2  3  .  5  .  7  .  .   .  11   .  13   .   .   .  17   .  19   .   .   .  23   .   .   .

Programmez la fonction Crible(n) retournant la liste des nombres premiers de [2,n].

Crible <- function(n) {
  x <- seq(3,n,2)
  y <- c(2)
  if(length(x) > 1) {
    while( x[1]*x[1] <= tail(x,1) ) {
      y <- append(y, x[1])
      x <- x[ x %% x[1] != 0]
    }
  }
  return(append(y, x))
}
print(Crible(100))

t <- 10
n <- 75000;
system.time(replicate(t,Crible(n)))
 [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
utilisateur     système      écoulé 
      0.173       0.002       0.174
crible_v2 <- function(n) {
  L <- c(2)
  C <- seq(3,n,2)
  while( length(C) > 0 && C[1]*C[1] <= n ) {
    L <- c(L,C[1])
    C <- C[ C %% C[1] != 0];
  }
  L <- c(L,C)
  return(L)
}
print(crible_v2(100))

t <- 10
n <- 75000;
system.time(replicate(t,crible_v2(n)))
 [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
utilisateur     système      écoulé 
      0.180       0.001       0.181
  1. Analyse des performances

    Comparez sur machine les temps de calcul pour construire la liste des nombres premiers jusqu’à 50000 :

    1. En utilisant le crible d’Eratosthènes vu ci-dessus.
    2. En utilisant la fonction estPremier du TP2.
    ####################  
    ## Voir TP 2
    ppdiv <- function(n) {
        ## avec au plus sqrt(n) operations
      if (n %% 2 == 0 ) return(2)
      d = 3
      ## acc to compute the d*d
      acc = 9 
      while (n %% d != 0 ) {
        ## <==> d > sqrt(n) mais plus rapide !
        if(acc > n) return(n)    
        d = d + 2
        acc = acc + 4*(d-1)
      }
      return(d)
    }
    
    estPremier <- function(n) (n >= 2) && (ppdiv(n) == n)      
    
    ####################  
    crible_naif <- function(n) {
        L=c(2);
        for (i in seq(3,n,2)) {
          if (estPremier(i)) L[length(L)+1] <- i;
        }
        return(L)
      } 
    crible_naif(100)
    
    t=10
    n=50000;
    system.time(replicate(t,crible_naif(n)))
    
     [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
    utilisateur     système      écoulé 
          6.336       0.002       6.332 
    

4 Dichotomie

Reprenons la recherche par dichotomie du cours 6. Le programme donné en cours utilise beaucoup de constructions de tranches qui sont des recopies coûteuses de portions de listes.

Pour les éviter et accélérer la recherche, vous allez travailler sur place (sans utiliser de listes supplémentaires) en gérant deux indices a et b. Au départ a=0 et b=length(L)-1. Vous allez faire se rapprocher les curseurs a et b, en abandonnant la moitié de l’intervalle [a,b] chaque fois. À vous !

## La liste triée
LT <- sort(runif(50000));
n <- length(LT)%/%20
## Les valeurs à chercher
xt <- c(sample(LT, n),runif(n))

micro_bench <- function(x, L, f_dicho) {
  res <- numeric(length(x));
  for (i in seq_along(x)) {
    res[i] <- f_dicho(x[i],L)     
  }
  return(res)
}
   rech_dicho <- function(x, L) {
     ## recursive, retourne un i tel que L[i] == x, ou bien NA
     m = (length(L)+1) %/% 2
     ## any pour la liste vide
     if (any(L[m] == x)) return(m)
     else if(length(L) > 1) {
         if(L[m] > x)  return(rech_dicho(x, L[1:(m-1)])) 
         else return( m + rech_dicho(x, L[(m+1):length(L)]))
     }
     return(NA)
   }
print(paste('Liste vide :', rech_dicho(0,c())))
x <- c(10,2,12)
L <- c(2,5,6,6,8,10,15,18,21,22,22,25,30)
print(sprintf('Liste trié [%s]', paste(L, collapse=', ')))     
for (i in x) {
  print(sprintf('Indice de %s : %d', i, rech_dicho(i,L)))     
} 
 system.time(micro_bench(xt, LT, rech_dicho))
[1] "Liste vide : NA"
[1] "Liste trié [2, 5, 6, 6, 8, 10, 15, 18, 21, 22, 22, 25, 30]"
[1] "Indice de 10 : 6"
[1] "Indice de 2 : 1"
Erreur dans sprintf("Indice de %s : %d", i, rech_dicho(i, L)) (depuis #2) : 
  format incorrect '%d' ; utilisez les formats %f, %e, %g ou %a pour les objets numériques
utilisateur     système      écoulé 
      2.487       0.000       2.484
rech_dicho <- function(x, L) {
  ## iterative, retourne un i tel que L[i] == x, ou bien -1
  a <- 1
  b <- length(L)
  while(a <= b) {
    m <- (a + b) %/% 2
    if (L[m] == x) return(m)
    else if(L[m] > x)  b <- m - 1
    else a <- m + 1
  }
  return(NA)
}
print(paste('Liste vide :', rech_dicho(0,c())))
print(sprintf('Liste trié [%s]', paste(L, collapse=', ')))     
for (i in x) {
  print(sprintf('Indice de %s : %d', i, rech_dicho(i,L)))     
} 
system.time(micro_bench(xt, LT, rech_dicho))
[1] "Liste vide : NA"
[1] "Liste trié [2, 5, 6, 6, 8, 10, 15, 18, 21, 22, 22, 25, 30]"
[1] "Indice de 10 : 6"
[1] "Indice de 2 : 1"
[1] "Indice de 12 : NA"
utilisateur     système      écoulé 
      0.116       0.000       0.116

5 Exceptions [d’après IUT Orsay].

  1. Définissez la fonction \(f(x) = sin(x)/x\).
  2. Calculez les valeurs de f pour les entiers compris entre -3 et 3
  3. Levez une erreur en cas de division par 0.
  4. Lancez un avertissement quand x devient très petit.
  5. Traitez les erreurs en renvoyant la limite de f(x) en 0.
  6. Traitez les avertissements en utilisant le developpement limité d'ordre 3 au voisinage de 0 de la fonction \(sin(x) \simeq x - \frac{x^3}{6}\).
  ## Evalute the desired series of expressions inside of tryCatch
  f <- function(x) {
    if(x == 0) stop("Division by 0");
    if(abs(x) <= 2**(-20)) warning("Number too small")
    return(sin(x) / x)      
}
## et tant pis si x == 0    :-)

  ## je prefere renvoyer une matrice plutot qu'une liste de listes :
  table_f <- function(x) {         
  ## tableau de points [(a,f(a)),(a+1,f(a+1)),...,(b,f(b))]
    M=numeric(length(x))
    for (i in seq_along(x)) {
      M[i] <- tryCatch({
        f(x[i])
      }, warning = function(war) {
        ## warning handler picks up where error was generated
        print(paste("MY_WARNING:  ",war))
        ## le dev. limité de sin(x) = x - x^3/6 => sin(x)/x ~= 1 - x^2/6
        return(1 - x[i]**2/6)
      }, error = function(err) {
        ## error handler picks up where error was generated
        print(paste("MY_ERROR:  ",err))
        return(1) # lim sin(x)/x = 1 lorsque x -> 0 (j'ai des maths...)
      }
      ) # END tryCatch
    }
    return(M)
  }
print('Tableau de points pour sin(x)/x :')
print(table_f(-3:3))
x <- 2**-16
print(sprintf("%.15f %.15f", table_f(x), sin(x)/x))
x <- 2**-27
print(sprintf("%.15f %.15f", table_f(x), sin(x)/x))
[1] "Tableau de points pour sin(x)/x :"
[1] "MY_ERROR:   Error in f(x[i]): Division by 0\n"
[1] 0.0470400 0.4546487 0.8414710 1.0000000 0.8414710 0.4546487 0.0470400
[1] "0.999999999961195 0.999999999961195"
[1] "MY_WARNING:   simpleWarning in f(x[i]): Number too small\n"
[1] "1.000000000000000 1.000000000000000"

Created: 2017-11-02 jeu. 13:59