R/mcompanion.R

In mcompanion: Objects and Methods for Multi-Companion Matrices

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

                                                             # 2014-06-11 substantially edited
mc_chain_scale <- function(ev, subset = NULL, fvec = NULL, fchain = NULL){
    if(is.null(subset)){
        n <- nrow(ev$eigvec)
        subset <- (n - ev$mo + 1) : n
    }else if(is.character(subset) && subset == "top")
        subset <- 1 : ev$mo
    ## else subset must be a vector of integers suitable for indexing

    if(is.null(fvec)){
        fvec <- function(v, ind){
            i <- which.max(abs( v[ind] ))
            v[ind][i]   # todo: check for positive length of ind? (in case v is all NA's)?
        }
    }

    if(is.null(fchain))
        fchain <- function(chain) chain / fvec(chain[ , 1], subset)

    new.chains <- lapply(chains_to_list(ev$eigvec, ev$len.block), fchain)
    ev$eigvec <- do.call("cbind", new.chains)
    ev
}

mc_chain_to_list <- function(ev){
    chains_to_list(ev$eigvec, ev$len.block)
}

mc_chain_subset <- function(ev, chainno){
    indx <- if(all(chainno > 0))
                chainno
            else
                (1:length(ev$len.block))[chainno]   # need positive indx for chain_ind()

    if(length(indx)==0)    # empty subset
        return( list() )

    len.block <- ev$len.block[indx]
    evindx <- chain_ind(indx, ev$len.block)
    eigvec    <- ev$eigvec[ , evindx, drop = FALSE]
    co <- if(is.null(ev$co))
              NULL
          else
              ev$co[, evindx, drop = FALSE]  # 2015-12-02 was: ev$co[, evindx]

    res <- list( mo        = ev$mo
               , mo.col    = ev$mo.col
               , eigval    = ev$eigval[indx]
               , len.block = len.block
               , eigvec    = eigvec
               , co        = co
                )
    res
}

mc_chain_merge <- function(ev1, ev2){                     # todo:??? merge more than 2 chains?
    if(identical(ev2, list()))
        return(ev1)
    if(identical(ev1, list()))
        return(ev2)

    mo        <- ev1$mo          # must be equal to ev2$mo
    mo.col    <- ev1$mo.col      # must be equal to ev2$mo.col
    eigval    <- c(ev1$eigval, ev2$eigval)
    len.block <- c(ev1$len.block, ev2$len.block)
    eigvec    <- cbind(ev1$eigvec, ev2$eigvec)


    co <- if(!is.null(ev1$co) && !is.null(ev2$co))    # merge the co's if both non-NULL
              cbind(ev1$co, ev2$co)
          else if(is.null(ev1$co) && is.null(ev2$co)) # keep co NULL if both NULL
              NULL
          else if(nrow(eigvec) < mo)               # co is  not usable here, so set it to NULL
              NULL                              # 2015-11-11 TODO: shouldn't this be an error?
          else  # 2015-12-02 was: eigvec[(nrow(eigvec)-mo+1):nrow(eigvec)]
                #                 (surely the intention is to get the bottom mo rows!)
              eigvec[(nrow(eigvec) - mo + 1):nrow(eigvec), , drop = FALSE]
                                                    # set co to the bottom, is this ok?
                                                    # maybe do this only if asked explicitly?
                                                    # and signal error otherwise?
    res <- list(  mo        = mo
                , mo.col    = mo.col
                , eigval    = eigval
                , len.block = len.block
                , eigvec    = eigvec
                , co        = co
                )
    res
}

                            # 2015-12-26 removed arguments mo (redundant, use ev$mo) and F0bot
