Marc Joiret 's web space Recreations and computer challenges in C++, Python and JavaScript e-mail icon. Opening mail in progress...Please wait.

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 :

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.
You can download here the complete JavaScript code I used ('sketch.js').

Computer challenge #2 : The Mandelbrot set: the classical fractal figure.

This computer challenge is about drawing the famous Mandelbrot set, the archetypical fractal object. The Mandelbrot set is the green figure you can see here in the background image. We want to combine Object Oriented Programming in C++ with the use of a graphic format called .ppm (portable pixmap format). This format is fully described on the web. Basically, a ppm file ("file.ppm") contains a three lines header and data triplets representing the color of each pixel to be drawn in RGB coding format. The three header lines are "P3" (first header line which is the format definition with end of line character), two integer numbers like 500 500 (second header line with the number of pixels in width and the number of pixels in height also with the end of line character), an integer number (third header line which is the maximum color value used in RGB coding also with the end of line character). The rest of the file contains all data triplets for the pixels. For instance the RGB (red green blue) triplet for pixel one could be 255 0 0 (the pixel would be red) or 0 0 0 (black pixel) or 255 255 255 (white pixel). A ppm file cannot be read by all graphical softwares. You will need, for instance, the GIMP software to read such files and most importantly to display the image. GIMP stands for GNU image manipulation program and is open source. You are advised to download the GIMP software before you start this computer challenge. Once you will have produced a ppm file with your C++ program, you will be able to display the image and to export it in other image formats (like .jpg or .png). You will have to prepare or use an existing header file (.h) in your C++ project implementing a Complex class (dealing with attributes, i.e. properties, of complex number ∈ ℂ, like their real and imaginary parts and their methods: for instance how the modulus of a complex number is calculated, or how two complex numbers are summed up or multiplied, how do you compute the power of a complex number, etc...). Finally, in the main C++ program, you will have to instanciate your complex number object. The procedure to determine if a complex number belongs or not to the Mandelbrot set is fully documented in Wikipedia. To spell it out simply: a complex number c ∈ ℂ belongs to the Mandelbrot set, if starting from z0 = 0, taking the second power of z, adding c to get the next z; then repeating this iteratively a number of times n to infinity, and finally checking if all z do not diverge. If they don't diverge or are all bounded, then c belongs to the Mandelbrot set. Equivalently, if the modulus of z remains smaller or equal to 2, whatever the number of iterations, then c belongs to the set. Of course, you should fix the number of iterations to a maximum integer value (256 for example in our code).

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.

Mandelbrot set figure

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++ :
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"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"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 }