R - Necesita ayuda para acelerar un bucle for

Question

Tengo dos dataframes; uno es de 48 filas de largo y tiene el siguiente aspecto:

name = Z31

 Est.Date Site Cultivar Planting 
1 24/07/2011 Birchip Axe   1 
2 08/08/2011 Birchip Bolac   1 
3 24/07/2011 Birchip Derrimut  1 
4 12/08/2011 Birchip Eaglehawk  1 
5 29/07/2011 Birchip Gregory  1 
6 29/07/2011 Birchip Lincoln  1 
7 23/07/2011 Birchip Mace   1 
8 29/07/2011 Birchip Scout   1 
9 17/09/2011 Birchip Axe   2 
10 19/09/2011 Birchip Bolac   2 

El otro es> 23000 filas y contiene la salida de un simulador. Se ve así:

name = pred

Date  maxt mint Cultivar Site Planting tt cum_tt 
1 5/05/2011 18  6.5 Axe  Birchip 1  12.25 0 
2 6/05/2011 17.5  2.5 Axe  Birchip 1  10  0 
3 7/05/2011 18  2.5 Axe  Birchip 1  10.25 0 
4 8/05/2011 19.5  2 Axe  Birchip 1  10.75 0 
5 9/05/2011 17  4.5 Axe  Birchip 1  10.75 0 
6 10/05/2011 15.5 -0.5 Axe  Birchip 1  7.5 0 
7 11/05/2011 14  5.5 Axe  Birchip 1  9.75 0 
8 12/05/2011 19   8 Axe  Birchip 1  13.5 0 
9 13/05/2011 18.5  7.5 Axe  Birchip 1  13  0 
10 14/05/2011 16  3.5 Axe  Birchip 1  9.75 0 

Lo que quiero hacer es tener la columna de la cum_tt empezar a añadir la columna de la tt de la fila actual a la cum_tt de la fila anterior (adición acumulativa) SOLAMENTE si la fecha en el DF pred es igual o mayor que el Z31 est.Date. He escrito el siguiente bucle:

for (i in 1:nrow(Z31)){ 
    for (j in 1:nrow(pred)){ 
    if (Z31[i,]$Site == pred[j,]$Site & Z31[i,]$Cultivar == pred[j,]$Cultivar & Z31[i,]$Planting == pred[j,]$Planting & 
     pred[j,]$Date >= Z31[i,]$Est.Date) 
    { 
     pred[j,]$cum_tt <- pred[j,]$tt + pred[j-1,]$cum_tt 
    } 
    } 
} 

Esto funciona, pero es tan lento que tomará alrededor de una hora para ejecutar todo el conjunto. Sé que los bucles no son el punto fuerte de R así que ¿alguien puede ayudarme a vectorizar esta operación?

Gracias de antemano.

ACTUALIZACIÓN

Aquí está la salida de dput (Z31):

