# Copyright 2006 Université de Liège, Nicolas Mazziotta, 
# sous les termes de la licence GPL2

# Tris
require(ade4)
require(vcd) # données catégorielles

tris <- function(df, plat=TRUE, croise=TRUE, pem.do=T, tdl.do=T,...) {
	tris<-NULL
	if (plat==TRUE)	{ tris$plat<-triplat(df) }
	if (croise==TRUE)	{ tris$croise<-tricroise(df, pem=pem.do, tdl=tdl.do) }
	return(tris)
}

triplat <- function(df) {
	colnames<-colnames(df)
	triplat<-NULL
	dim<-dim(df)
	if (length(dim)<2) {
			car<-c("");
			eff<-as.matrix(table(df))
				triplat[[car]]<-cbind(eff, round(eff/sum(eff)*100, digit=2))
				colnames(triplat[[car]])<-c("FAB","FRE")
				triplat[[car]]<-triplat[[car]][order(-triplat[[car]][,"FRE"],rownames(triplat[[car]])),]
				class(triplat[[car]])<-c("triplat")
	}
	else {
		for (i in 1:dim[2]) {
			car<-colnames[i];
			mod<-names(df[,i]);
			eff<-as.matrix(table(df[,i]))
				triplat[[car]]<-cbind(eff, round(eff/sum(eff)*100, digit=2))
				colnames(triplat[[car]])<-c("FAB","FRE")
				triplat[[car]]<-triplat[[car]][order(-triplat[[car]][,"FRE"]),rownames(triplat[[car]])]
				class(triplat[[car]])<-c("triplat")
		}
	}
	class(triplat)<-c("triplatcoll")
	return(triplat)
} 


tricroise <- function(x, ...) UseMethod("tricroise")

tricroise.default <- function(table, pem.do=T, tdl.do=F, fisher.do=2000000, rename=T, multiple=FALSE) {
	if (rename) {
		for (i in 1:length(dimnames(table))) {
			if (!is.null(names(dimnames(table))[i])) dimnames(table)[[i]]<-paste(names(dimnames(table))[i], dimnames(table)[[i]], sep=".")
		}
	}
	tricroise<-NULL
			if (length(table)<2) {return()}
			tricroise[['tab']]<-table
			chi2<-chisq.test(table)
			min.exp <- min(chi2$expected)
			tricroise$fisher<-NULL
			if ((min.exp < 5) & (as.logical(fisher.do)==T)) {tricroise$fisher<-fisher.test(table, workspace=fisher.do)}
			ecart<-chi2$observed-chi2$expected
			contributions<-chi2$resid**2
			tricroise[['chi2']]<-chi2
			tricroise[['multiple']]<-multiple
			tricroise[['sidak']]<-c(1-(1-.05)^(1/multiple), 1-(1-.01)^(1/multiple), 1-(1-.001)^(1/multiple))
			tricroise[['chi2minexp']]<-min.exp
			tricroise[['chi2contributions']]<-round(contributions, d=2)
			#if (length(dim(contributions)) > 1) {items<-rowSums(contributions)}
			#else { items<-contributions }
			items<-contributions
			signif<-round(pchisq(items, df=1, low=F), d=3)
			signif<-ifelse(signif <= 0.05, signif, NA);
			tricroise[['chi2signif']]<-signif;
			if (dim(table)[1]==2 & dim(table)[2]==2) {tricroise[['oddsratio']]<-oddsratio(table, log=F)}
			if (pem.do == TRUE) {
				pem<-pem(table, chi2=chi2, ecart=ecart)
				tricroise[['pem']]<-pem
			}
			if (tdl.do == TRUE) {
				tdl<-tdl(table, chi2=chi2, ecart=ecart)
				tricroise[['tdl']]<-tdl
			}
	class(tricroise)<-c("tricroise")
	return(tricroise)
}

tricroise.data.frame <- function(df, pem.do=T, tdl.do=T, ...) {
	colnames<-colnames(df)
	tricroise<-NULL
	if (length(dim(df))!=2) {
		return(NULL)
	}
	for (i in 1:dim(df)[2]) {
		for (j in i:dim(df)[2]) {
			if (i == j) {next}
			table<-table(df[,c(i,j)])
			current<-tricroise(table, pem.do=pem.do, tdl.do=tdl.do, ...)
			if (is.null(current)) {next}
			tricroise[[colnames[i]]][[colnames[j]]]<-current
			class(tricroise[[colnames[i]]][[colnames[j]]])<-c("tricroise")
		}
	}
	if(!is.null(tricroise)) {	attr(tricroise, "class")<-c("tricroisecoll") }
	return(tricroise)
}

#======================================================================================#
#  Coefficients
#======================================================================================#

