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 with 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.
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 }
