How to vectorize this process in R?

Question

How to vectorize this process in R without using too many loops ?

I have this function :

HM=function(CO,CS,CD,CSD){
  if(CO-CS)>1){
    return(2^(CS)/(2^(CO)-2^(CSD)))
  }
  else if(CO-CD)>1){
    return(1-2^(CD)/(2^(CO)-2^(CSD)))
  }
return(0)
}

Basically I need to get HM value for every combination of {CO,CS,CD,CSD} over thoses values :

CO  25.76031685 25.71126747 25.90163231
CS  24.40528297 24.09929848 23.51999092
CD  25.99405861 25.72906113 25.61374474
CSD 35.94195557 36.07263184 34.00024414

So I need to get those values :

HM(25.76031685,24.40528297,25.99405861,35.94195557)
HM(25.71126747,24.40528297,25.99405861,35.94195557)
HM(25.90163231,24.40528297,25.99405861,35.94195557)
HM(25.76031685,24.09929848,25.99405861,35.94195557)
HM(25.71126747,24.09929848,25.99405861,35.94195557)
HM(25.90163231,24.09929848,25.99405861,35.94195557)
HM(25.76031685,23.51999092,25.99405861,35.94195557)
HM(25.71126747,23.51999092,25.99405861,35.94195557)
HM(25.90163231,23.51999092,25.99405861,35.94195557)
etc...

Basically it's all the combination with 4 vectors of 3 elements :

Vectors :
a=c(1,2,3)
b=c(1,2,3)
c=c(1,2,3)
d=c(1,2,3)

Combinations :
1,1,1,1
2,1,1,1
1,2,1,1
1,1,2,1
1,1,1,2
3,1,1,1
1,3,1,1
etc...

I'm not sure how to count the number of combinations. Of course I could use 4 nested loops but I want to learn how to do it with vectorization since R is too slow for loops. I think we can use expand.grid but I don't know how. Also the table is in excel, I can export it in .csv but I'm not sure of the best way to implement this stuff so thank's for your help !

IRTFM · Answer 1 · 2013-06-21T22:52:39.820

The answer is rather uninteresting in this case because none of the conditions holds for these values and all zeros are returned:

> tdat  #  dataframe version of that data.
         CO       CS       CD      CSD
V2 25.76032 24.40528 25.99406 35.94196
V3 25.71127 24.09930 25.72906 36.07263
V4 25.90163 23.51999 25.61374 34.00024
> with( tdat, 
       ifelse( (CS-CO) > 1, 2^(CS)/(2^(CO)-2^(CSD)),  #1st consequent
                 ifelse ( (CD-CO) > 1, 1-2^(CD)/(2^(CO)-2^(CSD)), # 2nd
                                           0) ) )  # default
[1] 0 0 0

To do it on a matrix versioon of that data you need to first correct the mismatched parens in your code and then use apply while reference a single passed x value with the rownames:

mdat <- 
structure(c(25.76032, 24.40528, 25.99406, 35.94196, 25.71127, 
24.0993, 25.72906, 36.07263, 25.90163, 23.51999, 25.61374, 34.00024
), .Dim = c(4L, 3L), .Dimnames = list(c("CO", "CS", "CD", "CSD"
), NULL))

> apply(mdat, 2, function(x){
+   if( (x['CS']-x['CO'])>1){
+     return(2^(x['CS'])/(2^(x['CO'])-2^(x['CSD'])))
+   }
+   else if( (x['CD']-x['CO'])>1){
+     return(1-2^(x['CD'])/(2^(x['CO'])-2^(x['CSD'])))
+   }
+ return(0)
+ })
[1] 0 0 0

Thank you for your reply ! However I don't understand there is 3^4 combination possible for this list {3,3,3,3} so why are there only 3 returns ? By the way you were right I did a mistake, the condition is CO-CS|CD not the reverse. — Wicelo, Jun 22 '13 at 00:24
it's not 3^4 in fact I'm not sure how to count this but I added some examples of HM values I need. — Wicelo, Jun 22 '13 at 02:45
Sigh. It's probably going to involve `expand grid`, but without an example I don't see much point in continuing this guessing game. — IRTFM, Jun 22 '13 at 03:26

score 1 · Accepted Answer · answered Jun 22 '13 at 13:17

You can use expand.grid to get all combinations. But you need first to vectorize your function HM, using ifelse instead of if:

HM2 <- function(CO,CS,CD,CSD)
{
    den <- 2^CO-2^CSD

    ifelse(CO-CS>1, 2^CS/den,
        ifelse(CO-CD>1, 1-2^CD/den, 0))
}

Note that den is common to both results.

Now your data:

CO <- c(25.76031685, 25.71126747, 25.90163231)
CS <- c(24.40528297, 24.09929848, 23.51999092)
CD <- c(25.99405861, 25.72906113, 25.61374474)
CSD <- c(35.94195557, 36.07263184, 34.00024414)

The combinations:

cmbs <- expand.grid(CO, CS, CD, CSD)
names(cmbs) <- c("CO", "CS", "CD", "CSD")

Example:

> head(cmbs)
        CO       CS       CD      CSD
1 25.76032 24.40528 25.99406 35.94196
2 25.71127 24.40528 25.99406 35.94196
3 25.90163 24.40528 25.99406 35.94196
4 25.76032 24.09930 25.99406 35.94196
5 25.71127 24.09930 25.99406 35.94196
6 25.90163 24.09930 25.99406 35.94196

The final result can be obtained using within, to perform calculations inside the dataframe:

result <- within(cmbs, HM <- HM2(CO, CS, CD, CSD))

Example:

> head(result)
        CO       CS       CD      CSD            HM
1 25.76032 24.40528 25.99406 35.94196 -0.0003368911
2 25.71127 24.40528 25.99406 35.94196 -0.0003368814
3 25.90163 24.40528 25.99406 35.94196 -0.0003369210
4 25.76032 24.09930 25.99406 35.94196 -0.0002725079
5 25.71127 24.09930 25.99406 35.94196 -0.0002725000
6 25.90163 24.09930 25.99406 35.94196 -0.0002725321

Awesome I came to the same conclusion during this night. Thank's ! However when I add a third ifelse() to HM() it doesn't work anymore do you know why ? http://stackoverflow.com/questions/17252466/why-with-in-r-is-doing-vector-operation-in-one-case-and-not-in-the-other — Wicelo, Jun 22 '13 at 15:54
@Wicelo, it looks like you've figured it out. The problem was `&&`. Roland's answer is very nice, with a fully vectorized approach. See ya! — Ferdinand.kraft, Jun 22 '13 at 21:36

How to vectorize this process in R?

2 Answers2