#!/usr/bin/Rscript --silent
rem <- function(...) invisible(T)
rem( '
Rscript.exe "%~F0" "%1" "%2" "%3" "%4" "%5" "%~F6" "%~F7","%~F8","%~F9"
EXIT /B
rem ')
### above=BAT, below=R
pkgrequire <- function(pkg, preferred = 'http://cran.rstudio.com/') {
	if (!require(pkg, character.only = TRUE)) {
		cat("Cannot find", pkg, "package.\n")
		lib_path = Sys.getenv("R_LIBS_USER")
		if (!file.exists(lib_path)) dir.create(lib_path, rec = TRUE)
		install.packages(pkg, dep = TRUE, repos = preferred, lib = lib_path)
		if (!require(pkg, character.only = TRUE, lib.loc = lib_path)) stop("Package not found")
	}
}
pkgrequire("optparse")
pkgrequire("RColorBrewer")

NONE <- 0
CONVERT <- 1
OPTIMIZE <- 2
UNFOLD <- 4
PRINTSCORE <- 8
PLOT <- 16
ALL <- 1024
ACT <- NONE

getangle <- function(p1, p2, p3, p4, degs = TRUE) {
# Matlab reverse args
#	res <- atan2((p2$x - p1$x) * (p4$y - p3$y) - (p2$y - p1$y) * (p4$x - p3$x),
#			(p2$x - p1$x) * (p4$x - p3$x) + (p2$y - p1$y) * (p4$y -p3$y))
	a <- (p2$x - p1$x) * (p4$x - p3$x) + (p2$y - p1$y) * (p4$y - p3$y)
	b <- (p2$x - p1$x) * (p4$y - p3$y) - (p2$y - p1$y) * (p4$x - p3$x)
	res <- atan2(b, a)
	if (degs == TRUE) {
		res <- res * 180/pi
	}
}

option_list <- list(
	make_option(c("-v", "--verbose"), action="store_true", default=FALSE,
							help="Print extra output [default=false]"),
	make_option(c("-m", "--rmin"), type="double", default = 10,
							help = "rmin [default %default]"),
	make_option(c("-x", "--rmax"), type="double", default = 30, # 18.1321,
							help = "rmax [default %default]"),
	make_option(c("-s", "--scale"), type="double", default = 1,
							help = "scale [default %default]"),
	make_option(c("-b", "--breaks"), type="integer", default = 30,
							help="breaks [default %default]", metavar="number"),
	make_option(c("-t", "--smooth"), type="double", default = 2e-8,
							help="smooth [default %default]")
)
opt_parser <- OptionParser(usage = "usage: %prog [options]", option_list = option_list, description = "hmm3s", epilogue = "Send feedback to mackoel@gmail.com")

cmdline <- commandArgs(trailingOnly = TRUE)

flaggs <- parse_args(opt_parser, args = cmdline, print_help_and_exit = TRUE, positional_arguments = TRUE)
opts <- flaggs$options
args <- flaggs$args

USAGE_STRING = "============= Polarplot Script ========================
Usage:

$ polarplot <input>

*************************************************************************\n
"

if(length(args) < 1){
	cat(USAGE_STRING)
	stop("Please supply valid arguments\n")
}


flattabfile <- args[[1]]
outfiles <- unlist(strsplit(args[2], ",", fixed = TRUE))
flatresfile <-  outfiles[1]
rdatafile <- outfiles[2]
patternfile <- outfiles[3]

p <- read.csv(file = flattabfile, stringsAsFactors = F)
p <- p[p$nucmask_notempty > 0,]
N <- dim(p)[1]
xcenter <- mean(p$xmean)
ycenter <- mean(p$ymean)
xhat <- p$xmean - xcenter
yhat <- ycenter - p$ymean

phi <- xhat
rho <- sqrt(xhat^2 + yhat^2)
rho <- 100 * rho/max(rho)
for (t in 1:N) {
	p1 <- list(x=0, y=0)
	p2 <- list(x=1, y=0)
	p3 <- list(x=0, y=0)
	p4 <- list(x=xhat[t], y=yhat[t])
	phi[t] <- getangle(p1, p2, p3, p4) + 180
}

phihat <- abs(abs(phi - 180) - 90)
phihat <- phihat/90
rhohat <- opts$scale * (1 - phihat) * rho
maxrhohat <- max(rhohat)
minrhohat <- min(rhohat)
res <- cbind(xhat, yhat, rho, phi, rhohat, phihat, p)
#r <- res[opts$rmin * minrhohat < res$rhohat & res$rhohat < opts$rmax * maxrhohat,]
r <- res[opts$rmin < res$rho & res$rho < opts$rmax,]
r <- r[order(r$phi),]


conc <- ifelse(r$nucmask_notempty > 0, r$nucen_muc/r$nucmask_notempty, 0)
#runmed(x, k, endrule = c("median", "keep", "constant"),
#       algorithm = NULL, print.level = 0)
#sconc <- runmed(conc, opts$smooth)
#smooth.spline(x, y = NULL, w = NULL, df, spar = NULL, cv = FALSE,
#              all.knots = FALSE, nknots = .nknots.smspl,
#              keep.data = TRUE, df.offset = 0, penalty = 1,
#              control.spar = list(), tol = 1e-6 * IQR(x))
spconc <- smooth.spline(x = r$phi, y = conc, control.spar = list(eps = opts$smooth))
#predict(object, x, deriv = 0, ...)
pconc <- predict(spconc, r$phi)
write.csv(file=patternfile, data.frame(angle=r$phi, conc=pconc$y), row.names=F)

palette2 <- colorRampPalette(c("yellow", "red"))
palette <- palette2(opts$breaks)[as.numeric(cut(conc,breaks = opts$breaks))]

pdf(flatresfile)
#
# 1D plot of conc
#
plot(r$phi, conc, type='p', main = paste("Rmin =", opts$rmin, "%", ";", "Rmax =", opts$rmax, "%"))
lines(r$phi, pconc$y, type='l', col = "blue")
#
# 2D plot with threshold
#
plot(r$xmean, -r$ymean, pch=16, col = ifelse(conc > 90, "blue", "yellow"))
#
# 2D plot with posterized color
#
plot(r$xhat, r$yhat, pch=16, col = palette)
#
# 2D plot of phi
#
plot(r$xhat, r$yhat,
	pch=16,
	col = ifelse(r$phi > 90, ifelse(r$phi > 180, ifelse(r$phi > 270, "red", "yellow") , "black"), "blue") 
)
s <- dev.off()

save(res, spconc, file = rdatafile)