mc_chain_extend <- function(ev, newdim){  # function(ev, mo, newdim)
    indx0 <- which(ev$eigval == 0)       # need more careful test? even better maybe to
                                         # introduce an argument to suply a function for this.
    ev0    <- mc_chain_subset(ev, indx0)
    evnon0 <- mc_chain_subset(ev, if(length(indx0) > 0) -indx0 else 1:length(ev$eigval))

    chnon0 <- mC.non0chain.extend(evnon0, newdim)

    v0 <- chains_to_list(ev0$eigvec, ev0$len.block)
    ch0 <- mc_0chains(newdim, ev$mo, ev$mo.col, v0)

    mc_chain_merge(chnon0, ch0)
}

                                                                  # 2014-06-10 new arg. mo.col
                            # 2015-11-10 new arg. what.co for make_mcev() and make_mcgev(),
                            #            since '...'  is used also in the call of mc_0chains()
                          # 2015-11-10 dobavyam argument 'Mtop' i code za sluchaya mo.col < mo
                          # 2015-11-10 TODO: All this needs streamlining.
                          # 2015-12-25 removed arg. Mtop; use 'co' instead.
## make_mcchains() creates the chains using only make_mcev(), make_mcgev() and mc_0chains().
##    In particular, it doesn't call other low-level chain generating functions.
##    Also, if it calls  mc_0chains() with argument vec0, it always contains at least mo rows,
##      so argument F0bot is not passed on.
make_mcchains <- function(eigval, co, dim, len.block, eigval0 = FALSE, mo.col = NULL,
                          what.co = "bottom", ...){
    if(is(co, "SmallMultiCompanion")){# ignore args eigval, len.block and mo.col in this case
        Mch <- smc_chains(co)     # TODO: overwrite also what.co? (its value shouldn't matter)
        co <- Mch$co
        mo.col <- Mch$mo.col
        len.block <- Mch$len.block
        eigval <- Mch$eigval
    }else{
        if(is.vector(co))
            co <- as.matrix(co, ncol = 1)
        if(missing(len.block)){
            len.block <- if(length(eigval) == 1)
                             ncol(co)                # TODO: document this! Is this a feature?
                         else
                             rep(1, length(eigval))
        }
        if(length(eigval) != length(len.block))
            stop("Number of eigenvalues does not match number of chains.")
    }

    mo <- nrow(co)

    n.chain <- length(len.block)
    nc <- sum(len.block)
    if(nc > ncol(co))
        stop("Not enough seed coefficients.")
    else if(nc < ncol(co))
        warning("There are too many  seed coefficients, I ignore the excess.")

    ## generate the chains from the seeds using only make_mcev() and make_mcgev().
    ##     todo: for speed this chunk could be implemented in C/C++ (inlining make_xxx)
    xmat <- matrix(0, nrow = dim, ncol = nc)
    kcur <- 0
    for(i in 1:n.chain){
        kcur <- kcur + 1
        ev <- eigval[i]
                                          # 2015-11-10 remove '...' and use explicitly what.co
        xmat[ , kcur] <- make_mcev(ev, co[ , kcur], dim, what.co = what.co)
        if(len.block[i] > 1){
            for(j in 2:len.block[i]){
                v <- xmat[ , kcur]
                kcur <- kcur + 1
                                          # 2015-11-10 remove '...' and use explicitly what.co
                xmat[ , kcur] <- make_mcgev(ev, co[ , kcur], v, what.co = what.co)
            }
        }
    }

    if(is.null(mo.col))             # 2014-06-10 was unconditional
        mo.col <- ncol(xmat)        # 2015-11-11 was: if(missing(mo.col)) mo.col <- ncol(xmat)


    if(ncol(xmat) < dim  &&  eigval0){   # complete with structural  chains of zero eigenvalues
        eval0indx <- which(eigval == 0) # a crude test

        if(length(eval0indx) == 0){ # no zeroes in eigval
            if(mo.col < ncol(xmat))
                stop("Too many vectors (> mo.col) for non-zero eigenvalues.")
                                   # 2014-06-10 was: ... mc_0chains(dim, nrow(co), ncol(xmat))
            wrkvz <- mc_0chains(dim, mo, mo.col)
        }else{                      # zeroes in eigval
            if(sum(len.block[-eval0indx]) > mo.col)
                stop("Too many vectors (> mo.col) for non-zero eigenvalues.")

            vec0 <- vector("list", length(eval0indx))   # collect 0-chains in a list of chains
            for( i in seq(along = eval0indx)){ # 2015-11-10 was: for( i in eval0indx)
                indx0wrk <- chain_ind(eval0indx[i], len.block)
                vec0[[i]] <- xmat[ , indx0wrk, drop = FALSE]
            }
                 # 17/04/2007 was: flag0 <- check0chains(vec0,mo.col)
                 #                 if(!flag0)
                 #                   warning("echains for 0 eigenvalues are not good enough.")

            indx0all <- chain_ind(eval0indx, len.block)

            ## remove chains for 0 eval
            xmat <- xmat[ , -indx0all, drop = FALSE]   # 2015-11-11 xmat <- xmat[ , -indx0all]
            len.block <- len.block[-eval0indx]
            eigval <- eigval[-eval0indx]

                             # 2014-06-10 was: ... mc_0chains(dim, nrow(co), ncol(xmat), vec0)
                                                    # 2015-11-10 also pass "..." to mc_0chains
            wrkvz <- mc_0chains(dim, mo, mo.col, vec0, ...) # note arg. 'vec0' here!

        }
        ## merge the non0- and 0-chains
        xmat      <- cbind(xmat,  wrkvz$eigvec)
        eigval    <- c(eigval,    wrkvz$eigval)
        len.block <- c(len.block, wrkvz$len.block)
    }

    if(sum(len.block) != dim){     
        if(sum(len.block) > dim)     # todo: error, not warning?
            warning("The total length of the Jordan chains is greater than dim!")
        else if(eigval0) # here  sum(len.block) < dim
            warning("The total length of the Jordan chains is smaller than dim.")
    }

    list(eigval = eigval, len.block = len.block, mo = mo, eigvec = xmat, co = co,
         mo.col = mo.col )
}

                 # 18/04/2007 was: make_mcchains(eigval,co,dim,len.block,eigval0=TRUE,what.co)
                 # 2015-11-10 making eigval0 an argument to avoid error in the call to
                 #            make_mcchains() call when eigval0 is also in "..."

  # 2015-11-10 TODO: (..., type = "real", what.res = "matrix", eigval0)
  #                  need to check in pcts if and where this is called with unnamed arguments!
