Basic Question
I have a looping calculation that I need to perform on some large lists a very large number of times, and as I have written the functions the computational speed is prohibitively slow. I am looking for how to set the function up using Compile, or find other ways to speed up the calculation.
Context
I am doing Monte Carlo simulations of molecules on the surface of a spherical particle. Each molecule can have a charge of +1, -1, or 0. The probabilities are (P[0]=(1-(P[-1]+P[+1])) > P[-1] > P[+1], but need to be individually specified. Further, I want to shift probabilities based on adjacency; that is, a -1 molecule makes the adjacent molecules less likely to be negative, and more likely to be positive. We will assume that it shifts by a constant factor for each adjacent charge, and the shift is identical in magnitude regardless of direction.
Generating the particle
Thanks to an answer to another question here on StackExchange, I have an efficient way to make a particle:
<< IGraphM`
Needs["GraphUtilities`"];
size = 40; (* Modify this value to change the particle size *)
reg = BoundaryDiscretizeRegion[Ball[], PrecisionGoal -> 1, MaxCellMeasure -> 1/size];
g = IGMeshGraph[reg]
adj = AdjacencyMatrix[g];
Length[adj]
482
The above example, with size set to 40, gives a particle with 482 surface molecules. I need to be able to go up to 100,000 surface molecules. Generating the adjacency matrix only has to be done once for a given particle size, so this isn't the speed-prohibitive part.
Single Monte Carlo Calculation
I have a function, shown below, for calculating the surface charge density (total charge divided by number of surface molecules). The general idea is that I generate a random number from -0.5 to +0.5 for each molecule, and figure out cut-offs at each end that will correspond to the passed probabilities for positive and negative (e.g., a 5% chance of positive means that any value above +0.475 should become a positive charge). I then loop through my list of random numbers taking the most extreme values first. For each value I assign the charge, and then adjust the probabilities of the adjacent molecules by the shift value, and then reset the random number for the molecule in question to zero. None of the shifts will be big enough to take something that gets adjusted to zero and gets it back to the point that it would be reassigned as a charge. And I loop until I have no further charges outside the thresholds. Finally, I total the charges and divide by the number of molecules.
montecarlo[adjacency_, positivechance_, negativechance_, chanceshift_] :=
Module[
{len = Length[adjacency], charges, poscutoff = 0.5 -
(positivechance/2), negcutoff = -0.5 + (negativechance/2),
randomnumlist, topmargin, bottommargin, position, molstoadjust},
charges = Table[0, {i, 1, len}];
randomnumlist = RandomReal[{-0.5, 0.5}, len];
While[
(topmargin = Max[randomnumlist - poscutoff];
bottommargin = Min[randomnumlist - negcutoff];
(topmargin > 0) || (bottommargin < 0)),
If[
topmargin > -bottommargin,
(position = Position[randomnumlist, Max[randomnumlist]][[1, 1]];
molstoadjust =
adjacency[[position]]["NonzeroPositions"]\[Transpose][[1]];
randomnumlist[[molstoadjust]] = randomnumlist[[molstoadjust]] -
chanceshift/2;
randomnumlist[[position]] = 0;
charges[[position]] = 1;),
(position = Position[randomnumlist, Min[randomnumlist]][[1, 1]];
molstoadjust =
adjacency[[position]]["NonzeroPositions"]\[Transpose][[1]];
randomnumlist[[molstoadjust]] = randomnumlist[[molstoadjust]] +
chanceshift/2;
randomnumlist[[position]] = 0; charges[[position]] = -1;)]];
N[Total[charges]/len]]
Running Statistics
I need to run the above Monte Carlo simulation 10,000 or more times per particle. I will then generate histograms of the resulting data and extract means and standard deviations.
montecarloset[numparticles_, adjacency_, positivechance_, negativechance_,
chanceshift_] := Table[montecarlo[adjacency, positivechance,
negativechance, chanceshift],
{i, 1, numparticles}]
A sample run would be something like:
test = montecarloset[10000, adj, 0.03, 0.06, 0.03];
On my machine, the above example takes about two minutes, but that is for a particle with only 482 molecules (and I need to go up to 100,000) and 10,000 simulations (I would be happier going up to 100,000).
Is this something that can be Compiled? Are there other ways to speed this up that I may not have thought of? Thanks!