How to use the behseq package, v1.1

Copyright A . Grafen 2002,  2005

Introduction

Preliminaries

This file is intended to help you use "behseq", a Mathematica package that implements the optimisation model for sequences of behaviour introduced in the source paper (A. Grafen 2002 'A state-free optimisation model for sequences of behaviour', Animal Behaviour 63: 183-191). The source paper should be consulted for an introduction to the purpose and mathematical details of the model. This file discusses only the computer implementation.

You can execute the input cells in this notebook, and use the package as described.  So this is an active notebook.  You can  evaluate each input cell as you come to it (recognising them by the lack of an extra bar on the innermost bracket over to the right of the window).  (Or you can select 'Evaluate Notebook' from the menus, but the final three cells would take hours to complete on my machine, and anyway crash it if evaluated in the same session. So evaluating all the cells is not recommended for all users.)

Part 1 shows how to perform the simplest possible analysis, and should help a user make a start. Part 2 contains the code that produced the data for the figures in the source paper, and illustrates various extra tricks for using the package. For example, it will be usual to want to show how the optimum changes as a parameter changes.

The programming is rough and ready, and would win no prizes for elegance or efficiency). I offer it to encourage the use of the model, but would be delighted if someone were to write a better program for the purpose.

Caveats in relation to Mathematica before version 5, in 2002

The code has worked well for my situation, but I have not tried to test it with different kinds of problem. It may have new kinds of difficulties in new problems.

FindMinimum is the underlying Mathematica command used in behseq. It can operate in two modes, and behseq assumes FindMinimum is operating in the mode in which derivatives can be calculated. This worked fine for my situation, and the only time FindMinimum complained that it couldn't find derivatives was when some other error had occurred and the utilities were not properly specified (such as not setting one of the parameter values). If you really need to use the non-derivative mode, you will need to alter the underlying program as two sets of starting conditions are required.

It is not logical that tstart should be specified by type of activity rather than by place in the sequence. For sleeping, one might not want to start with a main sleep and an afternoon nap at the same initial durations. But it is simpler to specify tstart in that way, and this is how it is done in the program. For now.

Precision is not guaranteed. It is likely to be much better for the attained utility than for the durations. It is common to find that different starting points give slightly different durations, but almost exactly the same utility. This is a computational reflection of the zero slope at the maximum of a smoothly varying function.

Speed is OK on my machine for all except the final set of calculations, for Figure 5b. This will depend on your machine etc.

Caveats in relation to Mathematica 5 in 2005