make_mcmatrix <- function(type = "real", what.res = "matrix", ..., eigval0){
    x    <- make_mcchains(eigval0 = TRUE, ...)
    jmat <- Jordan_matrix(x$eigval, x$len.block)
    res  <- from_Jordan(x$eigvec, jmat)

    if(type == "real")
        res <- Re(res)

    if(what.res == "list"){
        x[["mat"]] <- res      # if used together with type="real" ev should be in conj pairs.
        res <- x
    }

    res
}

                            # 2015-12-26 removing argument F0bot and disallowing nrow(ev) < mo
                            #            removing argument mo (it is redundant, use ev$mo)
mC.non0chain.extend <- function(ev, newdim){  # 04/05/2007 x0 => F0bot
    eigvalshort <- ev$eigval # 02/05/2007 ev$values
    eigvecshort <- ev$eigvec #            ev$vectors
    mo          <- ev$mo
    mo.col      <- ev$mo.col
    len.block   <- ev$len.block


    if(is.null(len.block))
        len.block <- rep(1, length(eigvalshort))

    if(is.null(mo.col))
        mo.col <- sum(len.block)

    stopifnot(newdim >= nrow(eigvecshort), # only  extend, not shrink; TODO: allow shrink?
              nrow(eigvecshort) >= mo
              )

    v <- matrix(NA, nrow = newdim, ncol = ncol(eigvecshort))
    v[1:nrow(eigvecshort), 1:ncol(eigvecshort)] <- eigvecshort

    mfr <- max(mo.col+1, mo+1)
    mfill <- if(mfr <= nrow(v))
                 mfr:nrow(v)
             else
                 numeric(0)

    ## # 2015-12-26  commenting this chunk out; see also comments for dropping arg. F0bot
    ##
    ## if(mo.col < mo){      # as.matrix is redundant but %*% may not be defined for its class
    ##     wrk <- # the else part should work right for the "if" as well (less efficiently)
    ##         if(max(len.block) == 1)   # only simple eigenvectors here
    ##             (as.matrix(F0bot) %*% eigvecshort) /
    ##                 matrix(rep(eigvalshort, each = nrow(F0bot)), nrow = nrow(F0bot))
    ##         else
    ##             (as.matrix(F0bot) %*% eigvecshort) %*%
    ##                 solve(Jordan_matrix(eigvalshort, len.block))
    ##                                                         # assumes nrow(F0bot)=mo-mo.col
    ##     v[nrow(eigvecshort)+ 1:nrow(F0bot), 1:ncol(eigvecshort)] <- wrk
    ## }

    ## TODO: if mfill is empty, this if/else chunk does nothing  - check and streamline.
    if( max(len.block)==1 ){   # only simple eigenvectors here
        for( j in mfill ){
            v[j,] <- v[j-mo,]/eigvalshort
        }
    }else{                     # some chains have length greater than 1 here.
        kcur <- 0
        for(i in 1:length(len.block)){
            kcur <- kcur + 1
            for( j in mfill ){                          # parvo zapalvam evector
                v[j,kcur] <- v[j-mo,kcur]/eigvalshort[i]
            }
            if(len.block[i]>1){         # ... then the remaining vectors  in the chain, if any
                for(k in 2:len.block[i]){
                    for( j in mfill ){  # parvo zapalvam evector
                        kcur <- kcur + 1
                        v[j,kcur] <- (v[j-mo,kcur] - v[j,kcur-1] )/eigvalshort[i]
                    }
                }
            }
        }
    }
                                                        # 2015-12-02 added  drop = FALSE below
    list(eigval = eigvalshort, len.block = len.block, mo = mo, eigvec = v,
         co = v[(nrow(v)-mo+1):nrow(v), , drop = FALSE], mo.col = mo.col )
}

