I made this code for generating networks following the duplication-divergence model, in which we have to steps:

  • duplication: a node i is selected at random. A new node i', with a link to all neighbors of i is created, and with a probability p, a link betwenn i and i' is established;
  • divergence: for each of the nodes j linked to i and i' we choose randomly one of the two links (i,j) or (i',j) and remove it with probability q.

The code I have is the following:

  ddmodel[n_, p_, q_] /; n >= 3 := 
 Module[{g, vdeg, el, el2, node, newnode, links}, 
  g = ConstantArray[{}, n];
  g[[1 ;; 2]] = Table[AdjacencyList[CompleteGraph[2], j], {j, 2}];
  vdeg = ConstantArray[0, n];
  vdeg[[1 ;; 2]] = 1;

   node = RandomChoice[Range[j - 1]];
   g[[j]] = g[[node]];
   vdeg[[j]] = vdeg[[node]];
   If[RandomReal[] <= p,
    g[[node]] = Append[g[[node]], j];
    g[[j]] = Append[g[[j]], node];
    vdeg[[j]] += 1;
    vdeg[[node]] += 1;, Continue[]];

   If[ RandomReal[] <= q,
    Do[vdeg[[k]] += 1;
        g[[k]] = Append[g[[k]], j], {k, 
         Flatten[Table[Drop[g[[node]], {l}], {l, #}]]}] &@
    Do[vdeg[[k]] += 1; g[[k]] = Append[g[[k]], j], {k, g[[node]]}];];
   , {j, 3, n}];

    Map[Sort, Flatten[MapIndexed[Thread[{#2[[1]], #1}] &, g], 1]]]]]

The code works fine and quickly if i choose q=0.1 and p from 0.1 to 0.3. But since I have to test is for higher values, it is not working due to amount of time needed (in the order of days, for me). The duplication part is ok, because when I set q=0 I can generate networks for values of p from 0.1 to 1.0. But when I have to introduce the divergence in my models, it is not satisfatory. I already built the divergence code with some variations, but the results were the same.

  • $\begingroup$ The divergence step is not clear to me: Is it possible that {i,j} and {i,jnew} both remain? Or is it {i,j} with probability q and {i,jnew} with probability 1-q that remains? $\endgroup$ Mar 15, 2018 at 17:18
  • $\begingroup$ @HenrikSchumacher yes, it is possible that they remain. The probability q is the probability that one random link will be deleted or no. $\endgroup$ Mar 15, 2018 at 17:28

1 Answer 1


Well, this is my interpretation of the algorithm. It is still not very efficient, but it contains some minor improvements: First of all, removed the bookkeeping for vdeg since it hasn't been used. Second, I tried to generate as many random numbers at once as possible.

However, I could not get rid of any Append and Delete commands; they should be the main reasons for the slow performance.

f[n_, p_, q_] /; n >= 3 := 
 Module[{g, i, j, gi, ginew, prand, qrand, irand, coinflips},
  irand = 1 + RandomInteger /@ Range[0, n];
  prand = RandomReal[{0, 1}, n];
  g = ConstantArray[{}, n];
  g[[1]] = {2};
  g[[2]] = {1};
   i = irand[[inew-1]];
   gi = g[[i]];
   ginew = gi;
   qrand = RandomReal[{0, 1}, Length[gi]];
   coinflips = RandomInteger[{0, 1}, Length[gi]];
    j = gi[[k]];
    If[qrand[[k]] <= q,
     If[coinflips[[k]] == 1,
       g[[j]] = g[[j]] /. i -> inew;
       gi = Delete[gi, k];
       ginew = Delete[ginew, k];
     g[[j]] = Append[g[[j]], inew];
    , {k, Length[gi], 1, -1}];

   If[prand[[inew]] <= p,
    AppendTo[gi, inew];
    AppendTo[ginew, i];

   g[[i]] = gi;
   g[[inew]] = ginew;
   {inew, 3, n}];

   Sort /@ Join @@ Table[Thread[{k, g[[k]]}], {k, 1, n}]]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.