################################################################
# R source code for random  generation of multivariate data
# February 2005
# Modified August 2005
################################################################
library(stats)
library(MASS)

###########################################################
# Random generation of multivariate data normally distributed
###########################################################

# Random generation of normal distributed data with diagonal covariance matrix
# n : number of samples to be generated
# mean : vector of means
# sd : vector of standard deviations 
# seed : intialization of the random generator. If 0 (def.) no initialization is performed
# Output:
# a matrix of n columns with length (mean) rows. Row[i] has mean mean[i] and standard deviation sd[i].
rand.norm.generate <- function(n=5, mean=0, sd=1, seed=0){
  if (seed != 0)
	  set.seed(seed);
  len.mean <- length(mean);
  len.sd <- length(sd);
  if (len.mean!=len.sd)
    stop("rand.norm.generate: Length of mean and sd vector must match", call.=FALSE);
  len.vect <- n*len.mean;
  m <- rnorm(len.vect, mean, sd);
  m <- matrix(m, nrow=len.mean);
  m
}


# Random generation of normal distributed data with full covariance matrix
# n : number of samples to be generated
# mean : vector of means
# Sigma : Covariance matrix 
# seed : intialization of the random generator. If 0 (def.) no initialization is performed
# Output:
# a matrix m of n columns with length(mean) rows. Row[i] has mean mean[i].
rand.norm.generate.full <- function(n=5, mean=c(0,0), Sigma=matrix(c(0.1,0,0 ,0.1),2,2), seed=0){
  if (seed != 0)
	  set.seed(seed);
  len.mean <- length(mean);
  if (len.mean!=nrow(Sigma))
    stop("rand.norm.generate.full: Length of mean and number of rows of Sigma must match", call.=FALSE);
  if (len.mean!=ncol(Sigma))
    stop("rand.norm.generate.full: Length of mean and number of columns of Sigma must match", call.=FALSE);
  m <- t(mvrnorm(n, mean, Sigma));
  m
}



# Sample0 generator: multivariate normally  distributed data synthetic generator
# n examples for each from 3 classes are generated.
# All classes (each one of n examples) has dim components 
# The first class (first n examples) has its  components centered in  0 (of length dim)
# The second class (second n examples) has its   components centered in  m (of length dim)
# The third class (last n examples) has its   components centered in -m (of length dim)
# For all classes the covariance matrix is diagonal with values sigma
# Input:
# n : number of examples for each class
# m : mean value for the second class
# sigma : value of the diagonal elements of the covariance matrix
# dim : dimension of the examples
# Output:
# a matrix with dim rows (variables) and n*3 columns (examples)
generate.sample0 <- function(n=5, m=10, sigma=1, dim=2){
  
  sd <- rep(sigma,dim);
  mean <- rep(0,dim); 
  m1 <- rand.norm.generate(n, mean, sd);
  
  mean <- rep(m,dim); 
  m2 <- rand.norm.generate(n, mean, sd);
  
  mean <- rep(-m,dim); 
  m3 <- rand.norm.generate(n, mean, sd); 
  
  m <- cbind(m1,m2,m3);
  m  
}


# Sample1 generator: multivariate normally  distributed data synthetic generator
# n dim-dimensional examples for each class are generated.
# All classes (each one of n examples) have their last dim-500 variables centered in 0.
# The first class (first n examples) has its first 500  features centered in 0
# The second class (second n examples) has its first 500 features centered in  m
# The third class (last n examples) has its first 500 features centered in  -m
# For all classes the covariance matrix is diagonal with all values on the diagonal equal to sigma 
# Input:
# n : number of examples for each class
# m : center of the first 500 variables of the second class
# dim : number of variables (features)
# Output:
# a matrix with dim rows (variables) and n*3 columns (examples)
generate.sample1 <- function(n=2, m=6, sigma=1, dim=10000){
  
	if (dim<500) 
	  stop("generate.sample1 need at least 500 variables");
	
  mean <- rep(0,dim);
  sd <- rep(sigma,dim);
  m1 <- rand.norm.generate(n, mean, sd);
  
  mean <- rep(c(m,0),c(500,dim-500));
  m2 <- rand.norm.generate(n, mean, sd);
  
  mean <- rep(c(-m,0),c(500,dim-500));
  m3 <- rand.norm.generate(n, mean, sd); 
  
  m <- cbind(m1,m2,m3);
  m  
}


