# Copyright 2006 Université de Liège, Nicolas Mazziotta, 
# sous les termes de la licence GPL2
#
# Macros pour l'analyse des correspondances

require("ade4", quiet=T);

my.makebox <- function(x,y, main="", xlab="", ylab="") {
	par(font.main=1, font.lab=3, cex.axis=1, xpd=F);
	plot(
		x,y,
		type="n",
		asp=1,
		ax=T,
		main=main,
		xlab=xlab,
		ylab=ylab,
		las=2
	);
 	abline(h=0, col="gray", lty=1)
	abline(v=0, col="gray", lty=1)	
	#axis(1)
	#axis(2)
} 	

adjusteig <- function(coa, q) {
	# Ajustement des valeurs propres à partir de celles de B
	# Greenacre1994multipleandjoint: 156
	eig<-coa$eig^2 # Les valeurs propres de B valent le quarré de celles de Z, cf. Lebart Morineau Piron 1998: 126
	adj<-NULL
	for (i in 1:length(eig)) {
		lambda<-eig[i]
		if ( sqrt(lambda) > (1/q) ) {
			adj[i]<-( q/(q-1) )^2 * ( sqrt(lambda) - 1/q )^2
		} else {
			return(adj)
		}
	}
	return(adj)
}

my.ploteig <-function(coa, normal=FALSE, adjusted=TRUE) {
	# histogramme des valeurs propres (normales ou ajustées)
	if (normal | is.null(coa$adj))
	{
		eigRatio <- sum(coa$eig)/100;
			par(
				font.main=1, 
				font.lab=3, 
				cex.axis=.8,
				las=2
			);
			hist.labs<-c(rev(round(coa$eig/eigRatio, dig=2)));
			for (i in 1:length(hist.labs)) {
				hist.labs[i]<-paste(hist.labs[i],"%", sep="");
			}

			barplot(
				main="Valeurs propres",
					rev(coa$eig),
					names.arg=hist.labs,
					horiz=T,
					xlab=paste(""),
					ylab=paste(""),
			);
	}
	if (adjusted & !is.null(coa$adj)) {
		eigRatio <- sum(coa$adj)/100;
			par(
				font.main=1, 
				font.lab=3, 
				cex.axis=.8,
				las=2
			);
			hist.labs<-c(rev(round(coa$adj/eigRatio, dig=2)));
			for (i in 1:length(hist.labs)) {
				hist.labs[i]<-paste(hist.labs[i],"%", sep="");
			}

			barplot(
					rev(coa$adj),
					names.arg=hist.labs,
					horiz=T,
					xlab=paste(""),
					ylab=paste(""),
			);
	}
}

my.plotcoa <- function(coa, rows=T, cols=T, xax=1, yax=2, one=T, file="Rplots.pdf", font=c(1,3), col=c("black","black"), draw=draw.ellipse, ploteig=T,
adj=T, afc.lab="", main=NULL, xlab=NULL, ylab=NULL) {
	eigen<-coa$eig
	adjusted<-""
	if (adj & !is.null(coa$adj) & length(coa$adj) > 1) {
		eigen<-coa$adj
		adjusted<-" ajustée"
	}
	eigRatio <- sum(eigen)/100;
	rows.x<-coa$li[,xax];
	cols.x<-coa$co[,xax];
	rows.y<-coa$li[,yax];
	cols.y<-coa$co[,yax];
	

	# histogramme des valeurs propres
	if (ploteig == TRUE) {
		my.ploteig(coa)
	}

	x.all<-c()
	y.all<-c()
	if (rows == TRUE) {
		x.all<-c(x.all, rows.x)
		y.all<-c(y.all, rows.y)
	}
	if (cols == TRUE) {	
		x.all<-c(x.all, cols.x)
		y.all<-c(y.all, cols.y)
	}

	if (is.null(main)) main=paste(afc.lab, "\nContribution", adjusted, " du plan : ",round(eigen[xax]+eigen[yax], dig=4), " (",round((eigen[xax]+eigen[yax])/eigRatio, dig=2), " \%)", sep="")
	if (is.null(xlab)) xlab=paste("Contribution",adjusted," de f",xax," : ", round(eigen[xax], dig=4), " (", round(eigen[xax]/eigRatio, dig=2), " \%)", sep="")
	if (is.null(ylab)) ylab=paste("Contribution",adjusted," de f",yax," : ", round(eigen[yax], dig=4), " (", round(eigen[yax]/eigRatio, dig=2), " \%)", sep="")

	# scatterplot
	my.makebox(
		x=c(x.all),
		y=c(y.all),
		main=main,
		ylab=ylab,
		xlab=xlab
	)

		if (!is.null(coa$summary)) { 
			if (rows == TRUE) {	
				csq<-abs(coa$summary$rows[,c(paste("CS",xax,sep=""),paste("CS",yax,sep=""))])
				o<-order(-rowSums(csq))
					points(rows.x[o],rows.y[o], pch=21, cex=rev(sort(rowSums(csq)*10)), col="white", bg="dimgray")
			}
			if (cols == TRUE) {	
				csq<-abs(coa$summary$cols[,c(paste("CS",xax,sep=""),paste("CS",yax,sep=""))])
				o<-order(-rowSums(csq))
				points(cols.x[o],cols.y[o], pch=21, cex=rev(sort(rowSums(csq)*10)), col="white", bg="gray")
			}
			box()
		}
	if (rows == TRUE) {	
		text(rows.x,rows.y,lab=names(coa$lw),cex=.8,font=font[1], col=col[1]) 
	}
	if (cols == TRUE) {	
		text(cols.x,cols.y,lab=names(coa$cw),cex=.8,font=font[2], col=col[2]) 
	}

}

