Efficiency in calculation on graphs: Compile?

Question

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!

Answer 1

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 = Compile`GetElement[Ordering[rand, 1], 1];
      maxpos = Compile`GetElement[Ordering[rand, -1], 1];
      min = Compile`GetElement[rand, minpos];
      max = Compile`GetElement[rand, maxpos];
      topmargin = max - poscutoff;
      bottommargin = min - negcutoff;
      (topmargin > 0.) || (bottommargin < 0.)
      )
     ,
     If[topmargin > -bottommargin,
       Do[
        j = Compile`GetElement[adjacency, maxpos, i];
        If[j > 0, rand[[j]] -= shift];
        ,
        {i, 1, maxdegree}];
       rand[[maxpos]] = 0.;
       charges[[maxpos]] = 1.;
       ,
       Do[
        j = Compile`GetElement[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.