# Taux de liaison (cf. Cibois 1993: 45)
tdl <- function (df, chi2=NULL, ecart=NULL) {
	if (is.null(chi2)) {
		chi2<-chisq.test(df)
	}
	if (is.null(ecart)) {
		ecart<-chi2$observed-chi2$expected
	}
	tdl<-ecart/chi2$expected
	return(tdl)
}

pem <- function (df, chi2=NULL, ecart=NULL) {
# Pourcentage d'écart maximum (Cibois 1993)
	if (is.null(chi2)) {
		chi2<-chisq.test(df)
	}
	if (is.null(ecart)) {
		ecart<-chi2$observed-chi2$expected
	}
	n<-sum(df)
	dim<-dim(df)

	margin_c<-rep(colSums(df),times=dim[1])  
	margin_r<-rep(rowSums(df),times=dim[2])
	dim(margin_c)<-c(dim[2],dim[1])
	margin_c<-t(margin_c)
	dim(margin_r)<-c(dim)
	emax_pos<-pmin(margin_c,margin_r)-chi2$expected
	
	emax<-ecart

	if (length(dim(ecart))!=2) {return(NULL)}

	for (i in 1:dim(df)[1]) {
		for (j in 1:dim(df)[2]) {
			if ( ecart[i,j] >= 0 ) {
				emax[i,j]<-emax_pos[i,j]
			} else {
				compl<-margin_r[i,1]+margin_c[1,j]-n
				if (compl > 0) {
					emax[i,j]<-chi2$expected[i,j]-compl
				} 
				else {
					emax[i,j]<-chi2$expected[i,j]
				} 
			}
		}
	}
	pem<-ecart/emax*100
	return(round(pem, dig=2))
}

profil <- function(df) {
	# Peu utile si peu de caractères
	profils<-NULL
	disj<-acm.disjonctif(df)
	disj<-t(disj)
	profils$jaccard<-as.matrix(dist.binary(as.data.frame(disj), method=1))
	profils$pearson<-as.matrix(dist.binary(as.data.frame(disj), method=9))
	return(profils)
}

########################################################################################
#  Graphiques
########################################################################################

#======================================================================================#
#  Graphes de Bertin
#======================================================================================#

plot.tricroisecoll <- function(tri, device=pdf, ...) {
	if (is.null(tri)) {
		return(NULL)
	}
	device(...);
	for (i in 1:length(tri)) {
		for (j in 1:length(tri[[i]])) {
			if (1 %in% dim(tri[[i]][[j]]$tab)) { next }
			plot(tri[[i]][[j]])
		}
	}
	dev.off();
}

plot.tricroise <- function(tri, device=NULL, one=TRUE, file="Rplots.pdf", family="Times") {

	if (!is.null(device)) {device(file=file);}
			par(
				font=1,
				las=2
			);
			t<-tri$tab
			assocplot(t, col=c('white','lightgray'), space=0.1)
			p<-round(tri$chi2$p.value, d=3)
			min.exp<-signif(tri$chi2minexp)
			c<-round(tri$chi2$stat, d=2)
			d<-tri$chi2$para
			title(
				#	main=paste("","chi² =", )
				main=substitute(paste(chi, "² = ", x, ", ddl = ", z, ", p = ",y), list(x=c, y=p, z=d)),
			);
			if (min.exp < 5) {mtext(paste("Attention: certaines valeurs attendues sont inférieures à 5 ", "(min. ",round(min.exp, dig=2), ")", sep = ""), las=0)}
			else {
				if (min.exp < 10) {mtext(paste("Attention: certaines valeurs attendues sont inférieures à 10 ", "(min. ",round(min.exp, dig=2), ")", sep = ""), las=0)}
			}
	
	if (!is.null(device)) {dev.off();}
}

#======================================================================================#
#  Fonction de représentation des écarts
#======================================================================================#

