Marc Joiret 's web space Recreations and computer challenges in C++, Python and JavaScript

Below you will find recreational computer challenges. Their purpose is to introduce concepts in object oriented programming (OOP), famous numerical techniques, famous bioinformatics algorithms or simply help you dive in the current web technologies. The languages used will be C++, Python and possibly JavaScript associated to HTML5.
For the animations and interactivity that you will discover below, we used p5.js libraries.
p5.js is a recent interpretation of Processing written in JavaScript that makes it easy to interact with HTML5 objects, including text, input, video, webcam and sound. p5 is an interpreted programming language which promotes processing and sketching as an easy learning programming method.
Each computer challenge is divided into three sequential frames :

• The stated problem;
• The solution to the problem and results;
• A proposed code implementing the way to solve the problem.

Of course, it is recommended that you try to solve the problem on your own and find a better solution than the one you will find here :

# Computer challenge #1 : The Diffusion Limited Aggregation algorithm.

This problem is known as the DLA, « Diffusion Limited Aggregation », and is fully documented on Wikipedia. It was first proposed in 1981 by Witten and Sander in an article entitled « Diffusion limited aggregation, a kinetic critical phenomena » in Physical Review Letters, number 47. In DLA, model particles are added, one at a time, to a growing aggregate (cluster) of particules via random walk trajectories. The particules are imagined to have come from infinity but are actually lauched from a random position on a contour which encloses the seed aggregate. At some point in the random walk of the particle, if the particle comes at a close distance from the aggregate, it sticks to it. The aggregate has now grown larger and a new random particle is launched. The term « Diffusion » means the particules forming the structure wander around randomly before attaching (« Aggregating ») to the structure. It is called « Diffusion Limited » because the particules are considered to be at low concentrations so they don't interact which each other and the structure grows one particule at a time and not by chunks of pre-aggregated particules like it would be the case in coagulation and flocculation. Practical examples are found in coral growth, snow crystal growth, percolation phenomena, lightning paths, etc. There are a number of ways of simulating this process by computer, perhaps the most common is to start with a white image except for a single black pixel in the center. New points are introduced at the border contour and randomly (approximation of Brownian motion) walk until they are close enough to stick to an existing black pixel. In a first simple version, the random walk can be simulated in a 2-D lattice where a drunk ant has a clock and a coin. At each node position of the ant, and at each tick of the clock, the ant toss a fair coin twice. At the first toss, if it is heads or tails, the ant will move horizontally or vertically; at the second toss, the ant will move in one direction or the other. At subsequent ticks of the clock, the ant do it again and again until eventually, it reaches the aggregate. To speed up the algorithm a bit, if the ant finds itself in a position very far away from the aggregate, the ant is annihilated and a new one is lauched from the border contour. In a more general version of the algorithm, the 2D-lattice simplification is relaxed and the ant may now move in any direction : this is the off-lattice algorithm.

# Hints for a solution and results

What is the final result you have to achieve ? Bring your mouse cursor inside the left frame just below and you will see what Diffusion Limited Aggregation really is. The challenge is to reproduce this animation. Additional questions of interests might be the following :

• Find out the growth dynamic of the aggregate, i.e. what is the relation between the number of pixels added to the aggregate and the elapsed time or the relation between the weight of the aggregate and the vertical (horizontal) distance between the most distant two pixels belonging to the aggregate ?
• What is the fractal dimension of the aggregate ?
• What is the variation in entropy for the system (system = random walkers + aggregate) during the build-up of the aggregate ?
• On what assumptions depend the expression of the entropy ?
• Can you achieve a three dimensional extension of this algorithm and its graphical representation ?
To gain further insights on the dynamic geometrical properties of the aggregate while being built up, please click the mouse inside the animated frame. You will get elapsed time, aggregate box lengths, aggregate count and weight.
Collecting these data sequentially and plotting them in a log-log scale will provide the fractal dimension as the slope of curve (log of aggregate density versus log of aggregate length geometric mean).
A deeper treatment of DLA fractal properties can be found in Bunde, A. and Havlin, S. Fractals and Disordered Systems, Springer-Verlag, 1991.

