R/mixAR.R

In mixAR: Mixture Autoregressive Models

## Do not edit this file manually.
## It has been automatically generated from mixAR.org.

# setGeneric("coef")
# setMethod("coef", "mle", function(object) object@fullcoef )
# setMethod("coef", "summary.mle", function(object) object@coef )

setGeneric("lik_params", 
           function(model)
              standardGeneric("lik_params"), 
           useAsDefault=FALSE)

setGeneric("make_fcond_lik", 
           function(model, ts)
              standardGeneric("make_fcond_lik"), 
           useAsDefault=FALSE)

setGeneric("mix_hatk", 
           function(model, x, index, xcond)
              standardGeneric("mix_hatk"), 
           useAsDefault=FALSE)

setGeneric("mix_ek", 
           function(model, x, index, xcond, scale)
              standardGeneric("mix_ek"), 
           useAsDefault=FALSE)

setGeneric("mix_pdf", 
           function(model, x, index, xcond, scale)
              standardGeneric("mix_pdf"), 
           useAsDefault=FALSE)

setGeneric("mix_cdf", 
           function(model, x, index, xcond, scale)
              standardGeneric("mix_cdf"), 
           useAsDefault=FALSE)

setGeneric("mix_qf",  # 2021-04-30 new
           function(model, p, x, index, xcond)
              standardGeneric("mix_qf"), 
           useAsDefault=FALSE)

setClass("MixAR", 
         representation(prob = "numeric", 
                        order = "numeric", 
                        shift = "numeric", 
                        scale = "numeric", 
                        arcoef = "raggedCoef"
                        # todo: dist,  pdf, cdf, qdf
                        #       (This note is from 2012, is it still relevant?)
                        ), 
         contains="VIRTUAL"
         )

## zasega default stoynosti na prob sa NA, na shift 0, na scale 1.
##     todo: dali vsichki da sa NA?
setMethod("initialize", 
    signature(.Object = "MixAR"), 
    function (.Object, arcoef, order, prob, shift, scale, model, ...)#2012-10-30 new arg model
    {
        mis_order <- missing(order)
        if(missing(model)){    # 2012-10-30 lazy; but do not want to risk breaking this now.
            if(missing(arcoef))                   # new 2011-12-02
                arcoef <- raggedCoef(p = order)  # one of 'arcoef' and 'order' must be present
            if(is.list(arcoef) && any(names(arcoef) == "s")){
                ## 13-09-2018 Updated to create "raggedCoefS" objects
                ## Note: previous code not touched
                if(length(arcoef) != 3)
                    stop("Incorrect number of elements in the list for mixSAR model.")
                if(all(names(arcoef) != "a") && all(names(arcoef) != "as")){
                    names(arcoef)[which(names(arcoef) != "s")] <- c("a", "as")
                }else{
                    if(all(names(arcoef) != "a"))
                        names(arcoef)[ which(names(arcoef) == "")] <- "a"
                    if(all(names(arcoef) != "as"))
                        names(arcoef)[ which(names(arcoef) == "")] <- "as"
                }

                if(is.matrix(arcoef$a)){
                    arcoef$a <- split(arcoef$a, row(arcoef$a))
                }

                if(is.matrix(arcoef$as)){
                    arcoef$as <- split(arcoef$as, row(arcoef$as))
                }
                if(length(arcoef$a) != length(arcoef$as))
                    stop("Incorrect imputation of coefficient")

                arcoef <- new("raggedCoefS", a = arcoef$a, 
                              as = arcoef$as, s = arcoef$s)
            }else{
                if(!is(arcoef, "raggedCoef")){                 # 2012-09-04 dobavyam "matrix"
                    if(is.matrix(arcoef)){
                        arcoef <- split(arcoef, row(arcoef)) # turn into a list of rows
                        if(!mis_order){
                            arcoef <- lapply(1:length(arcoef), 
                                             function(x) arcoef[[x]][ seq_len(order[x]) ])
                        }
                    }
                    arcoef <- new("raggedCoef", arcoef) # hoping 'new' knows how to do this...
                }
            }
            wrkorder <- row_lengths(arcoef)
            if(!missing(order) &&
               (length(order) != length(wrkorder) || !all(order == wrkorder))){
                message("Arg's 'arcoef' and 'order' are not consistent, ignoring 'order'.")
            }
            order <- wrkorder

            ## adjust 'prob' if necessary
            if(missing(prob))
                prob <- rep(as.numeric(NA), length(order))
            else if(length(prob) == 1)
                ## 2020-03-25 TODO: shouldn't this be set to equal prob's?
                prob <- rep(prob, length(order))
            else if(length(prob) != length(order))
                stop("length of 'prob' should be 1 or the length of 'order'")

            ## adjust 'shift' if necessary
            if(missing(shift))
                shift <- rep(0, length(order))
            else if(length(shift) == 1)
                shift <- rep(shift, length(order))
            else if(length(shift) != length(order))
                stop("length of 'shift' should be 1 or the length of 'order'")

            ## adjust 'scale' if necessary
            if(missing(scale))
                scale <- rep(1, length(order))
            else if(length(scale) == 1)
                scale <- rep(scale, length(order))
            else if(length(scale) != length(order))
                stop("length of 'scale' should be 1 or the length of 'order'")
        }else{ # arg. 'model' supplied
            ncomp <- length(model@order)

            if(missing(order))
                order <- model@order
            else if(length(order) != ncomp || any(order != model@order))
                stop("Arg. 'order' (if present) must be the same as 'model@order'.")
            ## else order is same as model@order, nothing to do

            if(missing(prob))
                prob <- model@prob
            else if(length(prob) == 1)
                prob <- rep(prob, ncomp)
            else if(length(prob) != ncomp)
                stop("Arg. 'prob' (if present) must have the same length as 'model@prob'.")
            # else use 'prob' as is

            if(missing(shift))
                shift <- model@shift
            else if(length(shift) == 1)
                shift <- rep(shift, ncomp)
            else if(length(shift) != ncomp)
                stop("Arg. 'shift' (if present) must have the same length as 'model@shift'.")
            # else use 'shift' as is

            if(missing(scale))
                scale <- model@scale
            else if(length(scale) == 1)
                scale <- rep(scale, ncomp)
            else if(length(scale) != ncomp)
                stop("Arg. 'scale' (if present) must have the same length as 'model@scale'.")
            # else use 'scale' as is

            if(missing(arcoef))
                arcoef <- model@arcoef
            ## else{ ##arcoef supplied on top of model
            ##13-09-2018 Updated to create "raggedCoefS" objects
            ##Note: previous code not touched
            ##  if(is.list(arcoef) & any(names(arcoef) == "s")){
            ##   if(length(arcoef)!=3)stop(
            ##    "Incorrect number of elements in the list to build mixSAR model.
            ##   See help page for details")
            ##    if(all(names(arcoef)!="a") && all(names(arcoef)!="as")){
            ##     names(arcoef)[which(names(arcoef)!="s")] <- c("a", "as")
            ##    }else{
            ##      if(all(names(arcoef)!="a")) names(arcoef)[
            ##        which(names(arcoef) == "")] <- "a"
            ##  if(all(names(arcoef)!="as")) names(arcoef)[
            ##       which(names(arcoef) == "")] <- "as"
            ##   }
            
            ##   if(is.matrix(arcoef$a)){
            ##    arcoef$a <- split(arcoef$a, row(arcoef$a))
            ##   }
            
            
            ##   if(is.matrix(arcoef$as)){
            ##     arcoef$as <- split(arcoef$as, row(arcoef$as))
            
            ##  }
            ##   if(length(ar$a)!=length(ar$as))stop(
            ##    "Incorrect imputation of coefficient")
            ## TODO: This is wrong, I don't see variables a and s
            ##   arcoef <- new("raggedCoefS", a=a, as=as, s=s)
            ##If re-established, remember to add } before else
            else{
                if(!is(arcoef, "raggedCoef")){                  # 2012-09-04 dobavyam "matrix"
                    if(is.matrix(arcoef)){
                        arcoef <- split(arcoef, row(arcoef)) # turn into a list of rows
                        if(!mis_order){
                            arcoef <- lapply(1:length(arcoef), 
                                             function(x) arcoef[[x]][ seq_len(order[x]) ])
                        }
                    }
                    arcoef <- new("raggedCoef", arcoef) # hoping 'new' knows how to do this
                }
            }
            wrkorder <- row_lengths(arcoef)
            if(length(order) != length(wrkorder) || any(order != wrkorder))
                stop("Arg. 'arcoef' (if present) must be consistent with 'model@order'.")
        }
        callNextMethod(.Object, arcoef = arcoef, order = order, 
                       prob = prob, shift = shift, scale = scale, ...)
    }
)