plotchi <- function(t, colnames=c("Attractions","Répulsions (reste)"), do.plot=T) {
	
	poslist<-NULL
	neglist<-NULL

	chi<-chisq.test(t)
	if(is.null(dim(chi$res))) return(list(pos=poslist, neg=neglist))
	s<-sign(chi$res[,2])[rowSums(t)>=10]
	d<-apply(t,1, function(x) { 
		m<-rbind(x,colSums(t)-x)
		chi<-chisq.test(m)
		chi$p.value
	})[rowSums(t)>=10]
	d<-ifelse(d<=1e-300,1e-300,d)
	ratio<-cbind(t,((t[,"1"])/rowSums(t)*100))[rowSums(t)>=10,]
	pos<-d[s>0]
	neg<-d[s<=0]

	threshold<-1-(1-.05)**(1/length(d))
	maxlength<-max(length(pos),length(neg))

		par(mfrow=c(1,2))

	
	if (is.null(dim(ratio))) return(list(pos=poslist, neg=neglist))

	if (do.plot & (length(pos)!=0 | length(neg)!=0)) {
		# attractions
	  if (length(pos)!=0)	p_plotpoints(pos,rbind(ratio[s>0,]),threshold,paste("\n", colnames[1]), maxlength)
		# répulsions
		if (length(neg)!=0)	p_plotpoints(neg,rbind(ratio[s<=0,]),threshold,paste("\n", colnames[2]), maxlength)
	}

	if (length(d[s<0])>1) {neglist<-cbind(d[s<0],ratio[s<0,])} else {neglist<-c(d[s<0],ratio[s<0,])}
	if (length(d[s>0])>1) {poslist<-cbind(d[s>0],ratio[s>0,])} else {poslist<-c(d[s>0],ratio[s>0,])} 
	
	return(list(pos=poslist, neg=neglist))
}

plotbyMM<-function(d, select=c("arg","imm"), argheight=4) {

	argbox<-c(rep(0, 4-argheight), rep(1, argheight))

	d.bak<-d

	if ("arg" %in% select) {
		if ("imm" %in% select) {layout(matrix(c(argbox,argbox,2,2,2,3), 4,3))} 
		else {layout(matrix(c(argbox,argbox,0,0,0,2), 4,3))}
		par(las=1, xaxt="n",yaxt="n", mai=c(.28,1,0,0.2))
		d<-d[,c("phrase","pers","pers-arg","¬pers","¬pers-arg")]
		d<-d[rowSums(!is.na(d))!=0,]
		#d<-d[order(d[,4],d[,7],d[,5],apply(d,1,function (x) sum(na.omit(x)))),c(1:3,5,6,4,7)]   
		d<-d[order(apply(d,1,function (x) sum(na.omit(x))), rownames(d)), order(-apply(d,2,function (x) sum(abs(na.omit(x)))))]   


		#par(mfrow=c(1,2))
		plot(0,0, xlim=c(1,dim(d)[2]), ylim=c(1,dim(d)[1]), t="n",xlab="",ylab="")
		par(xaxt="s",yaxt="s")
		axis(2, at=1:dim(d)[1], rownames(d),cex.axis=.85,)
		axis(1, at=1:dim(d)[2], colnames(d))#,cex.axis=.5,)
		for (i in 1:dim(d)[1]) {
			col<-"white"
			if (i %% 2 == 0) col<-"lightgray" 
			rect(par()$usr[1],i-.5,par()$usr[2],i+.5, col=col, dens=NA)
			for (j in 1:dim(d)[2]) {
				
				if (is.na(d[i,j]) )	{ 
					# ne rien faire
				}
				else {	
					if (d[i,j]>0) {
						if (abs(d[i,j])==1000) {
							points(j,i, pch=22, cex=1.9, bg="black", lwd=.3)	
						}
						else if (abs(d[i,j])==100) {
							points(j,i, pch=24, cex=1.2, bg="black", lwd=.3)	
						}
						else {
							points(j,i, pch=21, cex=.5, bg="black", lwd=.3)	
						}
					} else {
						if (abs(d[i,j])==1000) {
							points(j,i, pch=22, cex=1.9, bg="white", lwd=.3)	
						} 
						else if (abs(d[i,j])==100) {
							points(j,i, pch=24, cex=1.2, bg="white", lwd=.3)	
						}
						else {
							points(j,i, pch=21, cex=.5, bg="white", lwd=.3)	
						}
					}
				}
			}
		}
		box()
	} else {
		layout(matrix(c(0,0,1,1,0,0,0,2,0,0,0,0), 4,3))
	}
	
	if ("imm" %in% select) {

		d<-d.bak
		d<-d[order(d[,"synt-arg"]),c("synt-arg","synt")]
		d<-d[rowSums(!is.na(d))!=0,]
		par(las=1, xaxt="n",yaxt="n", mai=c(.4,.8,0,0.2))
		plot(0,0, xlim=c(.6,dim(d)[2]+.3), ylim=c(1,dim(d)[1]), t="n",xlab="",ylab="")
		par(xaxt="s",yaxt="s")
		axis(2, at=1:dim(d)[1], rownames(d))#,cex.axis=.5,)
		par(cex.axis=1)
		axis(1, at=1:dim(d)[2], colnames(d))#,cex.axis=.5,)
		for (i in 1:dim(d)[1]) {
			col<-"white"
			if (i %% 2 == 0) col<-"lightgray" 
			rect(par()$usr[1],i-.5,par()$usr[2],i+.5, col=col, dens=NA)
			for (j in 1:dim(d)[2]) {
				
				if (is.na(d[i,j]) )	{ 
					# ne rien faire
				}
				else {	
					if (d[i,j]>0) {
						if (abs(d[i,j])==1000) {
							points(j,i, pch=22, cex=1.9, bg="black", lwd=.3)	
						}
						else if (abs(d[i,j])==100) {
							points(j,i, pch=24, cex=1.2, bg="black", lwd=.3)	
						}
						else {
							points(j,i, pch=21, cex=.5, bg="black", lwd=.3)	
						}
					} else {
						if (abs(d[i,j])==1000) {
							points(j,i, pch=22, cex=1.9, bg="white", lwd=.3)	
						} 
						else if (abs(d[i,j])==100) {
							points(j,i, pch=24, cex=1.2, bg="white", lwd=.3)	
						}
						else {
							points(j,i, pch=21, cex=.5, bg="white", lwd=.3)	
						}
					}
				}
			}
		}
		box()
	}
	par(las=1, xaxt="n",yaxt="n", mai=c(.3,0,0,0))
	plot(0,0, xlim=c(-1,10), ylim=c(-2,18), t="n",xlab="",ylab="", bty="n")
	points(00,10, pch=22, cex=1.9, bg="black", lwd=.3)	
	points(00,8, pch=24, cex=1.2, bg="black", lwd=.3)	
	points(00,6, pch=21, cex=.6, bg="black", lwd=.3)	
	points(00,4, pch=22, cex=1.9, bg="white", lwd=.3)	
	points(00,2, pch=24, cex=1.2, bg="white", lwd=.3)	
	points(00,0, pch=21, cex=.6, bg="white", lwd=.3)	
	text(00,10, "Attraction significative", pos=4)
	text(00,8, "Attraction significative (si non ajustée)", pos=4)
	text(00,6, "Attraction non significative (seuil 5%)", pos=4)
	text(00,4, "Répulsion significative", pos=4)
	text(00,2, "Répulsion significative (si non ajustée)", pos=4)
	text(00,0, "Répulsion non significative (seuil 5%)", pos=4)

}

