1D Random Walk with Excluded Sites and Variable Time-Step

Question

Here is the physical problem I would like to simulate:

Consider the following variations to the 1d random walk:

1. allow for certain, random sites to be excluded with probability 1/2.

2. allow for the walker to be able to hops either 1 or 2 units to left or right

3. whether hopping 1 or 2 units is weighted by hopping rate

4. time in between hops is taken from exponential distribution with width being total hopping rate (to all 4 possible sites)

In the end, I would like to find out how much slower a walker takes to reach, say x = 100, when the exclusions are present as to when they are not.

The following will be my attempt at the problem. I think I have the basics down but I am hoping there are areas to decrease the computation time. These problems were alluded to earlier in the question and answers of 1D Random Walk variant.

First the 4 different hop rates are defined:

kBack1[k0_] := k0 Exp[- 2 ];
kFor1[k0_] := k0 Exp[- 2 ];
kBack2[k0_] := k0 Exp[- 2  * 2];
kFor2[k0_] := k0 Exp[- 2 *2];

I know these are redundant since there are actually only two unique rates. However, I would like to eventually add a bias to the random walk which will make the 4 rates all different.

Now we need to choose which site the walker will hop to and how long it takes to make the hop. So the excluded sites are accounted for by If statements. If an excluded site is nearby, there is a probability that that site is unreachable. The new site (newCoords) is chosen by weighting the different hopping rates. The time-step is sampled from the exponential distribution and the simulation time is updated.

FindRate[k0_] := Module[{},
  kB1 = If[({Xi} \[Intersection] blockedSites + 1) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kBack1[k0]}], kBack1[k0]];
  kB2 = If[({Xi} \[Intersection] blockedSites + 2) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kBack2[k0]}], kBack2[k0]];
  kF1 = If[({Xi} \[Intersection] blockedSites - 1) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kFor1[k0]}], kFor1[k0]];
  kF2 = If[({Xi} \[Intersection] blockedSites - 2) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kFor2[k0]}], kFor2[k0]];

newCoords = 
RandomChoice[{kB1, kB2, kF1, kF2} -> {Xi - 1, Xi - 2, Xi + 1, 
  Xi + 2}];

dt = RandomReal[ExponentialDistribution[kB1 + kB2 + kF1 + kF2]];
t = simT;
{Xi, simT} = {newCoords, t + dt}
]

EDIT: More explanation on how the hop is chosen: Consider the walker to be at point 4. The next hop is chosen in the following way. Is the walker 1 or 2 spaces away from an excluded point? Let's say point 6 is excluded. The rate (kF2) to point 6 is then 0. The other points (2,3, and 5) can be hopped to with respective rate kB2,kB1, and kF1. Which one does it hop to? Well, compare the different rates. In the example we'll have kB2, kB1 = kF1, and kF2 = 0. Which rate is chosen is decided by RandomChoice where the choices are weighted by the rates. END EDIT

Now it's time to string the hops together

hopMod[k0_, Xf_] := Module[{},
Xi = 0; simT = 0;
blockedSites = RandomInteger[{-100, Xf-2}, 11];(*keep exclusion away from endpoint*)
NestWhileList[FindRate[k0] &, {0, 0}, Xi < Xf &]
]

where I define a random set of excluding sites. I would prefer to not have adjacent or next-adjacent points not both be exclusions (so separated by at least 2 normal sites). Is there a way to force RandomInteger to not sample consecutive integers?

Now calculate the time to reach a specified endpoint, Xf. And average this for n configurations.

With[{k0 = 1, Xf = 10, n = 20},(* n is # of runs to average over*)
 Print[AbsoluteTiming[
  Table[p[i] = hopMod[k0, Xf], {i, 1, n}];]];
 iTable = Flatten[Table[p[j][[-1]][[1]]/p[j][[-1]][[2]], {j, 1, n}]]]

This took me 1.33 seconds. The time is highly varying. The time increases dramatically as Xf or increase. For some reason, when I used ParallelTable, I got identical results every time I calculated another iTable. Must be something with the seed going wrong. To get convergence, n needs to be increased.

So how can this be sped up? Some different schemes for the hopMod module were discussed at 1D Random Walk variant. However I was not sure how to generalized all those approaches to my situation where I also need to keep track of the simulation time.

Answer 1

Here are some tweaks, plus a way to generate non-adjacent blocked sites. (I wasn't sure whether you wanted two free sites between blocked sites or just one. There's just one below, and comment where it could be changed to two.) The run time (per hop) is shorter by about 2/3, only a modest improvement.

kBack1[k0_] := k0 Exp[-2];
kFor1[k0_] := k0 Exp[-2];
kBack2[k0_] := k0 Exp[-2*2];
kFor2[k0_] := k0 Exp[-2*2];

FindRate[k0_, {Xi_, simT_}] :=
   ({kB1, kB2, kF1, kF2} =
      {kBack1[k0], kBack2[k0], kFor1[k0], kFor2[k0]} ({Xi - 1, Xi - 2, Xi + 1, Xi + 2} /. blockedSites);
    {RandomChoice[{kB1, kB2, kF1, kF2} -> {Xi - 1, Xi - 2, Xi + 1, Xi + 2}], 
     simT + RandomReal[ExponentialDistribution[kB1 + kB2 + kF1 + kF2]]});

generateBlockedSites[min_, max_, nSites_] := Reap[Nest[
     DeleteCases[#, 
       With[{site = RandomChoice[#]}, Sow[site :> RandomInteger[]]; 
        Alternatives @@ Range[site - 1, site + 1]]] &, (* Or Range[site - 2, site + 2] ?? *)
     Range[min, max],
     nSites
     ]; Sow[_Integer -> 1]][[-1, -1]];

hopMod[k0_, Xf_] := Module[{}, Xi = 0; simT = 0;
   blockedSites = generateBlockedSites[-100, Xf - 2, 11];(*keep exclusion away from endpoint*)
   NestWhileList[FindRate[k0, #] &, {0, 0}, #[[1]] < Xf &]];

Clear[p];
With[{k0 = 1, Xf = 10, n = 20},(*n is # of runs to average over*)
 Print[AbsoluteTiming[Table[p[i] = hopMod[k0, Xf], {i, 1, n}];]];
 iTable = Table[p[j][[-1]][[1]]/p[j][[-1]][[2]], {j, 1, n}]] (* Flatten unnecessary *)

The variable blockedSites is a list of replacement rules of the form

In[216]:= blockedSites

Out[216]= {-68 :> RandomInteger[], -24 :> RandomInteger[], -62 :> RandomInteger[],
  -88 :> RandomInteger[], -82 :> RandomInteger[], -39 :> RandomInteger[],
  -53 :> RandomInteger[], -27 :> RandomInteger[], -65 :> RandomInteger[],
  -60 :> RandomInteger[], -97 :> RandomInteger[], _Integer -> 1}

Basically it takes a list of sites {Xi - 1, Xi - 2, Xi + 1, Xi + 2} and changes them to 1s or 0s, according to whether the current site can hop to the corresponding sites. This is then multiplied by the rates {kBack1[k0], kBack2[k0], kFor1[k0], kFor2[k0]} to get the rates {kB1, kB2, kF1, kF2} for the hop. I hope that preserves the model for the simulation.