# [R] How to get rid of loop?

Martin Morgan mtmorgan at fhcrc.org
Mon Apr 27 20:53:41 CEST 2009

```Uwe Ligges <ligges at statistik.tu-dortmund.de> writes:

> Ken-JP wrote:
>> The code below shows what I'm trying to get rid of.
>> If there is no way to get rid of the loop, I will try to use
>> package( inline
>> ).
>> I'm just curious as to whether there is a "vector way" of doing this
>> algorithm.
>
>
> I don't see a vector way. If speed is an issue, I'd suggest to port

Parsing the code, it seems like a run of 1's starts with a value
greater than .75, and continues until a value less than .5. So this...

x1 <- x > .75
c(x1[[1]], x1[-1] & !x1[-n])

flags the start of each run, and this

y1 <- lapply(split(x > .5, cumsum(c(x1[[1]], x1[-1] & !x1[-n]))),
cumprod)

splits candidate runs x > .5 into subsets beginning at each start
point, and then processes the subset to extend the run as long as the
values are greater than .5. My 'vectorized' version, with an imperfect
solution to the first run and glossing over the point x==.5 in the
original code, is

encode.2 <- function(x) {
n <- length(x)
x0 <- x < .25
y0 <- lapply(split(x <= .5, cumsum(c(x0[[1]], x0[-1] & !x0[-n]))),
cumprod)
x1 <- x > .75
y1 <- lapply(split(x > .5, cumsum(c(x1[[1]], x1[-1] & !x1[-n]))),
cumprod)
as.vector((cumsum(abs(x-.5) > .25) != 0) *
(-unlist(y0, use.names=FALSE) + unlist(y1, use.names=FALSE)))
}

this seems to be 7-20x faster than encode.1

encode.1 <- function(x) {
n <- length( x )
y <- rep(NA, n)
yprev <- 0;
for ( i in (1:n)) {
if ( x[i]>0.75 ) {
y[i] <- 1;
} else if ( x[i]<0.25 ) {
y[i] <- -1;
} else
if ( yprev==1 & x[i]<0.5) {
y[i] <- 0;
} else if ( yprev==-1 & x[i]>0.5) {
y[i] <- 0;
} else {
y[i] <- yprev
}
yprev
<-
y[i];
}
y
}

> set.seed(1)
> mean(replicate(10, system.time(encode.1(runif(100000)), TRUE)[[3]]))
[1] 0.7256
> set.seed(1)
> mean(replicate(10, system.time(encode.2(runif(100000)), TRUE)[[3]]))
[1] 0.0983

> mean(replicate(10, system.time(encode.1(integer(100000)), TRUE)[[3]]))
[1] 0.6288
> mean(replicate(10, system.time(encode.2(integer(100000)), TRUE)[[3]]))
[1] 0.0338

and, modulo x==.5, even seems to produce the same result ;)

> set.seed(1); res.1 <- replicate(10, encode.1(runif(1000)))
> set.seed(1); res.2 <- replicate(10, encode.2(runif(1000)))
> identical(res.1, res.2)
[1] TRUE

Algorithm f7 of

https://stat.ethz.ch/pipermail/r-devel/2008-April/049111.html

might point to fast (comparable to C) R implementations.

Martin

> this small part to C, it's very simple and may yield a considerable
> performance boost (untested). It can probably be used as a textbook
> example for code where porting to C make sense.
>
> Best,
> Uwe Ligges
>
>
>
>> #
>> -----------------------------------------------------------------------------------------
>> set.seed(1)
>> x <- runif(100)
>> n <- length( x )
>> y <- rep(NA, n)
>> yprev <- 0;
>> for ( i in (1:n)) {
>> 	if ( x[i]>0.75 ) {
>> 		y[i] <- 1;
>> 	} else if ( x[i]<0.25 ) {
>> 		y[i] <- -1;
>> 	} else if ( yprev==1 & x[i]<0.5) {
>> 		y[i] <- 0;
>> 	} else if ( yprev==-1 & x[i]>0.5) {
>> 		y[i] <- 0;
>> 	} else {
>> 		y[i] <- yprev
>> 	}
>> 	yprev <- y[i];
>> }
>>
>>> y
>>   [1]  0  0  0  1 -1  1  1  1  1 -1 -1 -1  0  0  1  0  0  1  0  1  1 -1  0
>> -1 -1
>>  [26] -1 -1 -1  1  0  0  0  0 -1  1  1  1 -1  0  0  1  1  1  1  1  1
>> -1 -1
>> 0  0
>>  [51]  0  1  0 -1 -1 -1 -1  0  0  0  1  0  0  0  0  0  0  1 -1  1  0
>> 1  0 0  0
>>  [76]  1  1  0  1  1  0  0  0  0  1 -1  0 -1 -1 -1 -1 -1  0  1  1  1
>> 0  0 1  1
>>
>>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.

--
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793

```