library(lattice)

makeToyData<-function(n,
                      parameters,as.pvalues=TRUE,
                      correlated=TRUE,
                      correlatedInformative=TRUE,
                      pNoize=0,pNoizeBiMod=0
                      )
  {
    n1=n2=floor(n/2)
    
    class=c(rep('X',n1),rep('O',n2))
    data=data.frame(class)

    ## x1.true, x2.true > x.true are the real distribution from the parameters
    for(param in parameters)
      {
        if(!as.pvalues)
          {
            moy=param
            x1.true=rnorm(n1,mean=-moy, sd=1)
            x2.true=rnorm(n2,mean= moy, sd=1)
            x.true=c(x1.true,x2.true)
            preffix=paste('moy',param,sep='')
          }
        else
          {
            epsilon=param/10
            moy=getMean(p.val=param,n=(n1+n2))
            while(TRUE)
              {
                x1.true=rnorm(n1,mean=-moy)
                x2.true=rnorm(n2,mean=+moy)
                pval=t.test(x1.true, x2.true)$p.value
                if (abs(param-pval) < epsilon) break
              }
            x.true=c(x1.true,x2.true)
            preffix=paste('pval',param,sep='')
          }
        x = x.true+rnorm(length(x.true), sd=0.1*sd(x.true))
        x=data.frame(x)
        names(x)=preffix
        data=cbind(data,x)
        
        ## Correlated feature
        if(correlated)
          {
            xc = x.true+rnorm(length(x.true), sd=0.1*sd(x.true))
            xc=data.frame(xc)
            names(xc)=paste(preffix,'c',sep='')
            data=cbind(data,xc)
          }
        if(correlatedInformative)
          {
            ## Correlated but informative feature
            x1ci = x1.true+moy+rnorm(length(x1.true), sd=moy/2)##0.1*sd(x1.true))
            x2ci = x2.true-moy+rnorm(length(x2.true), sd=moy/2)##0.1*sd(x2.true))
            xci=c(x1ci,x2ci)
            xci=data.frame(xci)
            names(xci)=paste(preffix,'ci',sep='')
            data=cbind(data,xci)
          }
      }
    ## Make some noize
    if(pNoize!=0)
      data=makeNoize(data,n,pNoize)
    if(pNoizeBiMod!=0)
      data=makeNoize(data,n,pNoizeBiMod,bimod=TRUE)

    return( data )
  }

## Utils
makeNoize<-function(data,n,p,bimod=FALSE)
  {
    for(i in 1:p)
      {
        if(!bimod)
          {
            x=rnorm(n,sd=sqrt(20))
            preffix='normNoize'
          }
        else
          {
            x1=rnorm(n/2,mean=-3*sqrt(20),sd=sqrt(20))
            x2=rnorm(n/2,mean= 3*sqrt(20),sd=sqrt(20))
            x=sample(c(x1,x2))
            preffix='biNormNoize'
          }
        x=data.frame(x)
        names(x)=paste(preffix,i,sep='')
        data=cbind(data,x)
      }
    return(data)
  }

## Return the mean (mu) for distribution N(-mu,sd=1) N(mu,sd=1)
getMean<-function(p.val,n)
  {
    epsilon=p.val/10
    n1=n2=floor(n/2)
    moy=1
    while(TRUE)
      {
        sum.p.val=0
        for(i in 1:100)
          {
            sum.p.val=sum.p.val+
              t.test(rnorm(n1,mean=-moy), rnorm(n2,mean=moy))$p.value
          }
        p.val.esti=sum.p.val/100
        
        cat('moy=',moy,'\t p.val.estimated=',p.val.esti,'\n')

        if (abs(p.val.esti-p.val) < epsilon) break
        if(p.val.esti<p.val)moy=0.95*moy
        else                moy=1.05*moy
      }
    return(moy)
  }

## ----------------------------------------------------------------
## MAIN
##data=makeToyData(n=100,p=5,pNoize=10,by=0.1)


##getMean(p.val=1e-01,n=150)
##0.1840259
##0.1840259

##getMean(p.val=1e-02,n=150)
##0.2905309
##0.291989

##getMean(p.val=1e-03,n=150)
##0.357135
##0.3620782

##getMean(p.val=1e-04,n=150)
##0.4390264
##0.4249701
##0.487675

##getMean(p.val=1e-05,n=150)
##0.5403601
##0.4864558
##0.5403601
##0.5019067

##getMean(p.val=1e-06,n=150)
##0.5126893
##0.539673
##d=makeToyData(n=150,parameters=c(1.2e-6, 5.5e-5, 1.3e-4, 1.9e-4, 4.4e-4, 4.60e-4, 4.61e-4, 7e-4, 1.3e-3, 2.3e-3, 2.4e-3), as.pvalues=TRUE,correlatedInformative=FALSE)

##d=makeToyData(n=150,parameters=c(1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.98e-1, 0.99e-1, 1e-1, 1.01e-1, 1.02e-1), as.pvalues=TRUE, pNoize=30)
## d=makeToyData(n=50,parameters=c(1e-2),as.pvalues=TRUE,correlatedInformative=TRUE)
## gg(d,g='class')