my.plotsup <- function(sup, add=F, xax=1, yax=2) {
	x<-sup[[xax]][,"CRD"];
	y<-sup[[yax]][,"CRD"];
	cex<-.9
	srt<-0
	font<-1
	if ( add==FALSE ) { 
		my.makebox(
			x=x,
			y=y
		)
	}	else {
		mtext("Éléments supplémentaires en italiques", cex=.7, font=3, las=1)
		cex<-.7
		font<-3
		#srt<-45
	}
	text(x,y, lab=rownames(sup$"1"), font=font, cex=cex, srt=srt)
	box()
}


my.supvar <-function(coa, df, vars=list(), xax=1, yax=2, axis=c(1:2), threshold=.001) {

# FIXME: mal structuré, permettre plus d'axes

	# vars=list(name=c(values),...)
	#  'name' est le nom de la variable
	#  'values' sont les valeurs possibles de la variable (si vide,toutes)
	names<-names(vars);
	points<-list("1"=NULL,"2"=NULL);
	for(i in 1:length(names)) {
		if (length(vars[[i]])<=0) {
			print(levels(df[,names[i]]))
			vars[[i]]<-levels(df[,names[i]])
		}
		for(j in 1:length(vars[[i]])) {
			# indice des individus possédant la modalité
			r<-rownames(df[df[,names[i]] == vars[[i]][j],])
			for (k in 1:length(axis)) {
				# Calcul des coordonnées, cf. Lebart-Morineau-Piron 1998: 123-125
				l<-length(r)
				x<-sum(coa$li[r,axis[k]])/l/sqrt(coa$eig[axis[k]]);
				# vtest (disto en écarts-types)
				tx<-sqrt(l*(dim(df)[1]-1)/(dim(df)[1]-length(r)))*x
				tx<-ifelse(is.infinite(tx), NA, tx)
				# probabilité (selon la loi normale)
				px<-pnorm(sqrt(tx**2), low=F)
				#if (is.na(px) | px > threshold) {next}
				xline<-cbind(x,tx,px);
				rownames<-paste(names[[i]],vars[[i]][j],sep=".")
				rownames(xline)<-rownames
				points[[axis[k]]]<-rbind(points[[axis[k]]],xline)
				colnames(points[[axis[k]]])<-c("CRD","VTS","PRB")
				points[[axis[k]]]=na.omit(points[[axis[k]]])
			}
		}
	}
	points
} 

