args = commandArgs(trailingOnly=TRUE)
datPath=args[1]
cat(paste0("Simulating data and storing them in ",datPath,"/..."))
#setwd(datPath)
set.seed(1111)
M      <-  6000
minMAF <-  0.01
a      <-   1.1
p      <- minMAF + (1-2*minMAF)*rbeta(M,a,a)
dp     <- 2 * p
maf    <- ifelse(p>0.5,1-p,p)
h      <- 2*p*(1-p)
s      <- sqrt(h)
r2Tag  <- 1.0 * (4*p*(1-p))

simTag <- function(x,p,r2){
  D  <- sqrt(r2)*p*(1-p)
  l0 <- which(x==0); n0 <- length(l0)
  l1 <- which(x==1); n1 <- length(l1)
  l2 <- which(x==2); n2 <- length(l2)
  y  <- rep(NA,length(x))
  if(n0>0){
    y[l0] <- rbinom(n0,2,p-D/(1-p))
  }
  if(n1>0){
    y[l1] <- rbinom(n1,1,p-D/(1-p)) + rbinom(n1,1,p+D/p)
  }
  if(n2>0){
    y[l2] <- rbinom(n2,2,p+D/p)
  }
  ## print(cor(x,y)^2)
  return(y)
}

makeID <- function(i,base="IID"){
  if(i<10)           return(paste0(base,"000",i))
  if(i>=10  & i<100)   return(paste0(base,"00",i))
  if(i>=100 & i<1000)   return(paste0(base,"0",i))
  if(i>=1000) return(paste0(base,i))
}
makeSNP <- function(i){
  if(i<10)           return(paste0("rs000",i))
  if(i>=10   & i<100)   return(paste0("rs00",i))
  if(i>=100  & i<1000)   return(paste0("rs0",i))
  if(i>=1000 & i<10000) return(paste0("rs0",i))
  if(i>=10000) return(paste0("rs",i))
}

## Simulate founders
nFounders <- 5000
xFounders <- do.call("cbind",lapply(1:M,function(j) rbinom(nFounders,2,prob=p[j])))

## Simulate siblings
nsibpairs <- 500
nParents  <- 2*nsibpairs
iParents  <- sample(1:nFounders,nParents)
iParent1  <- iParents[+(1:nsibpairs)]
iParent2  <- iParents[-(1:nsibpairs)]
xParent1  <- xFounders[iParent1,]
xParent2  <- xFounders[iParent2,]

xp1_Sib1  <- do.call("rbind",lapply(1:nsibpairs,function(i) rbinom(M,1,prob=0.5*xParent1[i,])))
xp2_Sib1  <- do.call("rbind",lapply(1:nsibpairs,function(i) rbinom(M,1,prob=0.5*xParent2[i,])))

xp1_Sib2  <- do.call("rbind",lapply(1:nsibpairs,function(i) rbinom(M,1,prob=0.5*xParent1[i,])))
xp2_Sib2  <- do.call("rbind",lapply(1:nsibpairs,function(i) rbinom(M,1,prob=0.5*xParent2[i,])))
xSib1     <- xp1_Sib1 + xp2_Sib1
xSib2     <- xp1_Sib2 + xp2_Sib2

## Same for Tag

## ped/map file
xFam    <- do.call("rbind",lapply(1:nsibpairs,function(i) rbind(xParent1[i,],xParent2[i,],xSib1[i,],xSib2[i,])))
xUnrel  <- xFounders[-iParents,]
nTot    <- nrow(xFam) + nrow(xUnrel)
xTot    <- rbind(xFam,xUnrel)
xTags   <- do.call("cbind",lapply(1:M,function(j) simTag(xTot[,j],p[j],r2Tag[j])))
xGeno   <- do.call("cbind",lapply(1:M,function(j) cbind(xTot[,j],xTags[,j])))

nchr    <- 10
chr     <- rep(1:nchr,each=(2*M)/nchr)
pos     <- do.call("c",lapply(1:nchr,function(k) sort(sample(1e4:99999,(2*M)/nchr))))
a1a2    <- do.call("rbind",lapply(1:(2*M),function(j) sample(c("A","C","G","T"),2)))
snp     <- sapply(1:(2*M),makeSNP)

refGeno <- t(sapply(1:(2*M),function(j) c(paste0(a1a2[j,1],"\t",a1a2[j,1]),paste0(a1a2[j,1],"\t",a1a2[j,2]),paste0(a1a2[j,2],"\t",a1a2[j,2]))))
ped     <- do.call("cbind",lapply(1:(2*M),function(j) refGeno[j,1+xGeno[,j]]))

## fam file
iid    <- sapply(1:nTot,makeID,"IID")
fid    <- c(rep(sapply(1:nsibpairs,makeID,"FAM"),each=4),sapply(1:(nFounders-nParents),makeID,"UNREL"))
ipid   <- iid[seq(1,4*nsibpairs,by=4)]
imid   <- iid[seq(2,4*nsibpairs,by=4)]
pid    <- c(do.call("c",lapply(ipid,function(i) c(0,0,i,i))),rep(0,nFounders-nParents))
mid    <- c(do.call("c",lapply(imid,function(i) c(0,0,i,i))),rep(0,nFounders-nParents))
sex    <- c(do.call("c",lapply(1:nsibpairs,function(i) c(1,2,1,2))),sample(1:2,nFounders-nParents,replace=TRUE))
pheno  <- rep(-9,nTot)
fam    <- cbind.data.frame(fid,iid,pid,mid,sex,pheno)

## ped/geno
mapData <- cbind.data.frame(chr,snp,0,pos)
pedData <- cbind.data.frame(fam,ped)

write.table(mapData,paste0(datPath,"/mydata.map"),quote=FALSE,row.names=FALSE,col.names=FALSE,sep="\t")
write.table(pedData,paste0(datPath,"/mydata.ped"),quote=FALSE,row.names=FALSE,col.names=FALSE,sep="\t")

## Simulate phenotype
z  <- do.call("cbind",lapply(1:M,function(j) (xTot[,j]-dp[j])/s[j]))

## trait 1 - shared environt
vA <- 0.25
vC <- 0.50
vE <- 0.25
b1 <- rnorm(M,mean=0,sqrt(vA/M))
ce <- c(rep(rnorm(nsibpairs,mean=0,sqrt(vC)),each=4),rnorm(nFounders-nParents,mean=0,sqrt(vC)))
e1 <- rnorm(nTot,mean=0,sqrt(vE))
g1 <- c(z%*%b1)
y1 <- g1+ce+e1

## trait 2
vA1 <- 0.60
vA2 <- 0.00
vE  <- 1-vA1-vA2
l   <- which(maf<0.05)
M1  <- length(l)
M2  <- M-M1

b1 <- rnorm(M1,mean=0,sd=sqrt(vA1/M1))
b2 <- rnorm(M2,mean=0,sd=sqrt(vA2/M2))

g1 <- c(z[,+l]%*%b1)
g2 <- c(z[,-l]%*%b2)
e  <- rnorm(nTot,mean=0,sd=sqrt(vE))
y2 <- g1 + g2 + e

write.table(cbind(fam[,1:2],y1,y2),paste0(datPath,"/mydata.phen"),quote=FALSE,row.names=FALSE,col.names=FALSE,sep="\t")
cat("[done].\n")