From 31444abc970b7fe17463bcc916e95846272158db Mon Sep 17 00:00:00 2001
From: emilie <emilie@devijver.org>
Date: Fri, 17 Mar 2017 15:07:20 +0100
Subject: [PATCH] few last things

---
 pkg/R/selectiontotale.R              |  20 +-
 pkg/R/valse.R                        |   9 +-
 reports/.Rhistory                    | 512 +++++++++++++++++++++++++++
 reports/bazar Emilie/simulatedData.R |  10 +-
 reports/essai17mars.R                |   0
 reports/simulData_17mars.R           |  60 ++++
 6 files changed, 592 insertions(+), 19 deletions(-)
 create mode 100644 reports/.Rhistory
 create mode 100644 reports/essai17mars.R
 create mode 100644 reports/simulData_17mars.R

diff --git a/pkg/R/selectiontotale.R b/pkg/R/selectiontotale.R
index 7209fed..2cdac38 100644
--- a/pkg/R/selectiontotale.R
+++ b/pkg/R/selectiontotale.R
@@ -29,7 +29,7 @@ selectiontotale = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambd
     selectedVariables = list()
     Rho = list()
     Pi = list()
-    cpt = 0
+    cpt = 1
     #Pour chaque lambda de la grille, on calcule les coefficients
     for (lambdaIndex in 1:length(glambda)){
       print(lambdaIndex)
@@ -39,14 +39,16 @@ selectiontotale = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambd
       m = dim(phiInit)[2]
       #selectedVariables: list where element j contains vector of selected variables in [1,m]
       if (sum(params$phi) != 0){
-        cpt = cpt+1
-      selectedVariables[[cpt]] = sapply(1:p, function(j) {
-        #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector,
-        #and finally return the corresponding indices
-        c(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ], rep(0, m-length(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ] ) ))
-      })
-      Rho[[cpt]] = params$rho
-      Pi[[cpt]] = params$pi
+        selectedVariables[[cpt]] = sapply(1:p, function(j) {
+          #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector,
+          #and finally return the corresponding indices
+          c(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ], rep(0, m-length(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ] ) ))
+        })
+        if (length(unique(selectedVariables)) == length(selectedVariables)){
+          Rho[[cpt]] = params$rho
+          Pi[[cpt]] = params$pi
+          cpt = cpt+1
+        }
       }
     }
     list("selected"=selectedVariables,"Rho"=Rho,"Pi"=Pi)
diff --git a/pkg/R/valse.R b/pkg/R/valse.R
index 72e2d4d..9bc1f4a 100644
--- a/pkg/R/valse.R
+++ b/pkg/R/valse.R
@@ -16,7 +16,7 @@
 #' @export
 #-----------------------------------------------------------------------
 valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 10,
-                 maxi = 50,eps = 1e-4,kmin = 2,kmax = 3,
+                 maxi = 50,eps = 1e-4,kmin = 2,kmax = 2,
                  rang.min = 1,rang.max = 10) {
   ##################################
   #core workflow: compute all models
@@ -45,9 +45,9 @@ valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 1
     source('~/valse/pkg/R/gridLambda.R')
     grid_lambda <<- gridLambda(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, maxi, eps)
     
-    # if (length(grid_lambda)>50){
-    #   grid_lambda = grid_lambda[seq(1, length(grid_lambda), length.out = 50)]
-    # }
+    if (length(grid_lambda)>100){
+      grid_lambda = grid_lambda[seq(1, length(grid_lambda), length.out = 100)]
+    }
     print("Compute relevant parameters")
     #select variables according to each regularization parameter
     #from the grid: A1 corresponding to selected variables, and
@@ -56,7 +56,6 @@ valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 1
     params = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,grid_lambda,X,Y,1e-8,eps)
     #params2 = selectVariables(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,grid_lambda[seq(1,length(grid_lambda), by=3)],X,Y,1e-8,eps)
     ## etrange : params et params 2 sont différents ...
-    
     selected <<- params$selected
     Rho <<- params$Rho
     Pi <<- params$Pi