# Tableau avec les coordonnées, les contibutions et les cosinus carrés 
# des variables actives
my.summary <- function(coa, r=TRUE, c=TRUE) {
	inertia<-inertia.dudi(coa, r=r, c=c);

	summary<-NULL

	rlength<-dim(coa$li)[1]*(dim(coa$li)[2]*3+1)
	clength<-dim(coa$co)[1]*(dim(coa$co)[2]*3+1)

	rows<-c(1:rlength);
	dim(rows)<-c(dim(coa$li)[1],dim(coa$li)[2]*3+1)
	cols<-c(1:clength);
	dim(cols)<-c(dim(coa$co)[1],dim(coa$co)[2]*3+1)


	coacolnames<-c(1:((dim(coa$co)[2]*3)+1));

	for (i in 1:(dim(coa$co)[2]*3)) {
		mod<-i %% 3;
		div<-(i-mod)/3;
		if (mod == 1) {
			coacolnames[i]<-c(paste("CD",div+1, sep=""))
		} else if (mod == 2) {
			coacolnames[i]<-c(paste("CT",div+1, sep=""))
		} else {
			coacolnames[i]<-c(paste("CS",div, sep=""))
		}
	}

	coacolnames[dim(coa$co)[2]*3+1]<-c("R");

	if (r == TRUE) {
		for (i in 1:dim(coa$li)[2]*3) {
			if (i %% 3 > 0) {next}
			rows[,i-2]<-coa$li[,i/3];
			rows[,i-1]<-inertia$row.abs[,i/3]/10000 ;
			rows[,i]<-inertia$row.rel[,i/3]/10000 ;
		}
		rows[,dim(coa$li)[2]*3+1]<-round(inertia$row.cum[,"remain"]/10000,2);
		colnames(rows)<-coacolnames;
		rownames(rows)<-rownames(coa$li);
		summary$rows<-rows
	}
	
	if (c == TRUE) {
		for (i in 1:dim(coa$co)[2]*3) {
			if (i %% 3 > 0) {next}
			cols[,i-2]<-coa$co[,i/3];
			cols[,i-1]<-inertia$col.abs[,i/3]/10000 ;
			cols[,i]<-inertia$col.rel[,i/3]/10000 ;
		}
		cols[,dim(coa$co)[2]*3+1]<-round(inertia$col.cum[,"remain"]/10000,2);
		colnames(cols)<-coacolnames;
		rownames(cols)<-rownames(coa$co);
		summary$cols<-cols
	}

	class(summary)<-"coasummary"

	return(summary)

}

coa <- function(df, dudi=dudi.coa, print.row=T, q=NULL ,...) {
	coa<-dudi(df, ...);
	coa$summary<-my.summary(coa, r=print.row);
	if (!is.null(q)) {
		coa$adj<-adjusteig(coa, q)
	}
	return(coa)
}

my.bootstrap <- function(df, n=1) {
	length<-dim(df)
	results<-NULL
	for (i in 1:n) {
		s<-sample(1:length[1], length, rep=T);
		bdf<-df[s,]
		for (j in 1:length[2]) {
			bdf[,j]<-as.factor(bdf[,j])
		}
		results[[i]]<-bdf
		
	}
	return(results)
}


# TODO utiliser

# http://zoonek2.free.fr/UNIX/48_R/05.html [060816]
draw.ellipse <- function (
    x, y=NULL, N=100,
    method=lines, ...
  ) {
  if (is.null(y)) {
    y <- x[,2]
    x <- x[,1]
  }
  centre <- c(mean(x), mean(y))
  m <- matrix(c(var(x),cov(x,y),
                cov(x,y),var(y)),
              nr=2,nc=2)
  e <- eigen(m)
  r <- sqrt(e$values)
  v <- e$vectors
  theta <- seq(0,2*pi, length=N)
  x <- centre[1] + r[1]*v[1,1]*cos(theta) +
       r[2]*v[1,2]*sin(theta)
  y <- centre[2] + r[1]*v[2,1]*cos(theta) +
       r[2]*v[2,2]*sin(theta)
  method(x,y,...)
}
# http://zoonek2.free.fr/UNIX/48_R/05.html [060816]
draw.star <- function (x, y=NULL, ...) {
  if (is.null(y)) {
    y <- x[,2]
    x <- x[,1]
  }
  d <- cbind(x,y)
  m <- apply(d, 2, mean)
  segments(m[1],m[2],x,y,...)
}


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