Les petits bouts de codes suivants ont pour but de développer plusieurs points abordés au cours des présentations de lundi après-midi (organisation administrative française) et vendredi matin (postgis) :

Les données utilisées se trouvent dans le répertoire csv de ce répertoire :

qtlldir <- "/chemin/vers/le/répertoire/postgis/données/" 
setwd(qtlldir)

Commençons par importer le fichier des codes communes de l’Insee qui servira d’exemple :

##
commune2021 <- read.table(
  paste(
    "csv"
    , "commune2021.csv"
    , sep="/"
  )
  , sep = ","
  , header = T
  , quote = ''
  , stringsAsFactors = F
)

Génération d’une variable indicatrice des collectivités de Métropole

Reprenons l’exemple présenté lundi. On souhaite générer une variable indicatrice des DOM-COM d’une part et de la Métropole d’autre part en traitant les codes département comme des nombres malgré le fait que, depuis 1975, le code des départements est alphanumérique.

##
dcom <- as.integer(commune2021$DEP)>970
## Warning: NAs introduits lors de la conversion automatique

R nous avertit de l’échec de la conversion de certains codes. En effet,

##
table(dcom,useNA='always')
## dcom
## FALSE  TRUE  <NA> 
## 34521   129  3092

Les NA proviennent de la Corse mais aussi de communes pour lesquelles le département est manquant (il s’agit dans les faits de communes déléguées —voir après) :

##
table(commune2021$DEP[ is.na(dcom)])
## 
##        2A   2B 
## 2732  124  236

Donc, si on prend le complément de dcom :

##
table(!dcom,useNA='always')
## 
## FALSE  TRUE  <NA> 
##   129 34521  3092

3092 commnes sont manquantes. On pourrait évidemment combiner is.na() avec !dcom mais cela devient problématique si des valeurs étaient déjà NA avant la conversion qui se verraient donc située en Métropole alors que cette information n’est pas disponible. Il se trouve que, en l’espèce, tous les NA sont en Métropole mais ce n’est pas nécessairement le cas.

Une autre approche utilise le fait que les caractères ont un ordre

##
'2A' < '97'
## [1] TRUE

Notons toutefois que

##
'A' < '7'
## [1] FALSE

ce qui n’est pas un problème ici car le tri se fait de gauche à droite et

##
'2' < '9'
## [1] TRUE

Ce faisant,

##
dcom <- commune2021$DEP>'970'
table(dcom,useNA='always')
## dcom
## FALSE  TRUE  <NA> 
## 37613   129     0

Le problème est maintenant que

##
'' < '970'
## [1] TRUE

Donc les valeurs manquantes sont encore attribuées à la Métropole. Il faut donc remplacer les chaînes vides par NA (ce qui possible à l’importation en passant l’option na.strings ='' en argument).

##
commune2021 <- read.table(
  paste(
    "csv"
    , "commune2021.csv"
    , sep="/"
  )
  , sep = ","
  , header = T
  , quote = ''
  , na.strings =''
  , stringsAsFactors = F
)

On obtient enfin le résultat souhaité :

##
dcom <- commune2021$DEP>'970'
table(dcom,useNA='always')
## dcom
## FALSE  TRUE  <NA> 
## 34881   129  2732
table(!dcom,useNA='always')
## 
## FALSE  TRUE  <NA> 
##   129 34881  2732

Jointures

Opérations ensemblistes avec R

Commençons par générer des données

x <- data.frame(oid=1:6,lib='A')
y <- data.frame(oid=4:8,lib='B')

La colonne oid stocke la clefs du data.frame et c’est colonne qui servira à joindre les deux tables.

Appliquons maintenant les opérations ensemblistes vues dans la présentation aux clefs d’appariement :

Union des deux ensembles (toutes les clefs (univers)) :

## 
union(x$oid,y$oid)
## [1] 1 2 3 4 5 6 7 8

Intersection des deux ensembles (clefs en commun) :

intersect(x$oid,y$oid)
## [1] 4 5 6

Différence des deux ensembles :

## clefs de x absentes de y
setdiff(x$oid,y$oid)
## [1] 1 2 3
## clefs de y absentes de x
setdiff(y$oid,x$oid)
## [1] 7 8

On peut aussi calculer la différence symétrique (clefs qui sont soit absentes de x ou soit absentes de y càd le complément de l’union des clefs) :