# Sample2: multivariate normally  distributed data synthetic generator
# n 10000-dimensional examples for each class are generated.
# All classes (each one of n examples) has only no-noisy features but there is substantial overlapping between classes
# The first class (first n examples) has its features centered in 1 (first 5000 features) and 2 (last 5000 features)
# The second class (second n examples) has its features centered in -1 (first 5000 features) and -2 (last 5000 features)
# The diagonal of the covariance matrix of the first class has its first 2500 element equal to 0.5, the next 2500 equal to 1,
# the next 2500 to 0.5 and the last to 1.
# The diagonal of the covariance matrix of the second class has its first 5000 element equal to 1, the next 5000 equal to 2
# Input:
# n : number of examples for each class
# Output:
# a matrix with 10000 rows (variables) and n*2 columns (examples)
generate.sample2 <- function(n=2){
 
  mean <- rep(c(1,2),c(5000,5000));
  sd <- rep(c(0.5,1,0.5,1),rep(2500,4));
  m1 <- rand.norm.generate(n, mean, sd);
  
  sd <- rep(c(1,2),c(5000,5000));
  m2 <- rand.norm.generate(n, -mean, sd);
  
  m <- cbind(m1,m2);
  m  
}


# Sample3: multivariate normally  distributed data synthetic generator
# n 1000-dimensional examples for each class are generated.
# All classes (each one of n examples) has 300 no-noisy features and 700 noisy features. There is a certain overlap 
# between classes and a full covariance matrix (equal for all classes is used).
# The first class (first n examples) has its no-noisy features centered in 0.
# The second class (second n examples) has its no-noisy features centered in m 
# The third class (last n examples) has its no-noisy features centered in -m 
# Covariance matrix Sigma = (B, Zero; Zero', I) where B is a 300X300 matrix s.t. B[i,i]=1, B[i,i+1]=B[i,i-1]=0.5 and
# B[i,j]=0.1 j!=i-1,i,i+1; Zero is  a 300X700 zero matrix and Zero' its transpose; I is a 700X700 identity matrix
# Input:
# n : number of examples for each class
# Output:
# a matrix with 1000 rows (variables) and n*3 columns (examples)
generate.sample3 <- function(n=2, m=2){
  if (n < 2)
    stop("generate.sample3: You must generate at least 2 samples per class", call.=FALSE);
  Sigma <- matrix(rep(0,1000*1000), nrow=1000);
  B <- matrix(rep(1,300*300),nrow=300)*0.1;
  i<-1:300;
  B[cbind(i,i)]<-1;
  i<-1:299;
  B[cbind(i,i+1)]<-0.5;
  i <- 2:300;
  B[cbind(i,i-1)]<-0.5;
  Sigma[1:300,1:300] <- B;
  Sigma[cbind(301:1000,301:1000)] <- 1;
  
  mean <- rep(0,1000);
  m1<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- rep(c(m,0),c(300,700));
  m2<-rand.norm.generate.full(n, mean, Sigma);
  
  m3<-rand.norm.generate.full(n, -mean, Sigma);
  
  m <- cbind(m1,m2,m3);
  m
}


