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
)
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
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
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.
Lors la présentation de lundi nous avons vu que plusieurs type de communes existent :
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
)
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 :
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.
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
}
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
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.