###################################################
require("rrp")

set.seed(123)
nt <- 200
nc <- 200
n <- nt+nc

theta <- runif(nt)*2*pi
r <- runif(nt)*.8
x1 <- cos(theta)*r*1.3
y1 <- sin(theta)*r*.6

theta <- runif(nc)*2*pi
r <- runif(nc)*.8
x2 <- 0.5+cos(theta)*r*.5
y2 <- sin(theta)*r*1.4

dataA <- data.frame(X1 = c(x1, x2), X2 = c(y1, y2)) 


###################################################
tsubjects <- 1:nt
csubjects <- (nt+1):n
treated <- c(rep(TRUE,nt), rep(FALSE,nc))

maxY <- max(dataA[tsubjects,2])
minY <- min(dataA[tsubjects,2])
which( (dataA[csubjects,2] <= maxY) & (dataA[csubjects,2] >= minY) ) -> idxc
minX <- min(dataA[csubjects[idxc],1])
maxX <- max(dataA[csubjects[idxc],1])

which((dataA[tsubjects,1] <= maxX) & (dataA[tsubjects,1] >= minX))-> idxt
cat(paste("number of treated units in the common support:",length(idxt),"\n"))


###################################################
set.seed(123)
theta <- runif(nt)*2*pi
r <- runif(nt)*.8
x1 <- cos(theta)*r*1.3
y1 <- sin(theta)*r*.6

theta <- runif(nc)*2*pi
r <- runif(nc)*.8
x2 <- 1.5+cos(theta)*r*.5
y2 <- sin(theta)*r*1.4
dataB <- data.frame(X1 = c(x1, x2), X2 = c(y1, y2))


###################################################
MyCut <- 15
dati1 <- dataA
dati2 <- dataB

for(i in 1:2)
 dati1[,i] <- ordered(cut(dataA[, i], seq(min(dataA[, i], 
                    na.rm = TRUE), max(dataA[, i], na.rm = TRUE), 
                    length = MyCut), include.lowest = TRUE))

for(i in 1:2)
 dati2[,i] <- ordered(cut(dataB[, i], seq(min(dataB[, i], 
                    na.rm = TRUE), max(dataB[, i], na.rm = TRUE), 
                    length = MyCut), include.lowest = TRUE))

lambda <- seq(0,1,length=50)
nl <- length(lambda)

set.seed(123)
D1 <- rrp.dist(dati1, msplit=5, cut.in=0, asdist=TRUE)
P1 <- 1-as.matrix(D1)
px1 <- P1[tsubjects,csubjects]
S1 <- numeric(nl)
for(l in 1:nl)
	S1[l] <- sum(apply(px1,1,function(x) (length(which(x>=lambda[l]))>0)))/nt


set.seed(123)
D2 <- rrp.dist(dati2, msplit=5, cut.in=0, asdist=TRUE)
P2 <- 1-as.matrix(D2)
px2 <- P2[tsubjects,csubjects]
S2 <- numeric(nl)
for(l in 1:nl)
	S2[l] <- sum(apply(px2,1,function(x) (length(which(x>=lambda[l]))>0)))/nt 


###################################################
par(mar=c(5,4,1,1))
par(mfrow=c(2,2))
plot(dataA,col=treated+2,pch=ifelse(treated,18,17),xlim=c(-1.5,2.0),ylim=c(-1.5,1.5))
rect(minX, minY, maxX, maxY, lty=3)
plot(lambda,S1,type="l", ylim=c(0,1), xlab=expression(lambda), ylab=expression(S(lambda)))
plot(dataB,col=treated+2,pch=ifelse(treated,18,17),xlim=c(-1.5,2.0),ylim=c(-1.5,1.5))
plot(lambda,S2,type="l", ylim=c(0,1), xlab=expression(lambda), ylab=expression(S(lambda)))


###################################################
S <- function(lambda, P, group){
 g2 <- which(!group)
 g1 <- which(group)
 
 f <- function(x){
  if(lambda==0) return(TRUE)
  nm <- as.numeric(names(x))
  idx <- match(g2, nm)
  idx <- as.numeric(na.omit(idx))
  length(which(x[idx]>=lambda))>0
 }
 sum(as.numeric(lapply(P[g1], f))) / length(g1)
}

D3 <- rank.dist(dati1, asdist=FALSE)
gn <- function(x) S(x, D3, treated)
S3 <- sapply(lambda, gn)