# 2012-10-31 commenting out. It turns out that setClass now returns a generator which is
#                           naturally to be named after the class. Better to avoid confusion.
#                          (at the same time, assigning the return value of setClass to check)
# 2012-10-31 removing the comments. mixAR is virtual, so cannot create objects from it.
#            But it makes sense to use the name for something more useful.

                                            # this has diferent syntax from the "new"
                                            # function. Also, the defaults are NA, not zero.
mixAR <- function(template, coef, ..., filler = NA_real_){       # todo: make more versatile
    if(missing(coef))
        coef <- list()

    if(!missing(template)){
        coef$order <- template
    }

    coef <- .canonic_coef(coef, filler)

    g <- length(coef$order)
    if(g == 1)
        stop("A mixAR model must have two or more components.")

    new("MixARGaussian", arcoef = coef$arcoef, order = coef$order,
                         prob = coef$prob, shift = coef$shift, scale = coef$scale, ...)
}

setGeneric("mixAR")    # note: the name of the 1st arg. is appropriate for the default only.

setMethod("mixAR", signature(template = "MixAR"),
          function(template, coef, ..., filler = NA){
              model <- template
              if(missing(coef))
                  return(model)

              ## TODO: need to check that number of components in coef are consistent with
              ##       ncomp
              ncomp <- .nmix(model)
              for(s in c("prob", "shift", "scale", "arcoef", "order")){
                  val <- coef[[s]]
                  ## TODO: more care needed here
                  if(!is.null(val)){
                      if(length(val) == 1)
                          val <- rep(val, ncomp)
                      slot(model, s) <- val
                  }
              }

              model
          })

.canonic_coef <- function(coef, filler){
    wrk <- coef
    len <- sapply(c("prob", "shift", "scale", "arcoef", "order"),
                  function(x){ if(is.null(wrk[[x]]))
                                   wrk[[x]] <<- filler    # non-local assignment!
                               length(wrk[[x]])})
    len <- unique(len)
    g <- len[len > 1]
    if(length(g) > 1)
        stop("Conflicting lengths of arguments.")

    ## 2018-11-17: bug-fix: changed "<-" to "<<-"
    lapply(c("prob", "shift", "scale"),
           function(x) if(length(wrk[[x]]) == 1) wrk[[x]] <<- rep(wrk[[x]], g))

    wrk$arcoef <- if(is.null(wrk$order))
                        raggedCoef(value = wrk$arcoef)
                   else
                       raggedCoef(p = wrk$order, value = wrk$arcoef)
    wrk
}

setMethod("show", "MixAR", # 12-09-2018 "show" method adapted to handle class "raggedCoefS"
          function(object) { #note: previous content has not been changed, only "if" were added
              cl <- class(object)
              g <- length(object@prob)
              p <- max(object@order)
              mcoefs <- NULL
              if(inherits(object@arcoef, "raggedCoefS")){
                  ps <- max(object@arcoef@ps)
                  s  <- object@arcoef@s
              }
              ## cat("(To see the internal structure of the object, use function 'str'.)\n\n")
              cat("An object of class \"", cl, "\"\n", sep="")
              cat("Number of components:", length(object@prob), "\n")
              if(p > 0){
                  if(inherits(object@arcoef, "raggedCoefS")){
                      mcoefs <- ragged2char(object@arcoef@as)
                      colnames(mcoefs) <- paste("ar_", seq_len(ps)*s, sep="")
                      wrk <- lapply(1:g, function(i) c(object@prob[i], object@shift[i], 
                                                       object@scale[i], object@order[i], 
                                                       object@arcoef@a[[i]]))#), object@arcoef@ps[i]))
                      mcoef            <- cbind(ragged2char(wrk), as.character(object@arcoef@ps))
                      colnames(mcoef)  <- c("prob", "shift", "scale", "order", 
                                            paste("ar_", seq_len(p), sep=""), "s_order")
                  }else{
                      wrk <- lapply(1:g, function(i) c(object@prob[i], object@shift[i], 
                                                       object@scale[i], object@order[i], 
                                                       object@arcoef@a[[i]]))
                      mcoef <- ragged2char(wrk)
                      colnames(mcoef) <- c("prob", "shift", "scale", "order", 
                                           paste("ar_", seq_len(p), sep=""))
                  }
                  rownames(mcoef) <- paste("Comp_", 1:g, sep="")
                  print(cbind(mcoef, mcoefs), na.print="", quote=FALSE, justify="right")
              }else{
                  m <- cbind( prob = object@prob
                           , shift = object@shift
                           , scale = object@scale
                           , order = object@order
                             )
                  rownames(m) <- paste("Comp_", 1:g, sep="")
                  cat("The AR orders of all components are 0.\n")
                  print(m)
              }
              invisible(object)
          })

