Een tutorial over lussen in R - gebruik en alternatieven

Uitleg van for-lussen

Deze lusstructuur, bestaande uit het rechthoekige vak ‘init’ (of initialisatie), de ruitvormige beslissing en het rechthoekige vak i1, wordt een bekend aantal keren uitgevoerd.

In termen van stroomdiagrammen betekenen rechthoekige vakken zoiets als “doe iets dat geen beslissingen impliceert”. Ruiten of diamanten heten “beslissingssymbolen” en vertalen zich dus in vragen met slechts twee mogelijke logische antwoorden, namelijk True (T) of False (F).

Let op dat, om het simpel te houden, andere mogelijke symbolen uit de figuur zijn weggelaten.

Een of meer instructies binnen het initialisatievak worden gevolgd door de evaluatie van de voorwaarde op een variabele die waarden kan aannemen binnen een opgegeven reeks. In de figuur wordt dit voorgesteld door de ruit: de symbolen betekenen “behoort de huidige waarde van variabele v tot de reeks seq?”

Met andere woorden, je test of de huidige waarde van v binnen een opgegeven bereik ligt. Je definieert dit bereik normaal in de initialisatie, met iets als 1:100 om te zorgen dat de lus start.

Als niet aan de voorwaarde wordt voldaan en de uitkomst False is, wordt de lus nooit uitgevoerd. Dit wordt aangegeven door de losse pijl rechts van de for-lusstructuur. Het programma voert dan de eerste instructie uit na het lusblok.

Als aan de voorwaarde is voldaan, wordt een instructie—of een blok instructies—i1 uitgevoerd. En misschien is dit blok instructies weer een lus. In zulke gevallen spreek je van een geneste lus.

Daarna wordt de voorwaarde opnieuw geëvalueerd. Dit wordt aangegeven door de lijnen die van i1 terug naar boven lopen, direct na het initialisatievak. In R—en in Python—kun je dit in gewoon Engels uitdrukken door te vragen of onze variabele tot een bereik van waarden behoort of niet.

Let op dat in andere talen, bijvoorbeeld in C, de voorwaarde explicieter wordt gemaakt met behulp van een logische operator, zoals groter dan, kleiner dan, gelijk aan, …

Hier is een voorbeeld van een simpele for-lus:

# Create a vector filled with random normal values
u1 <- rnorm(30)
print("This loop calculates the square of the first 10 elements of vector u1")

# Initialize `usq`
usq <- 0. 

for (i in 1:10) {
    # i-th element of `u1` squared into `i`-th position of `usq`
    usq[i] <- u1[i] * u1[i]
    print(usq[i])
}

print(i)

Het for-blok staat tussen accolades. Die kun je direct na de testvoorwaarde plaatsen of eronder, bij voorkeur gevolgd door een inspringing. Niets hiervan is verplicht, maar de accolades verhogen zeker de leesbaarheid van je code en maken het makkelijk om het lusblok en eventuele fouten erin te spotten.

Let op dat de vector met de kwadraten, usq, is geïnitialiseerd. Dit zou niet nodig zijn in gewone RStudio-code, maar in de opgemaakte versie zou knitr niet compileren omdat er vóór gebruik in de lus geen verwijzing naar de vector wordt gevonden, wat in RStudio een fout oplevert. Ga voor meer informatie over knitr naar de R Markdown-pagina.

For-lussen nesten

Nu je weet dat for-lussen ook genest kunnen worden, vraag je je waarschijnlijk af waarom en wanneer je dit in je code zou gebruiken.

Stel dat je een tweedimensionale array wilt manipuleren door de elementen op specifieke waarden te zetten.

Dan kun je zoiets doen als dit:

# Insert your own integer here
my_int <- 42

nr <- as.integer(my_int)

# Create a `n` x `n` matrix with zeros
mymat <- matrix(0, nr, nr)

# For each row and for each column, assign values based on position
# These values are the product of two indexes
for (i in 1:dim(mymat)[1]) {
    for (j in 1:dim(mymat)[2]) {
        mymat[i, j] = i * j
    }
}

Tip: voor meer informatie over de functie matrix(), bezoek R Documentation.

Je hebt twee geneste for-lussen in de code hierboven en dus twee sets accolades, elk met een eigen blok en gestuurd door een eigen index. Dat wil zeggen, i loopt over de rijen en j loopt over de kolommen.

 

Wat is er geproduceerd?