p_plotpoints<-function(vals, ratio, threshold=.05, title="", maxlength=NULL, ignore) {
	vals<-rev(sort(vals))
	if (length(vals)>=2) ratio<-ratio[names(vals),]
	bottom<-5.9
	if (!is.null(maxlength)) { bottom <- bottom -.07*maxlength }
	if (bottom<=.3) {bottom<-.3}
  par(xaxt="n",yaxt="n", mai=c(bottom,.8,.5,1))
	margin<-((log10(max(vals))-log10(min(vals))))
	plot(c(vals),1:length(vals), 
		log="x",
		xlim=c(min(vals),1),#max(vals)*margin),
		ylim=c(0,length(vals)), 
		cex=1, adj=1, t="n", ylab="", xlab="", main=title, font.main=1,xaxs="r", )
	# seuil avec correction de Sidak
	vt<-(vals<threshold)
	if (FALSE %in% vt) {
		thresholdline<-length(vt[vt==FALSE])+.5
		rect(10^par()$usr[1],thresholdline,10^par()$usr[2],par()$usr[3], col="lightgray", dens=NA)
	}
	# Seuil sans correction
	lightthreshold<-.05
	vt<-(vals<lightthreshold)
	if (FALSE %in% vt) {
		lightthresholdline<-length(vt[vt==FALSE])+.5
		rect(10**par()$usr[1],lightthresholdline,10**par()$usr[2],par()$usr[3], col="dimgray", dens=NA)
	}

	for (i in 1:length(vals)) {
		#text(vals[i],i, lab=paste(, sep=""), cex=.7, adj=.5)
		points(vals[i],i, cex=.4, pch=19)
	}

	box()
	par(xaxt="s", yaxt="s", cex.axis=.8)
	par(cex.axis=.5)
	axis(2, at=c(1:length(vals)), names(vals), las=2)
	par(cex.axis=.5)
	axis(4,at=c(1:length(vals)), paste(ratio[,2],"/",ratio[,2]+ratio[,1], " (", format(ratio[,3],dig=2),"%)", sep=""), las=2, tick=T)
	par(cex.axis=.8)
	axis(1, at=c(min(vals),1), c(format(min(vals),dig=3),"1"), padj=-1)
	legend("bottomleft", c("significatif", paste("signif. à 5%, mais \n¬signif. à ",format(threshold*100, dig=3), "%",sep="") ,"¬signif."),
fill=c("white","lightgray","dimgray"), bg='white',cex=.45, merge=F)
}


#vim: encoding=utf-8 fileencoding=utf-8
