R/picard.R

In MarieEtienne/segTraj: Segmentation of Bivariate Time Series

picard <- function(x, ...)
{
    UseMethod("picard")
}


picard.default <- function(x,P,Kmax, pos=NULL, ...)
{
    if (ncol(x) !=2)
        stop("not yet implemented for more or less than two variables")
    if (!is.null(pos)) {
        if (length(pos)!=nrow(x))
            stop("pos should have the same length as x")
    }
    x <- as.matrix(x)
    xor <- x
    cons <- apply(x, 1, function(y) all(!is.na(y)))
    nna.places <- c(1:nrow(x))[cons]
    x <- x[nna.places,]

    Linc <- matrix(-Inf,nrow=Kmax,ncol=1)
    n <- dim(x)[2]
    param <- list()
    Kmin <- P

    if (P==1){
        rupt <- c(1, n)
        phi <- list(mu=matrix(rowMeans(x),ncol=1),
                    sigma <- matrix(apply(x,1,sd),ncol=1),
                    prop=1)
        Linc[Kmin:Kmax] <- logdens_simultanee(x,phi)
        param <- list(phi=phi,rupt=rupt, tau=NA)
        resu <- list(Linc=Linc, param=param)
    } else {
        for (i in 1:ncol(x))
            x[,i] <- as.double(x[,i])

        resu <- .Call("hybridSimultanee", x, as.integer(P), as.integer(Kmax), PACKAGE="segTraj")
        names(resu) <- c("Linc","param")
        nonnul <- c(1:length(resu$param))[!sapply(resu$param, is.null)]
        for (i in nonnul) {
            names(resu$param[[i]]) <- c("phi","rupt","tau")
            resu$param[[i]]$rupt <- resu$param[[i]]$rupt+1L
            names(resu$param[[i]]$phi) <- c("mu","sigma","prop")
        }
    }

    attr(resu, "P") <- P
    attr(resu, "Kmax") <- Kmax
    attr(resu, "seq.ori") <- xor
    attr(resu, "pos") <- pos
    attr(resu, "nna.places") <- nna.places
    class(resu) <- "picard"
    attr(resu, "typeseg") <- "default"
    return(resu)
}


print.picard <- function(x, ...)
{
    cat("=========\n\nSegmentation of a series using the method of Picard\n",
        "based on P =", attr(x, "P"),"different models",
        "with Kmax = ", attr(x, "Kmax"), ".\n\n=========\n",
        "This object is a list containing\n",
        "$Linc: the likelihood for the best segmentation with K-segment\n",
        "$param: a list containing the parameters of the segmentation\n",
        "        with one element per K; each element is itself a list containing:\n",
        "..$phi:  the parameters of the models. A list containing:\n",
        "......$mu:    the mean of each variable (column) for each model (row)\n",
        "......$sigma: the standard deviation of each variable (column) for each model (row)\n",
        "......$prop:  the proportion of observations generated by each model\n",
        "..$rupt: the indices of the limits (start, end) of the K segments (one per row)\n",
        "..$tau:  the likelihood of each model (row) for each segment (column)\n\n")
}


picard.ltraj <- function(x, P, Kmax, which=c("x","y"), ...)
{
    if (!inherits(x, "ltraj"))
        stop("ltraj should be of class \"ltraj\"")
    if (length(x)>1)
        stop("only implemented for one trajectory")
    if (length(which)!=2)
        stop("for the moment, only implemented for two variables")

    y <- x[[1]]
    if (!is.null(attr(y, "infolocs")))
        y <- cbind(y, attr(y, "infolocs"))

    if (!all(which%in%names(y)))
        stop("The variables in which were not found")

    coin <- y[,which]

    cons <- apply(coin, 1, function(x) all(!is.na(x)))
    nna.places <- c(1:nrow(coin))[cons]
    coin <- coin[nna.places,]
    pic <- picard(x=coin, P=P, Kmax=Kmax)
    attr(pic, "nna.places") <- nna.places
    attr(pic, "seq.ori") <- x
    attr(pic, "which") <- which
    attr(pic, "typeseg") <- "ltraj"
    return(pic)
}


