angulardensityrank<-function(x, y, k) 
{
	#x <- x-vector
	#y <- y-vextor
	#k <- no of upper order statistics to be used
	n <- min(length(x), length(y))
	rx <- n - rank(x[1:n])+1
	ry <- n - rank(y[1:n])+1
	theta <- atan2(rx, ry)
	rad <- k * sqrt((rx)^(-2) + (ry)^(-2))
	plot(density(theta[rad > 1]), type = "b", 
             xlab = "theta", ylab = "angular measure density")  
	abline(v = pi/4.)
}

angulardensityEstAlpha<-function(x, y, k, alpha1, alpha2)
{
	n <- min(length(x), length(y))
	x <- (x/(rev(sort(x)))[k])^(alpha1)
	y <- (y/(rev(sort(y)))[k])^(alpha2)
	r <- sqrt(x^2. + y^2.)
	theta <- atan2(x, y)
	plot(density(theta[r > 1.]), type = "b", 
		xlab = "theta", ylab = "angular measure density") 
	abline(v = pi/4.)
}

angularDFrank<-function(x, y, k)
{
	#x <- x-vector
	#y <- y-vector
	#k <- no of upper order statistics to be used
	if((n <- length(x)) != length(y)) {
		stop("The lengths of the data vectors do not match")
	}
	rx <- n - rank(x)+1
	ry <- n - rank(y)+1
	theta <- atan2(rx, ry)
	ord <- order(theta)
	l <- rx[ord] < k | ry[ord] < k
	stheta <- theta[ord]
	cl <- cumsum(l)
	sl <- sum(l)
	plot(stheta[l], cl[l]/sl, xlim = c(0., pi/2.), ylim = c(0., 1.), 
		type = "l", xlab = "theta\n\t\t", ylab 
		 = "S(theta)", main = "Spectral Distribution Function", 
		sub = paste(
		"number of upper order statistics used", 
			as.character(k)), font.sub = 3.) 
	abline(h = 0.5, lty = 3.)
	points(stheta[l], rep(0.5, sl), pch = "+")
}

Starica2dPlotnew<-function(x,y,k,PlotIt=TRUE,norm=L1norm) 
#function from Sid - two dimensional implementation;changed by Sid Jan12
#Inputs: x,y; choice of k, choice of whether to plot or not, choice of norm,
#Outputs: (scaled norms,scaled norms*j/k; scaled norms=R/R[k]
{
  if(length(x) != length(y)){
    stop("x and y different lengths\n")
  }
  n <- length(x)
  r <- apply(cbind(x,y),1,norm)
  s.r <- rev(sort(r))
  u <- s.r/s.r[k]
  ratio <- (u*(1:n))/k
  if(PlotIt){
    plot(u[u < 10], ratio[u < 10], xlim = c(0.1, 10), type = "l", 
	xlab = "scaling constant", ylab = "scaling ratio",col="blue")
    abline(h = 1,col="red")
    abline(v = 1,lty=2,lwd=0.5,col="red")
    title(paste("k =",k))
  }
  invisible(list(r=u,ratio=ratio,k=k,norm=norm))
}

ChooseK<-function (x, y, PlotIt = TRUE, MultByK=FALSE, norm = L1norm, Lower, Upper) 
{
    if (length(x) != length(y)) {
        stop("x and y different lengths\n")
    }
    n <- length(x)
rx<-ranktransform(x)
ry<-ranktransform(y)
 nk <- min(c(50, n/2))
    Kseq <- round(exp(seq(log(10), log(n/2), len = nk)))
    AbsDiff <- rep(0, nk)
    for (i in 1:nk) {
        Mult <- MultByK * Kseq[i] + (1 - MultByK) * 1
        est <- Starica2dPlotnew(Mult * rx, Mult * ry, Kseq[i], PlotIt = FALSE, 
            norm = norm)
        which.r <- est$r < Upper & est$r > Lower
        if (sum(est$r < Lower)) {
            ratio <- est$ratio[which.r]
            r <- est$r[which.r]
            diff <- (ratio[1:(sum(which.r) - 1)] - 1) * (r[1:(sum(which.r) - 
                1)] - r[2:sum(which.r)])/(1 + (r[-1] - 1)^2)
        }
        else {
            diff <- 10^6
        }
        AbsDiff[i] <- sum(abs(diff))
    }
    k <- Kseq[AbsDiff == min(AbsDiff)]
    if (PlotIt) {
        Mult <- MultByK * k + (1 - MultByK) * 1
        Starica2dPlotnew(Mult * rx, Mult * ry, k, PlotIt = TRUE, norm = norm)
        title(paste("k =", k))
    }
    k
}

