#### Example of using Picard et al.'s algorithm as implemented
#### in the tilingArray package

## simulate artificial chromosome with an aberration in the middle
LR <- rnorm(100, sd = 0.25)
LR[46:55] <- LR[46:55] + 0.75

## maxcp is the maximum number of segments to find in the profile
## maxk is the maximum size allowed for each segment
library(tilingArray)
obj <- findSegments(LR, maxcp = 5, maxk = 100)

selectionR <- function(J, Km, S = 0.5)
{
  Jtild <- (J[Km] - J) / (J[Km] - J[1]) * (Km - 1) + 1;
  temp <- which(diff(diff(Jtild)) > S);
  if( length(temp) > 0 )  {
    Kh <- max(temp) + 1;
  }
  else  {
    Kh <- 1;
  }  
  return(Kh)
}

Kh <- selectionR(obj$J, Km = 5, S = 0.5)
temp <- obj$th[Kh, ]

## isolate start and end coordinates for each segment
bp.start <- c(1, temp[1:(Kh-1)])
bp.end <- temp[1:Kh] - 1

## calculate means for each segment
segment.means <- numeric(Kh)
for( i in 1:Kh ) {
  segment.means[i] <- mean(LR[bp.start[i]:bp.end[i]])
}

## plot results for illustration purposes
##  * green horizontal lines denote segment means
##  * red vertical lines denote start coordinates
##  * blue vertical lines denote end coordinates
plot(LR, pch = 20, xlab = "Index", ylab = "Simulated Log2Ratio")
segments(x0 = bp.start, x1 = bp.end,
         y0 = segment.means, y1 = segment.means, col = "green")
abline(v = bp.start, col = "red")
abline(v = bp.end, col = "blue")