# Sample4: multivariate normally  distributed data synthetic generator
# n 6000-dimensional examples for each class (5 classes) are generated.
# All classes (each one of n examples) has 1000 no-noisy and 5000 noisy features but there is substantial overlapping 
# between classes 1 and 2 and 1 and 3, while class 4 and 5 are separated.
# The first class (first n examples) has its no noisy variables centered in 0.
# The second class (second n examples) has its no noisy variables centered in 1. 
# The third class (third n examples) has its no noisy variables centered in -1. 
# The fourth class (fourth n examples) has its no noisy variables centered in 5.
# The fifth class (fifth n examples) has its no noisy variables centered in -5.
# The diagonal of the covariance matrix for all classes has its elements equal to sigma (first 1000 variables) and equal to
# 2*sigma (last 5000 variables)
# # Input:
# n : number of examples for each class
# sigma : standard deviation of the first 1000 variables. The remaining variables have 2*sigma standard deviation
# Output:
# a matrix with 1000 rows (variables) and n*5 columns (examples)
generate.sample4 <- function(n=2, sigma=1){
 
  sd <- rep(c(sigma,2*sigma),c(1000,5000));
  
  mean <- rep(0,6000);
  m1 <- rand.norm.generate(n, mean, sd);
  
  mean <- rep(c(1,0),c(1000,5000));
  m2 <- rand.norm.generate(n, mean, sd);
  
  m3 <- rand.norm.generate(n, -mean, sd);
  
  mean <- rep(c(5,0),c(1000,5000));
  m4 <- rand.norm.generate(n, mean, sd);
  
  m5 <- rand.norm.generate(n, -mean, sd);
  
  m <- cbind(m1,m2,m3,m4,m5);
  m  
}


##########################################################
# Sample5: multivariate normally  distributed data synthetic generator
# n dim-dimensional examples for each class are generated.
# All classes (each one with n examples) has (1-ratio.noisy)*dim of no-noisy features and ratio.noisy*dim of noisy features. 
# Note that if the number on no-noisy feature is less than 2 the generation is aborted.
# A full covariance matrix (equal for all classes) is used.
# The first class (first n examples) has its no-noisy features centered in 0.
# The second class (second n examples) has its no-noisy features centered in m 
# The third class (third n examples) has its no-noisy features centered in -m 
# A fourth cluster (third n examples) has its no-noisy features centered in (m,-m) alternatively
# Covariance matrix Sigma = (B, Zero; Zero', I) where B is a (dim*(1-ratio.noisy))X(dim*(1-ratio.noisy)) matrix s.t. 
# B[i,i]=1, B[i,i+1]=B[i,i-1]=0.5 and
# B[i,j]=0.1 if j!=i-1,i,i+1; Zero is  a (dim*(1-ratio.noisy))X(dim*ratio.noisy) zero matrix and Zero' its transpose; 
# I is a (dim*ratio.noisy)X(dim*ratio.noisy) identity matrix
# Input:
# n : number of examples for each class
# dim : dimension of the examples (minimum=10).
# ratio.noisy : ratio of the noisy variables
# m : center of the II cluster (the third has center -m).
# Output:
# a matrix with dim rows (variables) and n*4 columns (examples)
generate.sample5 <- function(n=10, dim=10, ratio.noisy=0.8, m=2 ){
 
  if (n < 1)
    stop("generate.sample5: You must generate at least 1 sample per class", call.=FALSE);

	no.noisy.var <- round(dim*(1-ratio.noisy));
	noisy.var <- dim - no.noisy.var;
	if (no.noisy.var<2)
	  stop("generate.sample5: ratio.noisy parameter must be small enough to generate at least 2 no noisy variables", call.=FALSE);
	t <- no.noisy.var * no.noisy.var;
	x <- rep(0.1,times=t);
  B <- matrix(x,nrow=no.noisy.var);
  i<-1:no.noisy.var;
  B[cbind(i,i)]<-1;
  i<-(1:no.noisy.var-1);
  B[cbind(i,i+1)]<-0.5;
  i <- (2:no.noisy.var);
  B[cbind(i,i-1)]<-0.5;
	
	Sigma <- matrix(rep(0,dim*dim), nrow=dim);
	Sigma[1:no.noisy.var,1:no.noisy.var] <- B; 
  Sigma[cbind((no.noisy.var+1):dim,(no.noisy.var+1):dim)] <- 1;
  
  mean <- rep(0,dim);
  m1<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- rep(c(m,0),c(no.noisy.var,noisy.var));
  m2<-rand.norm.generate.full(n, mean, Sigma);
  		
  m3<-rand.norm.generate.full(n, -mean, Sigma);
	
	if (no.noisy.var%%2 == 0)	
	  mean <- c(rep(c(m,-m),floor(no.noisy.var/2)),rep(0,noisy.var))
	else
	  mean <- c(rep(c(m,-m),floor(no.noisy.var/2)),m,rep(0,noisy.var));
	m4<-rand.norm.generate.full(n, mean, Sigma);
  
  M <- cbind(m1,m2,m3,m4);
  M
}