## mc_0chains() is the main function here. It performs several steps which are implemented in
## separate functions for convenient development.

                                              # special Jordan chains for the zero eigenvalues
                                        # 2015-12-27 removing argument F0bot
mc_0chains <- function(dim, mo, mo.col, vec0, flagtriang = TRUE){ # , F0bot = NULL
    if(missing(vec0) || identical(vec0,list())){                     # only structural 0chains
        res <- mc.0chain.struct(dim, mo, mo.col)
    }else{                                                    # non-structural 0chains present
        nrvec0 <- nrow(vec0[[1]])
                                          # there is more to test here...  (TODO: this comment
                                          # is from ten years ago. Is it still relevant?)
                                          #
        ## 2015-12-27 commenting out since mc.0chain.dx() does not treat the case mo.col < mo
        ##            completely (for this, the non-0 eigenvectors may be needed).
        ##            argument F0bot is removed too, see above.
        ##
        ##                         # if (nrvec0 != mo.col) then the caller has set or
        ##                         #   calculated the e.v's for the mo x mo block, so we don't
        ##                         #   call mc.0chain.dx in this case even if mo.col < mo
        ## v0 <- if(mo.col < mo  &&  nrvec0 == mo.col)
        ##           mc.0chain.dx(mo, mo.col, chF0top = vec0, F0bot = F0bot)
        ##       else
        ##           vec0
        ##
        v0 <- vec0

        v0 <- if(flagtriang)
                  mc_chains_triangulate(v0, mo, mo.col)

        for(i in seq_along(v0))              # extend to dim the supplied chains, if necessary
            v0[[i]] <- mc.0chain.complete(dim, mo, v0[[i]]) # , F0bot = F0bot

        for( i in seq_along(v0))   # complete each chain with structural vectors, where needed
            v0[[i]] <- mc.0chain.structfill(mo, mo.col, v0[[i]])

                                                  # append whole structural 0chains
        res <- mc.0chain.struct(dim, mo, mo.col, v0)
    }

    res$eigvec <- do.call("cbind", res$chains)
    res
}

               # 2015-12-26 renamed from mc.0chain.triang() - not specific for eigenvalue zero