setAs("MixAR", "list",
      function(from){
              # nams <- slotNames(from)
              # structure(lapply(nams, function(x) slot(from, x)),
              #           names = nams)
          list(prob = from@prob, 
               order = from@order,
               shift = from@shift,
               scale = from@scale, 
               arcoef = from@arcoef@a ) # notice @a
      })

setMethod("lik_params", c(model = "MixAR"), 
          function(model){
              k <- .nmix(model)   # note: this is for estimation purposes, 
                                        # only k-1 of the mixing probs are passed on.
              c( model@prob[-k],  model@scale, model@shift, ragged2vec(model@arcoef)
                )
          })

setMethod("mix_ncomp", signature(x="MixAR"), 
          function(x){
              length(x@prob)
          })

setMethod("row_lengths", signature(x="MixAR"), 
          function(x){
              x@order
          })

## 2012-11-11 new
.nmix <- function(model) length(model@prob) # number of components in a mixAR model

setMethod("make_fcond_lik", signature(model="MixAR", ts="numeric"), 
          function(model, ts){            ## mixar order  is assumed constant in this function.
              locts <- ts
              mixar <- model
              k <- mix_ncomp(mixar)  # 2011-07-13: was length(mixar@prob)
              ## n <- length(lik_params(mixar))

              ## todo: is it worth making this function sufficiently general, so that
              ##       subclasses of mixAR do not need to provide their own version when they
              ##       have additional parameters?
              ind_prob   <- seq_len(k-1)            # these need to be in sync with lik_params
              ind_scale  <-      k:(k+k-1)
              ind_shift  <-  (2*k):(2*k+k-1)
              ind_arcoef <- - c( ind_prob, ind_scale, ind_shift)       ## note: negative index

              f <- function(x){
                  mixar@prob <- c(x[ind_prob], 1 - sum(x[ind_prob]))
                  mixar@scale <- x[ind_scale]
                  mixar@shift <- x[ind_shift]
                  mixar@arcoef <- rag_modify(mixar@arcoef, x[ind_arcoef])

                  res <- cond_loglik(mixar, locts)

                  ## 2020-06-12 was: cat("x = ", x, "\n")       
                  ##                 cat("f(x) = ", res, "\n\n")

                  res
              }
              f
          })

setMethod("mix_hatk", signature(model="MixAR", x="numeric", index="numeric", xcond="missing"), 
          function(model, x, index, xcond) {
              mixFilter(x, model@arcoef, index, shift=model@shift)
          })

setMethod("mix_ek", signature(model="MixAR", x="numeric", index="numeric", xcond="missing", 
                              scale="missing"), 
          function(model, x, index, xcond, scale) {
              mixFilter(x, model@arcoef, index, shift=model@shift, 
                        residual=TRUE)
          })

setMethod("mix_ek", signature(model="MixAR", x="numeric", index="numeric", xcond="missing", 
                              scale="logical"), 
          function(model, x, index, xcond, scale) {
              mixFilter(x, model@arcoef, index, shift=model@shift, 
                        residual=TRUE, scale=model@scale)
          })

## kogato xcond e dadeno, x se interpretira kato stoynosti na x_t za koito da se smyata
## Prognoz se smyata samo za xcond. Zasega samo edna stoynost na xcond.
##    Primerno prilozhenie - plot na pdf ili cdf pri dadeni xcond.
setMethod("mix_ek", signature(model="MixAR", x="numeric", index="missing", xcond="numeric", 
                              scale="missing"), 
          function(model, x, index, xcond) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              x - wrk      ## todo: check!
          })

setMethod("mix_ek", signature(model="MixAR", x="numeric", index="missing", xcond="numeric", 
                              scale="logical"), 
          function(model, x, index, xcond, scale) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              (x - wrk) / model@scale
          })

mixARgen <- setClass("MixARgen", representation(dist = "list" ), contains="MixAR")

setMethod("initialize", 
    signature(.Object = "MixARgen"), 
    function (.Object, dist, model, ...)   # 2012-10-30 new arg. 'model'
    {                                  # todo: investigate the following:
        # .Object <- callNextMethod()  # predava samo '...' args, seemingly contradicts doc.
        .Object <- if(missing(model)) callNextMethod()
                   else               callNextMethod(.Object, model = model, ...)
        if(!missing(dist)){
            if(is.list(dist))
                .Object@dist <- dist
            else
                stop("Argument 'dist' must be a list.")
        }else if(.hasSlot(model, "dist"))
            .Object@dist <- model@dist
        # else if(is(model, "MixARGaussian")) TODO: makes sense to set Gaussian components.
        else
            stop("Cannot set the distribution since argument 'dist' is missing.")

        .Object
    }
)

setMethod("show", "MixARgen",            # todo: this is only a quick fix, more work is needed
	  function(object) {
              callNextMethod()
                                       # cat("Parameters of the component distributions:\n\t")
              cat("\n")
              cat("Distributions of the error components:\n")
                                        # print(noise_params(object))
              b_show(get_edist(object))
              cat("\n")

              invisible(object)
	  })

setMethod("lik_params", c(model = "MixARgen"), 
          function(model){
              res <- callNextMethod()   # get the "usual" params
              ## 2017-04-21 was:  wrk <- mix_noiseparams(model) # get the noise params
              ##     it seems that noise_params() is the current name of the function:
              wrk <- noise_params(model) # get the noise params, TODO: option to name them?
              if(length(wrk)>0)
                  res <- c(res, wrk)
              res
          })

# 2011-07-17 pravya mix_pdf etc za "MixAR" (predi byacha samo za MixARGaussian)
## cdf
setMethod("mix_cdf", signature(model="MixARgen", x="numeric", index="numeric", 
                               xcond="missing", scale="ANY"), # note: scale is ignored here!
          function(model, x, index, xcond, scale) {
              ## 2020-06-12 was: print("gen!")
              cdf <- noise_dist(model, "cdf")

              wrk <- mix_ek(model, x, index, scale=TRUE)
              wrk <- cdf %of% wrk
              wrk <- inner(wrk, model@prob)
              wrk
          })