###############################################################################
# Sample6 generator: multivariate normally  distributed data synthetic generator
# n examples for each from 6 classes are generated.
# All classes (each one of n examples) has dim components 
# The clusters have a hierarchical structure: 2 or 6 clusters may be detected.
# Anyway note that the structure of the data depends on the parameters: two main clusters
# are centered in m and -m. Around each main cluster three other subclusters are generated using
# the displacement d.
# Input:
# n : number of examples for each class
# m : mean basic value 
# d : amount of the displacement from m
# dim : dimension of the examples
# s: value of the diagonal elements of the covariance matrix
# Output:
# a matrix with dim rows (variables) and n*6 columns (examples)
generate.sample6 <- function(n=20, m=10, dim=2, d=3, s=0.2){
  
  sd <- rep(s,dim);
  mean.base <- rep(m,dim); 
	range <- 1:floor(dim/2);	
	range2 <- (floor(dim/2)+1):dim;
	
	mean <- mean.base;
	mean[range] <- mean[range] - d;	
  m1 <- rand.norm.generate(n, mean, sd);
  
  mean <- mean.base;
	mean[range] <- mean[range] + d;	
  m2 <- rand.norm.generate(n, mean, sd);
  
  mean <- mean.base;
	mean[range2] <- mean[range2] - d;	
  m3 <- rand.norm.generate(n, mean, sd);
	
	mean <- -mean.base;
	mean[range] <- mean[range] - d;	
  m4 <- rand.norm.generate(n, mean, sd);
	
	mean <- -mean.base;
	mean[range] <- mean[range] + d;	
  m5 <- rand.norm.generate(n, mean, sd);
  
  mean <- -mean.base;
	mean[range2] <- mean[range2] - d;	
  m6 <- rand.norm.generate(n, mean, sd);
  
  M <- cbind(m1,m2,m3,m4,m5,m6);
	
	return(M);
}

########################################################



  generate.toy1 <- function(n=100, m=10) {
  mean <- c(0,0);
	Sigma <- matrix(c(2,0,0,1),nrow=2);
  m1<-rand.norm.generate.full(n*4, mean, Sigma);
  
  mean <- c(m,m);
	Sigma <- matrix(c(1,0,0,1),nrow=2);
  m2<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- c(-m,m);
	Sigma <- matrix(c(1,0,0,1),nrow=2);
  m3<-rand.norm.generate.full(n, mean, Sigma);
	
	mean <- c(-m,-m);
	Sigma <- matrix(c(1,0,0,1),nrow=2);
  m4<-rand.norm.generate.full(n, mean, Sigma);
	
	mean <- c(m,-m);
	Sigma <- matrix(c(1,0,0,1),nrow=2);
  m5<-rand.norm.generate.full(n, mean, Sigma);
	
	M <- cbind(m1,m2,m3,m4,m5);
	
	return(M);
	
}

########################################################
  generate.toy2 <- function(n=20, m=10, d=4, s=0.5) {
  mean <- c(m-d,m);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m1<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- c(m+d,m);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m2<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- c(m,m-d);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m3<-rand.norm.generate.full(n, mean, Sigma);
	
	mean <- c(-m-d,-m);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m4<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- c(-m+d,-m);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m5<-rand.norm.generate.full(n, mean, Sigma);
  
  mean <- c(-m,-m-d);
	Sigma <- matrix(c(s,0,0,s),nrow=2);
  m6<-rand.norm.generate.full(n, mean, Sigma);
	
	M <- cbind(m1,m2,m3,m4,m5,m6);
	
	return(M);
	
}



