Social Simulations

Real-World Computing

by Louis J. D'Andria

Introduction

The last few columns have discussed particular social and economic models in detail, with little attention given to real world programming concerns. For example, many simulations that you may want to run have memory requirements that exceed the memory capacity of your or any computer. This column focuses on these considerations, and offers solutions – in the section New Nesting Functions – which can be useful for programs having nothing to do with simulations. To illustrate these techniques, we look at the following new model.

Many social and economic situations can be viewed as bilateral interactions between two individuals. However, other social phenomena can be better described in terms of interactions between an individual and a group of other people. (This is not so-called group-level interaction between clusters of people but rather interactions between an individual and two or more other individuals). In this column, neighborhood formation resulting from individual-group interaction is discussed. This can lead to voluntary segregation and local conformism which can result in the development of social norms and the spread of fads and fashions.

Forming Neighborhoods

A model for self-forming neighborhoods based on the desire of people to live with their own kind was proposed by Schelling [1978, 147–53] in one of the first applications of computer simulation in social science. The Schelling model, which allows an individual who is not happy with the number of nearest neighbors who are like him, to move to the nearest empty site that has a sufficient number of similar nearest neighbors, demonstrated that spatial segregation, or ghettoization, occurs spontaneously, without being imposed by a central authority. It can result in the clustering of people by gender or age at a social gathering or in the clustering of people by ethnicity or race in society at large.

Because the Schelling model considers the movement of a person from where he is to be a result of his unhappiness with his neighbors, it is sometimes viewed as a model of an undesirable social phenomenon (fleeing). However, segregation, when it is voluntary, can have either positive or negative social and/or economic consequences [Durlauf 1997].

Here we recast the Schelling model so that an individual can decide to move, rather than stay where he is because of a desire to be with individuals who have similar attributes.

The System

Our model uses a square n by n lattice with wraparound boundary conditions. There is a population density p of individuals occupying lattice sites and the remaining sites are empty. The system evolves over a given number of time steps t.

Populating Society

The value of an empty site is 0.

The value of a site that is occupied by a person is an ordered pair. The first component of the ordered pair is a randomly chosen integer between 1 and 4 and indicates the direction the individual is facing (north, east, south, or west, respectively). The second component of the ordered pair is an attribute list of v elements, each of which has an integer value between 1 and w. The list elements represent various attributes of the individual and may include unchangeable traits such as race, gender, or ethnic identity and/or changeable beliefs such as political views, moral values, or personal interests.

Initially, the people are randomly placed on the lattice and face randomly chosen directions.

The code for this configuration is given by

RND := Random[Integer, {1, 4}]
society = 
      Table[Floor[p + Random[]], {n}, {n}] /. 
        1 :> {RND, 
            Table[Random[Integer, {1, w}], {v}]}

Executing a Time Step

During each time step, the following processes occur consecutively:

Deciding Whether to Stay or Move On

In this partial step, each individual decides whether he wants to stay in his present location or move. The criterion for this decision is based on the number of nearest neighbor sites containing individuals who share at least 50% of their attribute values with the individual.

For example, if the value of an individual is {2, {2, 2, 2, 1, 2}} and the eight sites in its Moore neighborhood have values {4, {1, 2, 2, 1, 2}}, {2, {1, 1, 1, 2, 1}}, 0, 0, {2, {1, 1, 2, 2, 1}}, {3, {2, 2, 2, 2, 1}}, 0, and {1, {1, 1, 2, 1, 1}}, then the individual shares 80%, 0%, 0%, 0%, 20%, 60%, 0%, and 40% of his attribute values with his nearest neighbor sites (an empty site has no attribute list and therefore shares no attribute values with a person). The individual shares more than 50% of his attribute values with only two of his eight nearest neighbor sites, or 25%, so this individual would decide to move.

The movestay rule for an empty site is simply to remain where it is.

movestay[0, __] := 0

The left-hand side of the movestay for a site occupied by an individual can be written as

movestay[{a_, b_}, res__]

where b is the attribute list of the individual and res is the sequence consisting of the values of the eight nearest neighbor sites in the individual's Moore neighborhood.

For a site occupied by an individual, we must determine if a majority of the nearest neighbor sites share at least 50% of their attribute values with that individual. This is done in several steps.