setMethod("mix_cdf", signature(model="MixARgen", x="numeric", index="missing", 
                               xcond="numeric", scale="ANY"), 
          function(model, x, index, xcond, scale) {
              ## 2020-06-12 was: print("gen!")
              cdf <- noise_dist(model, "cdf")

              wrk <- mix_hatk(model, xcond, length(xcond)+1)
                 # browser()
              wrk <- (x - wrk)/model@scale
              wrk <- cdf %of% wrk
              wrk <- inner(wrk, model@prob)
              wrk
          })

setMethod("mix_cdf", signature(model="MixARgen", x="missing", index="missing", 
                               xcond="numeric", scale="ANY"), 
          function(model, x, index, xcond, scale) {
              ## 2020-06-12 was: print("gen!")
              cdf <- noise_dist(model, "cdf")
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              f <- function(x){
                      wrk <- (x - wrk)/model@scale
                      wrk <- cdf %of% wrk
                      wrk <- inner(wrk, model@prob)
                      wrk
                  }
              f
          })

## pdf
setMethod("mix_pdf", signature(model="MixARgen", x="numeric", index="numeric", 
                               xcond="missing", scale="ANY"),  # note: scale is ignored here!
          function(model, x, index, xcond, scale) {
              pdf <- noise_dist(model, "pdf")

              wrk <- mix_ek(model, x, index, scale=TRUE)
              wrk <- pdf %of% wrk
              wrk <- inner(wrk, model@prob/model@scale)
              wrk
          })


setMethod("mix_pdf", signature(model="MixARgen", x="numeric", index="missing", 
                               xcond="numeric", scale="ANY"), 
          function(model, x, index, xcond, scale) {
              pdf <- noise_dist(model, "pdf")

              wrk <- mix_hatk(model, xcond, length(xcond)+1)
                 # browser()
              wrk <- (x - wrk)/model@scale
              wrk <- pdf %of% wrk
              wrk <- inner(wrk, model@prob/model@scale)
              wrk
          })


setMethod("mix_pdf", signature(model="MixARgen", x="missing", index="missing", 
                              xcond="numeric", scale="ANY"), 
          function(model, x, index, xcond, scale) {
              pdf <- noise_dist(model, "pdf")

              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              # browser()
              f <- function(x){
                      wrk <- (x - wrk)/model@scale
                      wrk <- pdf %of% wrk
                      wrk <- inner(wrk, model@prob/model@scale)
                      wrk
                  }
              f
          })

MixARGaussian <- setClass("MixARGaussian", 
         ## representation(), 
         ## prototype, 
         contains="MixAR"
         ## validity, access, where, version, sealed, package, 
         )

setMethod("show", "MixARGaussian", 
	  function(object) {
              callNextMethod()

              cat("\n")
              cat("Distributions of the error components:\n")
              cat("\tstandard Gaussian\n")
              cat("\n")

              invisible(object)
	  })

## cdf
setMethod("mix_cdf", signature(model="MixARGaussian", x="numeric", index="numeric", 
                               xcond="missing", scale="ANY"), # note: scale is ignored here!
          function(model, x, index, xcond, scale) {
              wrk <- mix_ek(model, x, index, scale=TRUE)
              wrk <- pnorm %of% wrk
              wrk <- inner(wrk, model@prob)
              wrk
          })


setMethod("mix_cdf", signature(model="MixARGaussian", x="numeric", index="missing", xcond="numeric", 
                              scale="ANY"), 
          function(model, x, index, xcond, scale) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              # browser()
              wrk <- (x - wrk)/model@scale
              wrk <- pnorm %of% wrk
              wrk <- inner(wrk, model@prob)
              wrk
          })

setMethod("mix_cdf", signature(model="MixARGaussian", x="missing", index="missing", xcond="numeric", 
                              scale="ANY"), 
          function(model, x, index, xcond, scale) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              f <- function(x){
                      wrk <- (x - wrk)/model@scale
                      wrk <- pnorm %of% wrk
                      wrk <- inner(wrk, model@prob)
                      wrk
                  }
              f
          })

## pdf
setMethod("mix_pdf", signature(model="MixARGaussian", x="numeric", index="numeric", xcond="missing", 
                              scale="ANY"),  # note: scale is ignored here!
          function(model, x, index, xcond, scale) {
              wrk <- mix_ek(model, x, index, scale=TRUE)
              wrk <- dnorm %of% wrk
              wrk <- inner(wrk, model@prob/model@scale)
              wrk
          })


setMethod("mix_pdf", signature(model="MixARGaussian", x="numeric", index="missing", xcond="numeric", 
                              scale="ANY"), 
          function(model, x, index, xcond, scale) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
              # browser()
              wrk <- (x - wrk)/model@scale
              wrk <- dnorm %of% wrk
              wrk <- inner(wrk, model@prob/model@scale)
              wrk
          })


setMethod("mix_pdf", signature(model="MixARGaussian", x="missing", index="missing", 
                              xcond="numeric", scale="ANY"), 
          function(model, x, index, xcond, scale) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1)
#browser()
              f <- function(x){
                      wrk <- (x - wrk)/model@scale
                      wrk <- dnorm %of% wrk
                      wrk <- inner(wrk, model@prob/model@scale)
                      wrk
                  }
              f
          })

## quantile function
setMethod("mix_qf", signature(model = "MixARGaussian", p = "numeric", x = "numeric",
                              index = "numeric", xcond = "missing"),
          function(model, p, x, index, xcond) { 
              ## index <- seq_along(x)[-(1:max(model@order))]
              ##
              ## TODO: verify that this is consistent with mix_pdf, mix_location, etc.
              ##       Is this  conditional quantile computed 'at time i' or 'for time i'?
              ##       (seems it should be 'for time i', i.e. at time i-1, 
              ##                                                     see ?mix_pdf-methods)
	      ## ind <- max(model@order) : 1

              res <- sapply(index, function(i) mix_qf(model, p, xcond = x[1:i]))
              res
          })

setMethod("mix_qf", signature(model = "MixARGaussian", p = "numeric", x = "missing", 
                              index = "missing", xcond = "numeric"), 
          function(model, p, x, index, xcond) {
              qf <- mix_qf(model, xcond = xcond)
              qf(p)
          })

