require(OpenMx)

fct2 <- function(vF1=NULL,vF2=NULL,corr=0)
 # vF1(ncond) et vF2(ncond) sont les variances de chacun des deux facteurs sous les ncond conditions
 # corr est la corrélation entre les deux facteurs, 0 par défaut
	{
	ncond=length(vF1)
	if (length(vF2)!=ncond) return(NULL)
	aF1<-sqrt(vF1)
	aF2<-sqrt(vF2)
	F<-matrix(c(.1,.2,.3,.7,.8,.9,.7,.7,.7,.3,.3,.3),ncol=2)  # définition pour des facteurs d'amplitude 1.0
	B<-mxMatrix("Diag",nrow=6,ncol=6,values = 1-rowSums(F^2)) # variance de bruit des variables identique pour toutes les conditions
	P<-NULL
	F<-F%*%vec2diag(1/sqrt(colSums(F^2)))
	nomsVar<-c("x1","x2","x3","x4","x5","x6")
	R<-matrix(c(1,corr,corr,1),nrow=2,ncol=2)
	for (g in 1:ncond)
		{
		A<-matrix(c(aF1[g],0,0,aF2[g]),nrow=2,ncol=2)
		C<-matrix(F%*%A%*%R%*%A%*%t(F)+B$values,nrow=6,ncol=6,dimnames=list(nomsVar,nomsVar))
		P<-c(P,mxMatrix(type="Full",nrow=6,ncol=6,values=C,name=paste("cov",g,sep="")))
		}
	return(P)
	}

PARFCT<- function(P=NULL,titre="PARAFAC à 2 facteurs")
 # P est produit par fct2
	{
	ncond=length(P)
	PF <- NULL
	typ <- "cor"
	fr <- c(FALSE,TRUE,FALSE)
 # définir les éléments à optimiser: 6 dans matPoids, 1 dans matCorr, 2*(ncond-1) dans matAmpl et 6 dans matResid
	matPoids <- mxMatrix(type="Full", nrow=6, ncol=2, values=c(8:3,3:8)/20, free=TRUE, name="A",
		labels=c("F11","F12","F13","F14","F15","F16","F21","F22","F23","F24","F25","F26"),
		lbound=c(0,NA,NA,NA,NA,NA,0,NA,NA,NA,NA,NA))
 	matCorr <- mxMatrix(type="Symm",nrow=2,ncol=2,values=c(1,0,1),labels=c("Un","corr","corr","un"),
		free=FALSE,name="R")
 #		free=c(FALSE,TRUE,TRUE,FALSE),name="R",lbound=c(NA,-.8,-.8,NA),ubound=c(NA,.8,.8,NA))
	matAmpl <- mxMatrix(type="Diag", nrow=2, ncol=2, values=c(1,1), free=FALSE, name="V", lbound=0)  #,labels=c(NA,"V1",NA,"V2"))
	matResid <- mxMatrix(type="Diag", nrow=6, ncol=6, values=1, free=TRUE, name="U",
		labels=c("E1","E2","E3","E4","E5","E6"))
 # contrainte pour empêcher d'interchanger les deux facteurs
	Contr <- mxConstraint(A.F11>A.F21, name = 'Con')
	grp <- NULL
	PFgr <-mxModel(paste(titre,", gr 1 a",ncond))
	for (cond in 1:ncond) {
		labl <- paste(titre,"_",cond,sep="")
		grp <- c(grp,labl)
		matAlg <- mxAlgebra(expression=A %*% V %*% R %*% V %*% t(A) + U, name="RR")
		Exp <- mxExpectationNormal(covariance="RR", dimnames=c("x1","x2","x3","x4","x5","x6"))
		PF1 <- mxModel(labl,matPoids,matCorr,matAmpl,matResid,matAlg,mxFitFunctionML(),Exp,
			mxData(observed=P[[cond]]$values, type=typ, numObs=500))
		PF <- c(PF,PF1)
		PFgr <- mxModel(PFgr,PF[[cond]])
		matAmpl$free<-c(TRUE,FALSE,FALSE,TRUE)
		typ <- TRUE
		typ <- "cov"}
	fun <- mxFitFunctionMultigroup(grp)
 #	PFgr <- mxModel(PFgr,fun,mxCI(c("A","corr","V")))
	PFgr <- mxModel(PFgr,fun,
	mxCI(c("F11","F12","F13","F14","F15","F16","F21","F22","F23","F24","F25","F26")))
	return(PFgr)
	}

vF1 <- c(2.08,3.08)
vF2 <- c(1.74,1.74)
P <- fct2(vF1,vF2)
P2 <- PARFCT(P,"PARAFAC_Example")
P2r<-mxRun(P2)
summary(P2r)