Here is the physical problem I would like to simulate:
Consider the following variations to the 1d random walk:
allow for certain, random sites to be excluded with probability 1/2.
allow for the walker to be able to hops either 1 or 2 units to left or right
whether hopping 1 or 2 units is weighted by hopping rate
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.
{time, position}
? $\endgroup$t
is localized in yourModule
and all the other variables assigned there are global. In fact, I don't thinkt
is needed at all. Why not just return{newCoords, simT + dt}
and make the assignment{Xi, simT} = FindRate[...]
at the point of call? $\endgroup$