structure(list(Est.Date = structure(c(15179, 15194, 15179, 15198, 15184, 15184, 15178, 15184, 15234, 15236, 15230, 15238, 15229, 15236, 15229, 15231, 15155, 15170, 15160, 15168, 15165, 15159, 15170, 15170, 15191, 15205, 15198, 15203, 15202, 15195, 15203, 15206, 15193, 15193, 15195, 15200, 15193, 15205, 15200, 15205, 15226, 15245, 15231, 15259, 15241, 15241, 15241, 15241), class = "Date"), Site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Birchip", "Gatton", "Tarlee"), class = "factor"), Cultivar = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), .Label = c("Axe", "Bolac", "Derrimut", "Eaglehawk", "Gregory", "Lincoln", "Mace", "Scout"), class = "factor"), Planting = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L)), .Names = c("Est.Date", "Site", "Cultivar", "Planting"), row.names = c(NA, -48L), class = "data.frame")

Aquí está Pred. Tenga en cuenta que hay algunas columnas adicionales aquí. Acabo de incluir las relevantes arriba para facilitar la lectura.

structure(list(Date = structure(c(15099, 15100, 15101, 15102, 
15103, 15104, 15105, 15106, 15107, 15108, 15109, 15110, 15111, 
15112, 15113, 15114, 15115, 15116, 15117, 15118), class = "Date"), 
    flowering_das = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Zadok = c(9, 9, 
    9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 11, 11.032, 11.085, 
    11.157), stagename = structure(c(8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 9L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 3L, 3L, 3L), .Label = c("emergence", 
    "end_grain_fill", "end_of_juvenil", "floral_initiat", "flowering", 
    "germination", "maturity", "out", "sowing", "start_grain_fi" 
    ), class = "factor"), node_no = c(0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 2, 2.032, 2.085, 2.157), maxt = c(18, 
    17.5, 18, 19.5, 17, 15.5, 14, 19, 18.5, 16, 16, 15, 16.5, 
    16.5, 20.5, 23, 25.5, 16.5, 16.5, 15), mint = c(6.5, 2.5, 
    2.5, 2, 4.5, -0.5, 5.5, 8, 7.5, 3.5, 6, 1, 5.5, 2, 7, 7, 
    9, 13.5, 11.5, 8.5), Cultivar = c("Axe", "Axe", "Axe", "Axe", 
    "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", 
    "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe"), Site = c("Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip"), Planting = c("1", "1", "1", "1", "1", "1", "1", 
    "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
    "1"), `NA` = c("Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out" 
    ), tt = c(12.25, 10, 10.25, 10.75, 10.75, 7.5, 9.75, 13.5, 
    13, 9.75, 11, 8, 11, 9.25, 13.75, 15, 17.25, 15, 14, 11.75 
    ), cum_tt = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0)), .Names = c("Date", "flowering_das", "Zadok", 
"stagename", "node_no", "maxt", "mint", "Cultivar", "Site", "Planting", 
NA, "tt", "cum_tt"), row.names = c(NA, 20L), class = "data.frame") 

ACTUALIZACIÓN

Gracias a todos por su ayuda. Todavía soy nuevo en la forma vectorial de hacer las cosas y no fui capaz de implementar algunas de las soluciones más complejas a tiempo. Tengo algunos momentos a continuación por la forma en que sugirió Subs. Es lo suficientemente rápido ahora para hacer lo que necesito. Estos son el número de segundos durante una iteración de Z sobre revés

mi camino: 59.77

Subs: 14,62

Subs utilizando fecha numérica: 11.12

Answer 1

seguro de que podemos hacer esto en pocos segundos ... mi primera respuesta de aquí a fin de ir suave

## first make sure we have dates in a suitable format for comparison 
## by using strptime, creating the columns estdate_tidy and date_tidy 
## in Z31 and pred respectively 

Z31$estdate_tidy = strptime(as.character(Z31$Est.Date), "%d/%m/%Y") 
pred$date_tidy = strptime(as.character(pred$Date), "%d/%m/%Y") 

## now map the estdate_tidy column over to pred using the match command - 
## Z31_m and pred_m are dummy variables that hopefully make this clear 

Z31_m = paste(Z31$Site, Z31$Cultivar, Z31$Planting) 
pred_m = paste(pred$Site, pred$Cultivar, pred$Planting) 
pred$estdate_tidy = Z31$estdate_tidy[match(pred_m, Z31_m)] 

## then define a ttfilter column that copies tt, except for being 0 when 
## estdate_tidy is after date_tidy (think this is what you described) 

pred$ttfilter = ifelse(pred$date_tidy >= pred$estdate_tidy, pred$tt, 0) 

## finally use cumsum function to sum this new column up (looks like you 
## wanted the answer in cum_tt so I've put it there) 

pred$cum_tt = cumsum(pred$ttfilter) 

Espero que esto ayude :)

ACTUALIZACIÓN (7 de junio):

El código vectorizado para hacer frente a la nueva especificación - es decir, que la cumSum debe realizarse por separado para cada conjunto de condiciones (sitio/Cultivar/plantación) - se muestra a continuación:

