# How can I improve this code?

The code below is working well, but it's very slow in the Manipulate box. I know that the Do and AppendTo parts aren't very efficient. Is there a better way in doing the same calculation? Take note that the order of the three Do/AppendTo is important, since the result isn't the same if the order is changed. I also need to know if these Do operations are well declared in that code, since Mathematica doesn't format the code in a nice way, while it doesn't give any error message. Please, take note that I'm using a very old version of Mathematica (7.0), so any suggestion should be compatible with old versions of Mma.

(* pts is initially an uniform distribution of "p" points in space, that I call "clusters". *)
(* Then I add randomly "q" points (galaxies) to that distribution of clusters. *)
(* The second and third Do add 5000 + 8000 more galaxies, to get irregular clusters of thousands of galaxies (the dots in the simulation). *)

galaxies[p_, q_, r_] := Module[
{pts = RandomReal[{-1, 1}, {p, 3}]},
Do[AppendTo[pts, RandomChoice[pts] + 1.00 r RandomReal[{-1, 1}, 3]],
{i, q}];pts
Do[AppendTo[pts, RandomChoice[pts] + 0.33 r RandomReal[{-1, 1}, 3]],
{i, 5000}];pts
Do[AppendTo[pts, RandomChoice[pts] + 0.11 r RandomReal[{-1, 1}, 3]],
{i, 8000}];pts
]

view[p_, q_, r_] := Show[
Graphics3D[{RGBColor[{0.5, 0.4, 1.0, 0.4}], PointSize[0.003], Point[galaxies[p, q, r]]}],
PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}},
Boxed -> True,
Background -> Black,
ImageSize -> {700, 700},
SphericalRegion -> True,
Method -> {"RotationControl" -> "Globe"}
]

Manipulate[
view[p, q, r],
{{p, 500, "p"}, 1, 1000, 1},
{{q, 10000, "q"}, 0, 10000, 1},
{{r, 0.05, "scale"}, 0.01, 0.2, 0.001}
]


Preview of what this code is doing:

