1 #Author: Pierre Ratinaud
2 #Copyright (c) 2008-2009 Pierre Ratinaud
6 #fonction pour la double classification
7 #cette fonction doit etre splitter en 4 ou 5 fonctions
9 AssignClasseToUce <- function(listuce, chd) {
10 print('assigne classe -> uce')
14 fille<-function(classe,classeuce) {
15 listfm<-unique(unlist(classeuce[classeuce[,classe%/%2]==classe,]))
16 listf<-listfm[listfm>=classe]
22 croiseeff <- function(croise, classeuce1, classeuce2) {
25 for (i in 1:ncol(classeuce1)) {
30 for (j in 1:ncol(classeuce2)) {
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))
42 addallfille <- function(lf) {
44 for (i in 1:length(lf)) {
45 if (! is.null(lf[[i]])) {
50 if (f1 > length(lf)) {
51 for (j in (length(lf) + 1):f2) {
62 getfille <- function(nlf, classe, pf) {
63 if (!length(nlf[[classe]])) {
66 for (cl in nlf[[classe]]) {
68 if (cl <= length(nlf)) {
69 pf <- getfille(nlf, cl, pf)
76 getmere <- function(list_mere, classe) {
77 i <- as.numeric(classe)
80 pf <- c(pf, list_mere[[i]])
86 getfillemere <- function(list_fille, list_mere, classe) {
87 return(c(getfille(list_fille, classe, NULL), getmere(list_mere, classe)))
90 getlength <- function(n1, clnb) {
92 tab <- table(n1[,colnb])
93 eff <- tab[which(names(tab) == as.character(clnb))]
98 find.terminales <- function(n1, list_mere, list_fille, mincl) {
99 tab <- table(n1[,ncol(n1)])
100 clnames <- rownames(tab)
101 terminales <- clnames[which(tab >= mincl)]
102 tocheck <- setdiff(clnames, terminales)
103 if ("0" %in% terminales) {
104 terminales <- terminales[which(terminales != 0)]
106 if (length(terminales) == 0) {
107 return('no clusters')
109 if ("0" %in% tocheck) {
110 tocheck <- tocheck[which(tocheck != "0")]
114 while (length(tocheck)!=0) {
115 for (val in tocheck) {
117 mere <- list_mere[[as.numeric(val)]]
121 ln.mere <- getlength(n1, mere)
124 filles.mere <- getfille(list_fille, mere, NULL)
127 filles.mere <- filles.mere[which(filles.mere != val)]
129 if ((ln.mere >= mincl) & (length(intersect(filles.mere, tocheck)) == 0) & (length(intersect(filles.mere, terminales)) == 0 )) {
131 terminales <- c(terminales, mere)
132 for (f in c(filles.mere, val, mere)) {
133 tocheck <- tocheck[which(tocheck != f)]
135 } else if ((ln.mere >= mincl) & (length(intersect(filles.mere, terminales)) == 0) & (length(intersect(filles.mere, tocheck))!=0)){
136 print('mere a checke cause fille ds tocheck')
137 tocheck <- tocheck[which(tocheck != val)]
138 tocheck <- c(mere, tocheck)
141 print('pas ok on vire du check')
142 tocheck <- tocheck[which(tocheck != val)]
146 tocheck <- tocheck[which(tocheck != val)]
156 make.classes <- function(terminales, n1, tree, lf) {
157 term.n1 <- unique(n1[,ncol(n1)])
158 tree.tip <- tree$tip.label
159 cl.n1 <- n1[,ncol(n1)]
160 classes <- rep(NA, nrow(n1))
161 cl.names <- 1:length(terminales)
163 for (i in 1:length(terminales)) {
164 if (terminales[i] %in% term.n1) {
165 classes[which(cl.n1==terminales[i])] <- cl.names[i]
166 new.cl[[terminales[i]]] <- cl.names[i]
167 tree.tip[which(tree.tip==terminales[i])] <- paste('a', cl.names[i], sep='')
169 filles <- getfille(lf, as.numeric(terminales[i]), NULL)
170 tochange <- intersect(filles, term.n1)
171 for (cl in tochange) {
172 classes[which(cl.n1==cl)] <- cl.names[i]
173 new.cl[[cl]] <- cl.names[i]
174 tree.tip[which(tree.tip==cl)] <- paste('a', cl.names[i], sep='')
178 make.tip <- function(x) {
179 if (substring(x,1,1)=='a') {
180 return(substring(x,2))
185 tree$tip.label <- tree.tip
187 tree$tip.label <- sapply(tree.tip, make.tip)
188 tovire <- sapply(tree.tip, function(x) {substring(x,1,1)!='a'})
189 tovire <- which(tovire)
190 ntree <- drop.tip(ntree, tip=tovire)
191 en.double <- ntree$tip.label[duplicated(ntree$tip.label)]
192 en.double <- unique(en.double)
193 tovire <- sapply(en.double, function(x) {which(ntree$tip.label == x)[1]})
194 ntree <- drop.tip(ntree, tip=tovire)
195 ntree$tip.label <- sapply(ntree$tip.label, function(x) {substring(x,2)})
196 classes[which(is.na(classes))] <- 0
197 res <- list(dendro_tot_cl = tree, tree.cl = ntree, n1=as.matrix(classes))
201 #nbt nbcl = nbt+1 tcl=((nbt+1) *2) - 2 n1[,ncol(n1)], nchd1[,ncol(nchd1)]
202 Rchdtxt<-function(uceout, chd1, chd2 = NULL, mincl=0, classif_mode=0, nbt = 9) {
203 #FIXME: le nombre de classe peut etre inferieur
205 tcl <- ((nbt+1) * 2) - 2
206 #Assignation des classes
207 classeuce1<-AssignClasseToUce(listuce1,chd1$n1)
208 if (classif_mode==0) {
209 classeuce2<-AssignClasseToUce(listuce2,chd2$n1)
212 # classeuce2<-classeuce1
215 #calcul des poids (effectifs)
217 makepoids <- function(classeuce, poids) {
223 poids[cl1 - 1] <- length(which(classeuce[,i] == cl1))
224 poids[cl2 - 1] <- length(which(classeuce[,i] == cl2))
229 # makepoids<-function(classeuce,poids) {
230 # for (classes in 2:(tcl + 1)){
231 # for (i in 1:ncol(classeuce)) {
232 # if (poids[(classes-1)]<length(classeuce[,i][classeuce[,i]==classes])) {
233 # poids[(classes-1)]<-length(classeuce[,i][classeuce[,i]==classes])
240 poids1<-vector(mode='integer',length = tcl)
241 poids1<-makepoids(classeuce1,poids1)
242 if (classif_mode==0) {
243 poids2<-vector(mode='integer',length = tcl)
244 poids2<-makepoids(classeuce2,poids2)
249 print('croisement classif')
251 # croise=matrix(ncol=tcl,nrow=tcl)
253 # docroise <- function(croise, classeuce1, classeuce2) {
254 # #production du tableau de contingence
255 # for (i in 1:ncol(classeuce1)) {
256 # #poids[i]<-length(classeuce1[,i][x==classes])
257 # for (j in 1:ncol(classeuce2)) {
258 # tablecroise<-table(classeuce1[,i],classeuce2[,j])
259 # tabcolnames<-as.numeric(colnames(tablecroise))
260 # #tabcolnames<-c(tabcolnames[(length(tabcolnames)-1)],tabcolnames[length(tabcolnames)])
261 # tabrownames<-as.numeric(rownames(tablecroise))
262 # #tabrownames<-c(tabrownames[(length(tabrownames)-1)],tabrownames[length(tabrownames)])
263 # for (k in (ncol(tablecroise)-1):ncol(tablecroise)) {
264 # for (l in (nrow(tablecroise)-1):nrow(tablecroise)) {
265 # croise[(tabrownames[l]-1),(tabcolnames[k]-1)]<-tablecroise[l,k]
272 if (classif_mode==0) {
273 croise <- croiseeff( matrix(ncol=tcl,nrow=tcl), classeuce1, classeuce2)
275 croise <- croiseeff( matrix(ncol=tcl,nrow=tcl), classeuce1, classeuce1)
278 if (classif_mode == 0) {ind <- (nbcl * 2)} else {ind <- nbcl}
280 mincl<-round(nrow(classeuce1)/ind)
288 #tableau des chi2 signes
292 # nr <- nrow(classeuce1)
293 # newchicroise <- function(croise, mincl, nr, poids1, poids2) {
294 # chicroise <- croise
295 # chicroise[which(croise < mincl)] <- 0
296 # tocompute <- which(chicroise > 0, arr.ind = TRUE)
297 # for (i in 1:nrow(tocompute)) {
298 # chitable <- matrix(ncol=2,nrow=2)
299 # chitable[1,1] <- croise[tocompute[i,1], tocompute[i,2]]
300 # chitable[1,2] <- poids1[tocompute[i,1]] - chitable[1,1]
301 # chitable[2,1] <- poids2[tocompute[i,2]] - chitable[1,1]
302 # chitable[2,2] <- nr - poids2[tocompute[i,2]] - chitable[1,2]
303 # chitest<-chisq.test(chitable,correct=FALSE)
304 # 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))
311 dochicroise <- function(croise, mincl) {
313 for (i in 1:nrow(croise)) {
314 for (j in 1:ncol(croise)) {
315 if (croise[i,j]==0) {
317 } else if (croise[i,j]<mincl) {
320 chitable<-matrix(ncol=2,nrow=2)
321 chitable[1,1]<-croise[i,j]
322 chitable[1,2]<-poids1[i]-chitable[1,1]
323 chitable[2,1]<-poids2[j]-chitable[1,1]
324 chitable[2,2]<-nrow(classeuce1)-poids2[j]-chitable[1,2]
325 chitest<-chisq.test(chitable,correct=FALSE)
326 if ((chitable[1,1]-chitest$expected[1,1])<0) {
327 chicroise[i,j]<--round(chitest$statistic,digits=7)
329 chicroise[i,j]<-round(chitest$statistic,digits=7)
338 dochicroisesimple <- function(croise, mincl) {
340 for (i in 1:nrow(croise)) {
341 for (j in 1:ncol(croise)) {
342 if (croise[i,j]==0) {
344 } else if (croise[i,j]<mincl) {
347 chitable<-matrix(ncol=2,nrow=2)
348 chitable[1,1]<-croise[i,j]
349 chitable[1,2]<-poids1[i]-chitable[1,1]
350 chitable[2,1]<-poids1[j]-chitable[1,1]
351 chitable[2,2]<-nrow(classeuce1)-poids1[j]-chitable[1,2]
352 chitest<-chisq.test(chitable,correct=FALSE)
353 if ((chitable[1,1]-chitest$expected[1,1])<0) {
354 chicroise[i,j]<--round(chitest$statistic,digits=7)
356 chicroise[i,j]<-round(chitest$statistic,digits=7)
364 if (classif_mode == 0) {
365 chicroise <- dochicroise(croise, mincl)
367 chicroise <- dochicroisesimple(croise, mincl)
372 #determination des chi2 les plus fort
373 chicroiseori<-chicroise
375 doxy <- function(chicroise) {
378 listxy <- which(chicroise > 3.84, arr.ind = TRUE)
380 val <- chicroise[which(chicroise > 3.84)]
381 ord <- order(val, decreasing = TRUE)
382 listxy <- listxy[ord,]
383 #for (i in 1:nrow(listxy)) {
384 # if ((!listxy[,2][i] %in% listx) & (!listxy[,1][i] %in% listy)) {
385 # listx <- c(listx, listxy[,2][i])
386 # listy <- c(listy, listxy[,1][i])
389 xy <- list(x = listxy[,2], y = listxy[,1])
392 xy <- doxy(chicroise)
399 # maxi[i]<-which.max(chicroise)
400 # chimax[i]<-chicroise[maxi[i]]
401 # chicroise[maxi[i]]<-0
403 # testpres<-function(x,listcoord) {
404 # for (i in 1:length(listcoord)) {
405 # if (x==listcoord[i]) {
413 # c.len=nrow(chicroise)
414 # r.len=ncol(chicroise)
419 # #on garde une valeur par ligne / colonne
420 # for (i in 1:length(maxi)) {
421 # #coordonnées de chi2 max
422 # #coord <- arrayInd(maxi[i], dim(chicroise))
423 # #x.co <- coord[1,2]
424 # #y.co <- coord[1,1]
425 # x.co<-ceiling(maxi[i]/c.len)
426 # y.co<-maxi[i]-(x.co-1)*c.len
429 # #print(arrayInd(maxi[i], dim(chicroise)))
430 # a<-testpres(x.co,listx)
431 # b<-testpres(y.co,listy)
443 #pour ecrire les resultats
444 for (i in 1:length(listx)) {
445 txt<-paste(listx[i]+1,listy[i]+1,sep=' ')
446 txt<-paste(txt,croise[listy[i],listx[i]],sep=' ')
447 txt<-paste(txt,chicroiseori[listy[i],listx[i]],sep=' ')
450 #colonne de la classe
451 #trouver les filles et les meres
452 trouvefillemere<-function(classe,chd) {
453 unique(unlist(chd[chd[,classe%/%2]==classe,]))
457 #----------------------------------------------------------------------
458 findbestcoord <- function(classeuce1, classeuce2, classif_mode, nbt) {
462 #fillemere1 <- unique(classeuce1)
463 #if (classif_mode == 0) {
464 # fillemere2 <- unique(classeuce2)
466 # fillemere2 <- fillemere1
470 listcoordok <- list()
473 lf1 <- addallfille(chd1$list_fille)
474 if (classif_mode == 0) {
475 lf2 <- addallfille(chd2$list_fille)
478 listx<-listx[1:((nbt+1)*2)]
479 listy<-listy[1:((nbt+1)*2)]
481 lme1 <- chd1$list_mere
482 if (classif_mode == 0) {
483 lme2 <- chd2$list_mere
487 print('length listx')
489 #if (classif_mode == 0) {
490 for (first in 1:length(listx)) {
497 #listxp<-listx[first:length(listx)]
498 #listxp<-c(listxp,listx[1:(first-1)])
499 #listyp<-listy[first:length(listy)]
500 #listyp<-c(listyp,listy[1:(first-1)])
501 listxp <- listxp[order(listx, decreasing = TRUE)]
502 listyp <- listyp[order(listx, decreasing = TRUE)]
503 #listxp<-c(listxp[first:length(listx)], listx[1:(first-1)])
504 #listyp<-c(listyp[first:length(listy)], listy[1:(first-1)])
505 for (i in 1:length(listx)) {
506 if( (!(listxp[i]+1) %in% f1) & (!(listyp[i]+1) %in% f2) ) {
510 coordok <- rbind(coordok, c(listyp[i] + 1,listxp[i] + 1))
511 #print(c(listyp[i] + 1,listxp[i] + 1))
512 un1 <- getfillemere(lf2, chd2$list_mere, listxp[i] + 1)
514 f1 <- c(f1, listxp[i] + 1)
515 un2 <- getfillemere(lf1, chd1$list_mere, listyp[i] + 1)
517 f2 <- c(f2, listyp[i] + 1)
521 #if (nrow(coordok) > maxcl) {
523 # listcoordok <- list()
524 listcoordok[[nb]] <- coordok
525 # maxcl <- nrow(coordok)
526 #} else if (nrow(coordok) == maxcl) {
528 # listcoordok[[nb]] <- coordok
532 # stopid <- ((nbt+1) * 2) - 2
533 # for (first in 1:stopid) {
540 # #listxp<-listx[first:length(listx)]
541 # #listxp<-c(listxp,listx[1:(first-1)])
542 # #listyp<-listy[first:length(listy)]
543 # #listyp<-c(listyp,listy[1:(first-1)])
544 # listxp <- listxp[order(listx, decreasing = TRUE)]
545 # listyp <- listyp[order(listx, decreasing = TRUE)]
546 # #listxp<-c(listxp[first:length(listx)], listx[1:(first-1)])
547 # #listyp<-c(listyp[first:length(listy)], listy[1:(first-1)])
548 # for (i in 1:stopid) {
549 # if( (!(listxp[i]+1) %in% f1) & (!(listyp[i]+1) %in% f2) ) {
550 # #print(listyp[i]+1)
553 # coordok <- rbind(coordok, c(listyp[i] + 1,listxp[i] + 1))
554 # #print(c(listyp[i] + 1,listxp[i] + 1))
555 # un1 <- getfillemere(lf2, chd2$list_mere, listxp[i] + 1)
557 # f1 <- c(f1, listxp[i] + 1)
558 # un2 <- getfillemere(lf1, chd1$list_mere, listyp[i] + 1)
560 # f2 <- c(f2, listyp[i] + 1)
564 # #if (nrow(coordok) > maxcl) {
566 # # listcoordok <- list()
567 # listcoordok[[nb]] <- coordok
568 # # maxcl <- nrow(coordok)
569 # #} else if (nrow(coordok) == maxcl) {
571 # # listcoordok[[nb]] <- coordok
576 listcoordok <- unique(listcoordok)
579 if (length(listcoordok) > 1) {
581 for (i in 1:length(listcoordok)) {
584 for (j in 1:nrow(listcoordok[[i]])) {
585 chi<-c(chi,chicroiseori[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
586 uce<-c(uce,croise[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
588 if (maxchi < sum(chi)) {
594 print(suce/nrow(classeuce1))
598 #---------------------------------------------------------------------------------
599 #pour trouver une valeur dans une liste
600 #is.element(elem, list)
602 oldfindbestcoord <- function(listx, listy) {
604 trouvecoordok<-function(first) {
609 listxp<-listx[first:length(listx)]
610 listxp<-c(listxp,listx[1:(first-1)])
611 listyp<-listy[first:length(listy)]
612 listyp<-c(listyp,listy[1:(first-1)])
613 for (i in 1:length(listxp)) {
614 if (!(listxp[i]+1)%in%fillemere1) {
615 if (!(listyp[i]+1)%in%fillemere2) {
616 coordok<-rbind(coordok,c(listyp[i]+1,listxp[i]+1))
617 fillemere1<-c(fillemere1,trouvefillemere(listxp[i]+1,chd2$n1))
618 fillemere2<-c(fillemere2,trouvefillemere(listyp[i]+1,chd1$n1))
624 #fonction pour trouver le nombre maximum de classes
625 findmaxclasse<-function(listx,listy) {
629 for (i in 1:length(listy)) {
630 coordok<-trouvecoordok(i)
631 if (maxcl <= nrow(coordok)) {
633 listcoordok[[nb]]<-coordok
637 listcoordok<-unique(listcoordok)
639 #si plusieurs ensemble avec le meme nombre de classe, on conserve
640 #la liste avec le plus fort chi2
641 if (length(listcoordok)>1) {
644 for (i in 1:length(listcoordok)) {
647 if (nrow(listcoordok[[i]])==maxcl) {
648 for (j in 1:nrow(listcoordok[[i]])) {
649 chi<-c(chi,croise[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
650 uce<-c(uce,chicroiseori[(listcoordok[[i]][j,1]-1),(listcoordok[[i]][j,2]-1)])
652 if (maxchi < sum(chi)) {
660 print((maxchi/nrow(classeuce1)*100))
664 coordok<-findmaxclasse(listx,listy)
667 #findmaxclasse(listx,listy)
668 #coordok<-trouvecoordok(1)
669 #coordok <- oldfindbestcoord(listx, listy)
670 print('begin bestcoord')
671 coordok <- findbestcoord(listx, listy, classif_mode, nbt)
674 lfilletot<-function(classeuce,x) {
676 for (classe in 1:nrow(coordok)) {
677 listfille<-unique(c(listfille,fille(coordok[classe,x],classeuce)))
682 listfille1<-lfilletot(classeuce1,1)
683 if (classif_mode == 0) {
684 listfille2<-lfilletot(classeuce2,2)
687 #utiliser rownames comme coordonnees dans un tableau de 0
688 Assignclasse<-function(classeuce,x) {
689 nchd<-matrix(0,ncol=ncol(classeuce),nrow=nrow(classeuce))
690 for (classe in 1:nrow(coordok)) {
691 clnb<-coordok[classe,x]
693 nchd[which(classeuce[,colnb]==clnb), colnb:ncol(nchd)] <- classe
697 print('commence assigne new classe')
698 nchd1<-Assignclasse(classeuce1,1)
699 if (classif_mode==0) {
700 nchd2<-Assignclasse(classeuce2,2)
702 print('fini assign new classe')
703 #croisep<-matrix(ncol=nrow(coordok),nrow=nrow(coordok))
704 if (classif_mode==0) {
705 nchd2[which(nchd1[,ncol(nchd1)]==0),] <- 0
706 nchd2[which(nchd1[,ncol(nchd1)]!=nchd2[,ncol(nchd2)]),] <- 0
707 nchd1[which(nchd2[,ncol(nchd2)]==0),] <- 0
711 elim<-which(nchd1[,ncol(nchd1)]==0)
712 keep<-which(nchd1[,ncol(nchd1)]!=0)
713 n1<-nchd1[nchd1[,ncol(nchd1)]!=0,]
714 if (classif_mode==0) {
715 n2<-nchd2[nchd2[,ncol(nchd2)]!=0,]
721 write.csv2(nchd1[,ncol(nchd1)],uceout)
722 res <- list(n1 = nchd1, coord_ok = coordok, cuce1 = classeuce1, cuce2 = classeuce2)