R/mean.ctmm.R

In ctmm: Continuous-Time Movement Modeling

# M %*% S %*% t(M) -> L %*% S[TRI]
quad2lin <- function(M,diag=FALSE)
{
  M <- rbind(M) # [n,m]
  n <- nrow(M)
  m <- ncol(M)

  TRI <- diag(1,m)
  TRI <- upper.tri(TRI,diag=TRUE)
  TRI <- which(TRI)
  p <- length(TRI)

  if(diag)
  { J <- matrix(0,n,p) }
  else
  { J <- matrix(0,p,p) }

  for(i in 1:p)
  {
    # E <- array(0,dim(M))
    E <- array(0,c(m,m))
    E[TRI[i]] <- 1
    E <- E + t(E)
    E <- E/max(E)

    if(diag) # n!=m and only care about n variances
    {
      # O(n^2) slow method like below
      # E <- M %*% E %*% t(M) # [n,n]
      # E <- diag(E)

      # O(n) fast method
      E <- sapply(1:n,function(i){M[i,] %*% E %*% M[i,]})
    }
    else # n==m and want full covariance
    {
      E <- M %*% E %*% t(M) # [n,n]
      E <- E[TRI]
    }
    J[,i] <- E
  }

  return(J)
}

#############
mean.features <- function(x,debias=TRUE,weights=NULL,trace=FALSE,IC="AICc",select="all",...)
{
  if(is.null(weights))
  { weights <- rep(1,length(x)) }

  N <- length(x)
  axes <- x[[1]]$axes

  ISO <- sapply(x,function(y){y$isotropic})
  ISO <- all(ISO)

  # only want biological parameters
  for(i in 1:length(x))
  {
    FEATURES <- x[[i]]$features
    ERRORS <- grepl('error',FEATURES)
    if(any(ERRORS))
    {
      x[[i]]$error <- FALSE
      x[[i]]$features <- FEATURES[!ERRORS]
      #FEATURES <- rownames(x[[i]]$COV)
      #ERRORS <- grepl('error',FEATURES) # how can this differ?
      x[[i]]$COV <- x[[i]]$COV[!ERRORS,!ERRORS,drop=FALSE]
    }
  }

  # all possible features
  FEATURES <- lapply(x,function(y){y$features})
  FEATURES <- unlist(FEATURES)
  FEATURES <- unique(FEATURES)

  TAU.EXPAND <- all(c("tau","tau position") %in% FEATURES)
  if(TAU.EXPAND)
  {
    FEATURES <- FEATURES[FEATURES!="tau"]
    if("tau velocity" %nin% FEATURES)
    { FEATURES <- c(FEATURES,"tau velocity") }
  }

  NFEAT <- FEATURES[FEATURES %nin% c('major','minor','angle')]

  M <- length(FEATURES)
  # all possible RSF beta
  BETA <- lapply(x,function(y){names(y$beta)})
  BETA <- unlist(BETA)
  BETA <- unique(BETA)

  MEANS <- VARS <- array(TRUE,M)
  names(MEANS) <- names(VARS) <- FEATURES

  MU <- array(0,c(N,M))
  SIGMA <- array(0,c(N,M,M))
  INF <- diag(Inf,M)
  colnames(MU) <- FEATURES
  dimnames(SIGMA) <- list(NULL,FEATURES,FEATURES)

  for(i in 1:N)
  {
    x[[i]] <- log_ctmm(x[[i]],debias=debias)
    features <- names(x[[i]]$par)

    if(TAU.EXPAND && "tau" %in% features)
    {
      # expand features
      features[features=="tau"] <- "tau position"
      features <- c(features,"tau velocity")

      # expand point estimates
      NAMES <- names(x[[i]]$par)
      names(x[[i]]$par)[NAMES=="tau"] <- "tau position"
      x[[i]]$par["tau velocity"] <- x[[i]]$par["tau position"]

      # expand covariance
      NAMES <- rownames(x[[i]]$COV)
      NAMES[NAMES=="tau"] <- "tau position"
      dimnames(x[[i]]$COV) <- list(NAMES,NAMES)

      ROW <- x[[i]]$COV['tau position',]
      x[[i]]$COV <- rbind(x[[i]]$COV,'tau velocity'=ROW)
      ROW['tau velocity'] <- ROW['tau position']
      x[[i]]$COV <- cbind(x[[i]]$COV,'tau velocity'=ROW)
    }

    MU[i,features] <- x[[i]]$par
    SIGMA[i,,] <- INF
    SIGMA[i,features,features] <- x[[i]]$COV

    P <- c("major","minor")
    if(!ISO && x[[i]]$isotropic)
    {
      MU[i,'minor'] <- MU[i,'major']
      # perfectly correlated observation
      SIGMA[i,P,P] <- matrix(SIGMA[i,'major','major'],2,2)
      # copy over other correlations as well
      SIGMA[i,'minor',NFEAT] <- SIGMA[i,'major',NFEAT]
      SIGMA[i,NFEAT,'minor'] <- SIGMA[i,NFEAT,'major']
    }
  } # for(i in 1:N)

  ###################
  J <- diag(1,ncol(MU))
  dimnames(J) <- list(FEATURES,FEATURES)

  # everything is already log transformed - diagonalize as much as possible
  if("minor" %in% FEATURES)
  {
    J["major",c("major","minor")] <- c(1/2,1/2) # average variance
    J["minor",c("major","minor")] <- c(1,-1) # difference / eccentricity

    if("tau" %in% FEATURES)
    { J["tau",c("major","minor","tau")] <- c(1/2,1/2,-1) } # D

    if("tau position" %in% FEATURES)
    { J["tau position",c("major","minor","tau position")] <- c(1/2,1/2,-1) } # D
  }
  else
  {
    if("tau" %in% FEATURES)
    { J["tau",c("major","tau")] <- c(1,-1) } # D

    if("tau position" %in% FEATURES)
    { J["tau position",c("major","tau position")] <- c(1,-1) } # D
  }

  # if("tau velocity" %in% FEATURES)
  # { J["tau velocity",c("major","tau position","tau velocity")] <- c(1,-1,-1) } # MSV correlates with D

  # if("omega" %in% FEATURES) # MSV - delta-MSV # sigma cancels out
  # {
  #   if("tau" %in% FEATURES) # OUf
  #   { J["omega",c("tau","omega")] <- c(1,-1) }
  #   else if(all(c("tau position","tau velocity") %in% FEATURES)) # OUF
  #   { J["omega",c("tau position","tau velocity","omega")] <- c(1,1,-2) }
  #   else if("tau velocity" %in% FEATURES) # IOU
  #   { J["omega",c("tau velocity","omega")] <- c(1,-1) }
  #   else if("tau position" %in% FEATURES) # OU
  #   { J["omega",c("tau position","omega")] <- c(1,-1) }
  # }

  # missing data / infinite uncertainties
  INF <- t(apply(SIGMA,1,function(M){diag(M)==Inf})) # [ind,par]

  tJ <- t(J)
  MU[INF] <- 0 # zero out infinite uncertainties (point estimates could be infinite after log transform)
  MU <- MU %*% tJ

  SIGMA <- SIGMA %.% tJ
  for(i in 1:nrow(SIGMA)) { for(j in 1:ncol(SIGMA)) { if(INF[i,j]) { SIGMA[i,j,] <- SIGMA[i,,j] <- 0; SIGMA[i,j,j] <- Inf } } }
  SIGMA <- aperm(SIGMA,c(1,3,2))
  SIGMA <- SIGMA %.% tJ
  for(i in 1:nrow(SIGMA)) { for(j in 1:ncol(SIGMA)) { if(INF[i,j]) { SIGMA[i,j,] <- SIGMA[i,,j] <- 0; SIGMA[i,j,j] <- Inf } } }
  dimnames(SIGMA) <- list(NULL,FEATURES,FEATURES)

  # can't make this infinite before or will mess up above calculations
  P <- c("minor","angle")
  for(i in 1:length(x))
  {
    if(!ISO && x[[i]]$isotropic)
    {
      INF[i,P] <- TRUE
      SIGMA[i,P,] <- 0
      SIGMA[i,,P] <- 0
      SIGMA[i,P,P] <- diag(Inf,2)
    }
  }

  # number of estimated features (amount of data)
  DIM <- length(FEATURES)

  make.names <- function(MEANS,VARS,OFF=FALSE)
  {
    MEANS <- rbind(MEANS)
    VARS <- rbind(VARS)
    OFF <- OFF && sum(VARS)>1 # no correlations present
    OFF <- ifelse(OFF,"COV","VAR")

    MEANS <- sapply(1:nrow(MEANS),function(i){ paste(FEATURES[MEANS[i,]],collapse=",") })
    VARS <- sapply(1:nrow(VARS),function(i){ paste(FEATURES[VARS[i,]],collapse=",") })

    NAMES <- paste0("E[",MEANS,"] ",OFF,"[",VARS,"]")

    return(NAMES)
  }

  # number of variance-covariance parameters modeled
  num.pars <- function(FIT,mean=FALSE)
  {
    V <- FIT$VARS
    V <- V[upper.tri(V,diag=TRUE)]
    V <- sum(V)
    if(mean) { V <- V + sum(FIT$MEANS) }
    return(V)
  }

  # data can potentially support these parameters only
  MEAN.MAX <- colSums(!INF)>0
  VAR.MAX <- colSums(!INF)>1
  # these mean parameters can't be turned off
  MEAN.MIN <- rep(TRUE,DIM)
  names(MEAN.MIN) <- FEATURES
  if("minor" %in% FEATURES) { MEAN.MIN[c("minor","angle")] <- FALSE } # can be mean zero
  MEAN.MIN[BETA] <- FALSE # can be mean-zero

  ##########
  ## build up phase
  ## minimal terms
  OLD <- list()
  MEANS <- MEAN.MIN
  VARS <- rep(FALSE,DIM)
  names(VARS) <- FEATURES

  S <- make.names(MEANS,VARS)
  if(trace) { message("Fitting covariance model ",S) }
  OLD[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=diag(VARS,DIM),MEANS=MEANS,WARN=FALSE,...)

  ## add terms one by one
  GUESS <- OLD[[1]]
  while(any(!VARS) || any(!MEANS))
  {
    NEW <- list()

    ########################
    ## add variance terms ##
    # which terms are presently off
    try <- which(!VARS & VAR.MAX & MEANS) # only add vars to means
    if(length(try))
    {
      TRYS <- array(VARS,c(DIM,length(try)))
      TRYS <- t(TRYS) # [off->on,all]
      colnames(TRYS) <- FEATURES
      for(i in 1:nrow(TRYS)) { TRYS[i,try[i]] <- TRUE }
      if("minor" %in% FEATURES) { TRYS[,"minor"] <- TRYS[,"angle"] <- TRYS[,"minor"] | TRYS[,"angle"] }
      TRYS <- unique(TRYS)

      # models that we will try by turning one term on
      NAMES <- make.names(MEANS,TRYS)
      # don't try what's been done
      SUB <- NAMES %nin% names(OLD)
      SUB <- which(SUB)

      TRYS <- TRYS[SUB,,drop=FALSE]
      NAMES <- NAMES[SUB]

      # fit all
      for(i in 1%:%length(NAMES))
      {
        S <- NAMES[i]
        if(trace) { message("Fitting covariance model ",S) }
        NEW[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=diag(TRYS[i,],DIM),MEANS=MEANS,GUESS=GUESS,WARN=FALSE,...)
      }
    }

    ####################
    ## add mean terms ##
    # which terms are presently off
    try <- which(!MEANS & MEAN.MAX)
    if(length(try))
    {
      TRYS <- array(MEANS,c(DIM,length(try)))
      TRYS <- t(TRYS) # [off->on,all]
      colnames(TRYS) <- FEATURES
      for(i in 1:nrow(TRYS)) { TRYS[i,try[i]] <- TRUE }
      if("minor" %in% FEATURES) { TRYS[,"minor"] <- TRYS[,"angle"] <- TRYS[,"minor"] | TRYS[,"angle"] }
      TRYS <- unique(TRYS)

      # models that we will try by turning one term on
      NAMES <- make.names(TRYS,VARS)
      # don't try what's been done
      SUB <- NAMES %nin% names(OLD)
      SUB <- which(SUB)

      TRYS <- TRYS[SUB,,drop=FALSE]
      NAMES <- NAMES[SUB]

      # fit all
      for(i in 1%:%length(NAMES))
      {
        S <- NAMES[i]
        if(trace) { message("Fitting covariance model ",S) }
        NEW[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=diag(VARS,DIM),MEANS=TRYS[i,],GUESS=GUESS,WARN=FALSE,...)
      }
    }

    ###############
    # wrap up
    if(!length(NEW)) { break } # no new models to fit

    ICS <- sapply(NEW,function(m){m[[IC]]})
    BEST <- which.min(ICS)
    BEST <- NEW[[BEST]]
    GUESS <- BEST
    MEANS <- BEST$MEANS
    VARS <- diag(BEST$VARS)

    OLD <- c(OLD,NEW)

    if(BEST[[IC]]==Inf) { break }
  } # finish build up variances

  ###########
  ## pair down terms phase
  # start with all variances or best (if IC==Inf)
  NAME.ALL <- make.names(MEAN.MAX,VAR.MAX)
  if(NAME.ALL %in% names(OLD))
  { BEST <- OLD[[NAME.ALL]] }
  else # limited data cannot support all variances
  {
    ICS <- sapply(OLD,function(m){m[[IC]]})
    BEST <- which.min(ICS)
    BEST <- OLD[[BEST]]
  }
  GUESS <- BEST
  MEANS <- BEST$MEANS
  VARS <- diag(BEST$VARS)

  ##############
  # pair down variances phase
  while(sum(VARS)>1 || sum(MEANS)>1)
  {
    NEW <- list()
    NAMES1 <- NAMES2 <- NULL

    #################
    ### var terms ###
    # which terms are presently off
    try <- which(VARS)
    if(length(try))
    {
      TRYS <- array(VARS,c(DIM,length(try)))
      TRYS <- t(TRYS) # [off->on,all]
      colnames(TRYS) <- FEATURES
      for(i in 1:nrow(TRYS)) { TRYS[i,try[i]] <- FALSE }
      if("minor" %in% FEATURES) { TRYS[,"minor"] <- TRYS[,"angle"] <- TRYS[,"minor"] & TRYS[,"angle"] }
      TRYS <- unique(TRYS)
      # models that we will try by turning one term off
      NAMES <- NAMES1 <- make.names(MEANS,TRYS)
      # paired down models that we haven't attempted yet
      SUB <- NAMES %nin% names(OLD)
      NEW.NAMES <- NAMES[SUB]
      TRYS <- TRYS[SUB,,drop=FALSE]

      # fit all
      for(i in 1%:%length(NEW.NAMES))
      {
        S <- NEW.NAMES[i]
        if(trace) { message("Fitting covariance model ",S) }
        NEW[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=diag(TRYS[i,],DIM),MEANS=MEANS,GUESS=GUESS,WARN=FALSE,...)
      }
    }

    #################
    ### mean terms ###
    # which terms are presently off
    try <- which(MEANS & !MEAN.MIN & !VARS)
    if(length(try))
    {
      TRYS <- array(MEANS,c(DIM,length(try)))
      TRYS <- t(TRYS) # [off->on,all]
      colnames(TRYS) <- FEATURES
      for(i in 1:nrow(TRYS)) { TRYS[i,try[i]] <- FALSE }
      if("minor" %in% FEATURES) { TRYS[,"minor"] <- TRYS[,"angle"] <- TRYS[,"minor"] & TRYS[,"angle"] }
      TRYS <- unique(TRYS)
      # models that we will try by turning one term off
      NAMES <- NAMES2 <- make.names(TRYS,VARS)
      # paired down models that we haven't attempted yet
      SUB <- NAMES %nin% names(OLD)
      NEW.NAMES <- NAMES[SUB]
      TRYS <- TRYS[SUB,,drop=FALSE]

      # fit all
      for(i in 1%:%length(NEW.NAMES))
      {
        S <- NEW.NAMES[i]
        if(trace) { message("Fitting covariance model ",S) }
        NEW[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=diag(VARS,DIM),MEANS=TRYS[i,],GUESS=GUESS,WARN=FALSE,...)
      }
    }

    ###############
    ## wrap up ##
    OLD <- c(OLD,NEW)
    NEW <- OLD[c(NAMES1,NAMES2)] # also consider models that we've already fitted (with the right number of parameters)
    if(!length(NEW)) { break }

    ICS <- sapply(NEW,function(m){m[[IC]]})
    BEST <- which.min(ICS)
    BEST <- NEW[[BEST]]
    GUESS <- BEST
    MEANS <- BEST$MEANS
    VARS <- diag(BEST$VARS)
  } # finish pair down variances

  # take best variance model
  ICS <- sapply(OLD,function(m){m[[IC]]})
  I <- which.min(ICS)
  BEST <- OLD[[I]]
  NP <- num.pars(BEST,mean=TRUE)
  NV <- num.pars(BEST,mean=FALSE)

  # now see if correlations are supported
  if(NV>1 && sum(!INF)>=NP)
  {
    GUESS <- BEST
    MEANS <- BEST$MEANS
    VARS <- diag(BEST$VARS)
    S <- make.names(MEANS,VARS,OFF=TRUE)
    if(S %nin% names(OLD))
    {
      if(trace) { message("Fitting covariance model ",S) }
      OLD[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,VARS=VARS,MEANS=MEANS,GUESS=GUESS,WARN=FALSE,...)
    }
  } # end correlation fit

  ICS <- sapply(OLD,function(m){m[[IC]]})
  names(ICS) <- names(OLD)
  I <- sort(ICS,index.return=TRUE)$ix
  ICS <- ICS[I]
  ICS <- ICS - ICS[1]
  I <- I[1]

  if(trace)
  {
    # report model selection
    ICS <- data.frame(ICS)
    colnames(ICS) <- paste0("\u0394",IC)
    message("* Model selection for autocovariance distribution.")
    print(ICS)
  }

  R <- OLD[[I]]
  MEANS <- R$MEANS
  VARS <- diag(R$VARS)
  R$name <- names(OLD)[I]
  R$isotropic <- !("minor" %in% FEATURES && MEANS["minor"]) # just the mean
  R$axes <- axes

  names(R)[ which(names(R)=="mu") ] <- "par" # population mean of features
  names(R)[ which(names(R)=="COV.mu") ] <- "COV" # uncertainty in population mean
  names(R)[ which(names(R)=="sigma") ] <- "POV" # population dispersion of features (under log)
  names(R)[ which(names(R)=="COV.sigma") ] <- "COV.POV" # uncertainty in population dispersion of features (under log)

  # transform from diagonal-ish basis for covariance model
  if(any(VARS))
  {
    J <- solve(J)
    tJ <- t(J)
    R$par <- (R$par %*% tJ)[1,]
    R$COV <- J %*% R$COV %*% tJ
    R$POV <- J %*% R$POV %*% tJ
    # SUB <- FEATURES[variance]
    # J <- J[SUB,SUB]
    J <- quad2lin(J)
    dimnames(J) <- dimnames(R$COV.POV)
    R$COV.POV <- J %*% R$COV.POV %*% t(J)
  }

  # information for fitting that we no longer use
  NIN <- !MEANS & !VARS
  if(any(NIN))
  {
    SUB <- !NIN
    NIN <- FEATURES[NIN]
    FEATURES <- FEATURES[SUB]
    R$par <- R$par[FEATURES]
    VARS <- VARS[SUB]
    # R$par <- NULL
    R$COV <- R$COV[FEATURES,FEATURES]
    R$POV <- R$POV[FEATURES,FEATURES]
    # COV.POV should be okay except when missing POV
    FEAT <- sapply(NIN,function(nin){grepl(nin,rownames(R$COV.POV),fixed=TRUE)})
    FEAT <- !apply(FEAT,1,any)
    R$COV.POV <- R$COV.POV[FEAT,FEAT]
  }
  R$features <- FEATURES

  # set beta structure (empty)
  BETA <- BETA[BETA %in% FEATURES]
  beta <- rep(0,length(BETA))
  names(beta) <- BETA
  if(length(BETA))
  {
    R$beta <- beta
    R$formula <- stats::as.formula(paste("~",paste(names(beta),collapse=" + ")))
  }

  # reverse log transformation
  # R <- set.parameters(R,par,linear.cov=TRUE)
  R <- exp_ctmm(R,debias=debias,variance=VARS)

  return(R)
}