## les deux expressions sont strictement équivalentes
symdiff <- function( x, y){
  sort(union(setdiff(x, y), setdiff( y, x)))
  ## ou : sort(setdiff( union(x, y), intersect(x, y)))
}
symdiff(y$oid,x$oid)
## [1] 1 2 3 7 8

Exemples des différents types de jointures

Join :

## 
merge(x,y, by='oid')
##   oid lib.x lib.y
## 1   4     A     B
## 2   5     A     B
## 3   6     A     B

Comme on peut le voir, R attribue une valeur manquante aux valeurs des clefs de x absentes de y.

Équivalent SQL :

select a.oid as xoid, a.lib as xlib, b.oid as yoid, b.lib as ylib
from x as a join y as b
on a.oid and b.oid
;

Left join :

##
merge(x,y, by='oid', all.x = T)
##   oid lib.x lib.y
## 1   1     A  <NA>
## 2   2     A  <NA>
## 3   3     A  <NA>
## 4   4     A     B
## 5   5     A     B
## 6   6     A     B

Équivalent SQL :

select a.oid as xoid, a.lib as xlib, b.oid as yoid, b.lib as ylib
from x as a left join y as b
on a.oid and b.oid
;

Right join :

## 
merge(x,y, by='oid', all.y = T)
##   oid lib.x lib.y
## 1   4     A     B
## 2   5     A     B
## 3   6     A     B
## 4   7  <NA>     B
## 5   8  <NA>     B

Équivalent SQL :

select a.oid as xoid, a.lib as xlib, b.oid as yoid, b.lib as ylib
from x as a right join y as b
on a.oid and b.oid
;

Full outer join :

## 
merge(x,y, by='oid', all = T)
##   oid lib.x lib.y
## 1   1     A  <NA>
## 2   2     A  <NA>
## 3   3     A  <NA>
## 4   4     A     B
## 5   5     A     B
## 6   6     A     B
## 7   7  <NA>     B
## 8   8  <NA>     B

Équivalent SQL :

select a.oid as xoid, a.lib as xlib, b.oid as yoid, b.lib as ylib
from x as a full outer join y as b
on a.oid and b.oid
;

Produit cartésien sur les clefs (toutes les combinaisons deux-à-deux des clefs)

## 
expand.grid(x$oid,y$oid)
##    Var1 Var2
## 1     1    4
## 2     2    4
## 3     3    4
## 4     4    4
## 5     5    4
## 6     6    4
## 7     1    5
## 8     2    5
## 9     3    5
## 10    4    5
## 11    5    5
## 12    6    5
## 13    1    6
## 14    2    6
## 15    3    6
## 16    4    6
## 17    5    6
## 18    6    6
## 19    1    7
## 20    2    7
## 21    3    7
## 22    4    7
## 23    5    7
## 24    6    7
## 25    1    8
## 26    2    8
## 27    3    8
## 28    4    8
## 29    5    8
## 30    6    8

Équivalent SQL :

select a.oid as xoid, a.lib as xlib, b.oid as yoid, b.lib as ylib
from x as a, y as b
;

Le nombre total de combinaisons vaut

nrow(x)*nrow(y)
## [1] 30

Ce nombre grandit donc très rapidement et conduit vite à des tables de grande taille.

Application au Code officiel géographique

Lors la présentation de lundi nous avons vu que plusieurs type de communes existent :

  • les associées et déléguées
  • les communes à arrondissements

Voyons maintenant comment cela se traduit concrètement dans le Code officiel géographique. Les fichiers Intercommunalite-Metropole_au_01-01-2021.csv et Intercommunalite-Metropole_au_01-01-2021–comp.csv du répertoire csv renseignent respectivement les EPCI et leur composition :

## Fichier des EPCI
epci2021 <- read.table(  
  paste(
    "csv"
    , "Intercommunalite-Metropole_au_01-01-2021.csv"
    , sep="/"
  )
  , sep=";"
  , header = T
  , stringsAsFactors = F
)
## Fichier de la composition des EPCI
epci2021.comp <- read.table(  
  paste(
    "csv"
    , "Intercommunalite-Metropole_au_01-01-2021--comp.csv"
    , sep="/"
  )
  , sep=";"
  , header = T
  , stringsAsFactors = F
)

Appariement des codes communes

Comparons maintenant les clefs des deux tables EPCI,

##
symdiff(epci2021$EPCI, epci2021.comp$EPCI)
## character(0)

on a bien les mêmes clefs dans les deux fichiers.

Comparons ensuite le fichier de composition des communes avec le fichier commune2021 :