plot.picard <- function(x, axes = FALSE, number=NULL,
                        addparam=FALSE, bg=TRUE,...)
{
    if (!inherits(x, "picard"))
        stop("x should be of class picard")

    typex <- attr(x, "typeseg")
    nna <- attr(x, "nna.places")

    ## prepare the two series and the x-coordinate
    if (typex == "ltraj") {

        xk <- attr(x, "seq.ori")
        y <- xk[[1]]
        pos <- y$date
        if (!is.null(attr(y, "infolocs")))
            y <- cbind(y, attr(y, "infolocs"))
        xk <- y[nna,attr(x,"which")]
        serie <- pos[nna]
    } else {
        xk <- attr(x, "seq.ori")
        pos <- attr(x, "pos")
        if (!is.null(pos)) {
            serie <- pos[nna]
        } else {
            serie <- c(1:nrow(xk))[nna]
        }
        xk <- xk[nna,]
    }

    ## ylim
    moxk <- apply(xk, 2, mean)
    sdxk <- apply(xk, 2, sd)
    xk <- scale(xk)
    ra <- apply(xk,2,range)
    ra <- c(min(ra[1,]),max(ra[2,]))

    ## Number of plots
    nonnul <- c(1:length(x$param))[!sapply(x$param, is.null)]
    if (axes) {
        opar <- par(mfrow = c(length(nonnul),1), mar=c(3,3,2,1))
    } else {
        opar <- par(mfrow = c(length(nonnul),1), mar=c(0,0,2,0))
    }
    on.exit(par(opar))

    ## Number of models
    P <- attr(x, "P")
    colo <- grey((1:P)/(P+1))

    for (i in nonnul) {
        plot(serie, xk[,1], type="n", axes=FALSE, main=paste("K =", i), xlab="", ylab="",
             ylim = ra)
        rupt <- x$param[[i]]$rupt

        ## the most likely model
        modt <- apply(x$param[[i]]$tau,1,which.max)
        if (bg) {
            tmp <- sapply(1:nrow(rupt), function(j) {
                ## the polygon
                pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
                               serie[rupt[j,2]],serie[rupt[j,1]]),
                             c(ra[1], ra[1], ra[2], ra[2]))
                m <- modt[j]
                polygon(pol[,1],pol[,2], col=colo[m])
            })
        }

        ## the two series
        if (!is.null(number)){
            if (addparam) {
                m <- sapply(1:nrow(rupt), function(j) {
                    mo <- (x$param[[i]]$phi$mu[modt[j],number]-moxk[number])/sdxk[number]
                    si <- x$param[[i]]$phi$sigma[modt[j],number]/sdxk[number]
                    segments(serie[rupt[j,1]], mo,
                             serie[rupt[j,2]], mo, col="black", lwd=2)
                    pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
                                   serie[rupt[j,2]],serie[rupt[j,1]]),
                                 c(mo-si, mo-si, mo+si,mo+si))
                    polygon(pol[,1],pol[,2], density=10)
                })
            }
            lines(serie, xk[,number], lwd=3, col="black")
            lines(serie, xk[,number], col="white")

        } else {
            if (addparam) {
                for (numberb in 1:2) {
                    m <- sapply(1:nrow(rupt), function(j) {
                        mo <- (x$param[[i]]$phi$mu[modt[j],numberb]-moxk[numberb])/sdxk[numberb]
                        si <- x$param[[i]]$phi$sigma[modt[j],numberb]/sdxk[numberb]
                        segments(serie[rupt[j,1]], mo,
                                 serie[rupt[j,2]], mo, col="black", lwd=2)
                        pol <- cbind(c(serie[rupt[j,1]],serie[rupt[j,2]],
                                       serie[rupt[j,2]],serie[rupt[j,1]]),
                                     c(mo-si, mo-si, mo+si,mo+si))
                        polygon(pol[,1],pol[,2], density=c(20,10)[numberb],
                                col=c("black","red")[numberb])
                    })
                }
            }
            lines(serie, xk[,1], lwd=3, col="black")
            lines(serie, xk[,1], col="white")
            lines(serie, xk[,2], lwd=3, col="red")
            lines(serie, xk[,2], col="white")
        }

        ## date:
        if (axes) {
            axis(1)
            axis(2)
        }
        box()

    }
    return(invisible(NULL))

}



testConvexHull <- function(y)
  {
    n <- length(y)
    print(y)
    print(n)
     resu <- .Call("testConvexHull", as.numeric(y), as.integer(n), PACKAGE="segTraj")
}