...
[iramuteq] / Rscripts / chdquest.R
1 #Author: Pierre Ratinaud
2 #Copyright (c) 2008-2009 Pierre Ratinaud
3 #License: GNU/GPL
4
5 fille<-function(classe,classeuce) {
6         listfm<-unique(unlist(classeuce[classeuce[,classe%/%2]==classe,]))
7         listf<-listfm[listfm>=classe]
8         listf<-unique(listf)
9         listf
10 }
11
12
13 croiseeff <- function(croise, classeuce1, classeuce2) {
14     cl1 <- 0
15     cl2 <- 1
16     for (i in 1:ncol(classeuce1)) {
17         cl1 <- cl1 + 2
18         cl2 <- cl2 + 2
19         clj1 <- 0
20         clj2 <- 1
21         for (j in 1:ncol(classeuce2)) {
22             clj1 <- clj1 + 2
23             clj2 <- clj2 + 2
24             croise[cl1 - 1, clj1 -1] <- length(which(classeuce1[,i] == cl1 & classeuce2[,j] == clj1))
25             croise[cl1 - 1, clj2 -1] <- length(which(classeuce1[,i] == cl1 & classeuce2[,j] == clj2))
26             croise[cl2 - 1, clj1 -1] <- length(which(classeuce1[,i] == cl2 & classeuce2[,j] == clj1))
27             croise[cl2 - 1, clj2 -1] <- length(which(classeuce1[,i] == cl2 & classeuce2[,j] == clj2))
28         }
29     }
30     croise
31 }
32
33
34 #fonction pour la double classification
35 #cette fonction doit etre splitter en 4 ou 5 fonctions
36
37 Rchdquest<-function(tableuc1,listeuce1,uceout ,nbt = 9, mincl = 2, mode.patate = FALSE, svd.method = 'irlba') {
38         #source('/home/pierre/workspace/iramuteq/Rscripts/CHD.R')
39     if (svd.method == 'irlba') {
40         library(irlba)
41     }
42         #lecture des tableaux
43         data1<-read.csv2(tableuc1)#,row.names=1)
44     cn.data1 <- colnames(data1)
45     data1 <- as.matrix(data1)
46     colnames(data1) <- cn.data1
47     rownames(data1) <- 1:nrow(data1)
48         data2<-data1
49         sc<-colSums(data2)
50         if (min(sc)<=4){
51                 data1<-data1[,-which(sc<=4)]
52                 sc<-sc[-which(sc<=4)]
53         }
54         #analyse des tableaux avec la fonction CHD qui doit etre sourcee avant
55         chd1<-CHD(data1, x = nbt, mode.patate = mode.patate, svd.method)
56         chd2<-chd1
57
58         #FIXME: le nombre de classe peut etre inferieur
59     nbcl <- nbt + 1
60     tcl <- ((nbt+1) * 2) - 2
61
62         #lecture des uce
63         listuce1<-read.csv2(listeuce1)
64         listuce2<-listuce1
65
66         #Une fonction pour assigner une classe a chaque UCE
67 #       AssignClasseToUce<-function(listuce,chd) {
68 #           out<-matrix(nrow=nrow(listuce),ncol=ncol(chd))
69 #           for (i in 1:nrow(listuce)) {
70 #                       for (j in 1:ncol(chd)) {
71 #                           out[i,j]<-chd[(listuce[i,2]+1),j]
72 #                       }
73 #           }
74 #           out
75 #       }
76
77     AssignClasseToUce <- function(listuce, chd) {
78         print('assigne classe -> uce')
79         chd[listuce[,2]+1,]
80     }
81         #Assignation des classes
82         classeuce1<-AssignClasseToUce(listuce1,chd1$n1)
83         classeuce2<-classeuce1
84         
85         #calcul des poids (effectifs)
86         poids1<-vector(mode='integer',length=tcl)
87 #       makepoids<-function(classeuce,poids) {
88 #           for (classes in 2:(tcl + 1)){
89 #               for (i in 1:ncol(classeuce)) {
90 #                   if (poids[(classes-1)]<length(classeuce[,i][classeuce[,i]==classes])) {
91 #                       poids[(classes-1)]<-length(classeuce[,i][classeuce[,i]==classes])
92 #                   }
93 #               }
94 #           }
95 #           poids
96 #       }
97     makepoids <- function(classeuce, poids) {
98         cl1 <- 0
99         cl2 <- 1
100         for (i in 1:nbt) {
101             cl1 <- cl1 + 2
102             cl2 <- cl2 + 2
103             poids[cl1 - 1] <- length(which(classeuce[,i] == cl1))
104             poids[cl2 - 1] <- length(which(classeuce[,i] == cl2))
105         }
106         poids
107     }
108     
109         poids1<-makepoids(classeuce1,poids1)
110         poids2<-poids1
111
112 #       croise=matrix(ncol=tcl,nrow=tcl)
113 #       #production du tableau de contingence
114 #       for (i in 1:ncol(classeuce1)) {
115 #           #poids[i]<-length(classeuce1[,i][x==classes])
116 #           for (j in 1:ncol(classeuce2)) {
117 #               tablecroise<-table(classeuce1[,i],classeuce2[,j])
118 #               tabcolnames<-as.numeric(colnames(tablecroise))
119 #               #tabcolnames<-c(tabcolnames[(length(tabcolnames)-1)],tabcolnames[length(tabcolnames)])
120 #               tabrownames<-as.numeric(rownames(tablecroise))
121 #               #tabrownames<-c(tabrownames[(length(tabrownames)-1)],tabrownames[length(tabrownames)])
122 #               for (k in (ncol(tablecroise)-1):ncol(tablecroise)) {
123 #                   for (l in (nrow(tablecroise)-1):nrow(tablecroise)) {
124 #                       croise[(tabrownames[l]-1),(tabcolnames[k]-1)]<-tablecroise[l,k]
125 #                   }
126 #               }
127 #           }
128 #           tablecroise
129 #       }
130     croise <- croiseeff( matrix(ncol=tcl,nrow=tcl), classeuce1, classeuce2)
131     if (mincl == 2) {
132             mincl<-round(nrow(classeuce1)/(nbt+1)) #valeur a calculer nbuce/nbt
133     }
134     #print('ATTENTION MINCL IMPOSE')
135         #mincl<-422
136         print('mincl')
137         print(mincl)
138         if (mincl < 3) {
139             mincl<-3
140         }
141     #print('ATTENTION : ON IMPOSE LA TAILLE DES CLASSES')
142     #mincl <- 15
143         print('table1')
144         print(croise)
145         #tableau des chi2 signes
146         chicroise<-croise
147         for (i in 1:nrow(croise)) {
148             for (j in 1:ncol(croise)) {
149                 if (croise[i,j]==0) {
150                     chicroise[i,j]<-0
151                 } else if (croise[i,j]<mincl) { 
152                     chicroise[i,j]<-0
153                 } else {
154                         chitable<-matrix(ncol=2,nrow=2)
155                         chitable[1,1]<-croise[i,j]
156                         chitable[1,2]<-poids1[i]-chitable[1,1]
157                         chitable[2,1]<-poids2[j]-chitable[1,1]
158                         chitable[2,2]<-nrow(classeuce1)-poids2[j]-chitable[1,2]
159                         chitest<-chisq.test(chitable,correct=FALSE)
160                         if ((chitable[1,1]-chitest$expected)<0) {
161                             chicroise[i,j]<--round(chitest$statistic,digits=7)
162                         } else {
163                             chicroise[i,j]<-round(chitest$statistic,digits=7)
164                 #print(chitest)
165                         }
166                 }
167             }   
168         }
169         print(chicroise)
170         #determination des chi2 les plus fort
171         chicroiseori<-chicroise
172     doxy <- function(chicroise) {
173         listx <- NULL
174         listy <- NULL
175         listxy <- which(chicroise > 3.84, arr.ind = TRUE)
176         #print(listxy)
177         val <- chicroise[which(chicroise > 3.84)]
178         ord <- order(val, decreasing = TRUE)
179         listxy <- listxy[ord,]
180         #for (i in 1:nrow(listxy)) {
181         #    if ((!listxy[,2][i] %in% listx) & (!listxy[,1][i] %in% listy)) {
182         #        listx <- c(listx, listxy[,2][i])
183         #        listy <- c(listy, listxy[,1][i])
184         #    }
185         #}
186         xy <- list(x = listxy[,2], y = listxy[,1])
187         xy
188     }
189     xy <- doxy(chicroise)
190     print(xy)
191     listx <- xy$x
192     listy <- xy$y
193     
194 #       maxi<-vector()
195 #       chimax<-vector()
196 #       for (i in 1:tcl) {
197 #           maxi[i]<-which.max(chicroise)
198 #           chimax[i]<-chicroise[maxi[i]]
199 #           chicroise[maxi[i]]<-0
200 #       }
201 #       testpres<-function(x,listcoord) {
202 #           for (i in 1:length(listcoord)) {
203 #               if (x==listcoord[i]) {
204 #                   return(-1)
205 #               } else {
206 #                   a<-1
207 #               }
208 #           }
209 #           a
210 #       }
211 #       c.len=nrow(chicroise)
212 #       r.len=ncol(chicroise)
213 #       listx<-c(0)
214 #       listy<-c(0)
215 #       rang<-0
216 #       cons<-list()
217 #       #on garde une valeur par ligne / colonne
218 #       for (i in 1:length(maxi)) {
219 #       #coordonnées de chi2 max
220 #           x.co<-ceiling(maxi[i]/c.len)
221 #           y.co<-maxi[i]-(x.co-1)*c.len
222 #           a<-testpres(x.co,listx)
223 #           b<-testpres(y.co,listy)
224 #           
225 #           if (a==1) {
226 #               if (b==1) {
227 #                   rang<-rang+1
228 #                   listx[rang]<-x.co
229 #                   listy[rang]<-y.co
230 #               }
231 #           }
232 #           cons[[1]]<-listx
233 #           cons[[2]]<-listy
234 #       }
235         #pour ecrire les resultats
236         for (i in 1:length(listx)) {
237             txt<-paste(listx[i]+1,listy[i]+1,sep=' ')
238             txt<-paste(txt,croise[listy[i],listx[i]],sep=' ')
239             txt<-paste(txt,chicroiseori[listy[i],listx[i]],sep=' ')
240             print(txt)
241         }
242
243         #colonne de la classe
244         #trouver les filles et les meres 
245         trouvefillemere<-function(classe,chd) {
246             unique(unlist(chd[chd[,classe%/%2]==classe,]))
247         }
248
249         #pour trouver une valeur dans une liste
250         #is.element(elem, list)
251         #== elem%in%list
252
253         coordok<-NULL
254         trouvecoordok<-function(first) {
255             fillemere1<-NULL
256             fillemere2<-NULL
257             listxp<-listx
258             listyp<-listy
259             listxp<-listx[first:length(listx)]
260             listxp<-c(listxp,listx[1:(first-1)])
261         #    listxp<-listxp[-first]
262             listyp<-listy[first:length(listy)]
263             listyp<-c(listyp,listy[1:(first-1)])
264         #    listyp<-listyp[-first]
265             for (i in 1:length(listxp)) {
266                if (!(listxp[i]+1)%in%fillemere1) {
267                   if (!(listyp[i]+1)%in%fillemere2) {
268                   coordok<-rbind(coordok,c(listyp[i]+1,listxp[i]+1))
269                   fillemere1<-c(fillemere1,trouvefillemere(listxp[i]+1,chd2$n1))
270                   fillemere2<-c(fillemere2,trouvefillemere(listyp[i]+1,chd1$n1))
271                   }
272                }
273             }
274             coordok
275         }
276 #fonction pour trouver le nombre maximum de classes
277         findmaxclasse<-function(listx,listy) {
278             listcoordok<-list()
279             maxcl<-0
280             nb<-1
281             for (i in 1:length(listy)) {
282
283                 coordok<-trouvecoordok(i)
284                 if (maxcl <= nrow(coordok)) {
285                     maxcl<-nrow(coordok)
286                     listcoordok[[nb]]<-coordok
287                     nb<-nb+1
288                 }
289             }
290             listcoordok<-unique(listcoordok)
291                 print('liste coord ok')
292 #               print('FIXME FIXME FIXME FIXME FIXME')
293 #               print(listcoordok)
294                 #si plusieurs ensemble avec le meme nombre de classe, on conserve
295                 #la liste avec le plus fort chi2
296             if (length(listcoordok)>1) {
297                 maxchi<-0
298                 best<-NULL
299                 for (i in 1:length(listcoordok)) {
300                     chi<-NULL
301                     uce<-NULL
302                     if (nrow(listcoordok[[i]])==maxcl) {
303                         for (j in 1:nrow(listcoordok[[i]])) {
304                             chi<-c(chi,croise[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
305                             uce<-c(uce,chicroiseori[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
306                         }
307                                 print(sum(uce))
308                         if (maxchi < sum(chi)) {
309                             maxchi<-sum(chi)
310                             suce<-sum(uce)
311                             best<-i
312                         } 
313                     }
314                 }
315             }
316             print((suce/nrow(classeuce1)*100))
317             listcoordok[[best]]
318         }
319         #findmaxclasse(listx,listy)
320         #coordok<-trouvecoordok(1)
321         coordok<-findmaxclasse(listx,listy)
322         print(coordok)
323
324         fille<-function(classe,classeuce) {
325             listfm<-unique(unlist(classeuce[classeuce[,classe%/%2]==classe,]))
326             listf<-listfm[listfm>=classe]
327             listf<-unique(listf)
328             listf
329         }
330
331
332         lfilletot<-function(classeuce) {
333             listfille<-NULL
334             for (classe in 1:nrow(coordok)) {
335                 listfille<-unique(c(listfille,fille(coordok[classe,1],classeuce)))
336                 listfille
337             }
338         }
339
340         listfille1<-lfilletot(classeuce1)
341         listfille2<-lfilletot(classeuce2)
342
343
344         #utiliser rownames comme coordonnees dans un tableau de 0
345         Assignclasse<-function(classeuce,x) {
346             nchd<-matrix(0,ncol=ncol(classeuce),nrow=nrow(classeuce))
347             for (classe in 1:nrow(coordok)) {
348                 
349                 clnb<-coordok[classe,x]
350                 colnb<-clnb%/%2
351                 tochange<-which(classeuce[,colnb]==clnb)
352             for (row in 1:length(tochange)) {
353                     nchd[tochange[row],colnb:ncol(nchd)]<-classe
354                 }
355             }
356             nchd
357         }
358         print('commence assigne new classe')
359         nchd1<-Assignclasse(classeuce1,1)
360         #nchd1<-Assignnewclasse(classeuce1,1)
361         nchd2<-Assignclasse(classeuce2,2)
362         print('fini assign new classe')
363         croisep<-matrix(ncol=nrow(coordok),nrow=nrow(coordok))
364         for (i in 1:nrow(nchd1)) {
365             if (nchd1[i,ncol(nchd1)]==0) {
366                 nchd2[i,]<-nchd2[i,]*0
367             }
368             if (nchd1[i,ncol(nchd1)]!=nchd2[i,ncol(nchd2)]) {
369                 nchd2[i,]<-nchd2[i,]*0
370             }
371             if (nchd2[i,ncol(nchd2)]==0) {
372                 nchd1[i,]<-nchd1[i,]*0
373             }
374         }
375         print('fini croise')
376         elim<-which(nchd1[,ncol(nchd1)]==0)
377         keep<-which(nchd1[,ncol(nchd1)]!=0)
378         n1<-nchd1[nchd1[,ncol(nchd1)]!=0,]
379         n2<-nchd2[nchd2[,ncol(nchd2)]!=0,]
380         print('debut graph')
381         clnb<-nrow(coordok)
382         print('fini')
383         write.csv2(nchd1[,ncol(nchd1)],uceout)
384         res <- list(n1 = nchd1, coord_ok = coordok, cuce1 = classeuce1, chd = chd1)
385         res
386 }
387 #n1<-Rchdtxt('/home/pierre/workspace/iramuteq/corpus/agir2sortie01.csv','/home/pierre/workspace/iramuteq/corpus/testuce.csv','/home/pierre/workspace/iramuteq/corpus/testuceout.csv')