##
nrow(epci2021.comp) - nrow(commune2021)
## [1] -2777

On a ici une première indication de différence entre les deux sources car la table commune2021 comporte 2777 communes en plus.

##
length(symdiff(epci2021.comp$CODGEO, commune2021$COM))
## [1] 2157
length(setdiff(epci2021.comp$CODGEO, commune2021$COM))
## [1] 0
length(setdiff(commune2021$COM,epci2021.comp$CODGEO))
## [1] 2157

Il manque donc 2157 clefs uniques dans le fichier commune2021. Par contre, aucune clefs ne manque dans le fichier EPCI.

Pour comprendre ce qui se passe, réalisons une jointure à gauche des deux fichiers (une jointure complète donnerait le même résultat) :

## 
commepci2021 <- merge(commune2021, epci2021.comp, by.x= 'COM', by.y='CODGEO', all.x = T, suffixes = c(".comm", ".epci"))
##
View(commepci2021)

Comme vu précédemment, les valeurs des communes présentes dans commune2021 mais absentes de epci2021.comp valent NA. On peut donc s’en servir comme indicatrice de la différence A - B :

## 
table( is.na(commepci2021$EPCI) )
## 
## FALSE  TRUE 
## 35585  2157

On retrouve le nombre de clefs absentes dans epci2021.comp. Du fait qu’il y a plusieurs façon de délimiter les communes, voyons si les différences proviennent de là :

## 
(frq<-table(commepci2021$TYPECOM, is.na(commepci2021$EPCI)))
##       
##        FALSE  TRUE
##   ARM      0    45
##   COM  34965     0
##   COMA     0   517
##   COMD   620  1595

Note :

  • COM : communes
  • COMA : communes associées
  • COMD : communes déléguées
  • ARM : arrondissements

En sommant la second colonne de frq

## 
sum(frq[,2])
## [1] 2157

on retrouve bien de clefs manquantes dans epci2021.comp. Les ( 2777 - 2157 ) = 620 lignes encore manquantes s’expliquent par la répétition du code commune pour certaines communes déléguées qui ont repris le code de la commune avec laquelle elles ont été fusionnées :

## 
table( table(commune2021$COM) >1)
## 
## FALSE  TRUE 
## 36502   620

La différence du nombre de ligne entre les deux tables s’explique donc bien par l’absence des arrondissements, des communes associées et des communes déléguées dans la table EPCI .

Comme souligné lors de la présentation de lundi, il faut donc toujours faire attention au découpage dans lequel les informations à la commune sont diffusées.

Fonctions de hachage

Cette partie vise à développer quelques aspects plus avancés des fonctions de hachage non abordés dans la partie cours et à illustrer leur fonctionnement en pratique.

Comme expliqué dans la présentation, les fonctions de hachage ont pour objet de transformer une séquence binaire en nombre dans une plage prédéterminée pour s’en servir d’index. L’exemple de fonction utilisé est le suivant :

h <- function(x, sz, init=5381, M=33){
  mod <- 2^32
  h = init
  for( l in  as.integer(charToRaw(x)) ){
    h <- ( (M * h)%%mod  + l )%%mod
  }
  as.integer(h %% sz)+1L
}

Remarques à propos de la fonction

Cette fonction de hachage est un classique du genre. Elle est notamment connue pour apparaître dans le Red Dragon Book qui est un incunable de la littérature sur la conception des compilateurs.

Les paramètres par défaut sont la version de Bernstein mais d’autres choix sont possibles comme celle de Kernighan & Ritchie, deux des concepteurs du langage C, (plus loin). Et, de leur aveu même, “This is not the best possible algorithm, but it has the merit of extreme simplicity”.

Au de-là de son caractère vénérable, elle a été avant tout choisie car elle peut être facilement implémentée avec R. De nombreuses autres fonctions existent ayant de bien meilleures propriétés mais elles auraient été beaucoup plus diffciles à implémenter ici.

Comme une fonction de hachage vise en transformer une séquence en nombre (quel que soit son type : chaîne de caractère mais aussi entier, floatant,…), l’appel de fonction transforme donc d’abord la chaîne de caractères en un vecteur d’octet qui n’est autre que la représentation sous-jacente de la chaîne de caractère en mémoire :

as.integer(charToRaw("L'Abergement-Clémenciat"))
##  [1]  76  39  65  98 101 114 103 101 109 101 110 116  45  67 108 195 169 109 101 110  99 105  97 116