# Code extracts in JavaScript with p5.js libraries that produced the solution above:

The animation above shows how powerful Javascript and the P5 associated libraries are in interacting with HTML5 to produce nice dynamical pictures which also can interact with the web end-user.
To set up the code for the DLA algorithm, it is recommended to build the random walker as an object. Then, you will probably want to create a bunch of such random walkers by setting up an array of this object. The random walker attributes are up to you but you certainly will have to include the walker position, its size (radius), a stuck status (is it or not stuck in an existing aggregate?), its color, the time elapsed before it stuck, etc... The methods in the object should include a display function, a random move function, a boolean function checking for the stuck status after a move. Make sure there is no overlapping between stuck random walkers : each pair of tangent stuck walkers may not overlap any existing stuck walker.
Now, once a walker has been stuck to the aggregate, it disappears (use the .split() method in JavaScript). Doing so, the aggregate is enriched by the stuck walker (use the .push() method in Javascript to append the array of the previously stuck walkers which is also built indeed as an array of walkers with special attributes).
The algorithm and code for this animation was freely inspired from nice P5.js tutorials shown on Youtube by Dan Shiffman.

# Hints for a solution and results

Prepare a C++ header file to create a complex class with properties (attributes) relevant to complex numbers (real, imaginary). Implement methods for this Complex class : modulus of a complex number, and more if you want... In your C++ main(), read from a .txt file, the limits of the Complex plane you want to investigate (see a figure of the Mandelbrot set in Wikipedia : x range from -1.5 to 0.7 and y range from -1.0 to 1.0), the maximum number of iterations you will carry out for the Mandelbrot set procedure (here we fixed 256). In this .txt file also get the relevant information on the size, i.e. number of pixels, you want to output for your final ppm file to be produced : for instance 512 times 512 = 262,144 pixels for the entire produced image. Actually, we advise you to have 7 pieces of information in your input file (.txt) : ImageWidth, ImageHeight, MaxValueofN in the Mandelbrot formula, minReal, maxReal, minImaginary, maxImaginary. Create appropriate functions to sum two complex numbers, to multiply two complex numbers, to take the power of a complex number if these operations are eventually useful. Map the real and imaginary part of c in the .ppm file for all c scanned (for loops) from the input section of the complex plane. The rgb triplet values of a pixel can be varied mathematically for different colors results depending on the observed number of iterations that were reached before z diverges (while loop in the function maxNumberOfIterationsMandelbrot). The header lines for the ppm format are written with all data triplets for all pixels in an 'output_image.ppm' file that you will eventually open with GIMP to display the result. # Code extracts that produced this drawing:

Part of the code was freely borrowed from a nice and fast youtube video about programming the Mandelbrot set in C++ : https://www.youtube.com/watch?v=entjSp3LIfQ.
The relevant input information you might need and you can read in an input file called "inputMandelbrot.txt" are :

```		512 512 256 -1.5 0.7 -1.0 1.0
There are 7 inputs : ImageWidth, ImageHeight,
MaxValueofN in Z_n (Mandelbrot formula),
minReal, maxReal, minImaginary, maxImaginary
```

