###########################################################
#    R code for ISOSIR (ISOSIR-Source-Rcode.R)            #
#     - [ ] ISOSIR-Source-Rcode.R                         #
#     - [x] ISOSIR-SimuData-Rcode.R                       #
#     - [ ] ISOSIR-Example-Rcode.R                        #
#                                                         #
#    Han-Ming Wu (hmwu@mail.tku.edu.tw)                   #
#    Department of Mathematics, Tamkang University        #
#                                                         #
#    For update or correctness,                           #
#    please visit at http://www.hmwu.idv.tw               #
#                                                         #
#    Version 2011-07-31, Initial Release                  #
###########################################################
###########################################################
## Generate some manifold data                            #
## Author: Han-Ming Wu                                    #
## E-mail: hmwu@mail.tku.edu.tw                           #
## Date: 2011/05/30                                       #
###########################################################

library(rgl)
library(scatterplot3d)
library(rpanel)

###########################################################
## swiss roll                                             #
###########################################################
swissroll <- function(n, sigma=0.05){
	
	 
	angle <- (3*pi/2)*(1+2*runif(n)); 
	height <- runif(n);
	xdata <-  cbind(angle*cos(angle), height, angle*sin(angle))
	
	#angle <-  seq((1.5*pi): (4.5*pi), length.out=n);
	#height <- sample(0:21, n, replace = TRUE);
	#xdata <-  cbind(angle*cos(angle), height, angle*sin(angle))
	
	xdata <- scale(xdata) + matrix(rnorm(n*3, 0, sigma), n, 3)

	order.angle <- order(angle)
	sort.angle <- sort(order.angle, index.return=TRUE)	
	col.id <- rainbow130(n)
	my.color <- col.id[sort.angle$ix]
	
	colnames(xdata) <- paste("x", 1:3, sep="")

    return(list(xdata=xdata, angle=angle, color=my.color))
}



isosir.swissroll <- function(n, sigma=0.05){
	
	 
	spread <- runif(n);
	angle <- (3*pi/2)*(1+2*spread); 
	xdata <-  cbind(angle*cos(angle), spread, angle*sin(angle))
	xdata <- scale(xdata) + matrix(rnorm(n*3, 0, sigma), n, 3)

	order.angle <- order(angle)
	sort.angle <- sort(order.angle, index.return=TRUE)	
	col.id <- rainbow130(n)
	my.color <- col.id[sort.angle$ix]
	
	colnames(xdata) <- paste("x", 1:3, sep="")

    return(list(xdata=xdata, angle=angle, color=my.color))
}

#n <- 500; 
#swissdata <- isosir.swissroll(n)
#xdata <- swissdata$xdata
#x.color <- swissdata$color
#plot3d(xdata[,1], xdata[,2], xdata[,3], col=x.color, size=5, xlab="", ylab="", zlab="", axes = FALSE)	         
#rgl.postscript("Fig3-swiss-spiral.eps", "eps", drawText=FALSE)



###########################################################
# 6. Noisy data (only swissroll, S-curve)
###########################################################
swissroll.addnoise <- function(sw.obj, m){
    
    n <- nrow(sw.obj$xdata)
    xndata <- matrix(0, ncol=3+m, nrow=n)    
	xndata[,1:3] <- sw.obj$xdata
    for(i in 1:m){
        xndata[, 3+i] <- as.matrix(rnorm(n, mean=0, sd=1))
    }	
	colnames(xndata) <- paste("x", 1:(3+m), sep="")
    
    return(list(xdata=xndata, angle=sw.obj$angle, color=sw.obj$color))

}