setMethod("mix_qf", signature(model = "MixARGaussian", p = "missing", x = "missing",
                              index = "missing", xcond = "numeric"), 
          function(model, p, x, index, xcond) {
              cdf <- mix_cdf(model, xcond = xcond)

              f <- function(p, ..., tol = .Machine$double.eps^0.5){
                  if(length(p) > 1)
                      # sapply(p, cdf2quantile, tol = tol, ..., MoreArgs = list(cdf = cdf))
                      sapply(p, cdf2quantile, cdf = cdf, tol = tol, ...)
                  else
                      cdf2quantile(p, cdf, tol = tol, ...)
              }

              f
          })

## todo: arg "1-B"
fit_mixAR <- function(x, model, init, fix, ...){
    stop("There is currently no default for this function")
}

setGeneric("fit_mixAR")

setMethod("fit_mixAR", signature(x = "ANY", model = "numeric", init = "missing"), 
          function(x, model, init, fix, ...){
              fit_mixAR(x, model, init = 1, fix = fix, ...)
          })

setMethod("fit_mixAR", signature(x = "ANY", model = "numeric", init = "numeric"), 
          function(x, model, init, fix, ...){
              model <- new("MixARGaussian", order = model)
              fit_mixAR(x, model, init = init, fix = fix, ...)
          })

setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "list"), 
          function(x, model, init, fix, ...){
              wrk <- list()
              for(i in init){
                  wrk <- c(wrk, list( fit_mixAR(x, model, init = i, fix = fix, ...) ))
              }
              wrk                     # todo: process the result, e.g. to select a best model?
          })

setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "numeric"), 
          function(x, model, init, fix, permute, select = TRUE, ...){
              if(missing(permute))
                 permute <- FALSE

              wrk <- list()
              est_shift <- missing(fix) ||
                          ! identical(fix, "shift")  # todo: tova e krapka, za da tragnat
                                                    # nestata.
              templ <- est_templ(model, shift = est_shift)
              for(i in 1:init){      # generate random mixAR models # todo: make 1:init safer!
                  ri <- em_rinit(as.numeric(x), model@order, templ)
                  wrk <- c(wrk, list( fit_mixAR(x, model, init = ri, fix = fix, ...) ))
                      # print(ri)
                      # print(wrk[[i]])
                      # browser()
              }
              so <- order(sapply(wrk, function(x) x$vallogf), decreasing = TRUE)
              wrk <- wrk[so]

                          # todo: for testing, reorganise later
              if(permute){      # todo: generiram otnovo vsichki permutatsii, ponezhe
                      # izglezhda, che inache izpuskam nesto. Obmisli i proveri pak!  =>
                      # Otkrich kakvo stava, neveroyatno! Dvoyna greshka!  permn(model@order)
                      # permutira elementite na order, napr. c(2, 2, 1) v rezultat
                      # mixAR_permute ne raboti kakto ochakvam, i ponezhe ne proveryava
                      # permutatsiite tichomalkom proizvezhda nevaliden model ('arcoef' ne
                      # saotvetsva na 'order'). Obache tova se kompensira po stechenie na
                      # obstoyatelstvata ot vtorata greshka, koyato se satoi v tova, che
                      # vmesto 'mixAR_permute' tryabva da se vika druga funktsiya (sega ste
                      # ya napisha). V rezultat izglezhda, che vse pak nestata se permutirat
                      # i rabotyat pone odnyakade kakto se ochakva kogato generiram vsichki
                      # permutatsii kakto po-dolu.
                  allperm <- permn(.nmix(model)) # todo: generate only unique ones.
                  # allperm <- unique( permn(model@order) )    # take only non-redundant perm.
                  cnt <- 0
                  f <- function(fm){
                      mo <- fm$model
                      locfits <- lapply(allperm, 
                                        function(p){
                                                               # locmo <- mixAR_permute(mo, p)
                                            locmo <- mixAR_switch(mo, p)
                                            fit_mixAR(x, model, init = locmo, fix = fix, ...)
                                        })
                      cnt <<- cnt + 1
                      ## 2020-06-12 was: cat("perm counter: ", cnt, "\n")
                      ##                 print( which.max( sapply(locfits, function(z) z$vallogf) ) )

                      # browser()   commented out on 2012-09-21

                      locfits[[ which.max( sapply(locfits, function(z) z$vallogf) )[1] ]]
                  }
                  newwrk <- lapply(wrk, f)
                                        # browser()
                  wrk <- newwrk
              }

                     # todo: needs more work; process the result, e.g. to select a best model;
              if(isTRUE(select)){
                                          # todo: more is needed here; the alg may stop before
                                          #       vallogf becomes NaN hence, krapka, for now.
                                          # todo: arbitrary constant below!
                  loli <- sapply(wrk, function(z) if(is.nan(z$vallogf)) -Inf else z$vallogf)
                  ololi <- order(loli, decreasing = TRUE)

                  goodsigma <- sapply(wrk[ololi], function(z) all(z$model@scale > 1e-7))

                  if(any(goodsigma))                      # todo: awful code!
                      res <- wrk[ololi][goodsigma][[1]]
                  else
                      res <- wrk[[ which.max(loli)[1] ]]
              }else{
                  res <- wrk
              }

              res
          })

## 12/09/2018 The following method now handles "raggedCoefS" objects
## note: previous content has not been changed, only added "if" conditions
setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "missing"), 
          function(x, model, init, fix, ...){
              if(inherits(model@arcoef, "raggedCoefS")){
                      # if(missing(fix)) {
                      #     mixSARfit(x, model)
                      # }else{
                      #     mixSARfixfit(x, model)
                      # }
                  fix <- if(missing(fix))
                             FALSE
                         else if(is.logical(fix))
                             fix
                         ## 2020-06-05 amended by Georgi
                         ##     you can't just use a completely different syntax for 
                         ##     an argument that works for other methods.
                         ##     An example in inst/slowtests/example.R actually uses 
                         ##     fix = "shift" and throws error (before this change).
                         ## TODO: error message below is still confusing since the user
                         ##     may not notice that it is about mixSARfit, not mixARfit
                         ##     and even if they did, that would be equally confusing.
                         else if(is.character(fix) && length(fix) == 1 && fix == "shift")
                             TRUE
                         else
                             stop("'fix' must be boolean or missing for 'mixSARfit()'")
                  mixSARfit(x, model, est_shift = fix)
              }else{# todo: process "..." for est_shift?
                  fit_mixAR(x, model, init = model, fix = fix, ...)
              }# todo: do more here?
          })

