Programmer en R/Programmation procédurale

Le langage R permet de faire de la programmation procédurale, ce qui donne accès à la mise en œuvre d'algorithmes complexes.

Structure générale d'un programme

Un programme est une suite d'instructions, qui peuvent être séparées par un point-virgule « ; » ou par une fin de ligne. Un point-virgule indique toujours la fin d'une instruction, alors qu'une fin de ligne est ignorée si la syntaxe de l'instruction n'est pas complète.

Les instructions peuvent être regroupées en blocs délimités par une paire d'accolades « {…} ».

Une instruction est évaluée lorsque sa syntaxe est complète et que l'on introduit une fin de ligne ; le programme R renvoie alors la valeur issue de l'évaluation de l'instruction. Un bloc est interprété comme une instruction, c'est-à-dire qu'il n'est évalué qu'à la fin de ligne suivant l'accolade fermante.

Un commentaire est introduit par le croisillon #.

Interactions avec l'utilisateur

On peut afficher des messages avec la commande print("Message").

Pour permettre à l'utilisateur de saisir une valeur, on utilise la commande scan(), sous la forme

a <- scan()

La fonction construit un vecteur : lorsque l'on appuie sur la touche [Entrée], il demande la valeur suivante. Pour terminer la saisie, on entre une valeur vide.

Par défaut, R considère que c'est une chaîne de caractères. On peut utiliser l'option what pour forcer un autre type :

what = numeric() : nombre réel ;
what = integer() : nombre entier ;
what = double() : nombre réel à double précision ;
what = logical() : booléen TRUE (ou T) et FALSE (ou F).

Par exemple

> a <- scan()
1: 1
2: 2
3: 3
4: 4
5: 
Read 4 items
> print(a)
[1] 1 2 3 4

On peut aussi avoir un éditeur de type tableur avec la commande data.entry :

> data.entry('a')

Il faut entrer les données dans la première colonne « a ». Il faut cliquer sur l'en-tête de la colonne pour définir si l'on veut entrer des nombres (numeric) ou des caractères (character). La variable est enregistrée lorsque l'on ferme la fenêtre du tableur.

Structures de contrôle

Évaluation conditionnelle

L'évaluation conditionnelle d'une instruction se fait avec la structure

if (instruction1)
    instruction2
else
    instruction3

où l'évaluation de instruction1 doit être un booléen :

instruction2 est évalué si instruction1 est vrai ;
instruction3 est évalué si instruction1 est faux.

La structure ci-dessus forme un bloc ; ainsi, on peut l'évaluer, et mettre la valeur dans une variable. Par exemple, les deux programmes ci-dessous sont équivalents :

if (a > b) max <- a else max <- b

et

max <- if (a > b) a else b

Si l'on veut écrire la structure sur plusieurs lignes, il faut la mettre dans un bloc :

{if (a > b)
   max <- a
else
   max <- b}

ou bien utiliser, pour instruction2, un bloc se fermant avant le else :

if (a > b)
   {max <- a
} else
   max <- b

On peut enchaîner les évaluations conditionnelles avec l'instruction elsif.

Évaluation à choix multiple

La structure switch permet de faire des choix multiples (l'équivalent de la structure case d'autres langages). La syntaxe générale est :

switch (instruction, liste)

où liste est une suite d'instructions séparées par une virgule.

Si l'évaluation de l'instruction est un nombre entier n, alors l'interpréteur exécute la n-ième instruction de la liste. Par exemple

n <- 5
switch (n, "a", "b", "c", "d", "e", "f", "g")

renvoie "e".

Si l'évaluation de l'instruction est une chaîne de caractères, alors la liste doit être une succession d'instructions de type nom = valeur. Si l'évaluation de l'instruction correspond exactement à un nom de la liste, alors l'évaluation de la structure renvoie la valeur correspondante. La liste peut contenir une valeur sans nom qui sert de valeur par défaut.

Par exemple

plat <- "fruit"
switch(plat, fruit = "banane", legume = "broccoli", "aucun")

renvoie "banane" ; avec plat <- "viande", l'évaluation renvoie "aucun". Notons que la valeur peut être quelconque : scalaire, vecteur, instruction, …

Tri, recherche et sélection

L'instruction sort() permet de classer des vecteurs. Par exemple, le code suivant classe par ordre croissant :

> x <- c(5, 8, 2, 1, 9, 6)
> sort(x)
[1] 1 2 5 6 8 9

Pour classer par ordre décroissant, on ajoute l'option decreasing = TRUE :

> x <- c(5, 8, 2, 1, 9, 6)
> sort(x, decreasing = TRUE)
[1] 9 8 6 5 2 1

On peut indiquer l'algorithme de tri utilisé : method = "shell" (cas général) pour la méthode de Shell (cas général), ou bien method = "quick" pour le tri rapide (uniquement pour les valeurs numériques).