###############################################################################
# Sample7 generator: multivariate normally  distributed data synthetic generator
# n examples for each from 6 classes are generated.
# All classes (each one of n examples) has dim components 
# The clusters have a hierarchical structure: 2 or 6 clusters may be detected.
# Anyway note that the structure of the data depends on the parameters: two main clusters
# are centered in m and -m. Around each main cluster two other subclusters are generated using
# the displacement d.
# Input:
# n : number of examples for each class
# m : mean basic value 
# d : amount of the displacement from m
# dim : dimension of the examples
# s: value of the diagonal elements of the covariance matrix
# Output:
# a matrix with dim rows (variables) and n*6 columns (examples)
generate.sample7 <- function(n=20, m=10, dim=1000, d=3, s=1){
  
  sd <- rep(s,dim);
  mean.base <- rep(m,dim); 
	range <- 1:floor(dim/2);	
	range2 <- (floor(dim/2)+1):dim;
	
	mean <- mean.base;	
  m1 <- rand.norm.generate(n, mean, sd);
  
  mean <- mean.base;
	mean[range] <- mean[range] - d;	
	mean[range2] <- mean[range2] + d;
  m2 <- rand.norm.generate(n, mean, sd);
  
  mean <- mean.base;
	mean[range] <- mean[range] - d;	
	mean[range2] <- mean[range2] - d;	
  m3 <- rand.norm.generate(n, mean, sd);
	
	mean <- -mean.base;
  m4 <- rand.norm.generate(n, mean, sd);
	
	mean <- -mean.base;
	mean[range] <- mean[range] - d;	
	mean[range2] <- mean[range2] + d;	
  m5 <- rand.norm.generate(n, mean, sd);
  
  mean <- -mean.base;
	mean[range] <- mean[range] - d;	
	mean[range2] <- mean[range2] - d;		
  m6 <- rand.norm.generate(n, mean, sd);
  
  M <- cbind(m1,m2,m3,m4,m5,m6);
	
	return(M);
}


#########################################
# Hierarchical cluster generator.
# A 2-dimensional two-level hierarchical cluster structure is generated.
# At a first level 3 distinct clusters at the vertices of an equilater tringle are generated.
# At a second level two other clusters at the left and right of the three "primary" clusters are generated.
# Input:
# n : number of examples for each cluster
# l: half length of the edge of the equilater triangle
# Delta.h : half of the "abscissa" distance between each pair of clusters inside the three major clusters
# sd: standard deviation
# with.I.level.examples: if true data centered at the vertices of the triangle are generated, otherwise only the secondary clusters are generated.
generate.sample.h1 <- function(n=20, l=5, Delta.h=1, sd=0.1, with.I.level.examples=FALSE){

	Center<-5;

	A <- numeric(2);
	B <- numeric(2);
	C <- numeric(2);
	sd <- rep(sd,2);
	sdI <- sd * 4;

	A[1] <- Center + l;
	B[1] <- Center - l;

	d <- l/cos(pi/6);
	y <- sin(pi/6) * d;

	A[2] <- B[2] <- Center + y;

	C <- c(Center,5 - y);
	
	A1 <- A2 <- A;
	B1 <- B2 <- B;
	C1 <- C2 <- C;
	A1[1] = A1[1] - Delta.h;
	B1[1] = B1[1] - Delta.h;
	C1[1] = C1[1] - Delta.h;
	A2[1] = A2[1] + Delta.h;
	B2[1] = B2[1] + Delta.h;
	C2[1] = C2[1] + Delta.h;
  
	if (with.I.level.examples==TRUE) {
	  mA <- rand.norm.generate(n/2, A, sdI);
	  mB <- rand.norm.generate(n/2, B, sdI);
	  mC <- rand.norm.generate(n/2, C, sdI);
	}
	
	mA1 <- rand.norm.generate(n, A1, sd);
  mB1 <- rand.norm.generate(n, B1, sd);
  mC1 <- rand.norm.generate(n, C1, sd);
	mA2 <- rand.norm.generate(n, A2, sd);
  mB2 <- rand.norm.generate(n, B2, sd);
  mC2 <- rand.norm.generate(n, C2, sd);
	
	if (with.I.level.examples==TRUE) 
	  M <- cbind(mA1,mA2,mB1,mB2,mC1,mC2,mA,mB,mC)
	else	
	  M <- cbind(mA1,mA2,mB1,mB2,mC1,mC2);
	return (M);
}