# Show the first 10x10 chunk or the first `nr` x `nr` chunk
if (nr > 10) {
    mymat[1:10, 1:10]
} else
    mymat[1:nr, 1:nr]

Let op dat je, om te voorkomen dat de gebruiker het scherm overspoelt met enorme tabellen, aan het eind een voorwaarde zet om alleen de eerste 10 x 10 sectie af te drukken als de gebruiker om een geheel getal groter dan 10 heeft gevraagd. Anders wordt een n x n-gedeelte afgedrukt.

De volledige code ziet er zo uit:

# Insert your own integer here
my_int <- 42

nr <- as.integer(my_int)

# Create a `n` x `n` matrix with zeroes
mymat <- matrix(0, nr, nr)

# For each row and for each column, assign values based on position
# These values are the product of two indexes
for (i in 1:dim(mymat)[1]) {
    for (j in 1:dim(mymat)[2]) {
        mymat[i, j] = i * j
    }
}

# Show the first 10x10 chunk or the first `nr` x `nr` chunk
if (nr > 10) {
    mymat[1:10, 1:10]
} else
    mymat[1:nr, 1:nr]

For-lussen: populair maar niet altijd perfect voor jouw gebruik

De for-lus is veruit het populairst, en de constructie impliceert dat het aantal iteraties vaststaat en vooraf bekend is, zoals in gevallen als “genereer de eerste 200 priemgetallen” of “zet de 10 belangrijkste klanten op een rij”.

Maar wat als je het aantal iteraties niet kent of niet stuurt, en er een of meer voorwaarden kunnen optreden die je niet vooraf kunt voorspellen?

Je wilt bijvoorbeeld het aantal klanten tellen die wonen in een gebied met een bepaalde postcode, of het aantal kliks op een banner op een webpagina in de laatste twee dagen, of soortgelijke onvoorziene gebeurtenissen.

In zulke gevallen kunnen de while-lus en diens familielid repeat uitkomst bieden…

Onderbreek je lussen met break

Wanneer de R-interpreter een break tegenkomt, geeft hij de controle door aan de instructie direct na het einde van de lus (als die er is). In het geval van geneste lussen zorgt break er alleen voor dat je de binnenste lus verlaat.

Hier is een voorbeeld.

Deze code definieert een m x n-matrix van nullen en gaat dan een geneste for-lus in om de posities van de matrix te vullen, maar alleen als de twee indexen verschillen. Het doel is een onderdriehoeksmatrix te maken, dat is een matrix waarvan de elementen onder de hoofddiagonaal ongelijk aan nul zijn. De andere blijven onaangeroerd op hun geïnitialiseerde nulwaarde.

Wanneer de indexen gelijk zijn en dus aan de voorwaarde in de binnenste lus, die loopt over de kolomindex j, wordt voldaan, wordt een break uitgevoerd en wordt de binnenste lus onderbroken met een directe sprong naar de instructie na de binnenste lus. Deze instructie is een print()-instructie. Vervolgens gaat de controle naar de buitenste for-voorwaarde (over de rijen, index i), die opnieuw wordt geëvalueerd.

Als de indexen verschillen, wordt de toekenning uitgevoerd en wordt de teller met 1 verhoogd. Aan het eind drukt het programma de teller ctr af, die het aantal toegewezen elementen bevat.

# Make a lower triangular matrix (zeroes in upper right corner)
m = 10
n = 10

# A counter to count the assignment
ctr = 0

# Create a 10 x 10 matrix with zeroes
mymat = matrix(0, m, n)

for (i in 1:m) {
    for (j in 1:n) {
        if (i == j) {
            break
        } else {
            # you assign the values only when i<>j
            mymat[i, j] = i * j
            ctr = ctr + 1
        }
    }
    print(i * j)
}

# Print how many matrix cells were assigned
print(ctr)

Let op dat het misschien wat overdreven voorzichtig is om ook accolades te zetten wanneer ze niet strikt noodzakelijk zijn. Meestal doe je dit om zeker te zijn dat alles wat met een { wordt geopend, ook met een } wordt gesloten. Dus als je een ongelijk aantal { of } ziet, weet je dat er een fout is, al geldt het omgekeerde niet per se!

Het gebruik van next in lussen

next stopt een specifieke iteratie en springt naar de volgende cyclus. In feite springt het naar de evaluatie van de voorwaarde van de huidige lus.

In andere talen vind je misschien het (licht verwarrende) equivalent “continue”, dat hetzelfde betekent: waar je je ook bevindt, spring bij verificatie van de voorwaarde naar de evaluatie van de lus.