###########################################################
## S-curve                                                #
###########################################################
Scurve <- function(n){

    ##############################
    #     upper half S curve     #
    ##############################
    theta1 <- runif((n/2), min=-0.9*(pi), max=(pi)/2)
    z1 <- runif((n/2), min=0, max=5)
    x1 <- -cos(theta1)
    y1 <- 2-sin(theta1)
    side.upper <- matrix(0, ncol=4, nrow=(n/2))
    for(i in 1:(n/2)){
        side.upper[i,1] <- x1[i]
        side.upper[i,2] <- y1[i]
        side.upper[i,3] <- z1[i]
        side.upper[i,4] <- theta1[i]
    }
    order.theta1 <- order(theta1)
    sort.theta1 <- sort(order.theta1, method="qu", index=TRUE)
    upper.color <- sort.theta1$ix

    ##############################
    #     lower half S curve     #
    ##############################
    theta2 <- runif((n/2), min=(pi)/2, max=1.9*(pi))
    z2 <- runif((n/2), min=0, max=5)
    x2 <- cos((pi)-theta2)
    y2 <- sin((pi)-theta2)
    side.lower <- matrix(0, ncol=4, nrow=(n/2))
    for(i in 1:(n/2)){
        side.lower[i,1] <- x2[i]
        side.lower[i,2] <- y2[i]
        side.lower[i,3] <- z2[i]
        side.lower[i,4] <- theta2[i]
    }
    order.theta2 <- order(theta2)
    sort.theta2 <- sort(order.theta2, method="qu", index=TRUE)
    lower.color <- sort.theta2$ix

    ###################
    #     S curve     #
    ###################
    xdata <- rbind(side.upper[,c(1,3,2)], side.lower[,c(1,3,2)])
    xdata <- scale(xdata) + matrix(rnorm(n*3,0, 0.05), n, 3)

    angle <- c(side.upper[,4], side.lower[,4])
    S.color <- c(upper.color, (n/2 + lower.color))
    my.color <- (rainbow130((1.2)*n)[S.color])[1:n]


	scatterplot3d(xdata, color=my.color, 
	              xlab="x", ylab="y", zlab="z", 	              
	              pch=20, angle=30)

	open3d()
	plot3d(xdata[,1], xdata[,2], xdata[,3], 
	       col=my.color, size=5, 
	       xlab="x", ylab="y", zlab="z")

    return(list(xdata=xdata, angle=angle, color=my.color))

}

###########################################################
## trig function on circle                                #
###########################################################
trigcircle <- function(n){
	
	angle <- seq(-pi, pi, length=n)
	n <- length(angle)
	xdata <- matrix(0, ncol=3, nrow=n)
	xdata[,1] <- cos(angle)
	xdata[,2] <- sin(angle)	
	xdata[,3] <- cos(3*angle) + rnorm(length(angle),0, 0.2)	
	my.color <- rainbow130(n)
	scatterplot3d(xdata[,1], xdata[,2], xdata[,3], pch=20, 
	              main="Cosine function supported on circle",
	              color=my.color)

	open3d()
	plot3d(xdata[,1], xdata[,2], xdata[,3],
	       col=my.color, size=5, 
	       xlab="x", ylab="y", zlab="z")

    return(list(xdata=xdata, angle=angle, color=my.color))

}


###########################################################
## noisy spiral                                          #
###########################################################
# from library(diffusionMap)
spiral <- function(n){
	
	xdata <- matrix(0, ncol=2, nrow=n)
	angle <- runif(n)^0.7*10
	beta1 <- 0.15
	beta0 <- 0.5
	xdata[,1] <- beta0*exp(beta1*angle)*cos(angle)+rnorm(n, 0, 0.1)
	xdata[,2] <- beta0*exp(beta1*angle)*sin(angle)+rnorm(n, 0, 0.1)

	order.angle <- order(angle)
	sort.angle <- sort(order.angle, index.return=TRUE)

	col.id <- rainbow130(n)
	my.color <- col.id[sort.angle$ix]
	
	plot(xdata[,1], xdata[,2], pch=20, col=my.color, cex=1, main="Noisy spiral")
	return(list(xdata=xdata, angle=angle, color=my.color))

}

###########################################################
## annulus data set                                       #
###########################################################
annulusData <- function(){
	library(diffusionMap)
	data(annulus)
	n <- nrow(annulus)
	my.color <- rainbow130(n)
	plot(annulus, main="Annulus Data", pch=20, cex=1, col=my.color)
	return(list(xdata=annulus, color=my.color))
}

###########################################################
## Chainlink data set                                      #
###########################################################
ChainlinkData <- function(){    
    
    library(diffusionMap)
    data(Chainlink)
    n <- nrow(Chainlink)
    h <- n/2
    angle1 <- atan2(Chainlink[1:500,1], Chainlink[1:500,2])
    angle2 <- atan2(Chainlink[501:1000,2], Chainlink[501:1000,3])
    angle <- c(angle1, angle2)
	order.angle1 <- order(angle1)
	sort.angle1 <- sort(order.angle1, index.return=TRUE)	
	col.id <- rainbow130(h)
	col.id1 <- col.id[sort.angle1$ix]
	
	order.angle2 <- order(angle2)
	sort.angle2 <- sort(order.angle2, index.return=TRUE)		
    col.id <- grey((1:h)/h) 
	col.id2 <- col.id[sort.angle2$ix]
	
	#plot(Chainlink[1:500,1], Chainlink[1:500,2], col=col.id1)
	#plot(Chainlink[501:1000,2], Chainlink[501:1000,3], col=col.id2)

    my.color <- c(col.id1, col.id2)
    scatterplot3d(Chainlink[,1], Chainlink[,2], Chainlink[,3], 
              color= my.color, main="Chainlink Data")
	return(list(xdata=Chainlink, angle=angle, color=my.color))
}