``` The Complex header file implements the class Complex, is saved in a "complex.h" header file to be included eventually in the C++ main program :
``` 1  #ifndef COMPLEX_H
2  #define COMPLEX_H
3
4  class Complex
5      {
6          private :
7              double a, b; // these are member variables for a complex number
8          public : //prototypes for the member functions
9              void setReal(double);
10              void setImag(double);
11              double getNorm(); // returns the norm (module) of the complex number
12              double getReal(); // give the real part
13              double getImag(); // give the imaginary part
14              double getTheta(); // angle phase (argument of complex)
15      };
16
17  //-------------start-----------
18      //ctor
19
20  void Complex::setReal(double re)
21              {
22                  a = re;
23              }
24
25  void Complex::setImag(double im)
26              {
27                  b = im;
28              }
29  double Complex::getReal()
30  {
31      return a;
32  }
33
34  double Complex::getImag()
35  {
36      return b;
37  }
38  double Complex::getNorm()
39  {
40      return sqrt(a*a + b*b);
41  }
42  double Complex::getTheta() // angle phase (argument of complex)
43              {
44                  const double PI = 3.14159265359;
45                  double theta;
46                  if (a == 0 && b == 0){theta = 0;}
47                  if (a == 0 && b > 0){theta = PI/2;}
48                  if (a == 0 && b < 0){theta = PI + PI/2;}
49                  if (a > 0 && b == 0){theta = 0;}
50                  if (a < 0 && b == 0){theta = PI;}
51                  if (a >0 && b >=0){ // premier quadrant :
52                          theta = acos(a/getNorm());}
53                  if (a < 0 && b > 0){ // deuxieme quadrant :
54                          theta = PI - acos(-a/getNorm());}
55                  if (a < 0 && b < 0){ // troisieme quadrant :
56                          theta = PI + acos(-a/getNorm());}
57                  if (a >=0 && b < 0){ // quatrieme quadrant :
58                          theta = 2*PI- acos(a/getNorm());}
59                  return theta;
60              }
61
62  //Complex::~Complex()
63  //{
64      //dtor
65  //}
66  //----------------------- end ---------------------------
67
68  #endif // COMPLEX_H
``````
```  The main() part of the program is displayed hereafter and uses the above complex.h header file :
```  1  #include <cmath>
2  #include <vector>
3  #include <iomanip>
4  #include <iostream>
5  #include <fstream>
6  #include <string>
7  #include <ctime>
8  #include "complex.h"
9
10  using namespace std;
11
12  const double PI = 3.14159265359;
13
14  // Function Prototype declaration :
15  double mapToReal(int x, int imageWidth, double minR, double maxR);
16  double mapToImaginary(int y, int imageHeight, double minI, double maxI);
17  Complex power(Complex c, int n);          // power of a complex number prototype declaration
18  Complex sum(Complex, Complex);            // operations defined on two complex numbers (sum and
19  Complex prod(Complex, Complex);           // product of 2 complex numbers).
20  int maxNumberOfIterationsMandelbrot(Complex c, int Nmax);
21
22  //main
23  int main()
24  {
25      double cr; // real part of c.
26      double ci; // imaginary part of c.
27      Complex c; // c, is a complex numbers.
28      int n;
29
30      // get the required values from an input file called 'inputMandelbrot.txt'
31      ifstream inFile;
32      inFile.open("C:\\Users\\Marc\\CloudStation\\AQUATION\\STATISTICS\\C++ Code Blocks\\PROJECTS\\LEARNINGc++\\inputMandelbrot.txt");
33      int imageWidth, imageHeight, maxN;
34      double minReal, maxReal, minImaginary, maxImaginary;
35      if (!inFile)
36      {
37          cout << "Could not open file. Make sure the file 'inputMandelbrot.txt' does exist." << endl;
38          cin.ignore();
39          return 0;
40      }
41      inFile >> imageWidth >> imageHeight >> maxN >> minReal >> maxReal >> minImaginary >> maxImaginary;
42      inFile.close(); // close input file when the required values have been read.
43      // cout << "The max number of iteration is : " << maxN << endl;
44
45      // open the PPM output file (ppm format means 'portable pixmap format')
46      // .ppm format has a header with 3 lines ended by a new line character :
47      // P3 endl;
48      // 3 2 endl;
49      // 255 endl;
50      // This header means : P3 means it is a RGB (Red Green Blue) color image in ASCII.
51      // 3 2 means width and height of the image in pixels.
52      // 255 is the maximum value for each color.
53      // After the header, image data are triplets of integer numbers. All triplets are separated by a single space character :
54      //   255 0 0 will produce a red pixel.
55      //   0 255 0 will produce a green pixel, 255 255 0 is a mix of red and green which is yellow, etc...
56      //   0 0 255 will produce a blue pixel.
57      //   0 0 0 a black pixel and 255 255 255 a white pixel.
58      // To read an image file in .ppm format you will need a special open source software like GIMP (GNU image manipulation program).
59      ofstream outFile;
60      outFile.open("C:\\Users\\Marc\\CloudStation\\AQUATION\\STATISTICS\\C++ Code Blocks\\PROJECTS\\LEARNINGc++\\output_image.ppm");
61      outFile << "P3" << endl; // header of ppm file meaning that the pixel data are triplets of RGB color in ASCII code.
62      outFile << imageWidth << " " << imageHeight << endl; // Width and Height of the image in pixels.
63      outFile << "256" << endl; // maximum value of the color scale for each of the three colors.
64
65      // For every pixel...
66      for (int y = 0; y < imageHeight; y++) // for all rows...
67      {
68          for (int x = 0; x < imageWidth; x++) // for all pixels in the row...
69          {
70              // ...Find the real and imaginary values for c, corresponding to that x, y pixel in the image.
71              //use mapTo functions for cr and ci and complex c *******************************************************************************
72              cr = mapToReal(x, imageWidth, minReal, maxReal);
73              ci = mapToImaginary(y, imageHeight, minImaginary, maxImaginary);
74              c.setReal(cr);
75              c.setImag(ci);
76
77              // ...Find the number of iterations in the Mandelbrot formula using said c :
78              n = maxNumberOfIterationsMandelbrot(c, maxN);
79
80              // ...Map the resulting value to an RGB value :
81              int r = ((int)(n*sinf(n))%256); //% 255;//(n*3) % 255; // change for more interesting colors...
82              int g = (n * 3) % 256;//255 - (n % 255);
83              int b = n % 256; //% 255;//255 - (n*2) % 255;
84
85              // ... output it to the ppm file :
86              outFile << r << " " << g << " " << b << " ";
87          }
88          outFile << endl;
89      }
90
91      outFile.close(); // close the created image output file
92      cout << "Mandelbrot set is completed in the PPM format file.\n";
93      cout << "Please open the file 'output_image.ppm' with GIMP tool." << endl;
94      return 0;
95  }
96
97  // Function implementations :
98  // Map the integer pixels to real and imaginary parts of the complex plane:
99  double mapToReal(int x, int imageWidth, double minR, double maxR)
100  {
101      double range = maxR - minR;
102      return x * (range / imageWidth) + minR;
103  }
104  double mapToImaginary(int y, int imageHeight, double minI, double maxI)
105  {
106      double range = maxI - minI;
107      return y * (range / imageHeight) + minI;
108  }
109  // Compute the power of a complex number :
110  Complex power(Complex c, int n) // returns the n power of complex number c
111              {
112                  double rho, arg, re, im;
113                  Complex out;
114                  rho = pow(c.getNorm(), n);
115                  arg = fmod(n*c.getTheta(), 2*PI);
116                  re = rho * cos(arg);
117                  im = rho * sin(arg);
118                  out.setReal(re);
119                  out.setImag(im);
120                  return out;
121              }
122  Complex sum(Complex c1, Complex c2) // returns the sum of two complex numbers
123      {
124          double re, im;
125          Complex out;
126          re = c1.getReal() + c2.getReal();
127          im = c1.getImag() + c2.getImag();
128          out.setReal(re);
129          out.setImag(im);
130          return out;
131      }
132
133  Complex prod(Complex c1, Complex c2) // returns the product of two complex numbers
134      {
135          double re, im;
136          Complex out;
137          re = c1.getReal() * c2.getReal() - c1.getImag() * c2.getImag();
138          im = c1.getReal() * c2.getImag() + c1.getImag() * c2.getReal();
139          out.setReal(re);
140          out.setImag(im);
141          return out;
142      }
143
144  int maxNumberOfIterationsMandelbrot(Complex c, int Nmax) // returns the number of iterations at which z diverges
145  // or the max allowed number of iterations
146      {
147          int iter;
148          Complex zNext;
149          double module;
150          zNext.setReal(0);
151          zNext.setImag(0);
152          iter = 0;
153          do
154          {
155              zNext = sum(power(zNext, 2), c);
156              module = zNext.getNorm();
157              iter += 1;
158          }
159          while (module <= 2.0 && iter <= Nmax);
160          return iter;
161      }
``````