##########
mean.mu <- function(x,debias=TRUE,weights=NULL,trace=FALSE,IC="AICc",...)
{
  if(is.null(weights))
  { weights <- rep(1,length(x)) }

  axes <- x[[1]]$axes
  AXES <- length(axes)
  N <- length(x)

  # Gaussian-Gaussian in all cases
  MU <- array(0,c(N,AXES))
  SIGMA <- array(0,c(N,AXES,AXES))

  colnames(MU) <- axes
  dimnames(SIGMA) <- list(NULL,axes,axes)

  for(i in 1:N)
  {
    MU[i,] <- x[[i]]$mu
    SIGMA[i,,] <- x[[i]]$COV.mu

    # TODO !!!
    # fill in with zeroes for non-stationary means
    # TODO !!!
  }

  # list of candidate models
  MM <- list()

  # no variance in mean location
  S <- "Dirac-\u03B4(\u03BC)"
  if(trace) { message("Fitting location-mean model ",S) }
  MM[[S]] <- meta.normal(MU,SIGMA,VARS=FALSE,isotropic=TRUE,debias=debias,weights=weights,WARN=FALSE,...)
  GUESS <- MM[[S]]

  # symmetric mean distribution
  if(2*N >= 2+1)
  {
    S <- "isotropic-\u03BC"
    if(trace) { message("Fitting location-mean model ",S) }
    MM[[S]] <- meta.normal(MU,SIGMA,isotropic=TRUE,debias=debias,weights=weights,GUESS=GUESS,WARN=FALSE,...)
    GUESS <- MM[[S]]

    # general distribution
    if(2*N >= 2+3)
    {
      S <- "anisotropic-\u03BC"
      if(trace) { message("Fitting location-mean model ",S) }
      MM[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,GUESS=GUESS,WARN=FALSE,...)
    }
  }

  ICS <- sapply(MM,function(m){m[[IC]]})
  names(ICS) <- names(MM)
  i <- which.min(ICS)
  i <- names(MM)[i]
  MM <- MM[[i]]
  MM$name <- i

  # report model selection
  if(trace)
  {
    i <- sort(ICS,index.return=TRUE)$ix
    ICS <- ICS[i] # sorted
    ICS <- ICS - ICS[1] # start at zero
    ICS <- data.frame(ICS)
    colnames(ICS) <- paste0("\u0394",IC)
    message("* Model selection for location-mean \u03BC distribution.")
    print(ICS)
  }

  # R$mu # population mean of mean locations
  # R$COV.mu # uncertainty in mean of means estimate
  names(MM)[ which(names(MM)=="sigma") ] <- "POV.mu" # dispersion of means
  names(MM)[ which(names(MM)=="COV.sigma") ] <- "COV.POV.mu" # uncertainty in dispersion of means

  return(MM)
}