diff --git a/reports/.Rhistory b/reports/.Rhistory
new file mode 100644
index 0000000..dc970cc
--- /dev/null
+++ b/reports/.Rhistory
@@ -0,0 +1,512 @@
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+Gamma
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) ))
+plot(wi, type='l')
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+# ydata = ydata.noise
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red")
+matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE)
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+eer
+mean(eer)
+source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R')
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+setwd("~/Dropbox/projet spectres/cluster")
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+sims= c(63:67,69,71:74,76:85,88:90,93:100)
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(9, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+sims= 1:100
+colnames(sims) = "isim"
+sims= 1:100
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(10, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",")
+load("~/valse/data/data.RData")
+load("~/valse/data/data.RData")
+X
+Y
+library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+library(roxygen2)
+setwd("~/valse")
+document()
+document()
+setwd("~/")
+install("valse")
+c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+notes = c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+hist(notes)
+hist(notes, nclass = 20)
+notesBis = c(7.375,19.375,10.125,9,6,12.25,11,18.75,10.5,11.5,15.5,13,2.375,12.75,17.75,15.375,8.125,16,8.5,14.5,11.625,11.25,9.625,14,
+6.5,11.375,11.875,16.25,10.125,12,6.25,9.75,8.75,3.5,0,0,0,5.75,2,3.75,6.625,5.25)
+hist(notesBis)
+notes = c(7,
+19,
+10,
+9,
+6,
+12,
+11,
+19,
+11,
+12,
+16,
+13,
+2,
+13,
+18,
+15,
+8,
+16,
+9,
+15,
+12,
+11,
+10,
+14,
+7,
+11,
+12,
+16,
+13,
+15,
+14,
+19,
+14,
+4,  9,
+8,
+7,
+8,
+8
+)
+hist(notes)
+shiny::runApp('Dropbox/cranview-master')
+packages = 'shock'
+min_date = Sys.Date() - 1
+for (pkg in packages)
+{
+# api data for package. we want the initial release - the first element of the "timeline"
+pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all"))
+pkg_data = httr::content(pkg_data)
+initial_release = pkg_data$timeline[[1]]
+min_date = min(min_date, as.Date(initial_release))
+}
+min_date
+library("capushe", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+A = matrix(rnorm(25), ncol = 5)
+A = rbind(A, matrix(0, ncol = 5, nrow = 2))
+A
+A[rowSums(A)==0,]=c()
+A[rowSums(A)==0,]=[]
+A = A[rowSums(A)!=0,]
+A
+source('~/valse/pkg/R/valse.R')
+setwd("~/valse/reports")
+XY = cbind(X,Y)
+X
+p = 10
+q = 8
+k = 2
+D = 20
+meanX = rep(0,p)
+covX = 0.1*diag(p)
+covY = array(dim = c(q,q,k))
+covY[,,1] = 0.1*diag(q)
+covY[,,2] = 0.2*diag(q)
+beta = array(dim = c(p,q,2))
+beta[,,2] = matrix(c(rep(2,(D)),rep(0, p*q-D)))
+beta[,,1] = matrix(c(rep(1,D),rep(0, p*q-D)))
+n = 100
+pi = c(0.4,0.6)
+data = generateXY(meanX,covX,covY, pi, beta, n)
+source('~/valse/pkg/R/generateSampleInputs.R')
+data = generateXY(meanX,covX,covY, pi, beta, n)
+X = data$X
+Y = data$Y
+XY = cbind(X,Y)
+data
+data$class
+affec = data$class
+XY = cbind(X,Y)
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+K= 2
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+XY_class= list()
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+}
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+meanPerClass= matrix()
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+apply(XY_class[[r]], 2, mean)
+p
+q
+meanPerClass[r,] = apply(XY_class[[r]], 2, mean)
+dim(XY)
+meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2])
+for (r in 1:K){
+XY_class[[r]] = XY[affec == r, ]
+meanPerClass[,r] = apply(XY_class[[r]], 2, mean)
+}
+matplot(meanPerClass)
+matplot(meanPerClass, type='l')
diff --git a/reports/bazar Emilie/simulatedData.R b/reports/bazar Emilie/simulatedData.R
index 7eaefb9..6031822 100644
--- a/reports/bazar Emilie/simulatedData.R	
+++ b/reports/bazar Emilie/simulatedData.R	
@@ -10,12 +10,12 @@ library("mclust")
   T = seq(0,1.5,length.out = p)
   T2 = seq(0,3, length.out = 2*p)
   n = 100