There is an upwards incompatibility between Mathematica 4.2 and Mathematica 5, in relation to FindMinimum. Apparently, its "execution semantics" have changed. This appeared to break the package, and this version 1.1 presents one solution; version A presents another (see website http://users.ox.ac.uk/~grafen/behseq/).

It turns out that the examples in this notebook will all work with the behseq package after all, provided the definition of f is changed, so that when an If is used, the value part is surrounded by "N[ ]". I presume this has something to do with embedded If statements and Mathematica's willingness to evaluate the conditions and pick an outcome rather than retain the whole If statement unevaluated. (I infer this from investigations of what goes wrong without the "N[ ]".) However, the fix is pretty much magic as far as I'm concerned, and I can offer no guarantee that other forms of f won't also cause problems.

Hence the package itself is unchanged, but the way to use it has changed, in version 1.1 for Mathematica 5.

So my advice for version 1.1 is (i) always surround the value parts of an If with "N[ ]" in the definitions of f and g (ii) if other properties of f and g cause problems, switch to version A, which seems slower but *should* work.

The speed of all the calculations is now very reasonable, and I have altered the notebook accordingly.

A note on the output

The output in this version 1.1 of howtousebehseq.nb is from Mathematica 5.2, though it still discusses Mathematica 4 as well. To see Mathematica 4 output, you need version 1.0.

Part 1 - A simple worked example - optseq

Loading the package

The first task is to load the package. You should place the folder "behseq" in the Extra Packages folder, which is inside the Add-Ons folder, which you will find in the Mathematica folder. Then you can load the file as follows. Don't forget the `. If Mathematica says it can't load the file, then it is probably in the wrong place. It may be necessary to right-click or control-click on the Mathematica icon and choose "Show Package Contents" to access these folders.

<< behseq`

Preparing to use the package

Defining the model involves specifying a number of functions and parameters. These have pre-defined name in the package, and so you must load the package first and then use those names.

You need to specify

1. nacts. The number of different activities. In our example, animals will be able to eat, drink and rest, so nacts=3.

2. f. These are the f-functions in the source paper. They are indexed functions so that f[1][t] should give the benefit derived from a bout of activity 1 of duration t. f[2][t] and so on up to f[nacts][t] should do the same for every other activity. See introduction for advice on defining f in Mathematica 5.

3. g. These are the g-functions in the source paper. Like f, except they are the benefit derived from not doing an activity for duration t. See introduction for advice on defining g in Mathematica 5.

4. u. The variable u has two types of component associated with it. First, a set of utility fragments need to be defined. Thus u[1][b] is the utility fragment for activity 1 when the net benefit rate is b, and u[2][b] up to u[nacts][b] must be defined. Second, the function u[{uf1,uf2...ufnacts}] must be defined to say how the utility fragments are combined into a single final utility.

These are all fundamental to the theory. There are two further variables that need to be defined before the program will work at all.

5. tstart. Mathematica starts searching for a maximum of utility, but it must be told where to start. tstart[1] gives the starting duration for all bouts of activity 1, and so on up to tstart[nacts]. Times are logged internally to avoid trouble at zero, so tstart values must be greater than and not equal to zero.

6. paramlist. The point of an analysis will usually be to show how the optimal behaviour changes as some parameter varies. Thus, the rate of food requirement per unit time will affect f[1] and g[1] if activity 1 is feeding. It is possible to include parameters in the definitions of f, g and u representing these quantities. But when Mathematica does the calculations, it needs to know the numerical value. If there are parameters foodreq and waterreq, then paramlist must set them in rules, e.g. paramlist={foodreq->3.1, waterreq->5.2}.

There is one further symbol that needs to be defined, called 'follow'. Its use is illustrated in Part 2.

This section shows sample definitions of the required symbols, and conducts an analysis.

nacts

Setting nacts is easy. We use := and want it to be 3 in our case.

nacts := 3

f and g

We now set the f-functions and g-functions. There is one f and one g for each of the three activities.

f[1][x_] := 
        If[x <= tmin[1], 0, N[1 - Exp[-lambda[1]*(x - tmin[1])]]] - 
          mu[1]*x;
g[1][x_] :=  - mu[1]*x + 0;
f[2][x_] := 
        If[x <= tmin[2], 0, N[1 - Exp[-lambda[2]*(x - tmin[2])]]] -  
          mu[2]*x;
g[2][x_] :=  - mu[2]*x;
f[3][x_] := x;
g[3][x_] := 0;

Notice that parameters are used. In f[1], there are tmin[1], lambda[1] and mu[1]. These must be set in paramlist before the optimum is found. These are the functions used in the source paper.  Changes between versions 4.2 and 5 of Mathematica make it necessary to enclose the value part of the If in "N[ ]", for reasons I don't understand. See introduction for discussion.

u

First we specify the utility  fragments

\(\(\(\(u[1]\)[b_] := 
        If[b <= 0, 
          0, \((b/omega[1])\)^\((omega[1]*kappa[1])\)*
            Exp[\(-kappa[1]\)*\((b - omega[1])\)]];\)\(
\)
    \)\)
\(\(\(\(u[2]\)[b_] := 
        If[b <= 0, 
          0, \((b/omega[2])\)^\((omega[2]*kappa[2])\)*
            Exp[\(-kappa[2]\)*\((b - omega[2])\)]];\)\(
\)
    \)\)
u[3][b_] := 10^3 *4*b*(1 - b);

Again there are parameters, explained in the source paper. It is also convenient to scale the utility. In this case the maximum utility is 1000. u[1] and u[2] have maximum possible values of 1, and u[3] has a maximum possible value of 1000. This helps to interpret the results, but is also important for numerical reasons. If the utilities drift to very small values, then differences between them become too small for Mathematica's FindMinimum command to worry about.

Next we specify how the utility fragments combine to produce a final utility

u[{uf1_, uf2_, uf3_}] := uf1*uf2*uf3;

Notice the curly brackets enclosing the list of values as well as the square brackets. Yes, they are needed.

tstart

The starting values for the search must be specified. If you start in the wrong place, you can have difficulties. One obvious difficulty is that if you start within a region where utility is at some minimal value, then all the derivatives with respect to the durations are zero, and the procedure stops. Thus you may need to use trial and error to find suitable starting points.

There are two ways to specify tstart. One is to specify it as follows

tstart[1] := 0.3;
tstart[2] := 0.3;
tstart[3] := 0.5;

The other way is to include these values in paramlist, as shown next

paramlist

When Mathematica is asked to maximise utility, it must be able to calculate the utilities as real numbers, so all the parameters included in the f- and g-functions, and in the utility functions, must have been assigned values. This can be done as follows, continuing the example from earlier

paramlist := {lambda[1] -> 1.0, lambda[2] -> 1.0, 
          mu[1] -> 0.01, mu[2] -> 0.01, tmin[1] -> 0.2, 
          tmin[2] -> 0.2, kappa[1] -> 10.0, kappa[2] -> 10.0, 
          tmin[3] -> 0.2, omega[i_] -> 0.19, 
          tstart[i_] -> 0.3};

This  sets lambda[1] to 1.0, etc. Notice that decimal values are used rather than integers. This is because Mathematica will attempt to operate in exact arithmetic when handling integers. This can be wasteful when, as in this case, the result will eventually be converted to real before being maximised. Notice also that patterns can be used, so that "omega[i_]->0.19" is equivalent to listing both "omega[1]->0.19" and "omega[2]->0.19". Finally, notice that tstart has been included in paramlist. This will override the setting made previously. We require only one of these methods to be used for tstart.

Obtaining the optimal durations for a given sequence

The command is optseq. It takes as argument the sequence of activities. Thus

optseq[{1, 2, 3}]

produces the optimal durations when the activities are performed in order eating, drinking, resting, eating, drinking resting etc. There is a warning about tolerances. The output looks like this

The first number of the utility obtained at the optimal durations. Then there is a list that specifies the duration of each type of bout.

A more complicated example would ask for the sequence {1,3,2,3}, as follows, in which there is a rest after feeding and after drinking.

optseq[{1, 3, 2, 3}]

Part 2 - Searching through sets of sequences - sortedseqs

To use optseq, you need to specify a given sequence. Often the question will instead be which sequence is optimal. Some tools are provided to help you search through sequences.Note that it can take a long time to process a lot of sequences. We will look at how you specify the sequences you want to search through, and then how to conduct the search.

Specifying the sequences to search

The facilities provided offer you two ways to determine the sequences to be searched. First, the length of the sequence, that is, the number of separate bouts.Second, you can specify which activities are allowed to follow which other activities. We specify the following rules by setting the variable 'follow'. For example, if eating must be followed by resting, and drinking must be followed by resting, but resting may be followed by either eating or drinking, then we would define follow as follows:

follow[1] := {3};
follow[2] := {3};
follow[3] := {1, 2};

Conducting the search

Let's look first at a request to check all sequences of lengths from 2 to 4, with the following rules defined earlier. This command is not evaluatable, to preserve the output, but there is an evaluatable copy given just afterwards that you can execute yourself.

sortedseqs[4]

sortedseqs[4]

This list produced at the end has five elements, corresponding to the five possible sequences. The program has eliminated other sequences because (i) some are not consistent with the settings of 'follow' and (ii) some are equivalent to an existing sequence because they are cyclically the same. That is, {3,2} has been eliminated because it is the same as {2,3}. It is the same because the sequence is considered to repeat indefinitely.

Each element gives (i) the sequence of activities (ii) the utility gained and (iii) the durations of the bouts. The list is sorted in decreasing order of utility, so that the optimal sequence comes first. It is useful to see the whole list because sometimes there are a number of equivalently good sequences

The warning from FindMinimum appears for those sequences that don't have some of each activity. The utility function we have chosen puts total utility to zero unless all activities are present, and so the derivatives of utility with respect to durations are all zero. Thus in this case the warnings are expected and nothing to worry about.

Let us now look for longer sequences. We avoid retesting shorter ones by giving two integers in a list, which are the minimum and maximum length of the sequences to be considered. Again this command is not evaluatable, to preserve the output, but an evaluatable copy follows it.

sortedseqs[{5, 6}]

sortedseqs[{5, 6}]

There are only three eligible sequences, and none matches the utility achieved by {1,3,2,3} in the previous output.

This section has shown how to search for optimal sequences of certain lengths. The command optseq can also be used in more tailored ways, as the next section shows.

Part 3 - Generating the figures for the source paper

This section shows how to use the package to generate the data for the figures in the source paper. Various methods of using optseq are illustrated. Variations on the model can be easily investigated by tweaking the code.

Tables 2 to 5 in the source paper were generated using the package, and we take them in turn. sortedseqs is *not* used, but instead other methods of harnessing optseq are.

Figure 2 - automating collection of results while a parameter is varied

The basic ingredients of the model are set up. These are repeats of the definitions in Part 1.

nacts := 3;

    
f[1][x_] := 
        If[x <= tmin[1], 0, N[1 - Exp[-lambda[1]*(x - tmin[1])]]] - 
          mu[1]*x;
g[1][x_] :=  - mu[1]*x;
f[2][x_] := 
        If[x <= tmin[2], 0, N[1 - Exp[-lambda[2]*(x - tmin[2])]]] -  
          mu[2]*x;
g[2][x_] :=  - mu[2]*x;
f[3][x_] := x;
\(\(\(g[3]\)[x_] := 0;\)
\)
\(\(\(u[1]\)[b_] := 
        If[b <= 0, 
          0, \((b/omega[1])\)^\((omega[1]*kappa[1])\)*
            Exp[\(-kappa[1]\)*\((b - omega[1])\)]];\)

\)
\(\(\(u[2]\)[b_] := 
        If[b <= 0, 
          0, \((b/omega[2])\)^\((omega[2]*kappa[2])\)*
            Exp[\(-kappa[2]\)*\((b - omega[2])\)]];\)

\)
u[3][b_] := 
        10^3 *4*
          b*(1 - b);

\(\(u[{uf1_, uf2_, uf3_}] := uf1*uf2*uf3;\)
\)
tstart[1] := 0.22;
tstart[2] := 0.22;
tstart[3] := 0.44;

The idea is to vary omega[1] and omega[2] while investigating the fixed sequence {1,2,3}. We set up 'baserules' to define all the other parameters. Then getresults is defined, which takes a value of omega as an argument, sets up 'paramrules' by adding the rule for omega, and applies optseq to {1,2,3}.

baserules = {lambda[1] -> 1, lambda[2] -> 1, 
          mu[1] -> 0.01`, mu[2] -> 0.01`, tmin[1] -> 0.2`, 
          tmin[2] -> 0.2`, kappa[1] -> 10.0, 
          kappa[2] -> 10.0`, 
          tmin[3] -> 0.2`};
getresults[
          myomega_] := 
(paramlist := 
            Join[baserules, {omega[i_] -> 
                  myomega}]; 
optseq[{1, 2, 
              3}]);

Now we apply getresults to values of omega from 0.1 to 3.0, and include in the results the value of omega, and then the optseq output. The results are held in the variable omtab, and include all the information needed to plot Figure 2 in the source paper.

omtab = 
      Map[{#, getresults[#]} &, Table[0.1*ii/10, {ii, 1, 30}]]

Figure 3a - resolving problems with tstart [but no problem with Mathematica 5]

We use the same basic definitions as Figure 2, and assume they have already been defined, and set out to investigate varying mu[2]. 'baserules' now defines everything except mu[2], and the output is stored in 'mu2tab'. But the results contain an unwelcome surprise. The mu2tab cell is not evaluatable, to preserve the output, but an evaluatable copy follows it.

baserules = {lambda[1] -> 1, lambda[2] -> 1, 
          mu[1] -> 0.01`, tmin[1] -> 0.2`, tmin[2] -> 0.2`, 
          kappa[1] -> 10.0, kappa[2] -> 10.0`, 
          tmin[3] -> 0.2`, omega[1] -> 0.2, 
          omega[2] -> 0.2};