###########
mean.ctmm <- function(x,weights=NULL,sample=TRUE,debias=TRUE,IC="AIC",trace=TRUE,...)
{
  select <- "all"

  # if(length(x)==1)
  # {
  #   warning("Only one individual in list.")
  #   return(x)
  # }

  n <- length(x)
  info <- mean_info(x)
  axes <- x[[1]]$axes

  if(is.null(weights))
  { weights <- rep(1,n) }
  else
  { weights <- weights/max(weights) }
  names(weights) <- names(x)

  isotropic <- c(TRUE,TRUE)
  names(isotropic) <- c("sigma","mu")

  #######################
  # average of mean locations
  if(sample)
  {
    MM <- mean.mu(x,debias=debias,weights=weights,IC=IC,trace=trace,...)
    isotropic['mu'] <- MM$isotropic

    FM <- mean.features(x,debias=debias,weights=weights,select=select,IC=IC,trace=trace,...)
    isotropic['sigma'] <- FM$isotropic

    x <- copy(from=FM,to=MM)
    x$isotropic <- isotropic

    # STORE EVERYTHING
    x$name <- paste(MM$name,FM$name)
    x$AIC <- MM$AIC + FM$AIC
    x$AICc <- MM$AICc + FM$AICc
    x$BIC <- MM$BIC + FM$BIC
  } # end sample
  else # !sample
  {
    # COV[mu^mu] diagonal-term
    cov.diag <- function(cov)
    {
      COV <- cov^2
      VEC <- c( diag(cov)*cov[1,2] , cov[1,1]*cov[2,2]+cov[1,2]^2 )
      COV <- rbind(COV,VEC[1:2])
      COV <- cbind(COV,VEC)
      COV <- 2*COV
      return(COV)
    }

    cov.off <- function(cov1,cov2)
    {
      COV <- 4 * cov1 * cov2
      VEC <- c(2*cov1[1,1]*cov2[1,2]+2*cov1[1,2]*cov2[1,1],2*cov1[2,2]*cov2[1,2]+2*cov1[1,2]*cov2[2,2],cov1[1,1]*cov2[2,2]+cov1[2,2]*cov2[1,1]+2*cov1[1,2]*cov2[1,2])
      COV <- rbind(COV,VEC[1:2])
      COV <- cbind(COV,VEC)
      return(COV)
    }

    MU <- COV.MU <- COV <- COV.COV <- 0
    w <- weights/sum(weights)
    for(i in 1:n)
    {
      MU <- MU + w[i]*x[[i]]$mu
      COV.MU <- COV.MU + w[i]^2*x[[i]]$COV.mu
      COV <- COV + w[i]*( x[[i]]$sigma + outer(c(x[[i]]$mu)) ) # M2
      COV.COV <- COV.COV + w[i]^2*sigma.COV(x[[i]]) + (w[i]-w[i]^2)^2*cov.diag(x[[i]]$COV.mu)
      for(j in (i+1)%:%n) { COV.COV <- COV.COV - w[i]^2*w[j]^2*cov.off(x[[i]]$COV.mu,x[[j]]$COV.mu) }
    }
    COV <- COV - outer(c(MU)) # M2 - M1^2
    rownames(COV.MU) <- colnames(COV.MU) <- axes

    isotropic['mu'] <- FALSE
    if(!all(sapply(x,function(y){y$isotropic[1]}))) { isotropic['sigma'] <- FALSE }

    sigma <- covm(COV,axes=axes,isotropic=isotropic['sigma'])
    if(isotropic['sigma'])
    {
      COV <- COV.COV[1,1,drop=FALSE]
      rownames(COV) <- colnames(COV) <- 'major'
      features <- 'major'
    }
    else
    {
      J <- J.par.sigma(sigma[c(1,4,2)])
      COV <- J %*% COV.COV %*% t(J)
      features <- c('major','minor','angle')
    }

    x <- ctmm(mu=MU,COV.mu=COV.MU,sigma=sigma,COV=COV,features=features)
    # TODO non-sigma features
  } # end !sample

  x$info <- info
  x <- do.call(ctmm,x)

  x$isotropic <- isotropic
  x$weights <- weights

  # return final result
  # x <- new.ctmm(x,info)
  return(x)
}