Starica2dplotrank<-function (x, y, k, PlotIt = TRUE, norm = L1norm) 
{
    if (length(x) != length(y)) {
        stop("x and y different lengths\n")
    }
    n <- length(x)
rx<-ranktransform(x)
ry<-ranktransform(y)
R <- apply(cbind(k*rx, k*ry), 1, norm)
u <- rev(sort(R))
    ratio <- (u * (1:n))/(length(R[R > 1]))
    if (PlotIt) {
        plot(u[u < 5], ratio[u < 5], xlim = c(0.1, 5), type = "l", 
            xlab = "scaling constant", ylab = "scaling ratio", 
            col = "blue")
        abline(h = 1, col = "red")
        abline(v = 1, lty = 2, lwd = 0.5, col = "red")
        title(paste("k =", k))
    }
    invisible(list(r = u, ratio = ratio, k = k, norm = norm))
}

ChooseKrank<-function (x, y, PlotIt = TRUE,  norm = L2norm, Lower,Upper)
{
    if (length(x) != length(y)) {
        stop("x and y different lengths\n")
    }
    n <- length(x)
rx<-ranktransform(x)
ry<-ranktransform(y)
R<-apply(cbind(rx,ry),1,norm) #norms of pairs after rank transform
R <- rev(sort(R)) #ordered norms; biggest first
nk <- min(c(500, ceiling(.5*n)))
Kseq <- (1:nk) #round(exp(seq(log(10), log(n/2), len = nk)))
dist<-rep(0,nk) #vector of length nk of zeros
for (i in 1:nk) {
dist[i]<-L2norm(       
(  i*(1:n)*R/length(i*R[i*R>=1])  )*(i*R>Lower & i*R <=Upper)-
(i*R>Lower & i*R <=Upper) 
)
}
 k <- (Kseq[dist == min(dist)])
if (PlotIt) {
 u <- k*R
  ratio <- k*R*(1:n)/length(k*R[k*R>=1])
 plot(u[u <= Upper & u>Lower], ratio[u <= Upper & u>Lower], 
	type = "l", xlab = "scaling constant", ylab = "scaling ratio",col="blue")
    abline(h = 1,col="red")
    abline(v = 1,lty=2,lwd=0.5,col="red")
    title(paste("k =",k))
  }
k
}

ranktransform<-function (x)
{
tr<-(length(x)-rank(x)+1)^(-1)
#invisible(list(tr=tr))
}

langmeas<-
function(x,y,k)
{      sum<-abs(x)+abs(y)
    sortsum<-sort(sum,decreasing=TRUE)
    theta1<- x/sum
    theta2<- y/sum
    theta1<-theta1[sum>sortsum[k]]
    theta2<-theta2[sum>sortsum[k]]
    cat("range negative theta1=", range(theta1[theta1<0]))
    cat("   range theta1=", range(theta1[theta1>0]),
        "\n")
    cat("   summary theta1", summary(theta1[theta1>0]))
    par(mfrow=c(1,2))
    plot(theta1,theta2,pch=20, xlim=c(-1,1),ylim=c(-1,1),main=expression("Diamond plot; L"[1] ),
         xlab=expression(theta[1]),ylab=expression(theta[2]))
#    plot(density(theta1),xlab="theta1",ylab="angular measure density")
    hist(theta1,nclass=10,xlim=c(-1,1),main=expression(paste("Histogram of ", theta[1])),xlab=expression(theta[1])    )
    par(mfrow=c(1,1))
   }