###########################################################
# 1. Li model 1
###########################################################
Li.model.1 <- function(n, dimension, i){
    set.seed(1234*i)
    n <- n
    dimension <- dimension
    x <- matrix(0, ncol=dimension, nrow=n)
    for(i in 1:dimension){
        x[,i] <- rnorm(n, mean=0, sd=1)
    }
    epsilon <- rnorm(n, mean=0, sd=1)
    e <- matrix(epsilon, ncol=1, nrow=n)
    y.inf <- x[,1] + x[,2] + x[,3] + x[,4] + 0 * x[,5] + e
    y <- matrix(y.inf, ncol=1, nrow=n)
    order.y <- order(y)
    sort.y <- sort(order.y, method="qu", index=TRUE)
    my.color <- (rainbow((1.2)*n)[sort.y$ix])[1:n]

    return(list(Y=y, X=x, color=my.color))
}

###########################################################
# 2. Li model 2
###########################################################
Li.model.2 <- function(n, dimension, sigma, i){
    set.seed(1234*i)
    n <- n
    dimension <- dimension
    sigma <- sigma
    x <- matrix(0, ncol=dimension, nrow=n)
    for(i in 1:dimension){
        x[,i] <- rnorm(n, mean=0, sd=1)
    }
    epsilon <- rnorm(n, mean=0, sd=1)
    e <- matrix(epsilon, ncol=1, nrow=n)
    y.inf <- x[,1] * (x[,1] + x[,2] + 1) + sigma * epsilon
    y <- matrix(y.inf, ncol=1, nrow=n)
    order.y <- order(y)
    sort.y <- sort(order.y, method="qu", index=TRUE)
    my.color <- (rainbow((1.2)*n)[sort.y$ix])[1:n]

    return(list(Y=y, X=x, color=my.color))
}

###########################################################
# 3. Li model 3
###########################################################
Li.model.3 <- function(n, dimension, sigma, i){
    set.seed(1234*i)
    n <- n
    dimension <- dimension
    sigma <- sigma
    x <- matrix(0, ncol=dimension, nrow=n)
    for(i in 1:dimension){
        x[,i] <- rnorm(n, mean=0, sd=1)
    }
    epsilon <- rnorm(n, mean=0, sd=1)
    e <- matrix(epsilon, ncol=1, nrow=n)
    y.inf <- x[,1] / (0.5 + (x[,2] + 1.5)^2 ) + sigma * epsilon
    y <- matrix(y.inf, ncol=1, nrow=n)
    order.y <- order(y)
    sort.y <- sort(order.y, method="qu", index=TRUE)
    my.color <- (rainbow((1.2)*n)[sort.y$ix])[1:n]

    return(list(Y=y, X=x, color=my.color))
}

###########################################################
# 6. Noisy data (only swissroll, S-curve)
###########################################################
superpan.noisydata <- function(datadim3, noisy){
    # datadim3 : n x 3 netrix
    # noisy : noise numbers
    n <- nrow(datadim3)
    NoisyMatrix <- matrix(0, nrow=n, ncol=(3+noisy))
    NoisyMatrix[, 1:3] <- datadim3
    for(i in 1:noisy){
        NoisyMatrix[, 3+i] <- as.matrix(rnorm(n))
    }
    NoisyMatrix
}

###########################################################
## GaussianData                                       #
###########################################################