Een eenvoudiger voorbeeld om de lus gaande te houden terwijl je een specifieke cyclus overslaat bij het optreden van een voorwaarde is:

m = 20

for (k in 1:m) {
    if (!k %% 2)
        next
    print(k)
}

Deze code drukt alle oneven getallen af binnen het interval 1:m (hier m=20). Met andere woorden, alle gehele getallen behalve die met een niet-nul rest bij deling door 2 (vandaar de modulus-operator %%), zoals gespecificeerd door de if-test, worden afgedrukt.

Getallen met rest nul worden niet afgedrukt, omdat het programma springt naar de evaluatie van de voorwaarde i in 1:m en eventuele volgende instructies negeert. In dit geval wordt print(k) genegeerd.

v1 <- c(3,4,5)
v2 <- c(3,4,6)
v3 <- c(0,0,0)
n = length(v1)

for (i in 1:n) {
    v3[i] <- v1[i] + v2[i]
}
v3

v3 = v1 + v2
v3

A <- matrix(1:6, nrow = 3, ncol = 2)
B <- matrix(c(2,4,3,1,5,7), nrow=3, ncol=2)

C= A + B
C

Vectorisatie in R uitgelegd

Waarom zou vectorisatie sneller zijn, gegeven dat het aantal elementaire bewerkingen ogenschijnlijk hetzelfde is?

Dit leg je het best uit door naar de interne werking van R te kijken, wat een apart artikel zou vergen, maar kort samengevat: R is om te beginnen een geïnterpreteerde taal en als zodanig worden alle details over variabeledefinitie afgehandeld door de interpreter. Je hoeft bijvoorbeeld niet te specificeren dat een getal een drijvende-kommawaarde is of geheugen te alloceren met een pointer.

De R-interpreter “begrijpt” deze zaken uit de context terwijl je je commando’s invoert, maar doet dit commando voor commando. Daarom moet hij met zulke definities omgaan elke keer dat je een commando geeft, zelfs als je het herhaalt.

Een compiler daarentegen lost letterlijk alle definities en declaraties op compileertijd op over de volledige code; die wordt vertaald naar compacte en geoptimaliseerde binaire code, vóór je iets uitvoert. Nu, omdat R-functies in een van deze lagere talen zijn geschreven, zijn ze efficiënter.

In de praktijk zou je op laag niveau aanroepen van C of C++ zien, meestal geïmplementeerd in wat een wrapper-code heet.

Ten tweede, in talen die vectorisatie ondersteunen (zoals R of Matlab), werkt elke instructie die een numeriek gegeven gebruikt op een object dat native is gedefinieerd als een vector, zelfs als het maar uit één element bestaat. Dit is de standaard wanneer je bijvoorbeeld een enkele numerieke variabele definieert: de interne representatie in R is altijd een vector, ook al bestaat die uit slechts één getal.

De lussen blijven onder de motorkap bestaan, maar dan op het lagere en veel snellere gecompileerde C/C++-niveau. Het voordeel van een vector is dat de definities door de interpreter slechts één keer worden opgelost, op de hele vector, ongeacht de grootte, in tegenstelling tot een lus op een scalaire waarde, waar definities en allocaties element-voor-element moeten worden gedaan, en dat is trager.

Tot slot maakt werken met het native vectorformaat het mogelijk om zeer efficiënte lineaire-algebaroutines te gebruiken (zoals BLAS of Basic Linear Algebra Subprograms), zodat R bij het uitvoeren van gevectoriseerde instructies op deze efficiënte numerieke routines leunt. De boodschap is dus: verwerk indien mogelijk hele datastructuren binnen één functieaanroep om alle kopieerbewerkingen te vermijden.

Maar genoeg uitweidingen, laten we een algemeen voorbeeld van vectorisatie maken, en in de volgende subsectie duik je in specifiekere en populaire gevectoriseerde R-functies ter vervanging van lussen.

Een voorbeeld van vectorisatie in R

Terug naar het idee van “data laten groeien”, een typisch inefficiënte manier om een data frame te “updaten”.

Eerst maak je een m x n-matrix met replicate(m, rnorm(n)) met m=10 kolomvectors van n=10 elementen elk, geconstrueerd met rnorm(n), dat willekeurige normale getallen maakt.

