Tenis Common opponents R I'm follow this article about "common opponent in tennis", my goal is to script it in the most efficient way. Below you can find my code but is so slow. To calculate the result of 1 match my laptop spent 120seconds more or less, and I have a dataset of 150k of rows to calculate.

The article: (PDF)

I need your help to clean and improve my code. Links to the data used:

tableA

tableB

data_tennis_co

data_tennis_co = read.table("test_co.csv", header=FALSE, sep=",", fill = TRUE)
A = read.table("tableA.csv", header=FALSE, sep=",", fill = TRUE)
B = read.table("tableB.csv", header=FALSE, sep=",", fill = TRUE)

#BASIC FUNCTIONS
G&lt;-function(p){res&lt;- p^4*(15-4*p-((10*p^2)/(1-2*p*(1-p))))}
d&lt;- function(p,q) {res&lt;- p*q*(1-(p*(1-q)+(1-p)*q))^-1}
TB &lt;- function(p,q) {res &lt;- foreach(i = seq_along(1:28), .combine = sum) %dopar% {tb&lt;-A[i,1]*p^A[i,2]*(1-p)^A[i,3]*q^A[i,4]*(1-q)^A[i,5]*d(p,q)^A[i,6]}}
S &lt;- function(p,q) {res &lt;- foreach(i = seq_along(1:21), .combine = rbind) %dopar% {s&lt;-B[i,1]*G(p)^B[i,2]*(1-G(p))^B[i,3]*G(q)^B[i,4]*(1-G(q))^B[i,5]*(G(p)*G(q)+(G(p)*(1-G(q))+(1-G(p))*G(q))*TB(p,q))^B[i,6]} sum(res)}
M3 &lt;- function(p,q) {res &lt;- S(p,q)^2*(1+2*(1-S(p,q)))}
DELTA_AB &lt;- function(spwAC,rpwAC,spwBC,rpwBC) {res &lt;- (spwAC-(1-rpwAC))-(spwBC-(1-rpwBC))}
PR&lt;- function(spwAC,rpwAC,spwBC,rpwBC) {res &lt;- (M3(0.6+DELTA_AB(spwAC,rpwAC,spwBC,rpwBC),(1-0.6))+M3(0.6,(1-(0.6-DELTA_AB(spwAC,rpwAC,spwBC,rpwBC)))))/2}

#COMMON OPPONENTS
MAL&lt;-function(id1,id2){
prova&lt;- subset(data_tennis_co, V3 == 1 &amp; V4==2)
  previous&lt;-subset(data_tennis_co, V2 &lt; prova$V2)
  s1 &lt;- subset(previous, V3 == 1 | V4==1)
  s2 &lt;- subset(previous, V3 ==2 | V4==2)
  s1$opp &lt;- ifelse(s1$V3==1, s1$V4, s1$V3)
  s2$opp &lt;- ifelse(s2$V3==2, s2$V4, s2$V3)
  inn&lt;- intersect(s1$opp,s2$opp)
  common1&lt;-s1[s1$opp %in% inn,]
  common2&lt;-s2[s2$opp %in% inn,]
  # fare media se id non unico
  COM &lt;- merge(common1, common2,by=c("opp"))
  COM$OMALLEY &lt;- unlist(mapply(PR, COM$V5.x, COM$V6.x, COM$V7.y, COM$V8.y))
  COM$OMALLEY[is.nan(COM$OMALLEY)] &lt;- 0.5
  return(tryCatch(sum(COM$OMALLEY)/nrow(COM), error=function(e) NaN))
}

RESA&lt;-MAL(1,2)
You could help us by providing the output of the profiler:

Rprof(tmp &lt;- tempfile())
RESA &lt;- MAL(1, 2)
Rprof()
summaryRprof(tmp)
unlink(tmp)


, but I suspect the main bottleneck in your code is the use of parallel loops in TB and S, for operations that can be done faster using vectorized R functions. Please give this a try and let us know if it is faster for you:

G  &lt;- function(p)    p^4*(15-4*p-((10*p^2)/(1-2*p*(1-p))))
d  &lt;- function(p, q) p*q*(1-(p*(1-q)+(1-p)*q))^-1
TB &lt;- function(p, q) sum(A[,1] * p^A[,2] * (1-p)^A[,3] *
                         q^A[,4] * (1-q)^A[,5] * d(p,q)^A[,6])
S  &lt;- function(p, q) {
   Gp &lt;- G(p)
   Gq &lt;- G(q)
   sum(B[,1] * Gp^B[,2] * (1-Gp)^B[,3] * Gq^B[,4] * (1-Gq)^B[,5] *
      (Gp*Gq+(Gp*(1-Gq)+(1-Gp)*Gq)*TB(p,q))^B[,6])
}
M3 &lt;- function(p, q) {
   s &lt;- S(p,q)
   s^2*(1+2*(1-s))
}
DELTA_AB &lt;- function(spwAC,rpwAC,spwBC,rpwBC) (spwAC-(1-rpwAC)) -
                                              (spwBC-(1-rpwBC))
PR &lt;- function(spwAC,rpwAC,spwBC,rpwBC) {
   D &lt;- DELTA_AB(spwAC, rpwAC, spwBC, rpwBC)
  (M3(p = 0.6 + D, q = (1 - 0.6)) +
   M3(p = 0.6, q = 1 - (0.6 - D))) / 2
}