R/proj.r

In dae: Functions Useful in the Design and ANOVA of Experiments

correct.degfree <- function(object)
{ 
  daeTolerance <- get("daeTolerance", envir=daeEnv)
  #check that degrees of freedom are correct
  check.df <- FALSE
  if (!is.projector(object))
    stop("Must supply a valid object of class projector")
  if (is.na(object@degfree))
    stop("Degrees of freedom are missing. Can use degfree<- to set.")
  e <- eigen(object@.Data, symmetric=T, only.values=T)
  nonzero.e <- e$values[e$values > daeTolerance[["eigen.tol"]]]
  dflen <- length(nonzero.e)
  if ((object@degfree - dflen) > daeTolerance[["eigen.tol"]])
    stop("Degrees of freedom of projector are not correct")   
  check.df <- TRUE
  check.df
}

validProjector <- function(object)
{ 
  daeTolerance <- get("daeTolerance", envir=daeEnv)
  #checks that a matrix is a square projection operator
  Q <- object@.Data
  isproj <- TRUE
  if (nrow(Q) != ncol(Q))
    isproj <- "Matrix is not square"
  else #if (!is.allzero(t(Q)-Q))
  {
    if (!isTRUE(all.equal(t(Q), Q, tolerance=daeTolerance[["element.tol"]], check.attributes = FALSE)))
      isproj <- "Matrix is not symmetric"
    else 
    { 
      #if (!is.allzero(Q%*%Q - Q))
      if (!isTRUE(all.equal(Q%*%Q, Q, tolerance=daeTolerance[["element.tol"]], check.attributes = FALSE)))
        isproj <- "Matrix is not idempotent"
    }
  }
  isproj
}

projector <- function(Q)
{ 
  daeTolerance <- get("daeTolerance", envir=daeEnv)
  p <- new("projector", .Data=Q)
  e <- eigen(Q, symmetric=T, only.values=T)
  nonzero.e <- e$values[abs(1 - e$values) < 0.9]
  dflen <- length(nonzero.e)
  p@degfree=dflen
  return(p)
}

setClass("projector", representation("matrix", degfree = "integer"), prototype(degfree=as.integer(NA)))
setAs(from="projector", to="matrix", 
      def=function(from){m <- from@.Data; m})
setValidity("projector", validProjector)

is.projector <- function(object)
{ 
  inherits(object, "projector") & validObject(object)
}

degfree <- function(object)
{ 
  if (!inherits(object, "projector"))
    stop("Must supply an object of class projector")
  object@degfree
}

"degfree<-" <- function(object, value)
  #A function to replace supplied or computed the degrees of freedom of the projector object
{ 
  if (!is.projector(object))
    stop("Must assign to a valid object of class projector")
  if (length(value) == 1)
    object@degfree <- as.integer(value)
  else
  { 
    daeTolerance <- get("daeTolerance", envir=daeEnv)
    e <- eigen(object@.Data, symmetric=T, only.values=T)
    nonzero.e <- e$values[e$values > daeTolerance[["eigen.tol"]]]
    object@degfree <- length(nonzero.e)
  }
  object
}

print.projector <- function(x, ...)
{ 
  if (!inherits(x, "projector"))
    stop("Must supply an object of class projector")
  print(as(x, "matrix"))
  cat("degfree: ",x@degfree,"\n")
  invisible(x)
}
setMethod("show", "projector", function(object) print.projector(object))

proj2.eigen <- function(Q1, Q2)
{ 
  #A procedure to compute the eigenvalues and eigenvectors for the decomposition 
  #of Q1 pertaining to Q2  i.e. the common eigenvalues of Q1Q2Q1 and Q2Q1Q2 and the 
  #eigenvectors of Q1 in their joint dcomposition.
  #They are stored in a list with elements named efficiencies and eigenvectors
  if (!is.projector(Q1) | !is.projector(Q2))
    stop("Must supply valid objects of class projector")
  if (nrow(Q1) != nrow(Q2))
    stop("Matrices not conformable.")
  daeTolerance <- get("daeTolerance", envir=daeEnv)
  Q121 <- Q1 %*% Q2 %*% Q1
	eff <- eigen(Q121, symmetric=T)
	nonzero.eff <- eff$values[eff$values > daeTolerance[["eigen.tol"]]]
	r <- length(nonzero.eff)
	if (r==0)
	{ 
	  nonzero.eff <- 0
	  nonzero.eigen <- NULL
	} else
	{ 
	  nonzero.eigen <- eff$vectors[,1:r]
	  nonzero.eigen <- (abs(nonzero.eigen) > daeTolerance[["eigen.tol"]])* nonzero.eigen
	}
	list(efficiencies = nonzero.eff, eigenvectors = nonzero.eigen)
}