GaussianData <- function(){
	mu.x <- 0
	mu.y <- 0
	sigma.x <- 1 
	sigma.y <- 1
	rho <- 0

	Q <- function(x, y){
	    s.x <- (x-mu.x)/sigma.x 
	    s.y <- (y-mu.y)/sigma.y
	    return((1/(1-rho^2))*(s.x^2 - 2 * rho * s.x * s.y + s.y^2))
	}
	f <- function(x, y){
	    a <- 1/(2 * pi * sigma.x * sigma.y * sqrt(1-rho^2))
	    return(a * exp(-0.5* Q(x, y)))
	}

	x.grid <- seq(-2, 2, length=25)
	y.grid <- seq(-2, 2, length=25)
	z.grid <- outer(x.grid, y.grid, FUN = f)

	fx <- function(x, y){   
	    return(x)
	}
	fy <- function(x, y){   
	    return(y)
	}
	xy1 <- outer(x.grid, y.grid, FUN = fx)
	xy2 <- outer(x.grid, y.grid, FUN = fy)
	xy <- cbind(c(xy1), c(xy2))

	xtemp <- cbind(xy, c(z.grid))
	rc <- rainbow130(nrow(xtemp))
	id <- order(xtemp[,3])
	x.color <- rc
	x <- xtemp[id, ]
	#library(scatterplot3d)
	#scatterplot3d(x[,1], x[,2], x[,3], main= "gaussian", color=rainbow130(nrow(x)))
   
    return(list(xdata=x, color=x.color))

}


###########################################################
## FishbowData                                            #
###########################################################

FishbowlData <- function(n){
 
	theta <- runif(n)*4/5*pi + 1/5*pi;
	fai <- runif(n)*2*pi;
	x0 <- sin(theta)*cos(fai) 
	y0 <- sin(theta)*sin(fai) 
	z0 <- cos(theta)
	xtemp <- cbind(x0,y0,z0)
	rc <- rainbow130(nrow(xtemp))
	id <- order(xtemp[,3])
	x.color <- rc
	x <- xtemp[id, ]
	#library(scatterplot3d)
	#scatterplot3d(x[,1], x[,2], x[,3], main= "FishbowData", color=rainbow130(nrow(x)))
	return(list(xdata=x, color=x.color))
}


###########################################################
## FishbowData                                            #
###########################################################
IncompleteTireData <- function(n){


	t <- pi*5*runif(n)/3
	s <- pi*5*runif(n)/3
	x0 <- (3+cos(s))*cos(t)
	y0 <- (3+cos(s))*sin(t)
	z0 <- sin(s)

	xtemp <- cbind(x0,y0,z0)
	rc <- rainbow130(nrow(xtemp))
	id <- order(xtemp[,2])
	x.color <- rc
	x <- xtemp[id, ]
	
	#scatterplot3d(x[,1], x[,2], x[,3], main= "FishbowData", color=rainbow130(nrow(x)))
    open3d()
    plot3d(x[,1], x[,2], x[,3], col=x.color, size=3, xlab="", ylab="", zlab="", axes = FALSE)	         
	return(list(xdata=x, color=x.color))

}



###########################################################
## color:rainbow130                                       #
###########################################################
rainbow130 <- function(n){
    rainbow130.rgb <- rgb(read.table("data-other\\color-rainbow130.rb.txt" , header=FALSE))
    a <- n%/%130
	b <- n-130*a
	id <- c(rep(c(1:b), each=(a+1)), rep(c((b+1):130), each=a))
	my.color <- rainbow130.rgb[id]
    return(color=my.color)
}


###########################################################
## use rpanel to plot 3D                                  #
###########################################################
rpanel3dplot <- function(xdata, x.color, pair){

	np <- nrow(pair)
	my.draw <- function(panel) {

        s3d <- scatterplot3d(xdata[,1], xdata[,2], xdata[,3], 
                             color=x.color, pch=20, angle= panel$angle, main = "")
	    for(i in 1:np){           
	        s3d$points3d(xdata[pair[i,],1], xdata[pair[i,],2], xdata[pair[i,],3],
	                     type="l", col="black")   
        }	    	   
	    panel
	}
	my.panel <- rp.control("Angle", angle = 10)
	rp.slider(panel = my.panel, var = angle, from = 10,  to = 360, action = my.draw)

}


###########################################################
## use rpanel to plot 3D                                  #
###########################################################
rpanel3dplot2 <- function(xdata, x.color){

	my.draw <- function(panel) {

        s3d <- scatterplot3d(xdata[,1], xdata[,2], xdata[,3], 
                             color=x.color, pch=20, angle= panel$angle, main = "")  	   
	    panel
	}
	my.panel <- rp.control("Angle", angle = 10)
	rp.slider(panel = my.panel, var = angle, from = 10,  to = 360, action = my.draw)

}