-  x1 = cos(2*pi*T) + 0.2*cos(4*2*pi*T) +2*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))
+  x1 = cos(2*base::pi*T) + 0.2*cos(4*2*base::pi*T) + 0.3*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))+1
   plot(T,x1)
   lines(T,x1)
-
-  sigmaX = 0.085
-  sigmaY = 0.1
+  
+  sigmaX = 0.12
+  sigmaY = 0.12
   beta = list()
   p1= 0.5
   beta[[1]] =diag(c(rep(p1,5),rep(1,5), rep(p1,5), rep(1, p-15)))
@@ -94,7 +94,7 @@ for (ite in c(1:ITE)){
   ###########
   ## k-means 1
   ###########
-  mod1 = Mclust(t(XY[ite,,]),G = 2, mode='VII')
+  mod1 = Mclust(t(XY[ite,,]),G = 1:2, mode='VII')
   ARI1[ite] = adjustedRandIndex(mod1$classification, affec[[ite]])
   Kmod1[ite] = mod1$G
   # ###########
diff --git a/reports/essai17mars.R b/reports/essai17mars.R
new file mode 100644
index 0000000..e69de29
diff --git a/reports/simulData_17mars.R b/reports/simulData_17mars.R
new file mode 100644
index 0000000..cdb476c
--- /dev/null
+++ b/reports/simulData_17mars.R
@@ -0,0 +1,60 @@
+simulData_17mars = function(ite){
+  set.seed = 22021989+ite
+  
+  ###########
+  ## Modele
+  ###########
+  K = 2
+  p = 48
+  T = seq(0,1.5,length.out = p)
+  T2 = seq(0,3, length.out = 2*p)
+  n = 100
+  x1 = cos(2*base::pi*T) + 0.2*cos(4*2*base::pi*T) + 0.3*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))+1
+  sigmaX = 0.12
+  sigmaY = 0.12
+  beta = list()
+  p1= 0.5
+  beta[[1]] =diag(c(rep(p1,5),rep(1,5), rep(p1,5), rep(1, p-15)))
+  p2 = 1
+  beta[[2]] = diag(c(rep(p2,5),rep(1,5), rep(p2,5), rep(1, p-15)))
+  ARI1 = ARI2 = ARI3 = 0
+  
+  ###########
+  ## Data + Projection
+  ###########
+  require(wavelets)
+  XY = array(0, dim = c(2*p,n))
+  Xproj = array(0, dim=c(48,n))
+  Yproj = array(0, dim=c(48,n))
+  XYproj = array(0, dim=c(96,n))
+  x = x1 + matrix(rnorm(n*p, 0, sigmaX), ncol = n)
+  affec = sample(c(1,2), n, replace = TRUE)
+  y  = x
+  xy = matrix(0,ncol=n, nrow= 2*p)
+  for (i in c(1:n)){
+    y[,i] = x[,i] %*% beta[[affec[i]]] + rnorm(p, 0, sigmaY)
+    xy[,i] = c(x[,i],y[,i])
+    XY[,i] = xy[,i] - mean(xy[,i])
+    Dx = dwt(x[,i], filter='haar')@W
+    Dx = rev(unlist(Dx))
+    Dx = Dx[2:(1+3+6+12+24)]
+    Ax = dwt(x[,i], filter='haar')@V
+    Ax = rev(unlist(Ax))
+    Ax = Ax[2:(1+3)]
+    Dy = dwt(y[,i], filter='haar')@W
+    Dy = rev(unlist(Dy))
+    Dy = Dy[2:(1+3+6+12+24)]
+    Ay = dwt(y[,i], filter='haar')@V
+    Ay = rev(unlist(Ay))
+    Ay = Ay[2:(1+3)]
+    Xproj[,i] = c(Ax,Dx)
+    Yproj[,i] = c(Ay,Dy)
+    XYproj[,i] = c(Ax,Dx,Ay,Dy)
+  }
+  
+  res_valse = valse(x,y)
+  res_valse_proj = valse(Xproj, Yproj)
+  
+  save(res_valse,file=paste("./Out/Res_",ite, ".RData",sep=""))
+  save(res_valse_proj,file=paste("./Out/ResProj_",ite, ".RData",sep=""))
+}
-- 
2.44.0