Dans cet article, je vais vous présenter les points suivants :
Définir les fractales de Koch.
Donner un programme C++ permettant de générer des fractales de type Koch de façon générique.
Expliquer les facteurs clés du code d'un point de vue mathématique, algorithmique et sémantique.
Note : Cette page sera lourde à charger car elle contient des fichiers SVG dont la taille peut atteindre 3,18 Mo au maximum.
Définition
Les fractales de Koch, ou d'une façon plus générale les fractales, sont des objets mathématiques dont la structure reste identique à toutes les échelles. C'est-à-dire que si l'on observe cette courbe (qui peut s'étendre à des objets mathématiques plus abstraits) sur une échelle locale bien définie, on constatera que l'objet initial se répète. Pour que cette propriété reste valable, il est nécessaire d'avoir une répétition à l'infini.
Propriétés
On dispose de propriétés importantes à connaître qui vont nous aider à comprendre et à implémenter le code :
Le nombre de segments :
`N = K*(P+1)^O`
Avec :
`N` : le nombre de segments.
`O` : l'ordre du polygone central.
`P` : l'ordre (la profondeur) de la fractale.
Cette valeur est calculée à l'aide de l'algorithme "Fast Exponentiation" (exponentiation rapide). Celui-ci calcule la puissance d'un nombre entier élevé à une puissance entière avec un coût algorithmique proportionnel au logarithme de la puissance.
Cette fractale a une structure linéaire fermée ou ouverte. Elle ne possède pas de branches comme dans les arbres ou les graphes. Sur cette base, on peut calculer le nombre de points qui est :
`Q = N + 1`
Le premier point et le dernier point sont identiques dans le cas d'une courbe fermée.
À chaque itération, on effectue un saut (dans le tableau contenant les points) calculé par :
`K*(P+1)^(O-D)`
Où D représente la profondeur actuelle.
Compilation de code
Pour compiler le code, on a besoin d'exécuter la commande suivante:
g++ code.cpp -fopenmp -o code -std==c++20
g++: est le compilateur pour le langage c++
-fopenmp: est un option qui ordonne au compilateur de respecter les standards openmp pour ce qui concerne la directive #pragma omp parallel for
-std=c++20: est utile seulement dans le contexte d'utilisation de fonctionnalités compatibles avec le standard c++ de 2020 comme pour le cas de std::format
Solution
#include<format>
#include<iostream>
#include<vector>
#include<math.h>
#include<fstream>
#include<iomanip>
#include<string>
#include<chrono>
using namespace std;
typedef unsigned int index_;
typedef unsigned int size_;
typedef size_ point_index_;
typedef size_ array_inc;
typedef unsigned char byte_;
template<typename T>
struct point
{
T x, y;
};
typedef double angle_type;
struct rotator
{
angle_type angle;
double ca, sa;
void set_angle(angle_type angle)
{
this->angle = angle;
this->ca = cos(angle);
this->sa = sin(angle);
}
template<typename T>
inline point<T> rotate(point<T>& a)
{
return {this->ca * a.x - this->sa * a.y, this->sa * a.x + this->ca * a.y};
}
template<typename T>
inline point<T> rotate(point<T>&& a)
{
return {this->ca * a.x - this->sa * a.y, this->sa * a.x + this->ca * a.y};
}
};
template<typename T>
point<T> operator+(const point<T>& a, const point<T>& b)
{
return {a.x + b.x, a.y + b.y};
}
template<typename T>
point<T> operator-(const point<T>& a, const point<T>& b)
{
return {a.x - b.x, a.y - b.y};
}
template<typename T>
point<T> operator*(T a, const point<T>& b)
{
return {a * b.x, a * b.y};
}
template<typename T>
point<T> operator*(const point<T>& b, T a)
{
return a * b;
}
template<typename T>
ostream& operator << (ostream& os, point<T>& b)
{
os << "{ " << setw(12) << setprecision(4) << b.x << " , " << setw(12) << setprecision(4) << b.y << " }" << endl;
return os;
}
long fast_power(long base, long exponent)
{
if (exponent < 0) return 0;
long result = 1;
while (exponent)
{
if (exponent & 1)
{
result *= base;
}
base *= base;
exponent = exponent >> 1;
}
return result;
}
class koch_fractal
{
private:
byte_ initial_segment_number, new_element_order, fractal_order;
size_ segment_number;
size_ point_number;
double initial_fractal_raduis, central_angle, factor;
vector<point<double>> points;
rotator rot;
void update_segment_number()
{
this->segment_number = this->initial_segment_number * fast_power(this->new_element_order + 1, this->fractal_order);
}
void create_initial_curve()
{
if (this->initial_segment_number == 1)
{
this->points[0] = {-this->initial_fractal_raduis, 0};
this->points[this->segment_number] = {this->initial_fractal_raduis, 0};
return ;
}
const size_ inc = this->segment_number / this->initial_segment_number;
const double angle_factor = 2 * M_PI / this->initial_segment_number;
for(int i = 0; i <= this->initial_segment_number; i++)
{
this->points[inc * i] = {this->initial_fractal_raduis * cos(angle_factor * i), this->initial_fractal_raduis * sin(angle_factor * i)};
}
}
inline void fractalize_segment(point_index_ begin, point_index_ end, array_inc inc)
{
point<double> A = this->points[begin], B = this->points[end];
begin += inc;
end -=inc;
this->points[begin] = (A + (B - A) * this->factor);
this->points[end] = (B + (A - B) * this->factor);
this->points[begin + inc] = (this->points[begin] + this->rot.rotate(this->points[end] - this->points[begin]));
begin += inc;
for (point_index_ i = begin; i < end; i+= inc)
{
this->points[i + inc] = (this->points[i] + this->rot.rotate(this->points[i - inc] - this->points[i]));
}
}
public:
koch_fractal(byte_ isn, byte_ neo, byte_ fo, double raduis, double factor)
{
this->initial_segment_number = isn;
this->new_element_order = neo;
this->fractal_order = fo;
this->factor = factor;
this->update_segment_number();
this->point_number = this->segment_number + 1;
this->initial_fractal_raduis = raduis;
this->central_angle = 2 * M_PI / this->new_element_order;
this->rot.set_angle(M_PI + this->central_angle);
this->points = vector<point<double>>(this->point_number);
}
void compute()
{
this->create_initial_curve();
size_ limit = this->initial_segment_number;
for(int i = 0; i < this->fractal_order; i++)
{
size_ inc = this->segment_number / limit;
#pragma omp parallel for
for (int j = 0; j < this->segment_number; j += inc)
{
this->fractalize_segment(j, j + inc, inc / (this->new_element_order + 1));
}
limit *= (this->new_element_order + 1);
}
}
friend ostream& operator << (ostream& os, koch_fractal frac);
void generate_svg(fstream& file_object)
{
string lines = "";
for (point_index_ i = 1; i < this->point_number; i++)
{
lines += format("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"blue\" stroke-width=\"0.25\" />\n", this->points[i-1].x, this->points[i-1].y, this->points[i].x, this->points[i].y);
}
file_object << format(R"(<svg width="500" height="500" viewBox="-250 -250 500 500" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
{}
</svg>)", lines);
}
};
ostream& operator << (ostream& os, koch_fractal frac)
{
for (point_index_ i = 0; i < frac.point_number; i++)
{
os << frac.points[i];
}
return os;
}
double mean(vector<int> v)
{
double m = 0;
for (int i =0; i < v.size(); i++)
{
m += v[i];
}
return m / v.size();
}
double stddev(vector<int> v, double mean)
{
double m = 0;
for (int i =0; i < v.size(); i++)
{
double temp = v[i] - mean;
m += (temp * temp);
}
return sqrt(m / v.size());
}
int main()
{
koch_fractal frac(1, 3, 7, 200.0, 1.0/3);
constexpr unsigned int limit = 1;
vector <int> times(limit);
for (int i =0; i < limit; i++)
{
auto start = chrono::system_clock::now();
frac.compute();
auto end = chrono::system_clock::now();
times[i] = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
}
double m = mean(times);
cout << "Le temps d'exécution moyen:\t" << m / fast_power(10,6) << " millisecondes" << endl;
cout << "Le variation standard du temps d'exécution:\t" << stddev(times, m) / fast_power(10,6) << " millisecondes" << endl;
// cout << frac;
fstream file("out.svg", ios_base::openmode::_S_out);
frac.generate_svg(file);
return 0;
}
Explication
Parallélisme du code
Un fait important que l'on peut déduire de la définition des fractales est leur propriété d'être embarrassingly parallel (parallélisme massif). Le fait que la structure se répète de la même façon à tous les niveaux nous permet de créer une seule fonction qui applique un seul ordre de fractalisation sur chaque segment. De ce fait, on peut exécuter le code sans crainte d'interdépendance.
La fonction qui s'occupe de cet aspect est la suivante :
inline void fractalize_segment(point_index_ begin, point_index_ end, array_inc inc)
{
point<double> A = this->points[begin], B = this->points[end];
begin += inc;
end -=inc;
this->points[begin] = (A + (B - A) * this->factor);
this->points[end] = (B + (A - B) * this->factor);
this->points[begin + inc] = (this->points[begin] + this->rot.rotate(this->points[end] - this->points[begin]));
begin += inc;
for (point_index_ i = begin; i < end; i+= inc)
{
this->points[i + inc] = (this->points[i] + this->rot.rotate(this->points[i - inc] - this->points[i]));
}
}
Cette fonction utilise la directive #pragma omp parallel for, qui permet d'exécuter la boucle for sur plusieurs threads physiques et plusieurs cœurs, en se basant sur le fait que cette fonction est embarrassingly parallel.
Un fait important à noter : les zones mémoire de calcul de chaque exécution de la boucle for n'ont pas d'intersection, ce qui nous permet d'éviter les data race (concurrence de données) et les race conditions (situations de compétition).
Objet de calcul
Comme je l'ai mentionné, la conception de ce programme est modulaire. Chaque partie est gérée par une classe. Par exemple, les points et le calcul de rotation sont gérés par les deux classes point et rotator. Puisque le calcul de rotation est constant tout au long de l'exécution, les valeurs de cos et de sin sont calculées une seule fois et réutilisées au besoin.
Le fait d'utiliser des classes est utile principalement si l'on souhaite calculer ou créer plusieurs fractales avec différents ordres du polygone central en parallèle. Dans le cas contraire, on pourrait utiliser des variables globales.
Ici, je n'ai pas accordé une importance majeure à l'aspect privé, public ou protégé des variables membres. De même, il n'y a pas eu recours à l'héritage ni à une nécessité de protection stricte des variables internes.
struct point
{
...........
};
struct rotator
{
...........
};
struct rotator
Cette structure permet de calculer les positions des points définissant le polygone central. Les variables qui restent constantes tout au long de l'exécution du code sont sauvegardées dans des variables internes à la structure. Comme dans les autres cas, la fonction de calcul est alignée (elle est insérée directement à l'endroit de son appel via inline).
Pour le passage de paramètres, on utilise const & (référence constante). Cela permet d'effectuer le calcul avec des rvalue et des lvalue.
Par exemple, rot.rotate(this->points[end] - this->points[begin]) est appelé avec une rvalue comme argument.
Celle-ci représente la notion de vecteurs au sein de l'algèbre linéaire avec des opérateurs surchargés permettant de réaliser le calcul vectoriel et la multiplication par un scalaire. Les opérateurs surchargés acceptent des références à des constantes, ce qui permet d'économiser la copie des données de points pour les types de données les plus lourds.
Calcul des points du polygone central
Figure 1 : Forme de base du flocon de Koch avec un triangle au centre.
Figure 2 : Forme de base du flocon de Koch avec un hexagone au centre.
En se basant sur la Figure 1, le calcul de chaque point du polygone repose sur une méthode récursive illustrée par la formule suivante :
`P_(n+1) = P_n + R(theta) * (P_(n-1) - P_n)`
Avec :
`theta` : l'angle défini par `/_P_(n-1)P_(n)P_(n+1)` et extérieur au polygone (par exemple dans la figure `/_ABC`).
Pour la première itération (quel que soit le polygone), on a `P_(n-1) = A, P_(n) = B, P_(n+1) = C`.
Dans le cas d'un polygone de degré 6, pour la deuxième itération, on a `P_(n-1) = B, P_(n) = C, P_(n+1) = F`.
L'inconvénient de cette méthode est qu'elle est récursive et que, pour des ordres plus élevés du polygone, on obtient un cumul d'erreurs au fur et à mesure des itérations.
Exportation dans un fichier SVG
L'exportation du code dans un fichier SVG n'est pas optimale en termes de temps d'exécution, mais elle remplit sa fonction d'exportation de données. L'exportation n'est pas la partie la plus critique du projet, sauf s'il y a une exportation massive de chaque courbe générée, ce qui deviendrait pénalisant.
Cependant, les points d'amélioration sont les suivants :
Le formatage des données avec un nombre à largeur fixe (ex: 4 chiffres avant et après la virgule selon l'ordre).
Une allocation statique de la chaîne de caractères finale, de la même manière que pour la génération des fractales.
Être le plus prévisible possible pour le compilateur.
Je vais lister seulement quelques images des fractales générées à partir du code. La dernière fractale s'auto-intersecte : les traits foncés correspondent aux segments redondants. La seule structure qui ne présente pas cet aspect est la fractale originelle de Koch.
Figure 3 : Fractale de Koch avec un polygone de base de degré 1, un polygone secondaire de degré 3 et un ordre de division 7.
Figure 4 : Fractale de Koch avec un polygone de base de degré 2, un polygone secondaire de degré 3 et un ordre de division 7.
Figure 5 : Fractale de Koch avec un polygone de base de degré 2, un polygone secondaire de degré 6 et un ordre de division 4.
Conclusion
Finalement, ce code a permis de générer les fractales de Koch et de les exporter dans un fichier SVG. De plus, nous avons vu que l'aspect embarrassingly parallel facilite grandement la parallélisation du code tout en évitant de nombreux problèmes de gestion de données.
Il faut toutefois s'assurer d'utiliser les bons mécanismes du langage, tels que le passage par référence et l'alignement des données. L'utilisation d'une allocation statique des données permet de ne pas surcharger l'allocateur dynamique. Il est également recommandé d'utiliser la spécification des patrons (templates) afin d'éviter d'avoir un code trop générique, laissant ainsi le compilateur choisir la fonction optimisée correspondante.
Dans un prochain article, je réaliserai le profilage de ce code en analysant ses performances sous plusieurs angles.