Vervolgens zet je die om in een data frame (dus 10 observaties van 10 variabelen) en voer je een algebraïsche bewerking uit op elk element met een geneste for-lus: bij elke iteratie verhoogt een sinusoïdale functie elk element waarnaar door de twee indexen wordt verwezen.

Het volgende voorbeeld is wat kunstmatig, maar het kan de toevoeging van een signaal aan wat willekeurige ruis voorstellen:

# This is a bad loop with 'growing' data
set.seed(42)
m = 10
n = 10

# Create matrix of normal random numbers
mymat <- replicate(m, rnorm(n))

# Transform into data frame
mydframe <- data.frame(mymat)

for (i in 1:m) {
  for (j in 1:n) {
    mydframe[i, j] <- mydframe[i, j] + 10 * sin(0.75 * pi)
    print(mydframe)
  }
}

Hier bestaat het grootste deel van de uitvoeringstijd uit kopiëren en het beheren van de lus.

Kijk eens hoe een gevectoriseerde oplossing eruitziet:

set.seed(42)
m = 10
n = 10
mymat <- replicate(m, rnorm(n))
mydframe <- data.frame(mymat)
mydframe <- mydframe + 10 * sin(0.75 * pi)
mydframe

Dit oogt eenvoudiger: de laatste regel vervangt de geneste for-lus. Let op het gebruik van set.seed() om te zorgen dat de twee implementaties exact hetzelfde resultaat geven.

Laten we nu de uitvoeringstijd voor de twee oplossingen kwantificeren.

Dat kan met R-systeemcommando’s, zoals system.time() waaraan je een codeblok kunt doorgeven zoals dit:

Tip: zet de code die je wilt evalueren tussen de haakjes van de functie system.time().

# Insert `system.time()` to measure loop execution
for (i in 1:m) {
  for (j in 1:n) {
    mydframe[i, j] <- mydframe[i, j] + 10 * sin(0.75 * pi)
  }
}

# Add `system.time()` to measure vectorized execution
mydframe <- mydframe + 10 * sin(0.75 * pi) 

In de bovenstaande code doe je het kiezen van m en n, het maken van de matrix en de omzetting naar een data frame slechts één keer aan het begin, en evalueer je het for-blok tegenover de “one-liner” van de gevectoriseerde versie met de twee afzonderlijke aanroepen van system.time().

Je ziet dat al met een minimale instelling van m=n=10 de gevectoriseerde versie 7 keer sneller is, hoewel dat bij zulke lage waarden voor de gebruiker nauwelijks uitmaakt.

Verschillen worden merkbaar (op menselijke schaal) als je m=n=100 zet, terwijl opschalen naar 1000 de for-lus doet lijken te hangen gedurende vele tientallen seconden, terwijl de gevectoriseerde vorm het nog steeds in een oogwenk doet.

Voor m=n=10000 hangt de for-lus meer dan een minuut, terwijl de gevectoriseerde 2,54 sec nodig heeft. Natuurlijk moet je deze metingen met een korrel zout nemen en hangen ze af van de hardware- en softwareconfiguratie; vermijd bij voorkeur het overladen van je laptop met tientallen open browsertabbladen en meerdere toepassingen op de achtergrond; maar de metingen illustreren de verschillen.

Eerlijkheidshalve: het vergroten van m en n beïnvloedt ook de matrixgeneratie sterk, zoals je gemakkelijk ziet door nog een system.time()-aanroep rond de replicatiefunctie te plaatsen.

Speel gerust met m en n om te zien hoe de uitvoeringstijd verandert, door de uitvoeringstijd te plotten als functie van het product m x n. Dit is de relevante indicator, omdat die de afmeting van de gemaakte matrices uitdrukt. Het kwantificeert dus ook het aantal iteraties dat nodig is om de taak via de nset-ed for-lus te voltooien.

Dit is dus een voorbeeld van vectorisatie. Maar er zijn er veel meer. In R News, een nieuwsbrief voor het R-project op pagina 46, staan zeer efficiënte functies voor het berekenen van sommen en gemiddelden voor bepaalde dimensies in arrays of matrices, zoals: rowSums(), colSums(), rowMeans() en colMeans().

Verder vermeldt de nieuwsbrief ook dat “...de functies in de ‘apply’-familie, genaamd [s,l,m,t]apply, bestaan om een andere functie toe te passen op de elementen/dimensies van objecten. Deze ‘apply’-functies bieden een compacte syntaxis voor soms nogal complexe taken die leesbaarder en sneller is dan slecht geschreven lussen.”