langmeaspos<-function(x,y,k)
{  #use when both components are positive
    sum<-abs(x)+abs(y)
    sortsum<-sort(sum,decreasing=TRUE)
    theta1<- x/sum
    theta2<- y/sum
    theta1<-theta1[sum>sortsum[k]]
    theta2<-theta2[sum>sortsum[k]]
#    cat("range negative theta1=", range(theta1[theta1<0]))
    cat("   range theta1=", range(theta1[theta1>0]))
    par(mfrow=c(1,2))
    plot(theta1,theta2,pch=20, xlim=c(0,1),ylim=c(0,1),main="Diamond plot; L1 norm")
#    plot(density(theta1),xlab="theta1",ylab="angular measure density")
    hist(theta1,nclass=10,xlim=c(0,1))    
    par(mfrow=c(1,1))
}

langmeas2<-function(x,k)
{#use to delete wedges from both pos and neg quadrants
      sum<-abs(x[,1])+abs(x[,2])
    sortsum<-sort(sum,decreasing=TRUE)
    theta<-cbind(x[,1]/sum,x[,2]/sum)
    theta<-subset(theta,sum>sortsum[k])
    cat("range negative theta1=", range(theta[,1][theta[,1]<0]))
    cat("   range pos theta1=", range(theta[,1][theta[,1]>0]),"\n")
    cat("   summary pos theta1", summary(theta[,1][theta[,1]>0]),"\n")
    cat("   10% and 90% quantile pos theta1", quantile(theta[,1][theta[,1]>0],probs=c(.1,.9)),"\n")
    cat("   summary neg theta1", summary(theta[,1][theta[,1]<0]),"\n")
    cat("   10% and 90$ quantile neg theta1", quantile(theta[,1][theta[,1]<0],probs=c(.1,.9)))
    par(mfrow=c(1,2))
    plot(theta,pch=20, xlim=c(-1,1),ylim=c(-1,1),main=expression("Diamond plot; L"[1] ),
         xlab=expression(theta[1]),ylab=expression(theta[2]))
#    plot(density(theta1),xlab="theta1",ylab="angular measure density")
    hist(theta[,1],nclass=10,xlim=c(-1,1),main=expression(paste("Histogram of ", theta[1])),xlab=expression(theta[1])    )
    par(mfrow=c(1,1))
}


hillish<-function(w,str="data")
{
  n=length(w[,1])

  ind=array(0,c(n,1))
  indneg=array(0,c(n,1))

  rx=array(0,c(n,1))
  rxneg=array(0,c(n,1))

  s=array(0,c(n,2))
  g=array(0,c(n,2))
  w1=array(0,c(n,2))

  w1<-w

  g<-cbind(runif(n),runif(n))

  w1<-w+cbind(min(w[,1])*g[,1],min(w[,2])*g[,2])
  s<-w1[order(w1[,2],decreasing=T),]


  #rx1=rank(s[1:n,1])
  #rx1[n+1]=n+1

  for (k in n:2){
    print(k)
    rx[1:k]=rank(s[1:k])
    rxneg[1:k]=rx[k:1]

    ind[k]=(1/k)*sum(log(k/(1:k))*log(k/(rx[1:k])))
    indneg[k]=(1/k)*sum(log(k/(1:k))*log(k/(rxneg[1:k])))

  }
   ind[1]=0
  indneg[1]=0


  plot(2:n,ind[2:n],type="l",col="blue",lwd=0.5,xlab="number of order statistics",ylab="Hillish statistic",main=paste("Hillish(R,theta): ",str,sep=" "))
  plot(2:n,indneg[2:n],type="l",col="blue",lwd=0.5,xlab="number of order statistics",ylab="Hillish statistic",main=paste("Hillish(R,-theta): ",str,sep=" "))
 }