setMethod("fit_mixAR", signature(x = "ANY", model = "MixAR", init = "MixAR"), 
          function(x, model, init, fix, ...){             # todo: process "..." for est_shift?
              ## 2012-11-02 premestvam processing na "fix" v mixARemFixedPoint
              ##    todo: da machna 'fix' ot spisaka na argumentite, za da se predava s '...'.
              ## est_shift <- missing(fix) ||
              ##             ! identical(fix, "shift")  # todo: tova e krapka, za da tragnat
                                                    # nestata.

              ## todo: in this case model is a template for the result. Currently it is not
              ##       used at all in this case but eventually should do some copying and
              ##       other housekeeping activities. In particular, 'init' and 'model' may be
              ##       from different subclasses of 'mixAR'.
              if(missing(fix))
                  mixARgenemFixedPoint(x, init, ...)
              else
                  mixARgenemFixedPoint(x, init, fix = fix, ...)
          })

setMethod("fit_mixAR", signature(x = "ANY", model = "MixARGaussian", init = "MixAR"), 
          function(x, model, init, fix, ...){             # todo: process "..." for est_shift?
              est_shift <- missing(fix) || ! identical(fix, "shift")  # tova beshe krapka, za
                                                                     # da tragnat nestata.

              ## todo: in this case model is a template for the result. Currently it is not
              ##       used at all in this case but eventually should do some copying and
              ##       other housekeeping activities. In particular, 'init' and 'model' may be
              ##       from different subclasses of 'mixAR'.
              mixARemFixedPoint(x, init, est_shift = est_shift, ...)
          })

# build_mixAR <- function(x, ncomp = 2, pmax = 2, ...){
#                                      # stop("There is currently no default for this funciton")
#     if(ncomp < 2)
#         ncomp <- 2
#
#     for(g in 2:ncomp){
#
#
#     }
# }

# setGeneric("build_mixAR")
#
# setMethod("build_mixAR", signature(x = "ANY", model = "numeric", init = "missing"), 
#           function(x, model, init, fix, ...){
#               build_mixAR(x, model, init = 1, fix = fix, ...)
#           })

setGeneric("mix_location",                         # 2012-11-08 new
          function(model, x, index, xcond)
              standardGeneric("mix_location"), 
           useAsDefault = FALSE)

setMethod("mix_location", 
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"), 
          function(model, x, index, xcond) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
              wrk <- inner(wrk, model@prob)
              wrk
          })

setMethod("mix_location", 
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "missing"), 
          function(model, x, index, xcond) {
              model <- model
              f <- function(xcond){
                  wrk <- mix_hatk(model, xcond, length(xcond)+1)
                  wrk <- inner(wrk, model@prob)
                  wrk
                  }
              f
          })

setMethod("mix_location", 
          signature(model = "MixAR", x = "numeric", index = "numeric", xcond = "missing"), 
          function(model, x, index, xcond) {
                                           # mix_pdf and others use mix_hatk for this purpose, 
                                           # mixFilter is cleaner (and the name is better)
              wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
              inner(wrk, model@prob)
          })

setMethod("mix_location", 
          signature(model = "MixAR", x = "numeric", index = "missing", xcond = "missing"), 
          function(model, x, index, xcond){
              index <- seq_along(x)[-(1:max(model@order))]
                                           # mix_pdf and others use mix_hatk for this purpose, 
                                           # mixFilter is cleaner (and the name is better)
              wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
              inner(wrk, model@prob)
          })


                                        # todo: needs to check existence of variance
setGeneric("mix_variance",                         # 2012-11-08 new
          function(model, x, index, xcond)
              standardGeneric("mix_variance"), 
           useAsDefault = FALSE)

setMethod("mix_variance", 
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"), 
          function(model, x, index, xcond) {
              wrk <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
              loc <- inner(wrk, model@prob)
                                     # wrk <- inner(wrk, model@prob, function(x, y) x^2 * y   )
              wrk <- drop( (wrk@m^2 + model@scale^2) %*% model@prob )
              wrk - loc^2
          })

setMethod("mix_variance", 
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "missing"), 
          function(model, x, index, xcond) {
              model <- model
              f <- function(xcond){
                  wrk <- mix_hatk(model, xcond, length(xcond)+1)
                  loc <- inner(wrk, model@prob)
                                     # smenyam, ponezhe "inner" ne e dostatachno gavkava
                                     #    wrk <- inner(wrk, model@prob, function(x, y) x^2 * y)
                  wrk <- drop(wrk@m^2 %*% model@prob)
                  wrk - loc^2
              }
              f
          })

setMethod("mix_variance", 
          signature(model = "MixAR", x = "numeric", index = "numeric", xcond = "missing"), 
          function(model, x, index, xcond) {
                                           # mix_pdf and others use mix_hatk for this purpose, 
                                           # mixFilter is cleaner (and the name is better)
              wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
              loc <- inner(wrk, model@prob)
                                     # wrk <- inner(wrk, model@prob, function(x, y) x^2 * y   )
              wrk <- drop(wrk@m^2 %*% model@prob)
              wrk - loc^2
          })

# setMethod("mix_variance", 
#           signature(model = "MixAR", x = "numeric", index = "missing", xcond = "missing"), 
#           function(model, x, index, xcond){
#               index <- seq_along(x)[-(1:max(model@order))]
#                                            # mix_pdf and others use mix_hatk for this purpose, 
#                                            # mixFilter is cleaner (and the name is better)
#               wrk <- mixFilter(x, model@arcoef, index, shift = model@shift)
#               inner(wrk, model@prob)
#           })




setGeneric("mix_moment",                         # 2012-11-08 new
          function(model, x, index, xcond, k)
              standardGeneric("mix_moment"), 
           useAsDefault = FALSE)

setMethod("mix_moment", 
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"), 
          function(model, x, index, xcond, k) {
              locall <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
              loc <- inner(locall, model@prob)

              gmom <- cbind(1, 0, model@scale^2, 0, 3*model@scale^4, 0) # todo: this is temp!

              facbin <- choose(k, k:0)
              facloc <- outer(drop(locall@m), k:0, "^")
              faceps <- gmom[ , 0:k + 1]

              # tuk sa obiknoveni matristi
              # sum( inner(facloc * faceps * model@prob, facbin) )
              sum((facloc * faceps * model@prob) %*% facbin)
          })