On peut aussi ne classer que certaines valeurs : on utilise l'option partial = vecteur où vecteur est un vecteur d'indices. Le tri n'est alors pas complet, mais on est sûr que pour chaque élément désigné de la liste triée, les éléments précédents ont une valeur inférieure, et les éléments suivants ont une valeur supérieure. Cela permet d'avoir un résultat plus rapidement. Par exemple, si l'on veut connaître la valeur de la médiane, il suffit de classer l'élément médian :

x <- rnorm(21)
indice <- round(length(x)/2)+1
y <- sort(x, partial = indice)
y[indice]

Le résultat peut être différent de median(x) si le nombre d'éléments est pair. Le script suivant donne la valeur conventionnelle de la médiane, et la compare avec l'évaluation de la fonction median().

n <- 20 # nombre d'éléments
x <- rnorm(n) # données

# ***** Calcul
indicebase <- round(length(x)/2)
pair <- (n%%2 == 0) # test de parité
if (pair) {indice <- c(indicebase, indicebase+1)
    } else indice <- indicebase+1
y <- sort(x, partial = indice)
mediane <- mean(y[indice])

# ***** Affichage du résultat
print(mediane)
mediane == median(x)

L'instruction rank() donne l'ordre des éléments d'un vecteur :

y <- rank(x)

si x est un vecteur de n éléments, alors y est un vecteur de n nombres entiers (ou éventuellement fractionnaire si certaines valeur de x sont en double) ; y(i) est le rang de la valeur x(i) dans la liste x triée.

> rank(c(1, 3, 2, 5))
[1] 1 3 2 4

indique que la valeur « 1 » est la 1^re, la valeur « 3 » est la 3^e, la valeur « 2 » est la 2^e et la valeur « 5 » est la 4^e de la liste triée. En cas de doublons (ties), on peut préciser la méthode de classement :

ties.method = "average" attribue à chaque élément la moyenne des rangs ; dans l'exemple suivant, on a trois fois la valeur « 1 », donc cette valeur a trois rangs (1, 2, 3) dont la moyenne vaut (1 + 2 + 3)/ = 2 ;
ties.method = "first" incrémente le rang à chaque occurrence de la valeur, soit dans notre exemple 1, 2 puis 3 ;
ties.method = "min" attribue à chaque élément le premier rang rencontré (ici, 1, soit trois premiers ex aequo) ; c'est typiquement la méthode du classement sportif ;
ties.method = "max" attribue à chaque élément le dernier rang rencontré ;
dans l'analyse de survie et de fiabilité, la fréquence cumulée par la méthode des modes est obtenue en prenant ce rang-là divisé par la taille de l'échantillon ;
ties.method = "random" est identique à <code"first" mais les rangs sont distribués aléatoirement à chaque valeur.

> rank(c(1, 1, 1), ties.method="average")
[1] 2 2 2
> rank(c(1, 1, 1), ties.method="first")
[1] 1 2 3
> rank(c(1, 1, 1), ties.method="min")
[1] 1 1 1
> rank(c(1, 1, 1), ties.method="max")
[1] 3 3 3
> rank(c(1, 1, 1), ties.method="random")
[1] 2 1 3

L'instruction order() permet de classer des matrices et data frames.

La recherche se fait avec une paire de double-crochet : l'expression

i <- x[[a]]

renvoie l'indice de la première occurrence de la valeur a dans le vecteur ou la matrice 'x. Si x est une matrice, i est un indice unique, obtenu en parcourant la matrice colonne par colonne, de haut en bas et de gauche à droite. Par exemple :

> x <- cbind(c(6, 5, 4), c(3, 2, 1))
> print(x)
     [,1] [,2]
[1,]    6    3
[2,]    5    2
[3,]    4    1
> x[[2]]
[1] 5
> a <- 4
> i <- x[[a]]
> x[i]
[1] 4

Si l'on veut trouver toutes les occurrences d'une valeur, on utilise la fonction which() :

> x <- c(1, 2, 3, 2, 1)
> i <- which(x == 2)
> print(i)
[1] 2 4
> x[i]
[1] 2 2

Cette commande fonctionne avec toute matrice booléenne, par exemple

>  x <- cbind(c(6, 5, 4), c(3, 2, 1))
> i <- which(x > 2)
> print(i)
[1] 1 2 3 4
> x[i]
[1] 6 5 4 3

On peut aussi utiliser directement une matrice booléenne pour extraire des éléments :

> x <- cbind(c(6, 5, 4), c(3, 2, 1))
> comparaison <- (x > 4)
> print(comparaison)
      [,1]  [,2]
[1,]  TRUE FALSE
[2,]  TRUE FALSE
[3,] FALSE FALSE
> x[comparaison]
[1] 6 5

Lecture et enregistrement de données

La commande scan() permet d'aller lire les valeurs dans un fichier texte. Si par exemple les données sont dans le fichier valeurs.txt sous la forme :

(une donnée par ligne), il faut entrer

a <- scan(file = "valeurs.txt")

mais si elles sont sous la forme

1, 2, 3, 4

il faut indiquer que le séparateur est la virgule