#########################################
# Three-level hierarchical cluster generator.
# A 2-dimensional three-level hierarchical cluster structure is generated.
# At a first level 3 distinct clusters at the vertices of an equilater triangle are generated.
# At a second level two other clusters at the left and right of the three "primary" clusters are generated (6 clusters)
# At a third level two other clusters above and below the secondaty clusters are generated (12 clusters)
# Input:
# n : number of examples for each cluster
# l: half length of the edge of the equilater triangle
# Delta.h : half of the "abscissa" distance between each pair of clusters inside the three major clusters
# Delta.v : half of the "ordinate" distance between each pair of clusters inside the three second order clusters
# sd: standard deviation
# with.I.II.level.examples: if TRUE data at the first and secondary level are generated (for a total of 21 clusters), otherwise only the third level
# clusters are generated.
generate.sample.h2 <- function(n=20, l=8, Delta.h=2,  Delta.v=1, sd=0.1, with.I.II.level.examples=FALSE){

	Center<-5;

	A <- numeric(2);
	B <- numeric(2);
	C <- numeric(2);
	sd <- rep(sd,2);
	sdI <- c(Delta.h/2,Delta.h/2);
	sdII <- c(Delta.v/2,Delta.v/2);

	A[1] <- Center + l;
	B[1] <- Center - l;

	d <- l/cos(pi/6);
	y <- sin(pi/6) * d;

	A[2] <- B[2] <- Center + y;

	C <- c(Center,5 - y);
	
	A1 <- A2 <- A;
	B1 <- B2 <- B;
	C1 <- C2 <- C;
	A1[1] = A1[1] - Delta.h;
	B1[1] = B1[1] - Delta.h;
	C1[1] = C1[1] - Delta.h;
	A2[1] = A2[1] + Delta.h;
	B2[1] = B2[1] + Delta.h;
	C2[1] = C2[1] + Delta.h;
	
	A1.1 <- A1.2 <- A1;
	A2.1 <- A2.2 <- A2;
	B1.1 <- B1.2 <- B1;
	B2.1 <- B2.2 <- B2;
	C1.1 <- C1.2 <- C1;
	C2.1 <- C2.2 <- C2;
	
	A1.1[2] = A1.1[2] - Delta.v;
	B1.1[2] = B1.1[2] - Delta.v;
	C1.1[2] = C1.1[2] - Delta.v;
	A2.1[2] = A2.1[2] - Delta.v;
	B2.1[2] = B2.1[2] - Delta.v;
	C2.1[2] = C2.1[2] - Delta.v;
	
	A1.2 [2]= A1.2[2] + Delta.v;
  B1.2 [2]= B1.2[2] + Delta.v;
  C1.2 [2]= C1.2[2] + Delta.v;
  A2.2 [2]= A2.2[2] + Delta.v;
  B2.2 [2]= B2.2[2] + Delta.v;
  C2.2 [2]= C2.2[2] + Delta.v;
	
	mA1.1 <- rand.norm.generate(n, A1.1, sd);
	mA1.2 <- rand.norm.generate(n, A1.2, sd);
	mA2.1 <- rand.norm.generate(n, A2.1, sd);
  mA2.2 <- rand.norm.generate(n, A2.2, sd);
  mB1.1 <- rand.norm.generate(n, B1.1, sd);
  mB1.2 <- rand.norm.generate(n, B1.2, sd);
  mB2.1 <- rand.norm.generate(n, B2.1, sd);
  mB2.2 <- rand.norm.generate(n, B2.2, sd);
	mC1.1 <- rand.norm.generate(n, C1.1, sd);
	mC1.2 <- rand.norm.generate(n, C1.2, sd);
	mC2.1 <- rand.norm.generate(n, C2.1, sd);
  mC2.2 <- rand.norm.generate(n, C2.2, sd);
  
	if (with.I.II.level.examples==TRUE) {
	  mA <- rand.norm.generate(n/4, A, sdI);
	  mB <- rand.norm.generate(n/4, B, sdI);
	  mC <- rand.norm.generate(n/4, C, sdI);
		mA1 <- rand.norm.generate(n/2, A1, sdII);
		mB1 <- rand.norm.generate(n/2, B1, sdII);
		mC1 <- rand.norm.generate(n/2, C1, sdII);
		mA2 <- rand.norm.generate(n/2, A2, sdII);
		mB2 <- rand.norm.generate(n/2, B2, sdII);
		mC2 <- rand.norm.generate(n/2, C2, sdII);    
	}
	
	if (with.I.II.level.examples==TRUE) 
	  M <- cbind(mA1.1,mA1.2,mA2.1,mA2.2,mB1.1,mB1.2,mB2.1,mB2.2,mC1.1,mC1.2,mC2.1,mC2.2,mA1,mA2,mB1,mB2,mC1,mC2,mA,mB,mC)
	else	
	  M <- cbind(mA1.1,mA1.2,mA2.1,mA2.2,mB1.1,mB1.2,mB2.1,mB2.2,mC1.1,mC1.2,mC2.1,mC2.2);
	return (M);
}