Notez que le résultat est susceptible de varier en fonction de l’encodage de la chaîne de caractère. Ce vecteur est celui que l’on obtient dans une session en UTF-8. Si on convertit la chaîne en WINDOWS-1252,

as.integer(charToRaw(iconv("L'Abergement-Clémenciat", to="WINDOWS-1252")))
##  [1]  76  39  65  98 101 114 103 101 109 101 110 116  45  67 108 233 109 101 110  99 105  97 116

le vecteur est plus court et le seizième octet est différent. Dans les deux jeux de caractères, “é” a le même point de code U+E9 (233) mais, en UTF-8, il est encodé avec deux octets ({0xC3, OxA9}) là où le code page Windows 1252 l’encode directement avec son point de code. Pour plus d’informations sur l’UTF-8, voir cette page ainsi que celle-ci sur encodage des jeux de caractères plus généralement.

En conséquence de quoi, les hachis seront différents en fonction de l’encodage utilisé ce qui ne pose pas de problèmes sauf si on utilise un encodage pour l’indexation et un autre pour la recherche.

On réalise de plus les calculs avec des flottants modulo (%%) 2^32 pour émuler les calculs avec des entiers non signé 32 bits. En effet, R ne propose que des entiers signés 32 bits qui donnerait des résultats différents. En outre, le débordement par le haut qui est normal ici conduirait à la production de NA ce qui n’est clairement pas le résultat attendu.

## débordement d'entier par le haut
print(as.integer(2L^31 - 1L)) +1L
## [1] 2147483647
## Warning in print(as.integer(2L^31 - 1L)) + 1L: NA produit par débordement d'entier par le haut
## [1] NA

Application de la fonction

Voyons maintenant ce que ça donne avec la première commune du Code officiel géographique dans l’ordre lexicographique :

(i<-h("L'Abergement-Clémenciat", 37742))
## [1] 24608

La fonction retourne 24608. La valeur (“01001”) sera donc stockée à la position 24608.

htab <- integer(37742)
htab[i] <- "01001"

Pour la retrouver, il suffit donc de hacher à nouveau la clef.

htab[h("L'Abergement-Clémenciat", 37742)]
## [1] "01001"

Mais ce serait sans compter sur le fait que, selon le théorème des tiroirs (ou le paradoxe des anniversaires), les clefs ne sont malheureusement pas nécessairement uniques. Ce qui peut être illustré avec une variante de la fonction, la version Kernighan & Ritchie :

hcom <- sapply(commune2021$LIBELLE, function(x)h(x, 37742, 0, 31))

Pas besoin de faire un test du Khi^2 pour voir que la distribution est n’est pas tout à fait uniforme.

(frq <- table(table(hcom)))
## 
##     1     2     3     4     5     6     7     8     9    10    11    12    13    14 
## 13329  6245  2277   671   230    84    42    21    14     6     5     2     1     1

De plus, des collisions apparaissent pour plus des 2/3 des clefs et .39 des indices ne sont pas occupées.

nrow(commune2021)-frq[1]
##     1 
## 24413
prop.table(table(match(1:nrow(commune2021), hcom,0)!=0))
## 
##    FALSE     TRUE 
## 0.392507 0.607493

En conséquence de quoi, L’Abergement-Clémenciat à la même indice que Gueugnon

which(hcom==hcom[1])
## L'Abergement-Clémenciat                Gueugnon 
##                       1                   30298

Comme la chaïne “Gueugnon” a été hachées après celle de “L’Abergement-Clémenciat”, la recherche

htab[h("L'Abergement-Clémenciat", 37742)]

retourne donc “71230” au lieu de “01001”.

Pour réduire le nombre de clefs ayant le même indice, on peut commencer par doubler la capacité de la table,

hcom <- sapply(commune2021$LIBELLE, function(x)h(x, 37742*2, 0, 31))
(frq <- table(table(hcom)))
## 
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 21150  5353  1166   278   111    38    24    19     7     5     2     1     2
nrow(commune2021)-frq[1]
##     1 
## 16592
prop.table(table(match(1:(nrow(commune2021)*2),hcom,0)!=0))
## 
##     FALSE      TRUE 
## 0.6269938 0.3730062

Le nombre de collisions se réduit mais reste conséquent et la proportion d’indices inutilisés augmente (.6).

Il conviendra donc d’utiliser une autre fonction tout en gardant à l’esprit qu’elle ne réduira pas la probabilité de collision à zéro (sauf cas particuiers) et qu’il faudra donc de toute façon rajouter une étape pour les gérer.