Z31$Lookup=with(Z31,paste(Site,Cultivar,Planting,sep="~")) 
Z31$LookupNum=match(Z31$Lookup,unique(Z31$Lookup)) 
pred$Lookup=with(pred,paste(Site,Cultivar,Planting,sep="~")) 
pred$LookupNum=match(pred$Lookup,unique(pred$Lookup)) 

pred$Est.Date = Z31$Est.Date[match(pred$Lookup,Z31$Lookup)] 
pred$ttValid = (pred$Date>=pred$Est.Date) 
pred$ttFiltered = ifelse(pred$ttValid, pred$tt, 0) 

### now fill in cumsum of ttFiltered separately for each LookupNum 
pred$cum_tt_Z31 = as.vector(unlist(tapply(pred$ttFiltered, 
              pred$LookupNum,cumsum))) 

Runtime es 0.16 segundos en mi máquina , y la columna final pred$cum_tt_Z31 coincide exactamente con la respuesta del código no vectorizado :)

Para completar, vale la pena señalar que la línea final complicada anterior puede ser reemplazada por el siguiente enfoque más simple con un ciclo corto sobre el 48 casos posibles:

pred$cum_tt_Z31 = rep(NA, nrow(pred)) 
for (lookup in unique(pred$Lookup)) { 
    subs = which(pred$Lookup==lookup) 
    pred$cum_tt_Z31[subs] = cumsum(pred$ttFiltered[subs]) 
} 

El tiempo de ejecución solo aumenta ligeramente a 0,25 segundos aproximadamente, porque el ciclo aquí es muy pequeño y el trabajo realizado dentro del ciclo se vectoriza.

¡Creo que lo hemos descifrado! :)

Algunas observaciones rápidas sobre la vectorización (8 de junio):

El proceso de vectorización los pasos del proceso de tiempo de ejecución presentadas abajo desde cerca de una hora a 0.16 segundos en total. Incluso teniendo en cuenta las diferentes velocidades de la máquina, esta es una aceleración de al menos un factor de 10.000, lo que empequeñece los factores de 2-5 que uno podría obtener al realizar ajustes menores, pero conservando una estructura de bucle.

La primera observación clave que se debe hacer: en la solución, cada línea crea, sin bucle, un nuevo vector completo con la misma longitud que las columnas en Z31 o pred. Para mayor pulcritud, a menudo me resulta útil crear estos nuevos vectores como nuevas columnas de marcos de datos, pero obviamente esto no es estrictamente necesario.

Segunda observación: la columna Est.Date requerida se transfiere correctamente de Z31 a pred utilizando una estrategia de "pegar 'en emparejamiento' '. Existen enfoques alternativos para este tipo de tarea (por ejemplo, usar merge), pero tomo esta ruta ya que es completamente a prueba de fallas y garantiza preservar el orden en pred (que es crucial). Esencialmente, la operación de pegar solo le permite hacer coincidir varios campos a la vez, porque si las cadenas pegadas coinciden, todas sus partes constituyentes coinciden. Utilizo ~ como separador (siempre que sepa que el símbolo no aparecerá en ninguno de los campos) para evitar que la operación de pegado cause ambigüedad. Si usa un separador de espacios, pegar juntos algo así como ("AB", "C", "D") daría el mismo resultado que pegar ("A", "BC", "D") - y queremos evitar dolores de cabeza!

Tercera observación: es fácil vectorizar operaciones lógicas como ver si un vector excede a otro (ver pred $ ttValid), o elegir un valor u otro basado en el valor de un vector (ver pred $ ttFiltered). En la situación actual, estos podrían combinarse en una sola línea, pero rompí las cosas un poco más como una demostración.