#########################################
# Two-levels hierarchical cluster generator.
# A 2-dimensional 2-levels hierarchical cluster structure is generated.
# At a first level 4 distinct clusters  are generated.
# At a second level two other clusters at the left and right of 2 of the 4 "primary" clusters are generated (6 clusters)
# Input:
# n : number of examples for each cluster
# DeltaA : vertical displacement of the the secondary clusters
# DeltaB : horizontal displacement of the the secondary clusters
# seed: seed for the random generator
generate.sample.h3 <- function(n=20, DeltaA=1, DeltaB=1, seed=0){
  if (seed != 0)
	  set.seed(seed);
  # I level cluster centers and standard dev.
	A <- c(2,8);
	sigmaA <- c(0.2,0.8);
	B <- c(7,7);
	sigmaB <- c(0.5,0.2);
	C <- c(6,3);
	SigmaC <- matrix(c(0.6,0.4,0.4,0.3),2,2);
	D <- c(4,5);
	sigmaD <- c(0.3,0.2);
	
  # I level cluster generation
	mA <- rand.norm.generate(n/2, mean=A, sd=sigmaA);
	mB <- rand.norm.generate(round(n/3), mean=B, sd=sigmaB);
	mC <- rand.norm.generate.full(2*n, mean=C, Sigma=SigmaC);
  mD <- rand.norm.generate(n/2, mean=D, sd=sigmaD);
	
	# II level cluster centers and stdev.
	A1 <- c(A[1],A[2]+DeltaA);
	A2 <- c(A[1],A[2]-DeltaA);
	sigmaA <- c(0.2,0.2);
	B1 <- c(B[1]+DeltaB,B[2]);
	B2 <- c(B[1]-DeltaB,B[2]);
	sigmaB <- c(0.2,0.2);
	
	 # I level cluster generation
	mA1 <- rand.norm.generate(n, mean=A1, sd=sigmaA);
	mA2 <- rand.norm.generate(n, mean=A2, sd=sigmaA);
	mB1 <- rand.norm.generate(n, mean=B1, sd=sigmaB);
  mB2 <- rand.norm.generate(n, mean=B2, sd=sigmaB);
	
	M <- cbind(mA,mB,mC,mD,mA1,mA2,mB1,mB2)
	return (M);
}


########################################
# Uniform bidimensional random data generator.
# Data are generated according to a uniform bidimensional random distribution. 
# Input:
# n : number of examples 
# range : vector with 2 values: min and max random uniform values
generate.uniform.random <- function(n=100, range=c(0,1)){
  x <- range[1] + runif(n) * (range[2]-range[1]);
	y <- range[1] + runif(n) * (range[2]-range[1]);
	M <- rbind(x,y);
	return(M);
}

########################################
# Uniform bidimensional  data generator
# Data are generated according to  bidimensional grid with equispatiated data.
# Input:
# n : number of examples 
# range : vector with 2 values: min and max  coordinates of the bidimensional grid
generate.uniform <- function(n=11, range=c(0,1)){
  x <- seq(from=range[1], to=range[2], length=n);
	M = matrix(numeric(n*n*2),nrow=2);
	for (i in 0:(n-1))
	  for (j in 1:n) {
		  M[1,i*n + j] <- x[j];
			M[2,i*n + j] <- x[i+1];
	  }
	return(M);
}