setGeneric("mix_central_moment",                         # 2012-11-08 new
          function(model, x, index, xcond, k)
              standardGeneric("mix_central_moment"), 
           useAsDefault = FALSE)

setMethod("mix_central_moment",                                # note: k is scalar (currently)
          signature(model = "MixAR", x = "missing", index = "missing", xcond = "numeric"), 
          function(model, x, index, xcond, k) {
              locall <- mix_hatk(model, xcond, length(xcond)+1) # todo: make simpler!
              loc <- inner(locall, model@prob)

                                   # gmom <- cbind(1, 0, model@scale^2, 0, 3*model@scale^4, 0)
              gmom <- outer(model@scale, 0:k, "^") * noise_moment(model, 0:k)

              facbin <- choose(k, k:0)
              facloc <- outer(drop(locall@m) - loc, k:0, "^")
              faceps <- gmom[ , 0:k + 1]

              # tuk sa obiknoveni matristi
              # sum( inner(facloc * faceps * model@prob, facbin) )
              sum( (facloc * faceps * model@prob) %*% facbin )
          })

mix_kurtosis <- function(...){
    mix_central_moment(..., k = 4) / mix_central_moment(..., k = 2)^2
}

mix_ekurtosis <- function(...){
    mix_kurtosis(...) - 3
}



################

noise_moment <- function(model, k){ #elements of k are >= 0; odd moments are zero, var is 1, 
    if(max(k) <= 3)
        c(1, 0, 1, 0)[k+1]
    else
        c(1, 0, 1, 0, rep(c(NA, 0), length.out = max(k-3)))[k+1]
}

setGeneric("noise_moment")

setMethod("noise_moment", signature(model = "MixARGaussian", k = "numeric"), 
          function(model, k){
              wrk <-
              if(max(k) <= 5)
                  c(1, 0, 1, 0, 3, 0)[k+1]
              else
                  stdnormmoment(k)
              res <- rep(wrk, .nmix(model))
              if(length(k) > 1)       # one row for each component
                  matrix(res, byrow = TRUE, nrow = .nmix(model))
              else
                  res
          })


setMethod("noise_moment", signature(model = "MixARgen", k = "numeric"), 
          function(model, k){
              dist <- get_edist(model)
              ## 2021-04-30 :TODO: 
              ##    the following throw error since there is no component 'moment':
              ## noise_moment(exampleModels$WL_At, 3)      # gives error
	      ## noise_moment(exampleModels$WL_ibm_gen, 2) # gives error
              wrk <- sapply(dist, function(x) x$moment(k))

              if(length(k) > 1)       # one row for each component
                  t(wrk)              # sapply returns one column for each, so transpose
              else
                  wrk
          })

setGeneric("noise_rand", 
           function(model, expand = FALSE)
              standardGeneric("noise_rand"), 
           useAsDefault = FALSE)

setMethod("noise_rand", signature(model = "MixAR"), 
          function(model, expand = FALSE){
              mess <- paste("No method for argument of class", class(model))
              stop(mess)
          })

setMethod("noise_rand", signature(model = "MixARGaussian"), 
          function(model, expand = FALSE){
              if(expand) rep(list("rnorm"), mix_ncomp(model))
              else       list("rnorm")
          })

setMethod("noise_rand", signature(model = "MixARgen"), 
          function(model, expand = FALSE){
              dist <- get_edist(model)
              res <- lapply(dist, function(x) x[["rand"]])

              if(expand) rep(res, mix_ncomp(model))
              else       res
          })



setGeneric("noise_dist", 
           function(model, what, expand = FALSE)
              standardGeneric("noise_dist"), 
           useAsDefault = FALSE)


setMethod("noise_dist", signature(model = "MixAR"), 
          function(model, what, expand = FALSE){
              mess <- paste("No method for argument of class", class(model))
              stop(mess)
          })

setMethod("noise_dist", signature(model = "MixARGaussian"), 
          function(model, what, expand = FALSE){
              res <- dist_norm[[what]]

              if(expand) rep(list(res), mix_ncomp(model))
              else       list(res)
          })



######################### new: 2011-11-16 caters for the new format of slot 'dist' in mixARgen
setGeneric("get_edist", 
           function(model)
              standardGeneric("get_edist"), 
           useAsDefault = FALSE)


setMethod("get_edist", signature(model = "MixAR"), 
          function(model){
              mess <- paste("No method for argument of class", class(model))
              stop(mess)
          })

setMethod("get_edist", signature(model = "MixARGaussian"), 
          function(model){
              list(dist_norm)
          })
#########################

setGeneric("noise_params", 
           function(model)
              standardGeneric("noise_params"), 
           useAsDefault = FALSE)


setMethod("noise_params", signature(model = "MixAR"), 
          function(model){
              list()    # the default is no parameters (but todo: should this be numeric(0)?
          })

setMethod("noise_params", signature(model = "MixARgen"),       # todo:  needs further thought !!
          function(model){
              if(!is.null(model@dist$param)){     # todo: krapka, to avoid breaking old code
                  ## for old code
                  model@dist$param  
              }else{
                  dist <- get_edist(model)
                  sapply(dist, function(x) x$get_param())
              }
          })



set_noise_params <- function(model, nu){
    if(!is.null(model@dist$param)){     # todo: krapka, to avoid breaking old code
        model@dist$param <- nu
    }else{
        for(i in seq_along(model@dist))
            model@dist[[i]]$set_param( nu[i] )  # very simplistic, assumes one param per dist.
    }
    model
}

setMethod("noise_dist", signature(model = "MixARgen"), 
          function(model, what, expand = FALSE){
              wrk <- get_edist(model)
              res <- lapply(wrk, function(x) x[[what]])

              # if(is.null(res))
              #     stop(paste("Property", what, "has not been specified for this model."))

              if(expand && length(res) == 1) rep(res, mix_ncomp(model))
              else       res
          })

setMethod("get_edist", signature(model = "MixARgen"), 
          function(model){ # todo: the else part of 'if' is a patch, to avoid breaking old code
              res <- if(!is.null(model@dist$generator)){
                         par <- model@dist$param
                         if(is.null(par) || length(par) == 0)
                             model@dist$generator()
                         else
                             model@dist$generator(par)
                     }else{ # tova li e za savmestimost? todo: Utochni!
                         model@dist
                     }
              res
          })

setGeneric("show_diff",
           function(model1, model2)
              standardGeneric("show_diff"),
           useAsDefault=FALSE)