• Please add comments to your galaxy code that explain what each line does. It would then be a lot easier for people to think of a new implementation, rather than trying to figure out what you did and fix it from there. Having to understand what you want from suboptimal code is a big barrier. Commented Aug 2, 2022 at 15:16
• Your code on Mathematica version 13.0.1.0 running on an Apple MacBook Pro M1 runs effortlessly, even with all three Manipulate controls in animation mode. Maybe it's time to upgrade Mathematica and/or your machine. Commented Aug 2, 2022 at 15:17
• @Jagra, yes, I agree. My version of Mma is very old (7.0), and can't change it until I upgrade the whole computer (maybe in october). So you say that the Manipulate box runs smooth on your setup? Please, can you post an animated gif?
– Cham
Commented Aug 2, 2022 at 15:21
• No time to create and post a gif. That said, your code might run in the [Wolfram Cloud[(wolframcloud.com). I think the processing would run on Wolfram's servers, you'd just see results in a browser (if memory serves). Basic Plan = free. Commented Aug 2, 2022 at 15:56
• Just tried it on Wolfram Cloud. It doesn't do animation but it does seem to work. Commented Aug 2, 2022 at 16:04

addPoints[pts_, num_, scale_] :=
Join @@ {pts,
RandomChoice[pts, num] + RandomReal[{-scale, scale}, {num, 3}]};

galaxiesNew[p_, q_, r_] := Module[
{pts = RandomReal[{-1, 1}, {p, 3}]},
pts = Nest[addPoints[#, 1, r] &, pts, q];
pts = Nest[addPoints[#, 1, 0.33*r] &, pts, 5000];
pts = Nest[addPoints[#, 1, 0.11*r] &, pts, 8000];
pts
]


Getting rid of the Do loops and AppendTo speeds up galaxy generation by about 30%.

First@AbsoluteTiming[galaxies[500, 10000, 0.05]]
First@AbsoluteTiming[galaxiesNew[500, 10000, 0.05]]


Gives me 0.47s v. 0.35s.

Edit: I missed the fact that the Do loop was updating pts before every call to RandomChoice meaning each subsequent call should have access to previously generated points. I accommodated this using Nest and only adding one point at a time, the increased number of calls to Join and RandomReal, and RandomChoice mean the speedup isn't so great.

• Apologies, I must have missed something. I'll take a look. Commented Aug 2, 2022 at 16:11
• Should provide the same functionality now, but the speed up isn't as great. Commented Aug 2, 2022 at 16:41
• Why the num variable in your definition of addPoints? Why it is set to 1 everywhere else?
– Cham
Commented Aug 2, 2022 at 16:52
• That's because RandomChoice[pts] returns {x,y,z}, whereas RandomChoice[pts,1] returns {{x,y,z}}, which matches dimensions with RandomReal[{-1,1},{1,3}]. If you want to use RandomChoice[pts] you need to use Join@@{pts,{RandomChoice[pts]+RandomReal[{-scale,scale},3]}}. Commented Aug 2, 2022 at 17:51
• Absolutely, you just change the second parameter of addPoints to 10. But be careful, this will change the behavior of the code as those ten points will all choose from the same set of points only the next batch of ten will have the last 10 integrated. You will also probably need to scale q, 5000, 8000, by 1/10 to get the same number of points. Commented Aug 2, 2022 at 20:09

First, get rid of AppendTo and replace it with a preallocation of the array because you know the size in advance.

Second, use RandomInteger instead of RandomSample.

galaxiesFaster[p_, q_, r_] := Module[{t1 = 5000, t2 = 8000, pts},
pts = ConstantArray[{0., 0., 0.}, p + q + t1 + t2];
pts[[1 ;; p]] = RandomReal[{-1., 1.}, {p, 3}];
Do[pts[[p + i]] =
pts[[RandomInteger[{1, p + i - 1}]]] +
1.00 r RandomReal[{-1., 1.}, 3], {i, q}];
Do[pts[[p + q + i]] =
pts[[RandomInteger[{1, p + q + i - 1}]]] +
0.33 r RandomReal[{-1., 1.}, 3], {i, t1}];
Do[pts[[p + q + t1 + i]] =
pts[[RandomInteger[{1, p + q + t1 + i - 1}]]] +
0.11 r RandomReal[{-1., 1.}, 3], {i, t2}];
pts];

galaxiesCompiled = Compile[{{p, _Integer}, {q, _Integer}, r},
Module[{t1 = 5000, t2 = 8000,
pts = ConstantArray[{0., 0., 0.}, p + q + t1 + t2]},
pts[[1 ;; p]] = RandomReal[{-1., 1.}, {p, 3}];
Do[pts[[p + i]] =
pts[[RandomInteger[{1, p + i - 1}]]] +
1.00 r RandomReal[{-1., 1.}, 3], {i, q}];
Do[pts[[p + q + i]] =
pts[[RandomInteger[{1, p + q + i - 1}]]] +
0.33 r RandomReal[{-1., 1.}, 3], {i, t1}];
Do[pts[[p + q + t1 + i]] =
pts[[RandomInteger[{1, p + q + t1 + i - 1}]]] +
0.11 r RandomReal[{-1., 1.}, 3], {i, t2}];
pts], RuntimeOptions -> "Speed"
];

galaxies[500, 10000, 0.05] // RepeatedTiming // First
(* 0.8353 *)

galaxiesFaster[500, 10000, 0.05] // RepeatedTiming // First
(* 0.2402 *)

galaxiesCompiled[500, 10000, 0.05] // RepeatedTiming // First
(* 0.0309 *)


The first two fix-ups (galaxiesFaster) produce approximately 3.5 times faster executions, and then compiling the code (galaxiesCompiled) is furthermore 7.5 times faster.

In version 7, there is no RuntimeOptions option for Compile, so you have to remove this.

• very nice solution. Commented Aug 2, 2022 at 17:08
• How do I insert your code into mine (in my question above)? I must admit that I don't understand your code at all! :-(
– Cham
Commented Aug 2, 2022 at 17:14
• @Cham, instead of using galaxies, use galaxiesCompiled. So for example: view[p_, q_, r_] := Show[ ..., Point[galaxiesCompiled[p, q, r]] ...]. Commented Aug 2, 2022 at 17:40
• @Domen, Even if you remove the compilation this is still about 5x faster than the original just using a pre-allocated array. Considerably smoother in Manipulate. Commented Aug 2, 2022 at 18:25
• @N.J.Evans, thanks! I have added this to my answer. Commented Aug 2, 2022 at 19:13