Nothing
# This class represents a matrix A = (1/rowSums(Toep)) * Toep
# where Toep is a Toeplitz matrix.
# This special structure allows fast matrix-vector and exponential computations.
#
# A Toeplitz matrix T of size NxN can be embedded into a circulant matrix of
# size 2Nx2N [ T S ]
# [ S T ]
# This circulant matrix is fully determined by its first column c.
# Then T*x = (I 0) A [ x ]
# [ 0 ]
# and A*X0 = IFFT(FFT(c) * FFT(X0))
# In order to improve FFT performance, we enlarge A such that its dimension
# is a power of two.
MPhiFFT <- setClass(
"MPhiFFT",
slots = c(
method= "character",
d0 = "vector",
firstColumn = "vector",
firstRow = "vector",
fftColumn = "vector",
n = "numeric",
zero = "numeric"
),
contains="MPhiAbstract"
)
setMethod(
f="initialize",
signature=c("MPhiFFT"),
definition=function(.Object, firstColumn, firstRow, ...){
callNextMethod()
args <- list(...)
if(is.null(args$method)) .Object@method="FFT"
else .Object@method=args$method
.Object@firstColumn <- firstColumn
.Object@firstRow <- firstRow
if(length(.Object@firstColumn) != length(.Object@firstRow)) {
stop("Error: first column and first row must have the same size")
}
if(.Object@firstColumn[1] != .Object@firstRow[1]) {
stop("Error: first row and first column must have the same first element")
}
.Object@n <- length(firstColumn)
.Object@d0 <- sapply(1:length(firstColumn), function(irow) {return(1/sum(c(firstColumn[irow:2], firstRow[1:(.Object@n+1-irow)])))})
.Object@d0[1] <- 1/sum(firstRow)
# Compute the next higher power of 2 for FFT, it is better to perform FFT on
# vectors of length power of 2.
N2 = 2 * .Object@n
size=2^ceiling(log(N2)/log(2))
circulantColumn=c(firstColumn, rep(0, size + 1 - 2*.Object@n), firstRow[length(firstRow):2])
.Object@zero <- rep(0, size - .Object@n)
.Object@fftColumn <- fft(circulantColumn)
return(.Object)
}
)
setMethod(
f="applyV",
signature=c("MPhiFFT","vector"),
definition=function(object,x){
return(object@d0 * Re(ifft_panda(object@fftColumn * fft(c(x,object@zero)))[1:object@n]))
}
)
# Compute exp(diag(diag1)+M*diag2)%*%x = exp(diag(diag1)+diag(diag2)%*%M)%*%x
# = exp(diag1*x + (diag2*d0)*(toeplitz%*%x))
setMethod(
f="expV",
signature=c("MPhiFFT","vector","vector","vector"),
definition=function(object,diag1,diag2,x){
mask_diag1 = diag1 < -20
x[mask_diag1] = 0
expGv=x
Gnv=x
i=1
epsnormv = 1e-10 * sum(abs(x))
MATVECT <- selectMethod(applyV, c("MPhiFFT","vector"))
while (sum(abs(Gnv)) > epsnormv){
Gnv = (diag1*Gnv + diag2*MATVECT(object,Gnv)) / i
Gnv[mask_diag1] = 0
expGv = expGv + Gnv
i=i+1
}
expGv[mask_diag1] = 2e-9
return(expGv)
}
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.