D4 <- rank.dist(dati2, asdist=FALSE)
gn <- function(x) S(x, D4, treated)
S4 <- sapply(lambda, gn)


###################################################
par(mar=c(5,4,1,1))
par(mfrow=c(2,2))
plot(lambda,S1,type="l",ylim=c(0,1),xlab=expression(lambda), ylab=expression(S(lambda)))
plot(lambda,S3,type="l",ylim=c(0,1),xlab=expression(lambda), ylab=expression(S(lambda)))
plot(lambda,S2,type="l",ylim=c(0,1),xlab=expression(lambda), ylab=expression(S(lambda)))
plot(lambda,S4,type="l",ylim=c(0,1),xlab=expression(lambda), ylab=expression(S(lambda)))


###################################################
require("rrp")
data("DWvsPSID")
n <- dim(DWvsPSID)[1]
ctr <- which(DWvsPSID$treated==0) 
trt <- which(DWvsPSID$treated==1) 
group <- (DWvsPSID$treated==1)
DWvsPSID$u74 <- factor(DWvsPSID$re74>0)
DWvsPSID$u75 <- factor(DWvsPSID$re75>0)
DWvsPSID$black <- factor(DWvsPSID$black)
DWvsPSID$married <- factor(DWvsPSID$married)
DWvsPSID$nodegree <- factor(DWvsPSID$nodegree)
DWvsPSID$hispanic <- factor(DWvsPSID$hispanic)

str(DWvsPSID)
DWvsPSID <- DWvsPSID[-c(1,9)]
px <- rank.dist(DWvsPSID,cut.in=20)


###################################################
thr <- 0.9
trt.m <- NULL
ctr.m <- NULL
for(i in trt){
 tmp <- as.numeric(names(which(px[[i]]>thr)))
 tmp <- tmp[tmp %in% ctr]
 if(length(tmp)>0){
  trt.m <- c(trt.m, i)
  ctr.m <- unique(c(ctr.m, tmp))
 }
}
ctr.m <- sort(ctr.m)
trt.m <- sort(trt.m)


###################################################
summary(DWvsPSID[trt,])
summary(DWvsPSID[ctr,])
summary(DWvsPSID[trt.m,])
summary(DWvsPSID[ctr.m,])


###################################################
cat(sprintf("\npre-match treated=%d, controls=%d",length(trt), length(ctr)))
cat(sprintf("\npost-match treated=%d, controls=%d",length(trt.m), length(ctr.m)))


###################################################
a <- newXPtr(10,5)
a
str(a)


###################################################
(XPtrToDist(a))
as.dist(matrix(5,10,10))


###################################################
M <- matrix(0,5,5)
d <- newXPtr(5,0)
x <- list(1:3, 4:5)
 
addXPtr(d,x,c(-1,+1))

M[1:3,1:3] <- M[1:3,1:3] - 1
M[4:5,4:5] <- M[4:5,4:5] + 1

(XPtrToDist(d))
as.dist(M)


###################################################
require("rrp")
data("iris")
str(iris,strict="cut",width=70)
set.seed(123)
sample(1:150, 10) -> test 
test
train <- (1:150)[-test]
rrp.dist(iris[,-5]) -> D 


###################################################
rrp.class(D, iris$Species[train], train, test) -> pred
table(pred, iris$Species[test])


###################################################
require("MASS")
data("birthwt")
attach(birthwt)
race <- factor(race, labels = c("white", "black", "other"))
ptd <- factor(ptl > 0)
ftv <- factor(ftv)
levels(ftv)[-(1:2)] <- "2+"
bwt <- data.frame(bwt, age, lwt, race, smoke = (smoke > 0), ptd, 
                   ht = (ht > 0), ui = (ui > 0), ftv)
detach()
rm(race, ptd, ftv)


###################################################
set.seed(123)
n <- dim(bwt)[1]
test <- sample(1:n, 15)
train <- (1:n)[-test]
D <- rrp.dist(bwt[, -1])


###################################################
true.wht <- bwt$bwt[test]
pred.wht <- rrp.predict(D, bwt$bwt[train], train, test)
mean(pred.wht - true.wht)


###################################################
data("iris")
X <- iris
n <- dim(X)[1]
set.seed(123)

miss <- sample(1:n, 10)
for(i in miss)
  X[i, sample(1:5, 2)] <- NA 
X[miss,] 

x <- rrp.impute(X)
x$new.data[miss,] 
iris[miss,]