mc_chains_triangulate <- function(chains, mo, mo.col){
    m <- length(chains)
    stopifnot(m <= mo)                                     # for mc matrix we must have m<=mo
    if(m==1)
        return(chains)

    dim <- NROW(chains[[1]])
    len <- sapply(chains, NCOL)

    wrk <- chains[ sort(len, method = "sh", index = TRUE, decreasing = TRUE)$ix ]# stable sort
    for(i in 1:(m - 1)){
        evec <- wrk[[i]][ , 1]
        # evec[j] ==0 for j<=dim-mo, no za da izbyagna iznenadi vklyuchvam testa po-dolu. !!!
                                        #testing equality to zero, may be problematic
        ii <- match(TRUE, evec != 0 & 1:dim > dim - mo)
        for(j in (i + 1):m){                    # todo:  look for the largest abs() or similar
            a <- wrk[[j]][ii,1] / evec[ii]          # for stability?
            wrk[[j]] <- wrk[[j]] - a * wrk[[i]][,1:ncol(wrk[[j]])]
            wrk[[j]][ii,1] <- 0                      # to ensure comparison with zero works ok
        }
    }
    wrk
}

                                                     # 2015-12-27 removed arg. F0bot
mc.0chain.complete <- function(dim, mo, chain, alt0){# 2014-06-07 some clean up and
    wrk <- chain                                            # bug fixing
    nrold <- nrow(wrk)
    ncold <- ncol(wrk)
    alt0skip <- 0

    ## 2015-12-27 removing the code below for the case nrold < mo. Now require nrold >= mo.
    stopifnot(nrold >= mo)

    ## if(nrold < mo ){                  # dopalvam do mo x mo;    assume nrold == mo.col here
    ##     wrk <- rbind(wrk,
    ##                  cbind(if(ncold > 1)   # no 'else' here; cbind() ignores NULL
    ##                           F0bot %*% chain[,2:ncold],
    ##                        numeric(nrow(F0bot)) )
    ##                  )                          # nrow(F0bot) should be equal to mo - nrold
    ##     if(!missing(alt0)){
    ##         wrk[(nrold+1):mo,ncold] <- alt0[1:(mo-nrold)]
    ##         alt0skip <- 1 : (mo - nrold)  # 2014-06-07 was: mo-nrold; not sure about the
    ##     }                                 #  intention but the old value cannot be correct!
    ##     tmp <- F0bot %*% chain[,1]
    ##     if(any(tmp!=0))                # todo: a better check here?
    ##         wrk <- cbind( c(numeric(nrold),tmp), wrk )
    ##
    ##     nrold <- nrow(wrk)                  # updates nrold and ncold !!!
    ##     ncold <- ncol(wrk)                       # update alt0 as well?  ???
    ## }


    if(dim > nrold){     # extend
        nrdiff <- dim - nrold
        ncdiff <- ceiling(nrdiff/mo)     # 10/05/2007: nrold replaces dim below
                # 2014-06-07 was wrong:
                #          if( nrdiff %% mo  > 0 && all( wrk[nrold-mo + 1:ncdiff , 1] == 0 ) )
        if( nrdiff %% mo  > 0 && all( wrk[nrold-mo + 1:(nrdiff %% mo) , 1] == 0 ) )
            ncdiff <- ncdiff - 1   # shorter chain in this case

        res <- cbind(  matrix(0, nrow=dim, ncol=ncdiff)
                     , rbind(wrk,  matrix(0, nrow=nrdiff, ncol=ncold) ) )

        ncnew <- ncold+ncdiff   # = ncol(res)
        if(!missing(alt0))                            # alternative init for the free elements
            res[(nrold+1):dim,ncnew] <- alt0[-alt0skip]

        ## 2015-12-26 TODO: Tova e krapka sled vavezhdaneto na small multi-companion for
        ##                  mo.col<mo.  Tozi klon predi ne se e izpalnyaval veroyatno poradi
        ##                  strukturni 0 v smc.
        ##            Proveri i vzh dali mozhe struturni nuli pak da ne idvat tuk.
        if(ncnew > 1){
            wrkrows <- nrold+ 1:nrdiff
            for(i in (ncnew-1):(ncdiff+1) ){
                res[wrkrows,i] <- res[wrkrows - mo,i+1]
            }

            for(i in ncdiff:1 ){ # vsastnost tozi for mozhe da e ot (dim-1):1, praveyki
                                 # predishniya for izlishen. Tova e vazmozhno ponezhe ako
                                 # chain e naistina 0ev veriga tya udovletvoryava tezi
                                 # usloviya i za dadenite elementi.
                res[(mo+1):dim,i] <- res[1:(dim-mo),i+1]
            }
        }
    }else if(dim < nrold){   # shrink
        res <- wrk[1:dim, , drop = FALSE]
        while( all(res[ , 1] == 0) )     # drop leading columns of zeroes, if present.
            res <- res[ , -1, drop = FALSE]   # 2014-06-07 was: res[,-1]
    }else{ # dim == nrold
        res <- wrk
    }

    res
}

                                # v0 may be different for different gen.ev's, not implemented.