step1  = 
      Map[Count[b - #[[2]], 0] &, {res} /. 0 -> {0, 0}]

Note: The value of an empty site is changed from 0 to {0, 0} during this calculation so that when its second component is used in (b - #[[2]]), a value of b is returned, and since b has no zero-valued elements, Count[b - #[[2]], 0]&[empty site] returns 0 indicating that the empty site has no shared attribute values with the person.

step2  =  Count[step1, _?(# < v/2 &)]/8

If this quantity is more than 50%, then the majority of the person's neighboring sites have sufficiently different attribute lists, and he will choose to move; otherwise, he will remain in place. This is done by rounding this percentage to 0 or 1, and multiplying the result by the person's direction. The walk rules will leave a person with zero direction stationary.

{a*Round[step2], b}

This returns {0, b} if the site is occupied by an individual with attribute list b who wants to stay where he is because at least 50% of his neighbor sites are occupied by individuals sharing at least 50% of their attribute values with the attribute values of the individual, or {a, b}, where a = 1, 2, 3, or 4 if the site is occupied by an individual who wants to move because not enough of his neighbors share enough of their attribute values with him.

We can therefore write the movestay rule for an individual (i.e., for a site occupied by an individual) as

movestay[{a_, b_}, 
        res__] := {a*
          Round[Count[
                Map[Count[b - #[[2]], 0] &, {res} /. 
                    0 -> {0, 0}], _?(# < v/2 &)]/8], b}

The movestay rules are applied to the lattice sites using the anonymous function

Moore[movestay, #] &

where

Moore[func_, lat_] := 
      MapThread[func, 
        Map[RotateRight[lat, #] &, {{0, 0}, {1, 0}, {0, -1}, {-1, 
              0}, {0, 1}, {1, -1}, {-1, -1}, {-1, 1}, {1, 
              1}}], 2]

Moving

Once individuals have decided whether to stay where they are or move, the 28 rules for the movement of an individual are given by the set of walk rules in [Gaylord 1998].

The walk rules are applied to the lattice sites using the anonymous function

GN[walk, Moore[movestay, #]] &

where

GN[func_, lat_] := 
      MapThread[func, 
        Map[RotateRight[lat, #] &, {{0, 0}, {1, 0}, {0, -1}, {-1, 
              0}, {0, 1}, {1, -1}, {-1, -1}, {-1, 1}, {1, 
              1}, {2, 0}, {0, -2}, {-2, 0}, {0, 2}}], 2]

Evolving the System

Here we prepare for memory-efficient code by replacing NestList with a user-defined function $NestingFunction, which is set to NestList for now. In the next section, $NestingFunction will be set to some other command. The system evolves over t time steps, starting with the initial lattice configuration, society, using the following nesting operation.

$NestingFunction = NestList;
$NestingFunction[GN[walk, Moore[movestay, #]] &, society, 
      t]

See the end of this article for the assembled program's code.

Running the Simulation

Run the neighborhood program on a 20 by 20 lattice having a 60% population density of individuals, each of whom has one attribute that can have two possible values (thereby reducing the model to the Schelling model) over 500 time steps. For later comparison with a more efficient solution, Print statements have been added.

Print["Memory in use before run: ", 
        MemoryInUse[]];
SeedRandom[9];
Print[
        First[Timing[
            results = 
              neighborhood[20, 0.60, 1, 2, 
                500]]]];
Print["Memory needed during run: ", 
        MaxMemoryUsed[]];
Print["Numbers contained in result: ", 
        Length[Flatten[results]]];

We can compare the initial and final system configurations, the latter showing definite neighborhood formation.

showLattice[lat_, opts___]  := 
      Show[Graphics[
          RasterArray[
            lat /. {0 -> RGBColor[0.7, 0.7, 0.7], {_, {1}} -> 
                  RGBColor[0, 1, 0], {_, {2}} -> RGBColor[0, 0, 1]}]], 
        AspectRatio -> Automatic, opts]

Show[
        GraphicsArray[
          Map[showLattice[#, DisplayFunction -> Identity] &, {First[
                results], Last[results]}]]];

New Nesting Functions

In the previous section, NestList was used to run a simulation without considering first what sort of analysis was going to be done with the data. By thinking of analysis first, one can greatly improve the real-world performance of both simulations and other nesting operations.

MapNestList

The reason that the results of a simulation are so large is simply that every possible piece of information is returned by NestList. Sometimes this will be necessary, but for some purposes, it is enough to summarize each matrix in a simulation run with some number or list of numbers. For instance, if results is a list of lattice configurations and matrixSummary is a function that computes some quantity from a matrix, then it may be that what's really desired is:

Map[matrixSummary,  results]

If this is the case, there is no need to create the entire list of matrices before reducing them to a summary. We can instead create the summary as we go, only keeping a complete matrix in memory as long as necessary, which usually is only long enough to generate the next step in the simulation run. Here is the definition for a new function, MapNestList, that does exactly this.

MapNestList[h_, f_, expr_, n_]  := 
      Block[{temp = expr, result  =  {h[expr]}}, 
        Do[temp = f[temp]; AppendTo[result, h[temp]], {i, 1, n}]; 
        result]

So MapNestList[h, f, expr, n] is equivalent to Map[h, NestList[f, expr, n]], and MapNestList[Identity, f, expr, n] is the same as NestList[f, expr, n].

Note: MapNestList could be defined in terms of Nest instead of Do, as well as building up the result list more efficiently – which is done in the Implementation section at the end of this article – but the above version makes the underlying procedure more transparent. This function does in fact do what is claimed.

h  /@  NestList[f, expr, 3]

MapNestList[h, f, expr, 3]

MapIndexedNestList

In some cases, even more flexibility is desired. It would be useful, for example, to return the entire initial and final states of the simulation, while summarizing the intermediate states. MapNestList cannot do this in a straight-forward way, because there is no way to keep track of which time step is the current time step. A small generalization of MapNestList is sufficient:

MapIndexedNestList[h_, f_, expr_, n_]  := 
      Block[{c = 0}, 
MapNestList[h[#, c++] &, f, expr, 
          n]]

While the first argument to MapNestList is a function of one variable, MapIndexedNestList uses a function of two variables, the second referring to the time step during which the matrix was created. Notice that the analogy to MapIndexed is somewhat loose, since we need only index items in the first level of the expression, and we start with the zero-th time step, corresponding to the initial matrix configuration.

MapIndexedNestList[h, f, expr, 3]

MapIndexedNestList[#1&, f, expr, n] is equivalent to NestList[f, expr, n]. MapIndexedNestList gives all of the flexibility as the NestListEffects function and related options from [Gaylord 1998], while being easier to code and understand. It is left as a (simple) exercise to code printCounter and printMemory from that reference in terms of MapIndexedNestList.

Nesting in Simulations

In order to use MapNestList or MapIndexedNestList effectively in a simulation, you must decide what information you wish to obtain, and which information can be discarded. For example, the previous section looked at a graphic generated from the first and last states of the simulation run. While this is a very simple example, one can cause the simulation to return just the first and last states via MapIndexedNestList. Recall the command that begins the simulation:

$NestingFunction[GN[walk, Moore[movestay, #]] &, society, t]

To change the behavior, we first define an auxiliary function which returns the lattice configuration for the initial and final configuration, and returns an empty sequence for all other time steps. Thus, the final output of this simulation will be a list of two matrices.

aux[lat_, 0] = lat;
aux[lat_, 500] = lat;
aux[lat_, _] = Sequence[];

Using this function as the first argument to MapIndexedNestList will provide the desired effect.

MapIndexedNestList[aux, GN[walk, Moore[movestay, #]] &, society, 
      t]

To achieve this in the intended way, with resetting $NestingFunction, make the following assignment:

$NestingFunction = MapIndexedNestList[aux, ##] &;

It's clear that this should be more memory efficient than the NestList solution. The following test supports that, while showing it to be just as fast as the original NestList version.

Print["Memory in use before run: ", 
        MemoryInUse[]];
SeedRandom[9];
$NestingFunction = 
        MapIndexedNestList[aux, ##] &;
Print[
        First[Timing[
            results2 = 
              neighborhood[20, 0.60, 1, 2, 
                500]]]];
Print["Memory needed during run: ", 
        MaxMemoryUsed[]];
Print["Numbers contained in result: ", 
        Length[Flatten[results2]]];

Again visualizing the result, we see that the same configurations are obtained using this procedure as NestList.

Show[
        GraphicsArray[
          Map[showLattice[#, DisplayFunction -> Identity] &, 
            results2]]];

Side Effects and Animations

MapNestList and MapIndexedNestList can be used to generate side effects as well as return values during a simulation run, by using a compound auxiliary function, or a CompoundExpression containing more than one auxiliary function. The following code snippets demonstrate, using commands new to version 3, how to open a new notebook window and display each lattice configuration as it is generated. At every step it replaces the previous graphic, which automatically prevents the front end from attempting to save hundreds or thousands of graphics. The values actually returned from this simulation run would be the same as those returned by using just aux.

pasteGraphic[nb_,  lat_, 
        counter_]  := (SelectionMove[nb, All, Notebook]; 
        NotebookWrite[nb, 
          Cell[GraphicsData["PostScript", 
              showLattice[lat, 
                DisplayFunction -> 
                  DisplayString]], "Graphics"]])

aux1[lat_, 0]  :=  (nb  =  NotebookCreate[]; 
          pasteGraphic[nb, lat, 0]);
aux1[lat_, n_]  :=  pasteGraphic[nb, lat, n];

$NestingFunction  = 
        MapIndexedNestList[(aux1[##]; aux[##]) &, ##] &;

Conclusion

One obstacle to running large simulations in Mathematica is relying on the result of NestList, filling the computer's memory with excessive amounts of irrelevant data. Two alternatives, MapNestList and MapIndexedNestList, are simple extensions of NestList which avoid memory problems common in simulations, while still allowing access to important information. These functions will not slow down a simulation run, and can be incorporated into existing code with minimal changes.

References

Axelrod, Robert. 1986. "An Evolutionary Approach to Norms." American Political Science Review 80: 1095–1111.
Durlauf, Steven N. 1997. "The Memberships Theory of Inequality: Ideas and Implications." Sante Fe Institute Working Papers No. 97-05-047. Abstract.
Frank, Robert H. 1987. "If Homo Economicus Could Choose His Own Utility Function, Would He Want One with a Conscience?" American Economic Review 77: 593–604.
Gaylord, Richard J. and Louis J. D'Andria. 1998. Simulating Society. New York: Springer-Verlag.
Schelling, Thomas C. 1978. Micromotives and Macrobehavior. New York: W. W. Norton.
Smith, John Maynard. 1982. Evolution and The Theory of Games. New York: Cambridge University Press, 173.
Tierney, John. 1997. "The Boor War." The New York Times Magazine, 5 January [includes statements attributed to Robert Axelrod].

Implementation

MapNestList::usage  = "MapNestList[h,f,expr,n] calculates \
Map[h,NestList[f,expr,n]] in a memory-efficient way. \
MapNestList[Identity,f,expr,n] is equivalent to NestList[f,expr,n].";

MapNestList1[h_, f_, expr_, n_]  := 
      Block[{temp = expr, result  =  {h[expr]}}, 
        Do[temp = f[temp]; AppendTo[result, h[temp]], {i, 1, n}]; 
        result]

MapNestList[h_, f_, expr_, n_]  := 
      Block[{temp = expr, wrap, result}, result = wrap[h[expr]]; 
        Do[temp = f[temp]; result  =  wrap[result, h[temp]], {i, 1, n}]; 
        List  @@  Flatten[result, Infinity, wrap]]

Note: rethinking the calls to AppendTo can result in significant increase in speed.

Timing[m1 = MapNestList1[a, b, 0, 5000];]

Timing[m2 = MapNestList[a, b, 0, 5000];]

m1 === m2

MapIndexedNestList::usage  =  "MapIndexedNestList[h,f,expr,n] \
calculates {h[expr,0], h[f[expr],1], h[f[f[expr]],2], ..., \
h[Nest[f,expr,n],n]} in a memory-efficient way. \
MapIndexedNestList[#1&,f,expr,n] is equivalent to \
NestList[f,expr,n]";

MapIndexedNestList[h_, f_, expr_, n_]  := 
      Block[{c = 0}, MapNestList[h[#, c++] &, f, expr, n]]

$NestingFunction  =  NestList;

neighborhood[n_, p_, v_, w_, t_] := 
      Module[{walk, movestay, society, RND, Moore, GN}, 
        RND := Random[Integer, {1, 4}]; 
society = 
          Table[Floor[p + Random[]], {n}, {n}] /. 
            1 :> {RND, 
                Table[Random[
                    Integer, {1, w}], {v}]}; 
movestay[
            0, __] := 0; 
movestay[{a_, b_}, 
            res__] := {a*
              Round[Count[
                    Map[Count[b - #[[2]], 0] &, {res} /. 
                        0 -> {0, 0}], _?(# < v/2 &)]/8], 
            b}; 
walk[{1, a___}, 
            0, _, _, _, {4, ___}, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[{1, a___}, 
            0, _, _, _, _, _, _, {2, ___}, _, _, _, _] := {RND, 
            a}; 
walk[{1, a___}, 
            0, _, _, _, _, _, _, _, {3, ___}, _, _, _] := {RND, 
            a}; 
walk[{1, a___}, 
            0, _, _, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[{2, a___}, _, 
            0, _, _, {3, ___}, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[{2, a___}, _, 
            0, _, _, _, {1, ___}, _, _, _, _, _, _] := {RND, 
            a}; 
walk[{2, a___}, _, 
            0, _, _, _, _, _, _, _, {4, ___}, _, _] := {RND, 
            a}; 
walk[{2, a___}, _, 
            0, _, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[{3, a___}, _, _, 
            0, _, _, {4, ___}, _, _, _, _, _, _] := {RND, 
            a}; 
walk[{3, a___}, _, _, 
            0, _, _, _, {2, ___}, _, _, _, _, _] := {RND, 
            a}; 
walk[{3, a___}, _, _, 
            0, _, _, _, _, _, _, _, {1, ___}, _] := {RND, 
            a}; 
walk[{3, a___}, _, _, 
            0, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[{4, a___}, _, _, _, 
            0, _, _, {1, ___}, _, _, _, _, _] := {RND, 
            a}; 
walk[{4, a___}, _, _, _, 
            0, _, _, _, {3, ___}, _, _, _, _] := {RND, 
            a}; 
walk[{4, a___}, _, _, _, 
            0, _, _, _, _, _, _, _, {2, ___}] := {RND, 
            a}; 
walk[{4, a___}, _, _, _, 
            0, _, _, _, _, _, _, _, _] := 
          0; 
walk[{_, 
              a___}, _, _, _, _, _, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[
            0, {3, ___}, {4, ___}, _, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, {3, ___}, _, {1, ___}, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, {3, ___}, _, _, {2, ___}, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, _, {4, ___}, {1, ___}, _, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, _, {4, ___}, _, {2, ___}, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, _, _, {1, ___}, {2, ___}, _, _, _, _, _, _, _, _] := 
          0; 
walk[
            0, {3, a___}, _, _, _, _, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[
            0, _, {4, a___}, _, _, _, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[
            0, _, _, {1, a___}, _, _, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[
            0, _, _, _, {2, a___}, _, _, _, _, _, _, _, _] := {RND, 
            a}; 
walk[
            0, _, _, _, _, _, _, _, _, _, _, _, _] := 
          0; 
Moore[func_, lat_] := 
          MapThread[func, 
            Map[RotateRight[lat, #] &, {{0, 0}, {1, 0}, {0, -1}, {-1, 
                  0}, {0, 1}, {1, -1}, {-1, -1}, {-1, 1}, {1, 
                  1}}], 2]; 
GN[func_, lat_] := 
          MapThread[func, 
            Map[RotateRight[lat, #] &, {{0, 0}, {1, 0}, {0, -1}, {-1, 
                  0}, {0, 1}, {1, -1}, {-1, -1}, {-1, 1}, {1, 
                  1}, {2, 0}, {0, -2}, {-2, 0}, {0, 2}}], 
            2]; 
$NestingFunction[
          GN[walk, Moore[movestay, #]] &, society, t]]

showLattice[lat_, opts___]  := 
      Show[Graphics[
          RasterArray[
            lat /. {0 -> RGBColor[0.7, 0.7, 0.7], {_, {1}} -> 
                  RGBColor[0, 1, 0], {_, {2}} -> RGBColor[0, 0, 1]}]], 
        AspectRatio -> Automatic, opts]