4

I have found a very subtle bug in my R code just now. The following code takes a list of objects as input and create new fields for each of the objects.

Each object originally has two fields (w, p, s, u), and then I create more, beta, phi, etc.. The normal variables are OK. However the dynamic functions (Q, K, K1, K2) are not correct. Suppose I have two nigs, nigs[[1]] and nigs[[2]], the functions Q, K, K1 and K2 for nigs[[1]] would be the same as nigs[[2]]!

I just found this bug and would consult on how to get this code correct (while keeping its elegance:) Thanks!

  D <- length(nigs)

  for (i in 1:D) {
    w <- nigs[[i]]$w
    p <- nigs[[i]]$p
    s <- nigs[[i]]$s
    u <- nigs[[i]]$u

    nigs[[i]]$beta <- beta <- w / s * p * (1-p^2)^(-1/2);
    nigs[[i]]$phi <- phi <- w^2 / s^2;

    nigs[[i]]$z <- z <- (x-u)/s;
    nigs[[i]]$alpha_bar <- alpha_bar <- w * (1-p^2)^(-1/2);
    nigs[[i]]$y_bar <- y_bar <- sqrt(1+z^2);

    nigs[[i]]$Q <- Q <- function(t) { sqrt(1 - (2*beta*t+t^2)/phi) }
    nigs[[i]]$K <- K <- function(t) { u*t - w*Q(t) + w }
    nigs[[i]]$K1 <- K1 <- function(t) { (u + w * (beta+t) / (Q(t)*phi)) }
    nigs[[i]]$K2 <- K2 <- function(t) { qt = Q(t); (w/(qt * phi) + w * (beta+t)^2 / (qt^3 * phi^2)); }
  }

EDIT

The primary error I made is that I assumed that for { } introduced new scopes, in that case, w,p,s,u are different w,p,s,u every time, actually not. Only functions in R introduce new scopes. And this scoping rule is different from C/Java.

Improve this question

2 Answers 2

Reset to default
6

That is normal behavior of the lexical scope. You can use closure instead.

f <- list()
g <- list()
for (i in 1:2) {
    j <- i * 2
    f[[i]] <- function() print(j)
    g[[i]] <- (function() {j <- j; function() print(j)}) ()
}

then,

> for (i in 1:2) f[[i]]()
[1] 4
[1] 4
> for (i in 1:2) g[[i]]()
[1] 2
[1] 4
Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

I think using local is a bit more elegant here.
6

In object oriented terminology each nigs[[i]] is an object and the functions Q, K, etc. are methods which act on the object's properties w, p, etc. Using the proto package we set each nigs[[i]] to a proto object and then update the object as indicated. Note that all methods take the object as the first argument so if p is a proto object containing method Q then p$Q(t) means to look in p for Q and then run it with the arguments p and t so p$Q(t) is the same as with(p, Q(p, t)). Thus we have added the extra first argument to each of the methods below. See proto home page for more.

library(proto)

# initialize
x <- 1
nigs <- lapply(1:2, function(i) proto(w = i/3, p = i/3, s = i/3, u = i/3))

for(p in nigs) with(p, {
    beta <- w / s * p * (1-p^2)^(-1/2)
    phi <- w^2 / s^2

    z <- (x-u)/s
    alpha_bar <- w * (1-p^2)^(-1/2)
    y_bar <- sqrt(1+z^2)

    Q <- function(., t) { sqrt(1 - (2*beta*t+t^2)/phi) }
    K <- function(., t) { u*t - w*.$Q(t) + w }
    K1 <- function(., t) { (u + w * (beta+t) / (.$Q(t)*phi)) }
    K2 <- function(., t) { 
        qt = .$Q(t)
        (w/(qt * phi) + w * (beta+t)^2 / (qt^3 * phi^2)) 
    }
  })

EDIT: A second possible design would be to create a parent object, meths to hold the methods instead of defining them over again in each separate proto object. In that case, within each method we must be sure that we use the properties of the object passed in the first argument since the methods and properties are now located in different objects:

meths <- proto(
        Q = function(., t) sqrt(1 - (2*.$beta*t+t^2)/.$phi),
        K = function(., t) .$u*t - .$w*.$Q(t) + .$w,
        K1 = function(., t) (.$u + .$w * (.$beta+t) / (.$Q(t)*.$phi)),
        K2 = function(., t) {
            qt = .$Q(t)
            (.$w/(qt * .$phi) + .$w * (.$beta+t)^2 / (qt^3 * .$phi^2)) 
        }
)

# initialize - meths$proto means define proto object with parent meths
x <- 1
nigs <- lapply(1:2, function(i) meths$proto(w = i/3, p = i/3, s = i/3, u = i/3))

for(p in nigs) with(p, {
    beta <- w / s * p * (1-p^2)^(-1/2)
    phi <- w^2 / s^2

    z <- (x-u)/s
    alpha_bar <- w * (1-p^2)^(-1/2)
    y_bar <- sqrt(1+z^2)
})

Now the following works by looking up Q in nigs[[1]] but not finding it there looking into its parent, meths, and running the Q found there. In nigs[[1]]$Q(.1) the call implicitly passes nigs[[1]] to Q as its first argument and we have defined all properties within the body of Q relative to the first argument so everything works:

> nigs[[1]]$Q(.1)
[1] 0.9587958
Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.