Efficiency in calculation on graphs: Compile?

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!

• I haven't looked at it in too much detail, but this seems like a candidate for Markov Chain Monte Carlo (e.g. Metropolis Hastings). This would make sampling a lot easier and probably faster. See also this demonstration in MMA. – Lukas Lang Jun 7 '18 at 7:05

1 Answer

The following should be 150 times faster. Too many micro changes to comment on all of them. Major speedup came from using Ordering to find position directly instead of searching the random number list for max/min entries. Moreover, the "AdjacencyLists" property of of the adjacency matrix was used (padded with zeroes since compiled functions cannot handle ragged lists) instead of the adjacency matrix itself.

montecarlo3 = Compile[{
{adjacency, _Integer, 2}, {positivechance, _Real}, {negativechance, _Real}, {chanceshift, _Real}
},
Module[{minpos, maxpos, min, max, len, maxdegree, charges,
poscutoff, negcutoff, rand, topmargin, bottommargin, shift, j},
len = Length[adjacency];
maxdegree = Dimensions[adjacency][[2]];
poscutoff = 0.5 - (positivechance/2);
negcutoff = -0.5 + (negativechance/2);
shift = chanceshift/2;

charges = Table[0., len];
rand = RandomReal[{-0.5, 0.5}, len];

While[(
minpos = CompileGetElement[Ordering[rand, 1], 1];
maxpos = CompileGetElement[Ordering[rand, -1], 1];
min = CompileGetElement[rand, minpos];
max = CompileGetElement[rand, maxpos];
topmargin = max - poscutoff;
bottommargin = min - negcutoff;
(topmargin > 0.) || (bottommargin < 0.)
)
,
If[topmargin > -bottommargin,
Do[
j = CompileGetElement[adjacency, maxpos, i];
If[j > 0, rand[[j]] -= shift];
,
{i, 1, maxdegree}];
rand[[maxpos]] = 0.;
charges[[maxpos]] = 1.;
,
Do[
j = CompileGetElement[adjacency, minpos, i];
If[j > 0, rand[[j]] += shift];
,
{i, 1, maxdegree}];
rand[[minpos]] = 0.;
charges[[minpos]] = -1.;
];

];
N[Total[charges]/len]],
CompilationTarget -> "C",
RuntimeAttributes -> Listable,
Parallelization -> True,
RuntimeOptions -> "Speed"
];

montecarloset3[numparticles_, adjacency_, positivechance_, negativechance_, chanceshift_] :=
montecarlo3[
PadRight[adjacency["AdjacencyLists"]],
positivechance,
negativechance,
ConstantArray[chanceshift, numparticles]
]


This graph has about 100000 vertices.

Needs["IGraphM"]
size = 14000;(*Modify this value to change the particle size*)
reg = BoundaryDiscretizeRegion[Ball[], MaxCellMeasure -> {2 -> 1/size}];
adj = IGMeshCellAdjacencyMatrix[reg, 0, 0];
Length[adj]
test = montecarloset3[48, adj, 0.03, 0.06, 0.03]; // AbsoluteTiming // First


108002

63.4267

48 Monte-Carlo simulations over it cost you about a minute on my quad core, so 100000 runs can be done in one and half days...

Alternative model

Note that everything could be made way faster if you would not require a certain ordering in the interactions. If all interactions with neighbors would happen simultaneously, this could be achieved with few matrix-vector products within just a few minutes.

For example, the following performs 10000 simulations with a slightly modified model for the graph above in under one minute:

montecarlo4[adjacency_, positivechance_, negativechance_,
chanceshift_] :=
Module[{len, poscutoff, negcutoff, shift, rand, charges, b, u},
len = Length[adjacency];
poscutoff = 0.5 - (positivechance/2);
negcutoff = -0.5 + (negativechance/2);
shift = chanceshift/2;

rand = RandomReal[{-0.5, 0.5}, {len}];
charges = Subtract[UnitStep[Subtract[negcutoff, rand]] - UnitStep[rand - poscutoff]];
b = Max[Abs[u]];
While[b != 0,
rand = ((1 - Unitize[u]) rand) + adjacency.(shift u);
With[{spu = SparseArray[u]},
charges[[Flatten[spu["NonzeroPositions"]]]] = spu["NonzeroValues"]
];
u = Subtract[UnitStep[Subtract[negcutoff, rand]] - UnitStep[rand - poscutoff]];
b = Max[Abs[u]]
];
Total[N[charges]]/len
];

data = ParallelTable[
montecarlo4[adj, 0.03, 0.06, 0.03], {i, 1, 10000},
Method -> "CoarsestGrained"]; // AbsoluteTiming // First
`

53.4714

Note that I did not check the code for correctness. I am almost sure that there is an odd number of sign errors in the code. But still, this might give you some hints how this can be implemented more efficiently.

• Incredible! It took me some time to figure out how to get a C compiler installed, but once I did your first approach worked beautifully (I guess my system better than I thought, b/c the 12 simulation example you gave took just under 7 s for me). I am lost on the mechanics of your alternate method. How does it ensure that an assigned charge that shifts an adjacent molecule below threshold doesn't result in a charge being assigned to that adjacent molecule? Regardless, the first solution you provided is brilliant, and does what I need in a tractable period of time. Thank you! – Kevin Ausman Jun 7 '18 at 18:40
• Hey Kevin. The alternate method guarantees nothing. That's the issue. I just tried to show how a small change in the model could open the door for more efficient treatments. A model is by definition a simplification that allows us to deal with it. So, no need to create a stone that you cannot lift, you know. ;) Of course, I am aware that you might not have a choice... I also checked my timings for the first method. Should have been 15 seconds for 12 simulations. I must have mixed something up. I am glad that I could help you! – Henrik Schumacher Jun 7 '18 at 19:13
• Gotcha, thanks. So this model is an outgrowth of one where I treated all molecules as independent, and thus could ignore adjacency. That original model predicted that the population standard deviation would scale as 1/Sqrt[n], which does not match my experimental data. The adjacency effect is an attempt to find a better-fitting model. And unfortunately I think trying to do the calculation all at once takes me back to the original model. – Kevin Ausman Jun 7 '18 at 21:32