memory
[iramuteq] / Rscripts / chdtxt.R
1 #Author: Pierre Ratinaud
2 #Copyright (c) 2008-2009 Pierre Ratinaud
3 #License: GNU/GPL
4
5
6 #fonction pour la double classification
7 #cette fonction doit etre splitter en 4 ou 5 fonctions
8
9 AssignClasseToUce <- function(listuce, chd) {
10     print('assigne classe -> uce')
11     chd[listuce[,2]+1,]
12 }
13
14 fille<-function(classe,classeuce) {
15         listfm<-unique(unlist(classeuce[classeuce[,classe%/%2]==classe,]))
16         listf<-listfm[listfm>=classe]
17         listf<-unique(listf)
18         listf
19 }
20
21
22 croiseeff <- function(croise, classeuce1, classeuce2) {
23     cl1 <- 0
24     cl2 <- 1
25     for (i in 1:ncol(classeuce1)) {
26         cl1 <- cl1 + 2
27         cl2 <- cl2 + 2
28         clj1 <- 0
29         clj2 <- 1
30         for (j in 1:ncol(classeuce2)) {
31             clj1 <- clj1 + 2
32             clj2 <- clj2 + 2
33             croise[cl1 - 1, clj1 -1] <- length(which(classeuce1[,i] == cl1 & classeuce2[,j] == clj1))
34             croise[cl1 - 1, clj2 -1] <- length(which(classeuce1[,i] == cl1 & classeuce2[,j] == clj2))
35             croise[cl2 - 1, clj1 -1] <- length(which(classeuce1[,i] == cl2 & classeuce2[,j] == clj1))
36             croise[cl2 - 1, clj2 -1] <- length(which(classeuce1[,i] == cl2 & classeuce2[,j] == clj2))
37         }
38     }
39     croise
40 }
41
42 addallfille <- function(lf) {
43     nlf <- list()
44     for (i in 1:length(lf)) {
45         if (! is.null(lf[[i]])) {
46             nlf[[i]] <- lf[[i]]
47             filles <- lf[[i]]
48             f1 <- filles[1]
49             f2 <- filles[2]
50             if (f1 > length(lf)) {
51                 for (j in (length(lf) + 1):f2) {
52                     nlf[[j]] <- 0
53                 }
54             }
55         } else {
56             nlf[[i]] <- 0
57         }
58     }
59 nlf
60 }
61
62 getfille <- function(nlf, classe, pf) {
63     if (length(nlf[[classe]]) == 1) {
64         return(pf)
65     } else {
66         pf <- c(pf, nlf[[classe]])
67         c1 <- nlf[[classe]][1]
68         c2 <- nlf[[classe]][2]
69         pf <- getfille(nlf, c1, pf)
70         pf <- getfille(nlf, c2, pf)
71     }
72     return(pf)
73 }
74
75 getmere <- function(list_mere, classe) {
76     i <- classe
77     pf <- NULL
78     while (i != 1 ) {
79         pf <- c(pf, list_mere[[i]])
80         i <- list_mere[[i]]
81     }
82     pf
83 }
84
85 getfillemere <- function(list_fille, list_mere, classe) {
86     return(c(getfille(list_fille, classe, NULL), getmere(list_mere, classe)))
87 }
88
89 #nbt nbcl = nbt+1 tcl=((nbt+1) *2) - 2  n1[,ncol(n1)], nchd1[,ncol(nchd1)]
90 Rchdtxt<-function(uceout, chd1, chd2 = NULL, mincl=0, classif_mode=0, nbt = 9) {
91         #FIXME: le nombre de classe peut etre inferieur
92     nbcl <- nbt + 1
93     tcl <- ((nbt+1) * 2) - 2
94     #Assignation des classes
95         classeuce1<-AssignClasseToUce(listuce1,chd1$n1)
96         if (classif_mode==0) {
97                 classeuce2<-AssignClasseToUce(listuce2,chd2$n1)
98         }
99         #} else {
100         #       classeuce2<-classeuce1
101     #}
102
103         #calcul des poids (effectifs)
104
105     makepoids <- function(classeuce, poids) {
106         cl1 <- 0
107         cl2 <- 1
108         for (i in 1:nbt) {
109             cl1 <- cl1 + 2
110             cl2 <- cl2 + 2
111             poids[cl1 - 1] <- length(which(classeuce[,i] == cl1))
112             poids[cl2 - 1] <- length(which(classeuce[,i] == cl2))
113         }
114         poids
115     }
116
117 #       makepoids<-function(classeuce,poids) {
118 #           for (classes in 2:(tcl + 1)){
119 #                   for (i in 1:ncol(classeuce)) {
120 #                       if (poids[(classes-1)]<length(classeuce[,i][classeuce[,i]==classes])) {
121 #                           poids[(classes-1)]<-length(classeuce[,i][classeuce[,i]==classes])
122 #                       }
123 #                   }
124 #           }
125 #           poids
126 #       }
127     print('make poids')
128         poids1<-vector(mode='integer',length = tcl)
129         poids1<-makepoids(classeuce1,poids1)
130         if (classif_mode==0) {
131                 poids2<-vector(mode='integer',length = tcl)
132                 poids2<-makepoids(classeuce2,poids2)
133         }# else {
134         #       poids2<-poids1
135         #}
136     
137     print('croisement classif')
138
139 #    croise=matrix(ncol=tcl,nrow=tcl)
140 #
141 #    docroise <- function(croise, classeuce1, classeuce2) {
142 #       #production du tableau de contingence
143 #       for (i in 1:ncol(classeuce1)) {
144 #           #poids[i]<-length(classeuce1[,i][x==classes])
145 #           for (j in 1:ncol(classeuce2)) {
146 #                   tablecroise<-table(classeuce1[,i],classeuce2[,j])
147 #                   tabcolnames<-as.numeric(colnames(tablecroise))
148 #                   #tabcolnames<-c(tabcolnames[(length(tabcolnames)-1)],tabcolnames[length(tabcolnames)])
149 #                   tabrownames<-as.numeric(rownames(tablecroise))
150 #                   #tabrownames<-c(tabrownames[(length(tabrownames)-1)],tabrownames[length(tabrownames)])
151 #                   for (k in (ncol(tablecroise)-1):ncol(tablecroise)) {
152 #                       for (l in (nrow(tablecroise)-1):nrow(tablecroise)) {
153 #                           croise[(tabrownames[l]-1),(tabcolnames[k]-1)]<-tablecroise[l,k]
154 #                       }
155 #                   }
156 #           }
157 #       }
158 #        croise
159 #    }
160         if (classif_mode==0) {
161         croise <- croiseeff( matrix(ncol=tcl,nrow=tcl), classeuce1, classeuce2)
162         } else {
163                 croise <- croiseeff( matrix(ncol=tcl,nrow=tcl), classeuce1, classeuce1)
164         }
165     print(croise)
166     if (classif_mode == 0) {ind <- (nbcl * 2)} else {ind <- nbcl}
167         if (mincl==0){
168                 mincl<-round(nrow(classeuce1)/ind)
169         }
170         if (mincl<3){
171                 mincl<-3
172         }
173     print(mincl)        
174         #print('table1')
175         #print(croise)
176         #tableau des chi2 signes
177     print('croise chi2')
178         #chicroise<-croise
179
180 #    nr <- nrow(classeuce1)
181 #    newchicroise <- function(croise, mincl, nr, poids1, poids2) {
182 #        chicroise <- croise
183 #        chicroise[which(croise < mincl)] <- 0
184 #        tocompute <- which(chicroise > 0, arr.ind = TRUE)
185 #        for (i in 1:nrow(tocompute)) {
186 #            chitable <- matrix(ncol=2,nrow=2)
187 #            chitable[1,1] <- croise[tocompute[i,1],  tocompute[i,2]]
188 #            chitable[1,2] <- poids1[tocompute[i,1]] - chitable[1,1]
189 #            chitable[2,1] <- poids2[tocompute[i,2]] - chitable[1,1]
190 #            chitable[2,2] <- nr - poids2[tocompute[i,2]] - chitable[1,2]
191 #            chitest<-chisq.test(chitable,correct=FALSE)
192 #            chicroise[tocompute[i,1],  tocompute[i,2]] <- ifelse(chitable[1,1] > chitest$expected[1,1], round(chitest$statistic,digits=7), -round(chitest$statistic,digits=7))
193 #        }
194 #        chicroise
195 #    }
196 #
197         
198
199         dochicroise <- function(croise, mincl) {
200         chicroise <- croise
201         for (i in 1:nrow(croise)) {
202                 for (j in 1:ncol(croise)) {
203                     if (croise[i,j]==0) {
204                         chicroise[i,j]<-0
205                     } else if (croise[i,j]<mincl) { 
206                         chicroise[i,j]<-0
207                     } else {
208                         chitable<-matrix(ncol=2,nrow=2)
209                         chitable[1,1]<-croise[i,j]
210                         chitable[1,2]<-poids1[i]-chitable[1,1]
211                         chitable[2,1]<-poids2[j]-chitable[1,1]
212                         chitable[2,2]<-nrow(classeuce1)-poids2[j]-chitable[1,2]
213                         chitest<-chisq.test(chitable,correct=FALSE)
214                         if ((chitable[1,1]-chitest$expected[1,1])<0) {
215                             chicroise[i,j]<--round(chitest$statistic,digits=7)
216                         } else {
217                             chicroise[i,j]<-round(chitest$statistic,digits=7)
218                 #print(chitest)
219                         }
220                     }
221                 }   
222             }
223         chicroise
224     }
225
226         dochicroisesimple <- function(croise, mincl) {
227                 chicroise <- croise
228                 for (i in 1:nrow(croise)) {
229                         for (j in 1:ncol(croise)) {
230                                 if (croise[i,j]==0) {
231                                         chicroise[i,j]<-0
232                                 } else if (croise[i,j]<mincl) { 
233                                         chicroise[i,j]<-0
234                                 } else {
235                                         chitable<-matrix(ncol=2,nrow=2)
236                                         chitable[1,1]<-croise[i,j]
237                                         chitable[1,2]<-poids1[i]-chitable[1,1]
238                                         chitable[2,1]<-poids1[j]-chitable[1,1]
239                                         chitable[2,2]<-nrow(classeuce1)-poids1[j]-chitable[1,2]
240                                         chitest<-chisq.test(chitable,correct=FALSE)
241                                         if ((chitable[1,1]-chitest$expected[1,1])<0) {
242                                                 chicroise[i,j]<--round(chitest$statistic,digits=7)
243                                         } else {
244                                                 chicroise[i,j]<-round(chitest$statistic,digits=7)
245                                                 #print(chitest)
246                                         }
247                                 }
248                         }   
249                 }
250                 chicroise
251         }
252         if (classif_mode == 0) {
253                 chicroise <- dochicroise(croise, mincl)
254         } else {
255                 chicroise <- dochicroisesimple(croise, mincl)
256         }
257     
258     print('fin croise')
259         #print(chicroise)
260         #determination des chi2 les plus fort
261         chicroiseori<-chicroise
262
263     doxy <- function(chicroise) {
264         listx <- NULL
265         listy <- NULL
266         listxy <- which(chicroise > 3.84, arr.ind = TRUE)
267         #print(listxy)
268         val <- chicroise[which(chicroise > 3.84)]
269         ord <- order(val, decreasing = TRUE)
270         listxy <- listxy[ord,]
271         #for (i in 1:nrow(listxy)) {
272         #    if ((!listxy[,2][i] %in% listx) & (!listxy[,1][i] %in% listy)) {
273         #        listx <- c(listx, listxy[,2][i])
274         #        listy <- c(listy, listxy[,1][i])
275         #    }
276         #}
277         xy <- list(x = listxy[,2], y = listxy[,1])
278         xy
279     }
280     xy <- doxy(chicroise)
281     print(xy)
282     listx <- xy$x
283     listy <- xy$y
284
285 #       maxi<-vector()
286 #       chimax<-vector()
287 #       for (i in 1:tcl) {
288 #           maxi[i]<-which.max(chicroise)
289 #           chimax[i]<-chicroise[maxi[i]]
290 #           chicroise[maxi[i]]<-0
291 #       }
292 #       testpres<-function(x,listcoord) {
293 #           for (i in 1:length(listcoord)) {
294 #                   if (x==listcoord[i]) {
295 #                       return(-1)
296 #                   } else {
297 #                       a<-1
298 #                   }
299 #           }
300 #           a
301 #       }
302 #       c.len=nrow(chicroise)
303 #       r.len=ncol(chicroise)
304 #       listx<-c(0)
305 #       listy<-c(0)
306 #       rang<-0
307 #       cons<-list()
308 #       #on garde une valeur par ligne / colonne
309 #       for (i in 1:length(maxi)) {
310 #       #coordonnées de chi2 max
311 #        #coord <- arrayInd(maxi[i], dim(chicroise))
312 #        #x.co <- coord[1,2]
313 #        #y.co <- coord[1,1]
314 #           x.co<-ceiling(maxi[i]/c.len)
315 #           y.co<-maxi[i]-(x.co-1)*c.len
316 #        #print(x.co)
317 #        #print(y.co)
318 #        #print(arrayInd(maxi[i], dim(chicroise)))
319 #           a<-testpres(x.co,listx)
320 #           b<-testpres(y.co,listy)
321 #           
322 #           if (a==1) {
323 #                       if (b==1) {
324 #                           rang<-rang+1
325 #                           listx[rang]<-x.co
326 #                           listy[rang]<-y.co
327 #                       }
328 #           }
329 #           cons[[1]]<-listx
330 #           cons[[2]]<-listy
331 #       }
332         #pour ecrire les resultats
333         for (i in 1:length(listx)) {
334             txt<-paste(listx[i]+1,listy[i]+1,sep=' ')
335             txt<-paste(txt,croise[listy[i],listx[i]],sep=' ')
336             txt<-paste(txt,chicroiseori[listy[i],listx[i]],sep=' ')
337             print(txt)
338         }
339
340         #colonne de la classe
341         #trouver les filles et les meres 
342         trouvefillemere<-function(classe,chd) {
343             unique(unlist(chd[chd[,classe%/%2]==classe,]))
344         }
345
346
347 #----------------------------------------------------------------------
348     findbestcoord <- function(classeuce1, classeuce2) {
349         #fillemere1 <- NULL
350         #fillemere2 <- NULL
351
352         #fillemere1 <- unique(classeuce1)
353         #if (classif_mode == 0) {
354         #    fillemere2 <- unique(classeuce2)
355         #} else {
356         #    fillemere2 <- fillemere1
357         #}
358
359         #
360         listcoordok <- list()
361         maxcl <- 0
362         nb <- 0
363         lf1 <- addallfille(chd1$list_fille) 
364         if (classif_mode == 0) {
365             lf2 <- addallfille(chd2$list_fille)
366         } else {
367             lf2 <- lf1
368         }
369         lme1 <- chd1$list_mere
370         if (classif_mode == 0) {
371             lme2 <- chd2$list_mere
372         } else {
373             lme2 <- lme1
374         }
375         for (first in 1:length(listx)) {
376             coordok <- NULL
377             f1 <- NULL
378             f2 <- NULL
379             listxp<-listx
380             listyp<-listy
381             
382             #listxp<-listx[first:length(listx)]
383             #listxp<-c(listxp,listx[1:(first-1)])
384             #listyp<-listy[first:length(listy)]
385             #listyp<-c(listyp,listy[1:(first-1)])
386             listxp <- listxp[order(listx, decreasing = TRUE)]
387             listyp <- listyp[order(listx, decreasing = TRUE)]
388             #listxp<-c(listxp[first:length(listx)], listx[1:(first-1)])
389             #listyp<-c(listyp[first:length(listy)], listy[1:(first-1)])
390             for (i in 1:length(listx)) {
391                 if( (!(listxp[i]+1) %in% f1) & (!(listyp[i]+1) %in% f2) ) {
392                     #print(listyp[i]+1)
393                     #print('not in')
394                     #print(f2)
395                     coordok <- rbind(coordok, c(listyp[i] + 1,listxp[i] + 1))
396                     #print(c(listyp[i] + 1,listxp[i] + 1))
397                     un1 <- getfillemere(lf2, chd2$list_mere, listxp[i] + 1)
398                     f1 <- c(f1, un1)
399                     f1 <- c(f1, listxp[i] + 1)
400                     un2 <- getfillemere(lf1, chd1$list_mere, listyp[i] + 1)
401                     f2 <- c(f2, un2)
402                     f2 <- c(f2, listyp[i] + 1)
403                 }
404                 #print(coordok)
405             }
406             #if (nrow(coordok) > maxcl) {
407                 nb <- 1
408             #    listcoordok <- list()
409                 listcoordok[[nb]] <- coordok
410             #    maxcl <- nrow(coordok)
411             #} else if (nrow(coordok) == maxcl) {
412                 nb <- nb + 1
413             #    listcoordok[[nb]] <- coordok
414             #}
415         }
416         listcoordok <- unique(listcoordok)
417         print(listcoordok)
418         best <- 1
419         if (length(listcoordok) > 1) {
420             maxchi <- 0
421             for (i in 1:length(listcoordok)) {
422                 chi <- NULL
423                 uce <- NULL
424                 for (j in 1:nrow(listcoordok[[i]])) {
425                     chi<-c(chi,chicroiseori[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
426                     uce<-c(uce,croise[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
427                 }
428                 if (maxchi < sum(chi)) {
429                     maxchi <- sum(chi)
430                     suce <- sum(uce)
431                     best <- i
432                 }
433             }
434         print(suce/nrow(classeuce1))
435         }
436         listcoordok[[best]]
437     }
438 #---------------------------------------------------------------------------------   
439         #pour trouver une valeur dans une liste
440         #is.element(elem, list)
441         #== elem%in%list
442     oldfindbestcoord <- function(listx, listy) {
443         coordok<-NULL
444         trouvecoordok<-function(first) {
445             fillemere1<-NULL
446             fillemere2<-NULL
447             listxp<-listx
448             listyp<-listy
449             listxp<-listx[first:length(listx)]
450             listxp<-c(listxp,listx[1:(first-1)])
451             listyp<-listy[first:length(listy)]
452             listyp<-c(listyp,listy[1:(first-1)])
453             for (i in 1:length(listxp)) {
454                 if (!(listxp[i]+1)%in%fillemere1) {
455                         if (!(listyp[i]+1)%in%fillemere2) {
456                             coordok<-rbind(coordok,c(listyp[i]+1,listxp[i]+1))
457                             fillemere1<-c(fillemere1,trouvefillemere(listxp[i]+1,chd2$n1))
458                             fillemere2<-c(fillemere2,trouvefillemere(listyp[i]+1,chd1$n1))
459                         }
460                }
461             }
462             coordok
463         }
464     #fonction pour trouver le nombre maximum de classes
465         findmaxclasse<-function(listx,listy) {
466             listcoordok<-list()
467             maxcl<-0
468             nb<-1
469             for (i in 1:length(listy)) {
470                         coordok<-trouvecoordok(i)
471                         if (maxcl <= nrow(coordok)) {
472                             maxcl<-nrow(coordok)
473                             listcoordok[[nb]]<-coordok
474                             nb<-nb+1
475                         }
476             }
477             listcoordok<-unique(listcoordok)
478             print(listcoordok)
479                 #si plusieurs ensemble avec le meme nombre de classe, on conserve
480                 #la liste avec le plus fort chi2
481             if (length(listcoordok)>1) {
482                     maxchi<-0
483                     best<-NULL
484                     for (i in 1:length(listcoordok)) {
485                         chi<-NULL
486                         uce<-NULL
487                         if (nrow(listcoordok[[i]])==maxcl) {
488                             for (j in 1:nrow(listcoordok[[i]])) {
489                                 chi<-c(chi,croise[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
490                                 uce<-c(uce,chicroiseori[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
491                             }
492                             if (maxchi < sum(chi)) {
493                                 maxchi <- sum(chi)
494                                 suce <- sum(uce)
495                                 best <- i
496                             } 
497                         }
498                     }
499             }
500             print((maxchi/nrow(classeuce1)*100))
501             listcoordok[[best]]
502         }
503         print('cherche max')
504             coordok<-findmaxclasse(listx,listy)
505             coordok
506     }
507         #findmaxclasse(listx,listy)
508         #coordok<-trouvecoordok(1)
509     #coordok <- oldfindbestcoord(listx, listy)
510     coordok <- findbestcoord(listx, listy)
511
512
513         lfilletot<-function(classeuce,x) {
514             listfille<-NULL
515             for (classe in 1:nrow(coordok)) {
516                         listfille<-unique(c(listfille,fille(coordok[classe,x],classeuce)))
517                         listfille
518             }
519         }
520     print('listfille')
521         listfille1<-lfilletot(classeuce1,1)
522         if (classif_mode == 0) {
523                 listfille2<-lfilletot(classeuce2,2)
524         }
525
526         #utiliser rownames comme coordonnees dans un tableau de 0
527         Assignclasse<-function(classeuce,x) {
528             nchd<-matrix(0,ncol=ncol(classeuce),nrow=nrow(classeuce))
529             for (classe in 1:nrow(coordok)) {
530                         clnb<-coordok[classe,x]
531                         colnb<-clnb%/%2
532             nchd[which(classeuce[,colnb]==clnb), colnb:ncol(nchd)] <- classe
533             }
534             nchd
535         }
536         print('commence assigne new classe')
537         nchd1<-Assignclasse(classeuce1,1)
538         if (classif_mode==0) {
539                 nchd2<-Assignclasse(classeuce2,2)
540         }
541         print('fini assign new classe')
542         #croisep<-matrix(ncol=nrow(coordok),nrow=nrow(coordok))
543         if (classif_mode==0) {
544         nchd2[which(nchd1[,ncol(nchd1)]==0),] <- 0
545         nchd2[which(nchd1[,ncol(nchd1)]!=nchd2[,ncol(nchd2)]),] <- 0
546         nchd1[which(nchd2[,ncol(nchd2)]==0),] <- 0
547         }
548
549         print('fini croise')
550         elim<-which(nchd1[,ncol(nchd1)]==0)
551         keep<-which(nchd1[,ncol(nchd1)]!=0)
552         n1<-nchd1[nchd1[,ncol(nchd1)]!=0,]
553         if (classif_mode==0) {
554                 n2<-nchd2[nchd2[,ncol(nchd2)]!=0,]
555         } else {
556                 classeuce2 <- NULL
557         }
558         #clnb<-nrow(coordok)
559         print('fini')
560         write.csv2(nchd1[,ncol(nchd1)],uceout)
561         res <- list(n1 = nchd1, coord_ok = coordok, cuce1 = classeuce1, cuce2 = classeuce2)
562         res
563 }