mc.0chain.structfill <- function(mo, mo.col, chain, v0 = rep(0,mo)){
    wrk <- cbind(chain)  # in case chain is a  vector
    if(nrow(wrk) <= mo)                  # dim <= mo:  if dim == mo, nothing to append;
        return(wrk)                      #             while dim < mo is probably an error.

    m <- ncol(wrk)  # m must be greater than 0,
    k <- 1:min(mo.col + mo, nrow(wrk))
    while(all(wrk[k, m] == 0) && any(wrk[-(1:mo), m] != 0)){
        wrk <- cbind(wrk, c(wrk[-(1:mo), m], v0))
        m <- m + 1
    }
    wrk
}

                                                       # append (structural) 0evecs, if needed
    # 04/05/2007 smenyam formata na rezultata, pravya go da e kato na mc.0chain.structObsolete
                                                                    # 2014-06-07 new arg. sort
mc.0chain.struct <- function(dim, mo, mo.col, chains = list(), sort = TRUE){
    m <- length(chains) # number of chains on input; if m == mo, then chains has the maximal
                        #     possible number of 0chains but we do not return as there may be
                        #     more 0vectors to be included.
    if(dim > mo.col){
        ## length(strupos) == min(mo, dim - mo.col), max number of structural 0chains
        if(m == 0){#patch; ensures largest blocks come first here, without changing other code
            strupos <- mo.col + (1:min(mo, dim - mo.col))
            wrkindx <- seq_along(strupos)   # 2014-05-31 was: strupos
        }else{
            strupos <- (max(mo.col, dim-mo) + 1) : dim

            moreindx <- numeric(0)
            for(i in 1:m){
                evec <- chains[[i]][ , 1]
                ix <- match(TRUE, evec != 0 ) # assumes triangulation; also. needs more care
                moreindx <- c(moreindx, ix)
            }

            wrk <- match(strupos, moreindx)
            wrk <- is.na(wrk)     # wrk[i] is TRUE  if strupos[i] is not found in moreindx.

            wrkindx <- which(wrk) # strupos[ wrkindx[i] ]  is not in moreindx.
        }

        ident <- diag(dim)
        moduli <- c(rep(-1, mo.col), ((mo.col+1):dim) %% mo)
        newchains <- lapply(wrkindx,
                            function(ind){
                                pos <- which(moduli  ==  strupos[ind] %% mo)
                                ident[ , rev(pos), drop = FALSE]
                            })
        chains <- c(chains, newchains)

        if(sort){            # 2014-06-07 sort the chains in descending order of their lengths
            indx <- order(sapply(chains, NCOL), decreasing = TRUE)
            chains <- chains[indx]
        }

    }else if(dim < mo.col){
        stop("dim must be >= mo.col")
    }#else dim == mo.col (nothing to add); # 2014-05-31 - return a list as in the other cases.

    eigval <- numeric(length(chains))
    len.block <- if(length(chains) > 0)        # 2014-10-31: inserted the 'if' clause
                     sapply(chains, NCOL)
                 else integer(0)   # sapply would give list() in this case

                                        # TODO: a check for the overall lengths of the chains?
                                        # TODO: insert mo and mo.col in the list?
    list(eigval = eigval, len.block = len.block, chains = chains)
}

                                             # parametric_gev_core