Cuarta observación: la última línea que crea pred $ cum_tt_Z31 básicamente hace una operación cumsum en las filas correspondientes a cada valor separado de pred $ LookupNum, usando tapply, que le permite aplicar la misma función sobre diferentes grupos de filas (aquí, estamos agrupando por pred $ LookupNum). La definición de pred $ LookupNum ayuda enormemente aquí: es un índice numérico con un bloque de 1, seguido de un bloque de 2, y así sucesivamente. Esto significa que la lista final de vectores cumsum que sale de la aplicación puede simplemente no ser incluida en la lista y poner en un vector y está automáticamente en el orden correcto. Si realiza un tapply y se divide por grupos que no están ordenados de esta manera, normalmente necesitará un par de líneas adicionales de código para hacer coincidir nuevamente las cosas correctamente (aunque no es complicado).

Observación final: si la aplicación final es demasiado aterradora, vale la pena destacar que un ciclo rápido en un pequeño número de casos (48 dicen) no siempre es desastroso si el trabajo dentro del ciclo está muy bien vectorizado. El "método alternativo" al final de la sección ACTUALIZACIÓN muestra que el paso cumsum-on-groups también se puede lograr pre-preparando una columna de respuesta (inicialmente todas las NA) y luego pasando por los 48 subconjuntos uno por uno y rellenando cada bloque con la cumsum apropiada. Como se señala en el texto, sin embargo, este único paso es aproximadamente la mitad de la velocidad que el enfoque astuto, y si se necesitaran más subconjuntos, esto sería considerablemente más lento.

Si alguien tiene alguna pregunta de seguimiento sobre este tipo de tarea, no dude en darme un grito.

Answer 2

una solución rápida es definir un vector fuera del bucle como:

temp_cumtt=c(rep(0,nrow(pred))) 

y luego usar este:

if (Z31[i,2] == pred[j,5] & Z31[i,3] == pred[j,4] & Z31[i,4] == pred[j,6] & pred[j,1] >= Z31[i,1]){ 
    temp_cumtt[j]=pred[j,7] + pred[j-1,8] 
} 

en lugar de actualizar la columna data.frame directamente.

después de salir del bucle, puede actualizar la columna:

pred$cum_tt = temp_cumtt 

La otra cosa es que hay que tener cuidado al usar j-1 con el índice de j a partir de las 1. En su ejemplo, no conduce a ese problema condicional.

EDIT:

Mira ahora en sus datos de formato, que tienen estas sugerencias.

1) No convierta en clase Date, en cambio consérvelo como un vector de valores.

2) Clasificar los data_frame Z31 acuerdo a la fecha en vectores: Z31=Z31[with(Z31, order(-Date)), ] (Nota en orden descendente, ya que desea comparar pred [, índice Fecha]> = Z31 [, índice Fecha]

3) Usar la primera loop como pred. Primero tome la fecha de pred ->pred[i,1] y trate de ordenar binariamente y encuentre qué índice en Z31 satisface y desde ese índice en adelante vaya a la lista. Si la condición Date satisface, verifique el resto de las condiciones y complete el temp_cumtt[i] como antes.

Esto debe ser rápida ardiente (ya tipo binario sólo en 48 filas de Z31 y se puede comparar el tiempo de ejecución .con la otra solución.

Answer 3

Hagámoslo con data.table, que debería acelerar mucho las cosas.

Z31 <- data.table(Z31,key="Site,Cultivar,Planting") 
pred <- data.table(pred) 

## First, let's create an extra column in `pred` to see the corresponding date from `Z31` 
## Note 1: The JT is necessary since both sets have the same column names 
## Note 2: I needed to use as.integer on Planting to make it work 

pred[,Z31Est.Date:={JT=J(Site,Cultivar,as.integer(Planting)); Z31[JT,Est.Date][[4]]}] 

## Now we can see for each row whether the date in `pred` is higher than or equal to that from `Z31`. 

pred[,DateTrue:=Date>=Z31Est.Date] 

## Finally, we only have to add up `pred[i,tt]` and `pred[i-1,cum_tt]` for each row where `DateTrue` equals `TRUE`. 

for (i in 1:nrow(pred)) set(pred,i,13L,if(pred[i,DateTrue]) pred[i-1,cum_tt]+pred[i,tt] else(0))