setMethod("show_diff", signature(model1="MixAR", model2="MixAR"),
          function(model1, model2){
              if(!identical(model1@order, model2@order)){
                  cat("The two models need to have the same number of components\n",
                      "  and the same AR orders.")
                  return("")
              }

              g <- length(model1@prob)
              p <- max(model1@order)
              f <- function(i){ wrkloc <- c(  model1@prob[i]  - model2@prob[i]
                                            , model1@shift[i] - model2@shift[i]
                                            , model1@scale[i] - model2@scale[i]
                                            , model1@order[i])
                                if(p > 0) c(wrkloc, model1@arcoef[[i]] - model2@arcoef[[i]])
                                else      wrkloc
                              }
              wrk <- lapply(1:g, f)

              mcoef <- ragged2char(wrk)

              colnames(mcoef) <- c("prob", "shift", "scale", "order",
                                   if(p > 0) paste("ar_", seq_len(p), sep="")
                                   else character(0))
              rownames(mcoef) <- paste("Comp_", 1:g, sep="")

              print(mcoef, na.print="", quote=FALSE, justify="right")


          })

                                     # todo: if all subclasses of mixAR can return a vector of
                                     # descriptions will not need so many additional methods.
setMethod("show_diff", signature(model1="MixARgen", model2="MixARgen"),
          function(model1, model2){
              callNextMethod()

              cat("\n")                                # todo: better formatting (in a table?)
              cat("Distributions of the error components:\n")

              cat("\n\tModel 1:\n"); b_show(get_edist(model1))
              cat("\n\tModel 2:\n"); b_show(get_edist(model2))

              invisible("")
	  })

setMethod("show_diff", signature(model1="MixARgen", model2="MixARGaussian"),
          function(model1, model2){
              callNextMethod()

              cat("\n")                                # todo: better formatting (in a table?)
              cat("Distributions of the error components:\n")

              cat("\n\tModel 1\n"); b_show(get_edist(model1))
              cat("\n\tModel 2\n"); cat("\t All components: Standard normal distribution\n")

              invisible("")
	  })

setMethod("show_diff", signature(model1="MixARGaussian", model2="MixARgen"),
          function(model1, model2){
              callNextMethod()

              cat("\n")                                # todo: better formatting (in a table?)
              cat("Distributions of the error components:\n")

              cat("\n\tModel 1\n"); cat("\t All components: Standard normal distribution\n")
              cat("\n\tModel 2\n"); b_show(get_edist(model2))

              invisible("")
	  })

                                                                  # 2012-11-02 new arg. 'drop'
parameters <- function(model, namesflag = FALSE, drop = character(0)){
    coef(model)
}
setGeneric("parameters")

setMethod("parameters", "MixAR",
          function(model, namesflag = FALSE, drop = character(0)){  # todo: return named list!
              nams <- c("order", "prob", "shift", "scale")
              nams <- nams[!(nams %in% drop)]  # todo: process arcoef similarly?

              res <- unlist(c(lapply(nams, function(x) slot(model, x)),
                              model@arcoef@a
                              ))
              if(namesflag)
                  names(res) <- c(sapply(nams,
                                         function(x) paste(x, 1:.nmix(model), sep="")),
                                  paste("ar_",
                                        rep(1:.nmix(model), times = model@order),
                                        unlist(lapply(model@order, function(x) seq_len(x))),
                                        sep=""))
              res
          })

"parameters<-" <- function(model, value){
    stop("'parameters<-' has no default.")
}
setGeneric("parameters<-")

setMethod("parameters<-", "MixAR",
          function(model, value){                                    # todo: return named list!
              nams <- c("order", "prob", "shift", "scale") # todo: put this in a variable?
              g <- .nmix(model)
              p <- model@order       # constant order

              # model@order <- value[1:g]     # krapka, za da raboti sas simulatsiite!
              model@prob  <- value[(g+1):(2*g)]
              model@shift <- value[(2*g+1):(3*g)]
              model@scale <- value[(3*g+1):(4*g)]
              model@arcoef <- rag_modify(model@arcoef, value[-(1:(4*g))])

              # browser()

              model
          })

# 2020-4-20
## companion_matrix and isStable now handle mixVAR objects.
companion_matrix <- function(v, ncol = length(v), nrow = ncol){
    mult <- ifelse(is.matrix(v), nrow(v), 1)
    if(ncol %% mult != 0) stop("Wrong dimensions for companion matrix")
    if(mult > 1){
        if(ncol(v) < ncol) v <- cbind(v, matrix(0, ncol = (ncol - ncol(v)), nrow = mult)) 
    }else{
        if(length(v) < ncol)
            v <- c(v, rep(0, ncol - length(v)))
    }
    rbind(v, diag(1, nrow = nrow-mult , ncol = ncol), deparse.level = 0 ) # no labels
}


isStable <- function(x){     # todo: make generic?
    cl <-inherits(x, "MixVAR")
    nc <- max(x@order)
    co <- if(cl) x@arcoef else x@arcoef[]    # a matrix
    prob <- x@prob
    if(ncol(co[]) == 0)   # i.i.d. mixture
        return(TRUE)

    f <- function(k){
        if(cl){ m <- companion_matrix(co[k, ], nc * nrow(co[]))
        }else m <- companion_matrix(co[k,], nc)
        prob[k] * kronecker(m,m)
    }

    # wrk <- lapply(seq_along(x@prob), f)
    wrk <- do.call(.mplus, lapply(seq_along(x@prob), f))

    abs(eigen(wrk)$values[1]) < 1                       # stable if  maximal eigenvalue is < 1
}

mixAR_permute <- function(model, perm){   # todo: make this generic; this is only a start.
    if(all(perm == seq_along(perm)))
        return(model)

    lapply(c("prob", "order", "shift", "scale"),
           function(nam) slot(model, nam) <<- slot(model, nam)[perm] )

    model@arcoef@a <- model@arcoef@a[perm]
    model@arcoef@p <- model@arcoef@p[perm]

    # print(model)

    model
}

mixAR_switch <- function(model, perm){  # todo: check the permutation?
    if(all(perm == seq_along(perm)))
        return(model)

    lapply(c("prob", "shift", "scale"),                        # note: "order" is omitted here
           function(nam) slot(model, nam) <<- slot(model, nam)[perm] )

    p <- model@order
    m <- model@arcoef[perm]
    for(i in 1:length(p)){
        if(p[i]>0)
            model@arcoef@a[[i]] <- m[i,1:p[i]]
    }
    model
}