As mentioned earlier, code that exists at a certain level of the environment hierarchy has at least read access to all the variables at the levels above it. On the other hand, direct write access to variables at higher levels via the standard <- operator is not possible.
If you wish to write to a global variable—or more generally, to any variable higher in the environment hierarchy than the level at which your write statement exists—you can use the superassignment operator, <<-, or the assign() function. Let’s discuss the superassignment operator first.
Consider the following code:
> two <- function(u) {
+ u <<- 2*u
+ z <- 2*z
+ }
> x <- 1
> z <- 3
> u
Error: object "u" not found
> two(x)
> x
[1] 1
> z
[1] 3
> u
[1] 2Let’s look at the impact (or not) on the three top-level variables x, z, and u:
x: Even though x was the actual argument to two() in the example, it retained the value 1 after the call. This is because its value 1 was copied to the formal argument u, which is treated as a local variable within the function. Thus, when u changed, x did not change with it.
z: The two z values are entirely unrelated to each other—one is top level, and the other is local to two(). The change in the local variable has no effect on the global variable. Of course, having two variables with the same name is probably not good programming practice.
u: The u value did not even exist as a top-level variable prior to our calling two(), hence the “not found” error message. However, it was created as a top-level variable by the superassignment operator within two(), as confirmed after the call.
Though <<- is typically used to write to top-level variables, as in our example, technically, it does something a bit different. Use of this operator to write to a variable w will result in a search up the environment hierarchy, stopping at the first level at which a variable of that name is encountered. If none is found, then the level selected will be global. Look what happens in this little example:
> f
function() {
inc <- function() {x <<- x + 1}
x <- 3
inc()
return(x)
}
> f()
[1] 4
> x
Error: object 'x' not foundHere, inc() is defined within f(). When inc() is executing, and the R interpreter sees a superassignment to x, it starts going up the hierarchy. At the first level up—the environment within f()—it does find an x, and so that x is the one that is written to, not x at the top level.
You can also use the assign() function to write to upper-level variables. Here’s an altered version of the previous example:
> two
function(u) {
assign("u",2*u,pos=.GlobalEnv)
z <- 2*z
}
> two(x)
> x
[1] 1
> u
[1] 2Here, we replaced the superassignment operator with a call to assign(). That call instructs R to assign the value 2*u (this is the local u) to a variable u further up the call stack, specifically in the top-level environment. In this case, that environment is only one call level higher, but if we had a chain of calls, it could be much further up.
The fact that you reference variables using character strings in assign() can come in handy. Recall the example in Chapter 5 concerning analysis of hiring patterns of various large corporations. We wanted to form a subdata frame for each firm, extracted from the overall data frame, all2006. For instance, consider this call:
makecorpdfs(c("MICROSOFT CORPORATION","ms","INTEL CORPORATION","intel","
SUN MICROSYSTEMS, INC.","sun","GOOGLE INC.","google")This would first extract all Microsoft records from the overall data frame, naming the resulting subdata frame ms2006. It would then create intel2006 for Intel, and so on. Here is the code (changed to function form, for clarity):
makecorpdfs <- function(corplist) {
for (i in 1:(length(corplist)/2)) {
corp <- corplist[2*i-1]
newdtf <- paste(corplist[2*i],"2006",sep="")
assign(newdtf,makecorp(corp),pos=.GlobalEnv)
}
}In the iteration i = 1, the code uses paste() to splice together the strings "ms" and "2006", resulting in "ms2006", the desired name.
Discrete-event simulation (DES) is widely used in business, industry, and government. The term discrete event refers to the fact that the state of the system changes only in discrete quantities, rather than changing continuously.
A typical example would involve a queuing system, say people lining up to use an ATM. Let’s define the state of our system at time t to be the number of people in the queue at that time. The state changes only by +1, when someone arrives, or by −1, when a person finishes an ATM transaction. This is in contrast to, for instance, a simulation of weather, in which temperature, barometric pressure, and so on change continuously.
This will be one of the longer, more involved examples in this book. But it exemplifies a number of important issues in R, especially concerning global variables, and will serve as an example when we discuss appropriate use global variables in the next section. Your patience will turn out to be a good investment of time. (It is not assumed here that the reader has any prior background in DES.)
Central to DES operation is maintenance of the event list, which is simply a list of scheduled events. This is a general DES term, so the word list here does not refer to the R data type. In fact, we’ll represent the event list by a data frame.
In the ATM example, for instance, the event list might at some point in the simulation look like this:
customer 1 arrives at time 23.12 customer 2 arrives at time 25.88 customer 3 arrives at time 25.97 customer 1 finishes service at time 26.02
Since the earliest event must always be handled next, the simplest form of coding the event list is to store it in time order, as in the example. (Readers with computer science background might notice that a more efficient approach might be to use some kind of binary tree for storage.) Here, we will implement it as a data frame, with the first row containing the earliest scheduled event, the second row containing the second earliest, and so on.
The main loop of the simulation repeatedly iterates. Each iteration pulls the earliest event off of the event list, updates the simulated time to reflect the occurrence of that event, and reacts to this event. The latter action will typically result in the creation of new events. For example, if a customer arrival occurs when the queue is empty, that customer’s service will begin—one event triggers setting up another. Our code must determine the customer’s service time, and then it will know the time at which service will be finished, which is another event that must be added to the event list.
One of the oldest approaches to writing DES code is the event-oriented paradigm. Here, the code to handle the occurrence of one event directly sets up another event, reflecting our preceding discussion.
As an example to guide your thinking, consider the ATM situation. At time 0, the queue is empty. The simulation code randomly generates the time of the first arrival, say 2.3. At this point, the event list is simply (2.3,“arrival”). This event is pulled off the list, simulated time is updated to 2.3, and we react to the arrival event as follows:
The queue for the ATM is empty, so we start the service by randomly generating the service time—say it is 1.2 time units. Then the completion of service will occur at simulated time 2.3 + 1.2 = 3.5.
We add the completion of service event to the event list, which will now consist of (3.5,“service done”)).
We also generate the time to the next arrival, say 0.6, which means the arrival will occur at time 2.9. Now the event list consists of (2.9,“arrival”) and (3.5,“service done”).
The code consists of a generally applicable library. We also have an example application, which simulates an M/M/1 queue, which is a single-server queue in which both interarrival time and service time are exponentially distributed.
The code in this example is hardly optimal, and the reader is invited to improve it, especially by rewriting some portions in C. (Chapter 15 shows how to interface C to R.) This example does, however, serve to illustrate a number of the issues we have discussed in this chapter.
Here is a summary of the library functions:
schedevnt(): Inserts a newly created event into the event list.
getnextevnt(): Pulls the earliest event off the event list.
dosim(): Includes the core loop of the simulation. Repeatedly calls getnextevnt() to get the earliest of the pending events; updates the current simulated time, sim$currtime, to reflect the occurrence of that event; and calls the application-specific function reactevnt() to process this newly occurred event.
The code uses the following application-specific functions:
initglbls(): Initializes the application-specific global variables.
reactevnt(): Takes the proper actions when an event occurs, typically generating new events as a result.
prntrslts(): Prints the application-specific results of the simulation.
Note that initglbls(), reactevnt(), and prntrslts() are written by the application programmer and then passed to dosim() as arguments. In the M/M/1 queue example included here, these functions are named mm1initglbls(), mm1reactevnt(), and mm1prntrslts(), respectively. Thus, in correspondence with the definition of dosim(),
dosim <- function(initglbls,reactevnt,prntrslts,maxsimtime,apppars=NULL,dbg=FALSE){our call is as follows:
dosim(mm1initglbls,mm1reactevnt,mm1prntrslts,10000.0, list(arrvrate=0.5,srvrate=1.0))
Here’s the library code:
1 # DES.R: R routines for discrete-event simulation (DES)
2
3 # each event will be represented by a data frame row consisting of the
4 # following components: evnttime, the time the event is to occur;
5 # evnttype, a character string for the programmer-defined event type;
6 # optional application-specific components, e.g.
7 # the job's arrival time in a queuing app
8
9 # a global list named "sim" holds the events data frame, evnts, and
10 # current simulated time, currtime; there is also a component dbg, which
11 # indicates debugging mode
12
13 # forms a row for an event of type evntty that will occur at time
14 # evnttm; see comments in schedevnt() regarding appin
15 evntrow <- function(evnttm,evntty,appin=NULL) {
16 rw <- c(list(evnttime=evnttm,evnttype=evntty),appin)
17 return(as.data.frame(rw))
18 }
19
20 # insert event with time evnttm and type evntty into event list;
21 # appin is an optional set of application-specific traits of this event,
22 # specified in the form a list with named components
23 schedevnt <- function(evnttm,evntty,appin=NULL) {
24 newevnt <- evntrow(evnttm,evntty,appin)
25 # if the event list is empty, set it to consist of evnt and return
26 if (is.null(sim$evnts)) {
27 sim$evnts <<- newevnt
28 return()
29 }
30 # otherwise, find insertion point
31 inspt <- binsearch((sim$evnts)$evnttime,evnttm)
32 # now "insert," by reconstructing the data frame; we find what
33 # portion of the current matrix should come before the new event and
34 # what portion should come after it, then string everything together
35 before <-
36 if (inspt == 1) NULL else sim$evnts[1:(inspt-1),]
37 nr <- nrow(sim$evnts)
38 after <- if (inspt <= nr) sim$evnts[inspt:nr,] else NULL
39 sim$evnts <<- rbind(before,newevnt,after)
40 }
41
42 # binary search of insertion point of y in the sorted vector x; returns
43 # the position in x before which y should be inserted, with the value
44 # length(x)+1 if y is larger than x[length(x)]; could be changed to C
45 # code for efficiency
46 binsearch <- function(x,y) {
47 n <- length(x)
48 lo <- 1
49 hi <- n
50 while(lo+1 < hi) {
51 mid <- floor((lo+hi)/2)
52 if (y == x[mid]) return(mid)
53 if (y < x[mid]) hi <- mid else lo <- mid
54 }
55 if (y <= x[lo]) return(lo)
56 if (y < x[hi]) return(hi)
57 return(hi+1)
58 }
59
60 # start to process next event (second half done by application
61 # programmer via call to reactevnt())
62 getnextevnt <- function() {
63 head <- sim$evnts[1,]
64 # delete head
65 if (nrow(sim$evnts) == 1) {
66 sim$evnts <<- NULL
67 } else sim$evnts <<- sim$evnts[-1,]
68 return(head)
69 }
70
71 # simulation body
72 # arguments:
73 # initglbls: application-specific initialization function; inits
74 # globals to statistical totals for the app, etc.; records apppars
75 # in globals; schedules the first event
76 # reactevnt: application-specific event handling function, coding the
77 # proper action for each type of event
78 # prntrslts: prints application-specific results, e.g. mean queue
79 # wait
80 # apppars: list of application-specific parameters, e.g.
81 # number of servers in a queuing app
82 # maxsimtime: simulation will be run until this simulated time
83 # dbg: debug flag; if TRUE, sim will be printed after each event
84 dosim <- function(initglbls,reactevnt,prntrslts,maxsimtime,apppars=NULL,
85 dbg=FALSE) {
86 sim <<- list()
87 sim$currtime <<- 0.0 # current simulated time
88 sim$evnts <<- NULL # events data frame
89 sim$dbg <<- dbg
90 initglbls(apppars)
91 while(sim$currtime < maxsimtime) {
92 head <- getnextevnt()
93 sim$currtime <<- head$evnttime # update current simulated time
94 reactevnt(head) # process this event
95 if (dbg) print(sim)
96 }
97 prntrslts()
98 }The following is an example application of the code. Again, the simulation models an M/M/1 queue, which is a single-server queuing system in which service times and times between job arrivals are exponentially distributed.
1 #DES application: M/M/1 queue, arrival rate 0.5, service rate 1.0
2
3 # the call
4 # dosim(mm1initglbls,mm1reactevnt,mm1prntrslts,10000.0,
5 # list(arrvrate=0.5,srvrate=1.0))
6 # should return a value of about 2 (may take a while)
7
8 # initializes global variables specific to this app
9 mm1initglbls <- function(apppars) {
10 mm1glbls <<- list()
11 # simulation parameters
12 mm1glbls$arrvrate <<- apppars$arrvrate
13 mm1glbls$srvrate <<- apppars$srvrate
14 # server queue, consisting of arrival times of queued jobs
15 mm1glbls$srvq <<- vector(length=0)
16 # statistics
17 mm1glbls$njobsdone <<- 0 # jobs done so far
18 mm1glbls$totwait <<- 0.0 # total wait time so far
19 # set up first event, an arrival; the application-specific data for
20 # each event will consist of its arrival time, which we need to
21 # record in order to later calculate the job's residence time in the
22 # system
23 arrvtime <- rexp(1,mm1glbls$arrvrate)
24 schedevnt(arrvtime,"arrv",list(arrvtime=arrvtime))
25 }
26
27 # application-specific event processing function called by dosim()
28 # in the general DES library
29 mm1reactevnt <- function(head) {
30 if (head$evnttype == "arrv") { # arrival
31 # if server free, start service, else add to queue (added to queue
32 # even if empty, for convenience)
33 if (length(mm1glbls$srvq) == 0) {
34 mm1glbls$srvq <<- head$arrvtime
35 srvdonetime <- sim$currtime + rexp(1,mm1glbls$srvrate)
36 schedevnt(srvdonetime,"srvdone",list(arrvtime=head$arrvtime))
37 } else mm1glbls$srvq <<- c(mm1glbls$srvq,head$arrvtime)
38 # generate next arrival
39 arrvtime <- sim$currtime + rexp(1,mm1glbls$arrvrate)
40 schedevnt(arrvtime,"arrv",list(arrvtime=arrvtime))
41 } else { # service done
42 # process job that just finished
43 # do accounting
44 mm1glbls$njobsdone <<- mm1glbls$njobsdone + 1
45 mm1glbls$totwait <<-
46 mm1glbls$totwait + sim$currtime - head$arrvtime
47 # remove from queue
48 mm1glbls$srvq <<- mm1glbls$srvq[-1]
49 # more still in the queue?
50 if (length(mm1glbls$srvq) > 0) {
51 # schedule new service
52 srvdonetime <- sim$currtime + rexp(1,mm1glbls$srvrate)
53 schedevnt(srvdonetime,"srvdone",list(arrvtime=mm1glbls$srvq[1]))
54 }
55 }
56 }
57
58 mm1prntrslts <- function() {
59 print("mean wait:")
60 print(mm1glbls$totwait/mm1glbls$njobsdone)
61 }To see how all this works, take a look at the M/M/1 application code. There, we have set up a global variable, mm1glbls, which contains variables relevant to the M/M/1 code, such as mm1glbls$totwait, the running total of the wait time of all jobs simulated so far. As you can see, the superassignment operator is used to write to such variables, as in this statement:
mm1glbls$srvq <<- mm1glbls$srvq[-1]
Let’s look at mm1reactevnt() to see how the simulation works, focusing on the code portion in which a “service done” event is handled.
} else { # service done
# process job that just finished
# do accounting
mm1glbls$njobsdone <<- mm1glbls$njobsdone + 1
mm1glbls$totwait <<-
mm1glbls$totwait + sim$currtime - head$arrvtime
# remove this job from queue
mm1glbls$srvq <<- mm1glbls$srvq[-1]
# more still in the queue?
if (length(mm1glbls$srvq) > 0) {
# schedule new service
srvdonetime <- sim$currtime + rexp(1,mm1glbls$srvrate)
schedevnt(srvdonetime,"srvdone",list(arrvtime=mm1glbls$srvq[1]))
}
}First, this code does some bookkeeping, updating the totals of number of jobs completed and wait time. It then removes this newly completed job from the server queue. Finally, it checks if there are still jobs in the queue and, if so, calls schedevnt() to arrange for the service of the one at the head.
What about the DES library code itself? First note that the simulation state, consisting of the current simulated time and the event list, has been placed in an R list structure, sim. This was done in order to encapsulate all the main information into one package, which in R, typically means using a list. The sim list has been made a global variable.
As mentioned, a key issue in writing a DES library is the event list. This code implements it as a data frame, sim$evnts. Each row of the data frame corresponds to one scheduled event, with information about the event time, a character string representing the event type (say arrival or service completion), and any application-specific data the programmer wishes to add. Since each row consists of both numeric and character data, it was natural to choose a data frame for representing this event list. The rows of the data frame are in ascending order of event time, which is contained in the first column.
The main loop of the simulation is in dosim() of the DES library code, beginning at line 91:
while(sim$currtime < maxsimtime) {
head <- getnextevnt()
sim$currtime <<- head$evnttime # update current simulated time
reactevnt(head) # process this event
if (dbg) print(sim)
}First, getnextevnt() is called to remove the head (the earliest event) from the event list. (Note the side effect: The event list changes.) Then the current simulated time is updated according to the scheduled time in the head event. Finally, the programmer-supplied function reactevnt() is called to process the event (as seen in the M/M/1 code discussed earlier).
The main potential advantage of using a data frame as our structure here is that it enables us to maintain the event list in ascending order by time via a binary search operation by event time. This is done in line 31 within schedevnt(), the function that inserts a newly created event into the event list:
inspt <- binsearch((sim$evnts)$evnttime,evnttm)
Here, we wish to insert a newly created event into the event list, and the fact that we are working with a vector enables the use of a fast binary search. (As noted in the comments in the code, though, this really should be implemented in C for good performance.)
A later line in schedevnt()is a good example of the use of rbind():
sim$evnts <<- rbind(before,newevnt,after)
Now, we have extracted the events in the event list whose times are earlier than that of evnt and stored them in before. We also constructed a similar set in after for the events whose times are later than that of newevnt. We then use rbind() to put all these together in the proper order.
Use of global variables is a subject of controversy in the programming community. Obviously, the question raised by the title of this section cannot be answered in any formulaic way, as it is a matter of personal taste and style. Nevertheless, most programmers would probably consider the outright banning of global variables, which is encouraged by many teachers of programming, to be overly rigid. In this section, we will explore the possible value of globals in the context of the structures of R. Here, the term global variable, or just global, will be used to include any variable located higher in the environment hierarchy than the level of the given code of interest.
The use of global variables in R is more common than you may have guessed. You might be surprised to learn that R itself makes very substantial use of globals internally, both in its C code and in its R routines. The superassignment operator <<-, for instance, is used in many of the R library functions (albeit typically in writing to a variable just one level up in the environment hierarchy). Threaded code and GPU code, which are used for writing fast programs (as described in Chapter 16), tend to make heavy use of global variables, which provide the main avenue of communication between parallel actors.
Now, to make our discussion concrete, let’s return to the earlier example from Section 7.7:
f <- function(lxxyy) { # lxxyy is a list containing x and y
...
lxxyy$x <- ...
lxxyy$y <- ...
return(lxxyy)
}
# set x and y
lxy$x <- ...
lxy$y <- ...
lxy <- f(lxy)
# use new x and y
... <- lxy$x
... <- lxy$yAs noted earlier, this code might be a bit unwieldy, especially if x and y are themselves lists.
By contrast, here is an alternate pattern that uses globals:
f <- function() {
...
x <<- ...
y <<- ...
}
# set x and y
x <- ...
y <- ...
f() # x and y are changed in here
# use new x and y
... <- x
... <- yArguably, this second version is much cleaner, being less cluttered and not requiring manipulation of lists. Cleaner code is usually easier to write, debug, and maintain.
It is for these reasons—avoiding clutter and simplifying the code—that we chose to use globals, rather than to return lists, in the DES code earlier in this chapter. Let’s explore that example further.
We had two global variables (both lists, encapsulating various information): sim, associated with the library code, and mm1glbls, associated with our M/M/1 application code. Let’s consider sim first.
Even many programmers who have reservations about using globals agree that such variables may be justified if they are truly global, in the sense that they are used broadly in the program. This is the case for sim in our DES example. It is used both in the library code (in schedevnt(), getnextevnt(), and dosim()) and in in our M/M/1 application code (in mm1reactevnt()). The latter access to sim is on a read-only basis in this particular instance, but it could involve writes in some applications. A common example of such writes is when an event needs to be canceled. This might arise in modeling a “whichever comes first” situation; two events are scheduled, and when one of them occurs, the other must be canceled.
So, using sim as a global seems justified. Nevertheless, if we were bound and determined to avoid using globals, we could have placed sim as a local within dosim(). This function would then pass sim as an argument to all of the functions mentioned in the previous paragraph (schedevnt(), getnextevnt(), and so on), and each of these functions would return the modified sim. Line 94 for example, would change from this:
reactevnt(head)
to this:
sim <- reactevnt(head)
We would then need to add a line like the following to our application-specific function mm1reactevnt():
return(sim)
We could do something similar with mm1glbls, placing a variable called, say, appvars as a local within dosim(). However, if we did this with sim as well, we would need to place them together in a list so that both would be returned, as in our earlier example function f(). We would then have the lists-within-lists clutter described earlier—well, lists within lists within lists in this case.
On the other hand, critics of the use of global variables counter that the simplicity of the code comes at a price. They worry that it may be difficult during debugging to track down locations at which a global variable changes value, since such a change could occur anywhere in the program. This seems to be less of a concern in view of our modern text editors and integrated development tools (the original article calling for avoiding use of globals was published in 1970!), which can be used to find all instances of a variable. However, it should be taken into consideration.
Another concern raised by critics involves situations in which a function is called in several unrelated parts of the overall program using different values. For example, consider using our example f() function in different parts of our program, each call with its own values of x and y, rather than just a single value of each, as assumed earlier. This could be solved by setting up vectors of x and y values, with one element for each instance of f() in your program. You would lose some of the simplicity of using globals, though.
The above issues apply generally, not just to R. However, for R there is an additional concern for globals at the top level, as the user will typically have lots of variables there. The danger is that code that uses globals may accidentally overwrite an unrelated variable with the same name.
This can be avoided easily, of course, by simply choosing long, very application-specific names for globals in your code. But a compromise is also available in the form of environments, such as the following for the DES example above.
sim <<- list()
would be replaced by
assign("simenv",new.env(),envir=.GlobalEnv)This would create a new environment, pointed to by simenv at the top level. It would serve as a package in which to encapsulate our globals. We would access them via get() and assign(). For instance, the lines
if (is.null(sim$evnts)) {
sim$evnts <<- newevntin schedevnt() would become
if (is.null(get("evnts",envir=simenv))) {
assign("evnts",newevnt,envir=simenv)Yes, this is cluttered too, but at least it is not complex like lists of lists of lists. And it does protect against unwittingly writing to an unrelated variable the user has at the top level. Using the superassignment operator still yields the least cluttered code, but this compromise is worth considering.
As usual, there is no single style of programming that produces the best solution in all applications. The globals approach is another option to consider for your programming tool kit.
Recall that an R closure consists of a function’s arguments and body together with its environment at the time of the call. The fact that the environment is included is exploited in a type of programming that uses a feature also known (in a slight overloading of terminology) as a closure.
A closure consists of a function that sets up a local variable and then creates another function that accesses that variable. This is a very abstract description, so let’s go right to an example.[2]
1 > counter
2 function () {
3 ctr <- 0
4 f <- function() {
5 ctr <<- ctr + 1
6 cat("this count currently has value",ctr,"\n")
7 }
8 return(f)
9 }Let’s try this out before going into the internal details:
> c1 <- counter()
> c2 <- counter()
> c1
function() {
ctr <<- ctr + 1
cat("this count currently has value",ctr,"\n")
}
<environment: 0x8d445c0>
> c2
function() {
ctr <<- ctr + 1
cat("this count currently has value",ctr,"\n")
}
<environment: 0x8d447d4>
> c1()
this count currently has value 1
> c1()
this count currently has value 2
> c2()
this count currently has value 1
> c2()
this count currently has value 2
> c2()
this count currently has value 3
> c1()
this count currently has value 3Here, we called counter() twice, assigning the results to c1 and c2. As expected, those two variables will consist of functions, specifically copies of f().
However, f() accesses a variable ctr through the superassignment operator, and that variable will be the one of that name that is local to counter(), as it is the first one up the environment hierarchy. It is part of the environment of f() and, as such, is packaged in what is returned to the caller of counter().
The key point is that each time counter() is called, the variable ctr will be in a different environment (in the example, the environments were at memory addresses 0x8d445c0 and 0x8d447d4). In other words, different calls to counter() will produce physically different ctrs.
The result, then, is that our functions c1() and c2() serve as completely independent counters, as seen in the example, where we invoke each of them a few times.
[2] Adapted from an example in “Top-level Task Callbacks in R,” by Duncan Temple Lang (2001), http://developer.r-project.org/TaskHandlers.pdf.