proj2.efficiency <- function(Q1, Q2)
{ 
  #A procedure to compute the canonical efficiency factors (eigenvalues) in the joint 
  #decomposition of Q1 and Q2 
  proj.Q1Q2 <- proj2.eigen(Q1, Q2)
  nonzero.eff <- proj.Q1Q2$efficiencies
	nonzero.eff
}

decomp.relate <- function(decomp1, decomp2)
{ 
  #A procedure to examine the the relationship between the eigenvectors in 
  #decomp1 and decomp2
  #decomp is a list produced by proj2.eigen or proj2.combine
  daeTolerance <- get("daeTolerance", envir=daeEnv)
  relmat <- crossprod(decomp1$eigenvectors, decomp2$eigenvectors)
  relmat <- (abs(relmat) > daeTolerance[["element.tol"]])* relmat
  dimnames(relmat) <- list(as.character(round(decomp1$efficiencies,4)), 
                           as.character(round(decomp2$efficiencies,4)))
  relmat
}

"proj2.combine" <- function(Q1, Q2)
{ #A procedure to compute the Residual operator for P remove Q when P and Q are nonorthogonal.
  
  #  Corresponding projection operator for Q in P is also obtained.
  n <- nrow(Q1)
  if (n != nrow(Q2))
    stop("Matrices not conformable.")
  isproj <- is.projector(Q1) & is.projector(Q2)
  Qconf <- projector(matrix(0, nrow = n, ncol = n))
  Qres <- Q1
  Eff.Q1.Q2 <- 0
  eigenvec <- NULL
  if (degfree(Q1) > 0)
  { 
    #compute efficiencies
    decomp <- proj2.eigen(Q1, Q2)
    Eff.Q1.Q2 <- decomp$efficiencies 
    if (length(Eff.Q1.Q2) == 1 & Eff.Q1.Q2[1]==0) #check matrices are orthogonal
    { 
      warning("Matrices are orthogonal.")
    } else
    { 
      daeTolerance <- get("daeTolerance", envir=daeEnv)
      EffUnique.Q1.Q2 <- remove.repeats(Eff.Q1.Q2)
      K <- length(EffUnique.Q1.Q2)
      #check for all confounded (i.e. eff = 1)
      if (K == 1 & EffUnique.Q1.Q2[1] == 1 & length(Eff.Q1.Q2) == degfree(Q2))
      { 
        Qconf <- projector(Q2)
        Qres <- projector(Q1 - Q2)
      } else if(length(Eff.Q1.Q2) == degfree(Q1)) # all of Q1 is confounded by Q2
      { 
        Qconf <- projector(Q1)
        Qres <- projector(matrix(0, nrow = nrow(Q1), ncol = ncol(Q1)))
      } else  
      { 
        if (K < 10) #compute projection operators for partially confounded case by recursion
        { 
          Qres <- Q1
          Q121 <- Q1 %*% Q2 %*% Q1
          for(eff in EffUnique.Q1.Q2[1:K])
          { 
            Qres <- Qres %*% (Q1 - (Q121/eff))
            Qres <- (Qres + t(Qres))/2  #force symmetry because know that it must be
          }
          #if have not got a projector, then compute the residual operator directly
          if (inherits(validProjector(Qres), what = "character"))
            Qres <- Q1 - Q1 %*% mat.ginv(Q2 %*% Q1 %*% Q2) %*% Q1
        } else #compute projection operators for partially confounded case by inversion
          Qres <- Q1 - Q1 %*% mat.ginv(Q2 %*% Q1 %*% Q2) %*% Q1
      
          #Form projectors
          Qres <- projector(Qres)
          if (degfree(Qres) > 0)
            Qconf <- projector(Q1 - Qres)
          else
            Qconf <- Q1
      } 
    }
    eigenvec <- decomp$eigenvectors
  }
  list(efficiencies = Eff.Q1.Q2, eigenvectors=eigenvec, Qconf = Qconf, Qres = Qres)
}