getresults[
          mymu2_] := 
(paramlist := 
            Join[baserules, {mu[2] -> 
                  mymu2}]; 
optseq[{1, 2, 
              3}]);

This is Mathematica 4 output

mu2tab = 
      Map[{#, getresults[#]} &, Table[0.01*ii/10, {ii, 1, 30}]]

This is Mathematica 5 output.

mu2tab = 
      Map[{#, getresults[#]} &, Table[0.01*ii/10, {ii, 1, 30}]]

The surprise is that for the final eight values of mu[2], the utility is zero, under Mathematica 4. (Users of Mathematica 5 will probably find that these eight zero values do not occur, and so it is not necessary or useful to alter tstart. The explanation is retained for the immediate benefit of users of earlier versions, and also to help users of Mathematica 5 who find it necessary in other examples to alter tstart.) The error messages tell us that there were zero derivatives in those cases. What has happened is that the starting durations achieved zero utility in those cases, while they attained positive utility in the rest. In a region of zero utility, with the current definitions of u, the derivatives are all zero, and FindMinimum doesn't know which direction to go in. This could in principle have meant that there was no viable sequence with such high values of mu[2]. But here it didn't, and by inspecting tstart and changing it, the data was all generated on the next attempt.

?? tstart

tstart[1] := 0.73;
tstart[2] := 0.81;
tstart[3] := 0.42;

Here is Mathematica 4 output with the new values of tstart.

mu2tab = 
      Map[{#, getresults[#]} &, Table[0.01*ii/10, {ii, 1, 30}]]

Figure 3b - checking out an analytical equivalence

Now we want to do the same as Figure 3a, except that there are two rest periods i.e. the sequence is {1,3,2,3}. Having just calculated for Figure 3a, this requires only changing the argument of optseq in the definition of getresults. We store the results in mudoubresttab. The key point here is that (because of the linearity of the g-functions and the fact that the benefit of rest depends only on the total over the cycle and not on individual durations) the utilities and durations of eating and drinking are the same, and the total duration of rest over the cycle is also the same. This acts as a general check on the whole procedure, as the program confirms an analytically known equality.

getresults[
          mymu2_] := 
(paramlist := 
            Join[baserules, {mu[2] -> 
                  mymu2}]; 
optseq[{1, 3, 2, 
              3}]);

mu2doubresttab = 
      Map[{#, getresults[#]} &, Table[0.01*ii/10, {ii, 1, 30}]]

Figure 4 - adding doors as activities

Figure 4 introduces doors. A door must be opened to get to the eating chamber, and another door must opened to get to the drinking chamber. Technically, we declare opening the two kinds of door as two further activities. So we redefine nacts as 5,  add the missing f- and g-functions, redefine the u[{uf1,...ufnacts}] function, and define the missing utility fragment u[i] functions.

The door-opening is deemed in the source paper to take no time. The program does not allow this. So we slightly fake it. By defining the f- and g-functions so that the time taken in opening doors has no effect on the utility fragments, and by defining tstart[4] and tstart[5] to be very small, we persuade the program not to change the durations of door-opening, and not to allow the durations to have any measurable impact on the durations of feeding and drinking or on the final utility.

nacts := 5;

f[4][x_] := beta[1];
g[4][x_] := 0;
f[5][x_] := beta[2];
g[5][x_] := 0;

u[4][b_] := Exp[b];
u[5][b_] := Exp[b];

    
u[{uf1_, uf2_, uf3_, uf4_, uf5_}] := 
        uf1*uf2*uf3*uf4*uf5;

tstart[4] := 1.0*^-10;
tstart[5] := 1.0*^-10;

To work with doors, some extra procedures are now defined that will be used in all the remaining figures.

doorexp[l_] := 
        l /. {1 -> Sequence[4, 1], 
            2 -> Sequence[5, 2]};
doorseqs[n_] := 
        doorexp /@ (Append[#, 3] & /@ 
              Table[Join[Array[1 &, ii], Array[2 &, n - ii]], {ii, 1, 
                  n - 1}]);
sorted[nacts_, 
        seqs_] := 
Sort[({#1, optseq[#1]} &) /@ 
          seqs, #1[[2, 
              1]] > #2[[2, 
              1]] &]

'sorted' applies optseq to a list of sequences and sorts the results by utility. 'doorexp' takes a list of sequences containing 1's, 2's and 3's, and substitutes '4,1' for every 1 and '5,2' for every 2. Thus it expands a sequence by adding the necessary door-openings. 'doorseqs[n]' generates all sequences of length n that have at least one 1, at least one 2 and exactly one 3; and then adds the door-openings to them. (Because of the linearity of the g-functions and the fact that the benefit rate of rest depends only on the total rest in one cycle, these are the only cases we need to consider. If the g-functions were non-linear, then the different orders of eating and drinking would need to be enumerated and tested.)

Now we define 'doorresults'. This takes the disutility of opening the eating door as first argument, the disutility of opening the drinking door as the second argument, and the number of eating-or-drinking bouts as the third argument. For the case of Figure 4, the betas are always equal and nbouts=2.

doorresults[beta1_, beta2_, 
          nbouts_] := 
(paramlist = 
            Join[baserules, {beta[1] -> beta1, 
                beta[2] -> beta2}]; 
Print[beta1, 
            beta2]; 
sorted[5, 
            doorseqs[nbouts]]);

We can now generate the data for Figure 4 and store it in the variable bothdoorstab, as follows (again defining baserules to set all parameters except the ones to be varied). doorresults prints out (beta1,beta2) each time just to confirm how far its got.

baserules = {lambda[1] -> 1, lambda[2] -> 1, 
          mu[1] -> 0.01`, mu[2] -> 0.01`, tmin[1] -> 0.2`, 
          tmin[2] -> 0.2`, tmin[3] -> 0.2`, omega[1] -> 0.2`, 
          omega[2] -> 0.2`, kappa[1] -> 10.0`, 
          kappa[2] -> 10.0`, tmin[4] -> 1.0*^-10, 
          tmin[5] -> 1.0*^-10};
bothdoorstab = 
      Map[{-#, doorresults[-#, -#, 2]} &, 
        Table[ii/10.0, {ii, 0, 30}]]

Figures 5a and 5b - further doors

In principle, all we need to do is to supply different arguments to 'doorresults'. For Figure 5a, we keep the eating door costless and allow the drinking door to vary, but keep nbouts=2. I repeat 'baserules' for the record, and store the results in onedoortab.

baserules = {lambda[1] -> 1.0, lambda[2] -> 1.0, 
          mu[1] -> 0.01`, mu[2] -> 0.01`, tmin[1] -> 0.2`, 
          tmin[2] -> 0.2`, tmin[3] -> 0.2`, omega[1] -> 0.2`, 
          omega[2] -> 0.2`, kappa[1] -> 10.0`, 
          kappa[2] -> 10.0`, tmin[4] -> 1.0*^-10, 
          tmin[5] -> 1.0*^-10};

    
onedoortab = 
      Map[{-#, doorresults[0.0, -#, 2]} &, 
        Table[ii/10.0, {ii, 0, 30}]]

This worked for me. Notice that there is an extra pair of {} around the optseq output. That's because doorresults delivers a list, but in this case, there is only one element in the list.

THE REST of this cell is about Mathematica 4 and based on 2002 technology. Following the commands, there is a cell about Mathematica 5 and 2005 technology, for calculating the Fig 5b data.

The difficulty with Figure 5b was that nbouts=12, to allow different numbers of eating and drinking bouts. This increased greatly the time taken by Mathematica, and also the memory required. For some reason or other, the program seized up part way through, and so I needed to do Figure 5b in three sections, and add them together afterwards. In fact, I wrote each out to a file and read them in again afterwards to join them together. If your machine is bigger and/or more efficient in some relevant sense, then you may be able to do all this in one.

This time, the first element of doorresults is extracted (with the [[1]]) to take just the optimal sequence. For exploration, one might want to look at the whole list. This is very slow on my machine, taking minutes per value of beta[2]. This is a good example of a common conflict in simulation. It is clear that the bouts of the same activity must have the same length in the optimal sequence, and a more specialised program could be written to take advantage of that and obtain very much speedier results. On the other hand, it would cost time and effort to write and debug that program, which would have only limited use.

The first command does three values of beta[2], and takes ten to fifteen minutes on my machine. Unless you need to, it probably isn't worth doing the rest. To prevent accidental waste of time, I've made the final three cells unevaluatable - you can swap that back, or copy and paste the contents to a new cell if you do want to execute them.

onedoor12bouts1 = 
      Map[{-#, doorresults[0.0, -#, 12][[1]]} &, 
        Table[ii/10.0, {ii, 0, 2}]]

onedoor12bouts2 = 
        Map[{-#, doorresults[0.0, -#, 12][[1]]} &, 
          Table[ii/10.0, {ii, 3, 10}]];

onedoor12bouts3 = 
        Map[{-#, doorresults[0.0, -#, 12][[1]]} &, 
          Table[ii/10.0, {ii, 11, 20}]];

onedoor12bouts4 = 
        Map[{-#, doorresults[0.0, -#, 12][[1]]} &, 
          Table[ii/10.0, {ii, 21, 30}]];

For Mathematica 5 users in 2005 an ordinary laptop has no difficulty with all the Fig 5b calculations at once, so there is no need to split the calculation up. The following command does all at once. I've made it unevaluatable but you can copy and paste it into a new cell. The whole calculation for ii from 0 to 30 takes under four minutes.

onedoor12boutsAll = 
      Map[{-#, doorresults[0.0, -#, 12][[1]]} &, 
        Table[ii/10.0, {ii, 0, 30}]]