a <- scan(file = "valeurs.txt", sep=",")

Si l'on indique un chemin d'accès (répertoire, dossier), il faut utiliser la barre oblique (slash) / même sous Microsoft Windows, par exemple

scan(file = "C:/Documents and Settings/foo.bar/valeurs.txt", sep = ",")

Si l'on a un fichier au format CSV (comma separated values), par exemple exporté depuis un tableur, on peut utiliser la commande

a <- read.csv(file = "valeurs.csv")

Bonnes pratiques de programmation

Nous rappelons ci-dessous quelques règles de bonne pratique. Ces règles sont destinées à réduire le risque d'erreur. En particulier, une autre personne que le programmeur doit pouvoir lire le code et le comprendre, et le programmeur doit être capable de se relire plusieurs mois ou années après avoir écrit son code. Cela permet le déverminage (débogage, recherche et correction des erreurs) et l'évolution du programme.

Le code doit être le plus général possible, ce qui permet sa transposition à d'autres applications (copier-coller). En particulier, on évite les valeurs numériques « en dur », on définit à la place des constantes en début de programme.

Le code doit être fractionné : il comprend typiquement trois grandes parties :

Initialisation et définition des constantes.
Définition des fonctions.
Programme principal.

Le programme principal doit être très court, et consister essentiellement à des appels de fonctions.

Une fonction ne doit pas faire plus d'une page de long (25 lignes), et ne doit pas contenir plus d'une dizaine de variables locales ; si ce n'est pas le cas, c'est qu'elle peut probablement être divisée en plusieurs fonctions. En particulier, chaque fois qu'une portion de code est dupliquée, elle doit être remplacée par une fonction.

Et une ligne ne doit pas faire plus de la largeur d'une page (80 caractères), et ne doit contenir qu'une seule instruction.

Le code doit être indenté : des blancs en début de ligne marquent l'imbrication des commandes, en particulier à l'intérieur d'un bloc délimitant une fonction (function() {…}), de boucles (for, while, repeat) et des exécutions conditionnelles (if … else …).

Le code doit être abondamment commenté. En particulier, pour chaque fonction, on indique

à quoi elle sert (son but, son objectif, sa fonctionnalité) ;
quelles sont les données d'entrée (arguments), et en particulier leur format (réel, complexe, booléen, chaîne de caractères) ;
quelles sont les données de sortie (idem).

Les grandes étapes de la fonction doivent être expliquées.

Exemples

Courbe de von Koch

Courbe de Koch

Le programme suivant trace la courbe du flocon de Koch de manière récursive.

#============================================================================
# nom : von_koch.sce
# auteur : Christophe Dang Ngoc Chan
# date de création : 2013-07-25
# dates de modification : 
#   2013-07-29 : calcul de la rotation par multiplication matricielle 
#----------------------------------------------------------------------------
# version de R : 2.15.2
# extension requises : aucune
#----------------------------------------------------------------------------
# Objectif : trace la courbe du "flocon de neige" de von Koch
# Entrées : aucun (paramètres codés en dur)
# Sorties : fenêtre graphique avec une courbe
#============================================================================


# **************
# * Constantes *
# **************

n <- 6 
# nombre d'étapes ;
# 9 donne 262 145 points
# 6 étapes (4 097 points) sont suffisantes pour un bon rendu

matrice_rotation = matrix(c(0.5, sqrt(3)/2, -sqrt(3)/2, 0.5), 2, 2)
# rotation de +60°
l <- 1
# longueur du segment initial (unité arbitraire)


# *************
# * Fonctions *
# *************
  
vonkoch <- function (pA, pB, i) {
  # Objectif : transforme le segment [AB] en la ligne brisée [ACDEB]
  #   et trace cette ligne si elle a une "longueur unité"
  # Entrées : pA et pB sont des vecteurs de nombres (coordonnées des points)
  #   i est un entier qui décroît à chaque étape
  # Sorties : tracé graphique
  u = (pB - pA)/3 # tiers du vecteur AB
  v = matrice_rotation%*%u
  # vecteur tourné de +60°
  pC = pA + u
  pD = pC + v
  pE = pB - u ;

  # points de la ligne
  if (i == 1) {
    # tracé des segments les plus petits
    x = c(pA[1], pC[1], pD[1], pE[1], pB[1])
    y = c(pA[2], pC[2], pD[2], pE[2], pB[2])
    lines(x, y)
  }
  else {
    # appel récursif
    j <- i-1 ;
    vonkoch(pA, pC, j)
    vonkoch(pC, pD, j)
    vonkoch(pD, pE, j)
    vonkoch(pE, pB, j)
 }
}


# ***********************
# * Programme principal *
# ***********************
  
# Extrémités de la ligne de base
debut <- c(0, 0)
fin <- c(l, 0)

# Tracé de la courbe
plot(c(0, l/2, l), c(0, sqrt(3)*l/6, 0), pch = ".", asp = 1,
     axes = F, xlab = "", ylab = "")
vonkoch(debut, fin, n)