## population stationary distribution
mean_pop <- function(CTMM)
{
  # spread (x,y)
  sigma <- CTMM$POV.mu + CTMM$sigma
  # uncertainty of spread (x-x,x-y,y-y)
  if(CTMM$isotropic['mu'])
  { COV <- c(CTMM$COV.POV.mu) * diag(c(1,1,0)) }
  else
  {
    P <- c(1,3,2) # xx,yy,xy
    COV <- CTMM$COV.POV.mu[P,P]
  }

  if(CTMM$isotropic['sigma'])
  {
    P <- 'major'
    COV <- COV + c(CTMM$COV[P,P]+CTMM$POV[P,P])*diag(c(1,1,0))
  }
  else
  {
    P <- c('major','minor','angle')
    J <- J.sigma.par(CTMM$sigma@par)
    COV <- COV + J%*%(CTMM$COV[P,P]+CTMM$POV[P,P])%*%t(J)
  }

  isotropic <- all(CTMM$isotropic)
  sigma <- covm(sigma,isotropic=isotropic)

  if(isotropic)
  {
    P <- 'major'
    COV <- COV[1,1,drop=FALSE]
  }
  else
  {
    P <- c('major','minor','angle')
    J <- J.par.sigma(sigma[c(1,4,2)])
    COV <- J %*% COV %*% t(J)
  }
  dimnames(COV) <- list(P,P)

  CTMM$isotropic <- isotropic
  CTMM$sigma <- sigma
  CTMM$COV <- COV
  CTMM$POV <- CTMM$COV.POV <- CTMM$COV.POV.mu <- NULL
  CTMM$features <- rownames(COV)
  CTMM$tau <- NULL
  CTMM$circle <- FALSE

  # BETA ARE TOTALLY UNFINISHED
  # will need eventually...

  return(CTMM)
}