Skip to content

Commit

Permalink
At last, shaped nD-SSA
Browse files Browse the repository at this point in the history 
New saga begins!
  • Loading branch information
@eodus
eodus committed Aug 16, 2014
1 parent 69d3472 commit b3b8840
Show file tree
Hide file tree
Showing 4 changed files with 60 additions and 51 deletions.
47 changes: 23 additions & 24 deletions R/hbhankel.R
View file
Expand Up @@ -53,39 +53,37 @@ fft2 <- function(X, inverse = FALSE) {
x.dim + y.dim - 1,
x.dim)

names(input.dim) <- names(output.dim) <- names(x.dim)

list(input.dim = input.dim,
output.dim = output.dim)
}

convolve2 <- function(x, y, conj = TRUE, type = "circular") {
if (length(type) > 2)
warning("Incorrect argument length: length(type) > 2, two leading values will be used")
if (length(type) != 2)
type <- rep(type, 2)[1:2]
type <- sapply(type, match.arg, choices = c("circular", "open", "filter"))
convolven <- function(x, y, conj = TRUE, type = "circular") {
if (is.null(dim(x))) dim(x) <- length(x)
if (is.null(dim(y))) dim(y) <- length(y)

io.dims <- .convolution.dims(dim(x), dim(y), type)
input.dim <- io.dims$input.dim
output.dim <- io.dims$output.dim

X <- Y <- matrix(0, input.dim[1], input.dim[2])
storage.mode(x) <- storage.mode(y) <- "double"
storage.mode(conj) <- "logical"
storage.mode(input.dim) <- storage.mode(output.dim) <- "integer"

X[seq_len(nrow(x)), seq_len(ncol(x))] <- x
Y[seq_len(nrow(y)), seq_len(ncol(y))] <- y
res <- .Call("convolveN", x, y, input.dim, output.dim, conj)
if (length(output.dim) == 1)
dim(res) <- NULL
else
dim(res) <- output.dim

if (is.null(dim(X))) dim(X) <- length(X)
if (is.null(dim(Y))) dim(Y) <- length(Y)
storage.mode(X) <- storage.mode(Y) <- "double"
storage.mode(conj) <- "logical"
# tmp <- Re(fft2(fft2(X) * if (conj) Conj(fft2(Y)) else fft2(Y), inverse = TRUE)) / prod(input.dim)
tmp <- .Call("convolveN", X, Y, conj)
tmp[seq_len(output.dim[1]), seq_len(output.dim[2]), drop = FALSE]
res
}

.factor.mask.2d <- function(field.mask, window.mask, circular = FALSE) {
field.mask[] <- as.numeric(field.mask)
window.mask[] <- as.numeric(window.mask)
tmp <- convolve2(field.mask, window.mask, conj = TRUE,
tmp <- convolven(field.mask, window.mask, conj = TRUE,
type = ifelse(circular, "circular", "filter"))

abs(tmp - sum(window.mask)) < 0.5 # ==0, but not exact in case of numeric error
Expand All @@ -94,13 +92,14 @@ convolve2 <- function(x, y, conj = TRUE, type = "circular") {
.field.weights.2d <- function(window.mask, factor.mask, circular = FALSE) {
window.mask[] <- as.numeric(window.mask)
factor.mask[] <- as.numeric(factor.mask)
res <- convolve2(factor.mask, window.mask, conj = FALSE,
res <- convolven(factor.mask, window.mask, conj = FALSE,
type = ifelse(circular, "circular", "open"))
res[] <- as.integer(round(res))

res
}

# TODO generalize to ndim case
circle.mask <- function(R) {
I <- matrix(seq_len(2*R - 1), 2*R - 1, 2*R - 1)
J <- t(I)
Expand All @@ -118,10 +117,11 @@ triangle.mask <- function(side) {
new.hbhmat <- function(F, L = (N + 1) %/% 2,
wmask = NULL, fmask = NULL, weights = NULL,
circular = FALSE) {
if (length(circular) > 2)
warning("Incorrect argument length: length(circular) > 2, two leading values will be used")
if (length(circular) != 2)
circular <- rep(circular, 2)[1:2]
rank <- length(dim(F))
if (length(circular) > rank)
warning("Incorrect argument length: length(circular) > rank, two leading values will be used")
if (length(circular) != rank)
circular <- circular[(seq_len(rank) - 1) %% length(circular) + 1]

N <- dim(F)

Expand All @@ -141,8 +141,7 @@ new.hbhmat <- function(F, L = (N + 1) %/% 2,
mask <- weights > 0
F[!mask] <- mean(F[mask]) # Improve FFT stability & remove NAs
} else {
weights <- tcrossprod(.hweights.default(N[1], L[1]),
.hweights.default(N[2], L[2]))
weights <- .hweightsn(N, L)
}
storage.mode(weights) <- "integer"

Expand Down
33 changes: 19 additions & 14 deletions R/ssa.R
View file
Expand Up @@ -72,7 +72,7 @@ ssa <- function(x,
kind <- "cssa"
else if (inherits(x, "mts") || inherits(x, "data.frame") || inherits(x, "list") || inherits(x, "series.list"))
kind <- "mssa"
else if (is.matrix(x))
else if (is.matrix(x) || is.array(x))
kind <- "2d-ssa"
}
kind <- match.arg(kind)
Expand Down Expand Up @@ -153,30 +153,35 @@ ssa <- function(x,
column.projector <- row.projector <- NULL
}
} else if (identical(kind, "2d-ssa")) {
if (length(circular) > 2)
warning("Incorrect argument length: length(circular) > 2, two leading values will be used")
if (length(circular) != 2)
circular <- rep(circular, 2)[1:2]

# Coerce input to matrix if necessary
if (!is.matrix(x))
x <- as.matrix(x)

N <- dim(x);

mask <- if (is.null(mask)) !is.na(x) else mask & !is.na(x)
# Coerce input to array if necessary
if (!is.array(x))
x <- as.array(x)
N <- dim(x)

wmask <- .fiface.eval(substitute(wmask),
envir = parent.frame(),
circle = circle.mask,
triangle = triangle.mask)
ecall$wmask <- wmask
if (is.null(wmask)) {
wmask <- matrix(TRUE, L[1], L[2])
wmask <- array(TRUE, dim = L)
} else {
L <- dim(wmask)
}

# Fix rank (ndims) of x
rank <- length(L)
if (length(dim(x)) < rank)
dim(x) <- c(dim(x), rep(1, rank - length(dim(x))))

mask <- if (is.null(mask)) !is.na(x) else mask & !is.na(x)

# Check `circular' argument
if (length(circular) > rank)
warning("Incorrect argument length: length(circular) > rank, the leading values will be used")
if (length(circular) != rank)
circular <- circular[(seq_len(rank) - 1) %% length(circular) + 1]

K <- ifelse(circular, N, N - L + 1)

if (is.null(neig))
Expand Down
4 changes: 2 additions & 2 deletions inst/tests/test-sh2dssa.R
View file
Expand Up @@ -5,11 +5,11 @@ source(system.file("extdata", "common.test.methods.R", package = "Rssa"))
context("Shaped 2dSSA")

convolve2.open <- function(X, Y, conj = FALSE) {
convolve2(X, Y, conj = conj, type = "open")
convolven(X, Y, conj = conj, type = "open")
}

convolve2.filter <- function(X, Y, conj = TRUE) {
convolve2(X, Y, conj = conj, type = "filter")
convolven(X, Y, conj = conj, type = "filter")
}

test_that("new.hbhmat returns matrix with proper dimension", {
Expand Down
27 changes: 16 additions & 11 deletions src/hbhankel.c
View file
Expand Up @@ -618,10 +618,13 @@ SEXP hbhankelize_one_fft(SEXP U, SEXP V, SEXP hmat) {
return F;
}

SEXP convolveN(SEXP X, SEXP Y, SEXP Conj) {
SEXP DIM = getAttrib(X, R_DimSymbol);
R_len_t rank = length(DIM);
R_len_t *N = INTEGER(DIM);
SEXP convolveN(SEXP x, SEXP y,
SEXP input_dim, SEXP output_dim,
SEXP Conj) {
SEXP x_dim = getAttrib(x, R_DimSymbol);
SEXP y_dim = getAttrib(y, R_DimSymbol);
R_len_t rank = length(input_dim);
R_len_t *N = INTEGER(input_dim);
R_len_t pN = prod(rank, N), phN = hprod(rank, N);
int conjugate = LOGICAL(Conj)[0];

Expand All @@ -645,13 +648,15 @@ SEXP convolveN(SEXP X, SEXP Y, SEXP Conj) {
Free(revN);

/* Fill input buffer by X values*/
memcpy(circ, REAL(X), pN * sizeof(double));
memset(circ, 0, pN * sizeof(double));
fill_subarray(circ, REAL(x), rank, N, INTEGER(x_dim), 1);

/* Run the plan on X-input data */
fftw_execute_dft_r2c(r2c_plan, circ, ox);

/* Fill input buffer by Y values*/
memcpy(circ, REAL(Y), pN * sizeof(double));
memset(circ, 0, pN * sizeof(double));
fill_subarray(circ, REAL(y), rank, N, INTEGER(y_dim), 1);

/* Run the plan on Y-input data */
fftw_execute_dft_r2c(r2c_plan, circ, oy);
Expand All @@ -661,19 +666,19 @@ SEXP convolveN(SEXP X, SEXP Y, SEXP Conj) {
for (i = 0; i < phN; ++i)
oy[i] = conj(oy[i]);


/* Dot-multiply ox and oy, and divide by Nx*...*Nz*/
for (i = 0; i < phN; ++i)
oy[i] *= ox[i] / pN;

/* Compute the reverse transform to obtain result */
fftw_execute_dft_c2r(c2r_plan, oy, circ);


SEXP res;
PROTECT(res = allocVector(REALSXP, pN));
memcpy(REAL(res), circ, pN * sizeof(double));
setAttrib(res, R_DimSymbol, DIM);
PROTECT(res = allocVector(REALSXP, prod(rank, INTEGER(output_dim))));
fill_subarray(circ, REAL(res), rank, N, INTEGER(output_dim), 0);
/* setAttrib(output_dim, R_NamesSymbol, R_NilValue); */
setAttrib(res, R_DimSymbol, output_dim);
/* setAttrib(res, R_DimNamesSymbol, R_NilValue); */

/* Cleanup */
fftw_free(ox);
Expand Down

0 comments on commit b3b8840

Please sign in to comment.