spec_core <- function(mo, evalue, heights,
                      ubasis = NULL,
                      uorth = NULL,
                      evspace = NULL
                      ){   # mo - mc-order

    heights <- sort(heights, decreasing = TRUE) # heights of chains, a.k.a. block lengths
    hmax <- heights[1]   # maximal height
    s <- length(heights)  # no. of chains corresponding to evalue same as number of e.vectors

    stopifnot(s <= mo) # s - dim of ev.space?
                                        # ni[i] is the number of chains with height at least i
                                        # in particular, ni[1] is the number of e.vec.
    ni <- sapply(1:hmax, function(x) sum(heights >= x) )

    talli <- c(ni[-1],0)     # talli[i] is the number of vectors in chain i that are part of
                             #          taller chains

    mi <- ni - talli  # mi[i] is the number of over-hanging vectors in chain i
                      #    (i.e. vectors that are the last in their chains)

    full.core.basis <- diag(nrow = mo) # todo: consider other spaces, e.g. Haan

                       # consider this fixed, usually ev.universe.basis is the standard basis.
                       #     these two are considered fixed
    ev.universe.basis <-
        if(is.null(ubasis))           # matrix(NA_real_, nrow = mo, ncol = s)
            full.core.basis
        else
            ubasis

                                                          # todo: is ev.uninverse.orth needed?
    ev.uninverse.orth <-        # matrix(NA_real_, nrow = mo, ncol = mo - s)
        if(is.null(uorth))
            Null(ev.universe.basis)
        else
            uorth

                                               # if ev.space coincides with ev.universe, then
                                               # ev.space.orth contains only the zero vector
    ev.space <-       # ncol(evspace)  == s?
        if(!is.null(evspace))
            evspace
        else if(s == ncol(ev.universe.basis)) # number of e.v. equal to dim of allowed space
            ev.universe.basis          # this covers the case s = mo
        else if(s <  ncol(ev.universe.basis))
            matrix(NA_real_, nrow = mo, ncol = s)
        else
            stop("Dimension of allowed e.v. space should be >= no. of e.vectors")


                      # this is orth. w.r.t. "full space", which maybe larger than ev.universe
    ev.space.orth <-
        if(anyNA(ev.space))  # todo: this could be more refined; use null Complement, i.e.
                             # replace with something like:
                             #     null_complement(ev.space, universe = full.core.basis)
            matrix(NA_real_, nrow = mo, ncol = ncol(full.core.basis) - s)
        else
            Null(ev.space)

    evso.dim <- ncol(ev.space.orth) # space for the cores of hanging vectors

    sp.tall <- vector(length = hmax, mode = "list")
    sp.hang <- vector(length = hmax, mode = "list")

    for(i in 1:hmax){
        sp.tall[[i]] <- if(talli[i] > 0)
                            matrix(NA_real_, nrow = mo, ncol = talli[i])
                        else
                            NA
        sp.hang[[i]] <- if(mi[i] > 0)
                            matrix(NA_real_, nrow = mo, ncol = mi[i])
                        else
                            NA
    }

    param.tall <- vector(length = hmax, mode = "list")
    param.hang <- vector(length = hmax, mode = "list")

    ## the cores of the hanging vectors are in ev.space.orth
    for(i in 1:hmax){
        ## for now the param tails are same as sp.tail
        param.tall[[i]] <- if(talli[i] > 0)
                               matrix(NA_real_, nrow = mo, ncol = talli[i])
                           else
                               NA
    }

    ev.space.locorth <- # todo: this whole if-else may be can be replaced with null_complement
        if(anyNA(ev.space)){  # todo: this could be more refined
            if(isNA(param.tall[[1]]))
                matrix(NA_real_, nrow = mo, ncol = ncol(ev.space))
            else
                matrix(NA_real_, nrow = mo, ncol = ncol(ev.space) - ncol(param.tall[[1]]))
        }else{
            if(isNA(param.tall[[1]]))
                ev.space
            else
                null_complement(param.tall[[1]], universe = ev.space)
        }


    for(i in 1:hmax){
        param.hang[[i]] <-
            if(mi[i] > 0){
                if(evso.dim == 0 && i >= 2) # if evso.dim = 0 these cores can be set to zero
                             # matrix(numeric(1), nrow = ncol(ev.space.locorth), ncol = mi[i])
                    matrix(numeric(1), nrow = 0, ncol = mi[i])
                else
                    matrix(NA_real_, nrow = ncol(ev.space.locorth), ncol = mi[i])
            }else
                NA

    }





    core.vectors <- vector(length = hmax, mode = "list")

    ## e.v.
    i <- 1
    if(!isNA(sp.tall[[i]])){
        v.tall <- sp.tall[[i]]
    }else
        v.tall <- NA

    if(!isNA(sp.hang[[i]])){
        if(nrow(param.hang[[i]]) == 0){
            if(ncol(ev.space.orth) == 0)
                v.hang <- ev.space.locorth # ????????? matrix( ev.space.orth
            else
                v.hang <- ev.space.orth
        }else
            v.hang <- sp.hang[[i]] %*% param.hang[[i]]
    }else
        v.hang <- NA


    ## krapka; TODO: make more refined!
    ## if the eigenvectors span the ev universe and all heights are the same,
    ##     the e.v. can be considered fixed
    if(all(heights == heights[1]) &&
       ncol(ev.space) == ncol(ev.universe.basis)){

        v.tall <- if(all(is.na(ev.space)))
                      ev.universe.basis
                  else
                      ev.space
    }


    v.both <-
        if(!isNA(v.tall)){
            if(!isNA(v.hang))
                cbind(v.tall, v.hang)
            else
                v.tall
        }else if(!isNA(v.hang))
            v.hang
        else # both are NA, should not happen
            stop("no vectors for height one!")


    generators <- vector(length = hmax, mode = "list")

    if(ncol(v.both) == ncol(ev.universe.basis) &&
       ncol(v.tall) < ncol(ev.universe.basis)
       ){ # e.v. span the space, so any of v.tall or
                                                 # v.hang determines the other, choose the
                                                 # smaller

        ## FOR TESTING ONLY: setting things to -Inf
        if(ncol(v.tall) <= ncol(v.hang)){ # v.tall are the parameters
            par <- "tall"
            v.hang[] <- Inf # for now, to allow xx.ss handle this case without much change.
            v.tall[] <- -Inf    # for now, see the comment above
                                # todo: derive from others?
        }else{
            par <- "hang"
            v.tall[] <- Inf    # for now, see the comment above
            v.hang[] <- -Inf    # for now, to allow xx.ss handle this case without too much
        }
        v.both <- cbind(v.tall, v.hang)

        generators[[1]] <- list(param = "tall", tall = v.tall, hang = v.hang,
                                universe = ev.universe.basis, method = "complement")

    }

    core.vectors[[i]] <- list(tall = v.tall, hang = v.hang, both = v.both)
                # i=1 here

    ## g.e.v. i >= 2
    for(i in (1:hmax)[-1]){
        if(!isNA(sp.tall[[i]])){ # TODO: why not in both cases? (to get logical NA?)
            v.tall <- sp.tall[[i]]
        }else
            v.tall <- NA

        if(!isNA(sp.hang[[i]])){
            if(nrow(param.hang[[i]]) == 0){
                if(i == 1) # ev  ; but TODO: i cannot be one here!
                    if(ncol(ev.space.orth) == 0)
                        v.hang <- ev.space.locorth # ????????? matrix( ev.space.orth
                    else
                        v.hang <- ev.space.orth
                else # gev
                    v.hang <- matrix(numeric(1), nrow = mo, ncol = ncol(param.hang[[i]]))
            }else
                ## 2014-11-24 was: v.hang <- sp.hang[[i]] %*% param.hang[[i]]
                ##  TODO: check  if the change is correct!
                v.hang <- ev.space.locorth %*% param.hang[[i]]
        }else
            v.hang <- NA

        v.both <-
            if(!isNA(v.tall)){
                if(!isNA(v.hang))
                    cbind(v.tall, v.hang)
                else
                    v.tall
            }else if(!isNA(v.hang))
                v.hang
            else # both are NA, should not happen
                stop("no vectors for height i!")

        core.vectors[[i]] <- list(tall = v.tall, hang = v.hang, both = v.both)
    }

    ## TODO: describe what is going on here.
    wrk <- vector(length = s, mode = "list")
    cur.col <- 1
    for(i in 1:s){
        li <- lapply(core.vectors[1:heights[i]],
                     function(x) x$both[ , i, drop = FALSE] )
        wrk[[i]] <- do.call("cbind", li)
    }

    co <- do.call("cbind", wrk)

    list(evalue = evalue, heights = heights,
         co = co, core.vectors = core.vectors,
         param.tall = param.tall,
         param.hang = param.hang,
         generators = generators
         )
}