diff --git a/cours/audio/Makefile b/cours/audio/Makefile new file mode 100644 index 0000000..7d29470 --- /dev/null +++ b/cours/audio/Makefile @@ -0,0 +1,44 @@ +HEPTC=heptc +HEPTLIB=$(shell heptc -where) +CC=gcc +CFLAGS=-g $(shell pkg-config --cflags --libs sdl2 sndfile) -I$(HEPTLIB)/c -lm +PYGMENTS=python -m pygments -x + +TARGET=audio +SOURCES=audio_c/audio_types.c \ + audio_c/audio.c \ + buffer.c \ + vcd_lib.c \ + mathext.c \ + vcd.c \ + main.c + +.PHONY: all clean test + +all: $(TARGET) audio.html main.html + +%.html: %.ept + $(PYGMENTS) -l ../../../notes/heptagon.py:HeptagonLexer -O full -o $@ $^ + +%.html: %.c + $(PYGMENTS) -O full -o $@ $^ + +clean: + rm -f $(TARGET) *.{epci,log,mls,obc,html} $(TARGET).{pdf,tex} + rm -rf audio_c + +test: $(TARGET) + ./$(TARGET) + +$(TARGET): $(SOURCES) + @pkg-config --exists sdl2 || \ + ( echo "La bibliothèque SDL2 est absente."; exit 1 ) + @pkg-config --exists sndfile || \ + ( echo "La bibliothèque sndfile est absente."; exit 1 ) + $(CC) $(CFLAGS) -o $@ -I audio_c -I. $^ + +audio_c/audio_types.c audio_c/audio.c: audio.ept main.c mathext.epci vcd.epci + $(HEPTC) -target c $< + +%.epci: %.epi + $(HEPTC) $< diff --git a/cours/audio/README.md b/cours/audio/README.md new file mode 100644 index 0000000..335e71b --- /dev/null +++ b/cours/audio/README.md @@ -0,0 +1,14 @@ +# Synthèse sonore élémentaire en Heptagon + +Ce dossier contient quelques noeuds Heptagon élementaires qui produisent du son. + +Il utilise les bibliothèques SDL2 et sndfile. Elles doivent être installées via +le gestionnaire de paquet de votre système d'exploitation (`apt-get` sous +GNU/Linux Debian ou Ubuntu, `pacman` sous Arch Linux, `brew` sous macOS, etc.). + +```shell +$ make +$ ./audio +``` + +Pour essayer différents codes, il faut éditer le noeud `main` dans `audio.ept`. diff --git a/cours/audio/audio.ept b/cours/audio/audio.ept new file mode 100644 index 0000000..f273bd6 --- /dev/null +++ b/cours/audio/audio.ept @@ -0,0 +1,515 @@ +(******************************************************************************) +(* SYNTHÈSE SONORE ÉLÉMENTAIRE EN HEPTAGON *) +(******************************************************************************) + +(* Le but de ce fichier est de démontrer quelques techniques élémentaires de + génération de son en Heptagon, à travers de modestes expérimentations. Il n'a + bien sûr pas vocation à se substituer à un cours de traitement du signal ou + d'acoustique. En revanche, il peut facilement servir d'illustration de + diverses techniques de programmation en Heptagon. *) + +(* Le flot d'échantillons sonores produits par ce programme synchrone est + branché à un petit bout de code C qui les envoie au système sonore de votre + système d'exploitation par le biais de la bibliothèque SDL2. Votre OS les + transmet à la carte son qui elle même les envoie à vos enceintes, casque ou + écouteurs. *) + +(* Les langages synchrones ont été utilisés pour la synthèse sonore. Si ce sujet + vous intéresse, vous pouvez par exemple consulter la page du langage Faust, à + la syntaxe rudimentaire mais aux bibliothèques acoustiques, sonores et + musicales très développées : http://faust.grame.fr *) + +(* On va utiliser une petite bibliothèque de composants mathématiques. Les + curieuses et curieux pourront aller voir mathext.epci. *) + +open Mathext + +(* Avant de commencer, on a besoin de quelques définitions et outils. *) + +(* Le nombre d'échantillons, c'est à dire ici de pas synchrones, que le système + sonore va consommer une seconde. *) + +const period : int = 44100 + +(* Un signal mono est un simple flot de nombres à virgule flottante. *) + +type mono = float + +(* Un signal stéréo fournit deux échantillons, gauche et droit, la carte son se + chargeant de les mixer pour donner l'impression d'un son 'surround'. *) + +type stereo = { l : float; r : float } + +(* Le signal constant silencieux. *) + +const silence : stereo = { l = 0.0; r = 0.0 } + +(* On peut dupliquer un signal mono pour obtenir un signal stéréo + inintéressant, les deux canaux portant la même valeur. *) + +fun stereo_of_mono(a : mono) returns (o : stereo) +let + o = { l = a; r = a } +tel + +(* On peut appliquer un gain à un signal stéréo, c'est à dire le multiplier par + un flottant pour l'amener à une amplitude différente. *) + +fun stereo_gain(g : float; s : stereo) returns (o : stereo) +let + o = { l = g *. s.l; r = g *. s.r }; +tel + +(* Étant donné deux signaux, on peut les combiner via leur somme. *) + +fun stereo_sum(s1, s2 : stereo) returns (o : stereo) +let + o = { l = s1.l +. s2.l; r = s1.r +. s2.r } +tel + +(* Quand on utilise les deux fonctions qu'on vient de définir, gare à + l'amplitude en sortie ! Une amplitude trop élevée risque de dépasser la + capacité de votre carte son, enceintes ou écouteurs, ce qui cause un + phénomène de saturation : tous les échantillons d'amplitude trop élevée sont + écrasés sur l'amplitude maximale. *) + +(* La fonction mix ci-dessous pallie le défaut de la fonction stereo_sum en + renormalisant le résultat. De plus, elle traite un tableau de signaux, et + donc moralement un nombre d'entrées arbitraires. *) + +fun stereo_mix<>(s : stereo^n) returns (o : stereo) +let + o = stereo_gain(1.0 /. float(n), fold<> stereo_sum(s, silence)); +tel + +(* On peut commencer à écouter un peu de son, par exemple celui du silence. *) + +node main0() returns (o : stereo) +let + o = silence; +tel + +(* Quid du noeud suivant ? *) + +node cracks() returns (o : stereo) +let + o = { l = 4200.0; r = 4200.0 }; +tel + +(* On entendu un craquement, puis plus rien, puis un craquement lorsqu'on + interromp le programme. Pourquoi ? + + Physiquement, le son est une vibration produit par une onde acoustique, c'est + à dire une oscillation de la pression de l'air. Autrement dit, il s'agit + d'une *variation*. Donc, le signal constant ne peut pas donner lieu à un son, + sauf au premier instant (passage de 0 à 4200) puis lorsqu'on interromp le + programme (passage de 4200 à 0). + + Et si on essayait un signal qui varie ? Par exemple, un signal carré qui + passe de 1 à 0 toutes les demi-secondes. *) + +node periodic(p : int) returns (o : int) +var n : int; +let + o = 0 fby (if n = p then 0 else n); + n = o + 1; +tel + +node beats_1() returns (o : stereo) +let + o = stereo_of_mono(if periodic(period) <= period / 2 then 1.0 else -. 1.0); +tel + +(* On obtient une série de battements simples. Faire en sorte que le canal droit + soit l'opposé du canal gauche produit un effet intéressant. *) + +node beats_2() returns (o : stereo) +var l : float; +let + l = if periodic(period) <= period / 2 then 1.0 else -. 1.0; + o = { l = l; r = -. l }; +tel + +(* Essayons maintenant de générer un signal qui croît indéfiniment. *) + +node fcnt(ini : float; step : float) returns (o : float) +let + o = ini fby (o +. step); +tel + +node sawtooth_1() returns (o : stereo) +let + o = stereo_of_mono(fcnt(0.0, 1.0)); +tel + +(* On entend quelques craquements, puis plus rien. Normal : ce signal n'oscille + pas vraiment, ou du moins pas avant d'atteindre l'overflow. Pourquoi ne pas + tester un signal périodique en dents de scie, dans ce cas ? *) + +node sawtooth_2() returns (o : stereo) +var t : float; +let + t = float(periodic(128)); + o = stereo_of_mono(t); +tel + +(* Tiens, un son à peu près constant ! Pas très harmonieux cependant. *) + +(* Est-ce qu'appliquer un gain ferait une différence ? Pour bien observer la + différence, on n'a qu'à faire passer le gain de 0 à 1 à chaque seconde. + C'est très facile à programmer en Heptagon. *) + +node sawtooth_3() returns (o : stereo) +var t : float; g : float; +let + t = float(periodic(128)); + g = float(periodic(period)) /. float(period); + o = stereo_gain(g, stereo_of_mono(t)); +tel + +(* On entend nettement le signal en dent de scie, avec un pic à la fin de la + seconde. De façon intéressante, si on augmente le gain, le son apparaît comme + plus pincé, un peu comme les notes d'une guitare. *) + +node sawtooth_4() returns (o : stereo) +var t : float; g : float; +let + t = float(periodic(128)); + g = 3.0 *. float(periodic(period)) /. float(period); + o = stereo_gain(g, stereo_of_mono(t)); +tel + +(* En augmentant la période, les pics s'éloignent, en la diminuant, les + pics se rapprochent. *) + +node period_per_sec(a : int) returns (o : float) +let + o = float(periodic(period / a)) /. float(period / a); +tel + +node sawtooth_5() returns (o : stereo) +var t : float; g : float; +let + t = float(periodic(128)); + g = period_per_sec(2); + o = stereo_gain(g, stereo_of_mono(t)); +tel + +(* On peut aussi appliquer des gains différents sur le canal mono et stéréo. *) + +node every_sec(s : int) returns (c : bool) +let + c = periodic(period * s) = ((- 1) fby 0); +tel + +node sawtooth_6() returns (o : stereo) +var t : float; g1, g2 : float; +let + t = float(periodic(128)); + o = { l = g1 *. t; r = g2 *. t }; + automaton + state FastLeftSlowRight + do g1 = period_per_sec(1); + g2 = period_per_sec(8); + until every_sec(5) then SlowLeftFastRight + + state SlowLeftFastRight + do g1 = period_per_sec(5); + g2 = period_per_sec(1); + until every_sec(5) then FastLeftSlowRight + end +tel + +(* Tous ces sons ne sont pas très harmonieux. Peut-on en obtenir de plus purs ? + + Le traitement du signal nous enseigne, via la théorie de la transformée de + Fourier, que tout signal raisonnablement régulier peut se décomposer en une + somme (infinie) de sinusoïde. Autrement dit, les signaux sinusoïdaux peuvent + servir de briques de base élémentaires mais universelles. Considérés comme + des signaux audio, ils forment des tons purs, élémentaires. + *) + +node pure_tone(p : float) returns (o : float) +var t : float; +let + t = fcnt(0.0, 1.0); + o = sin(t *. (p /. float(period)) *. 2.0 *. Mathext.pi); +tel + +(* Par exemple, la sinusoïde de fréquence 440.1 Hz, communément désignée sous le + nom de La 440, devrait vous être familière. *) + +node main_pure_1() returns (o : stereo) +let + o = stereo_of_mono(pure_tone(440.0)); +tel + +(* En plus d'être la tonalité du téléphone, elle sert de référence pour + l'accordage des pianos, violons et d'autres instruments. + + https://fr.wikipedia.org/wiki/La_440 *) + +(* En mélangeant plusieurs sinusoïdes ensembles, on peut obtenir des effets + rétro assez amusants. *) + +node some_pure_tone(p : float; i : int) returns (s : stereo) +let + s = stereo_gain(period_per_sec(i + 1), stereo_of_mono(pure_tone(p))); +tel + +node oscillating_counter<>(i : int) returns (last o : int = 0) +var step : int; +let + step = if every_sec(1) then 1 else 0; + automaton + state Init + do o = i + until true then Increase + + state Increase + do o = last o + step + until o >= m then Decrease + + state Decrease + do o = last o - step + until o <= 0 then Increase + end +tel + +node main_pure_2() returns (o : stereo) +var periods : float^3; speeds : int^3; +let + periods = [440.0, 261.6256, 4186.009]; + speeds = map<<3>>(oscillating_counter<<10>>)([1, 3, 7]); + o = stereo_mix<<3>>(map<<3>> some_pure_tone(periods, speeds)); +tel + +(* Enfin, les amatrices et amateurs de piano pourront trouver sur la page + + https://en.wikipedia.org/wiki/Piano_key_frequencies + + une formule associant une fréquence de sinusoïde à une note de piano. On peut + l'utiliser comme suit. *) + +fun piano_freq_of_key(k : int) returns (f : float) +let + f = Mathext.pow(2.0, (Mathext.float(k) -. 49.0) /. 12.0) *. 440.0; +tel + +node tone_of_piano_key(k : int) returns (o : stereo) +let + o = stereo_of_mono(pure_tone(piano_freq_of_key(k))); +tel + +node maintain(c : bool; x : int on c; ini : int) returns (o : int) +let + o = merge c x ((ini fby o) whenot c); +tel + +node main_pure_3() returns (o : stereo) +var k : int; c : bool; +let + o = tone_of_piano_key(k); + k = maintain(c, 40 + periodic(53 - 40), 40); + c = periodic(period) = 0; +tel + +(* On peut essayer de programmer un piano midi. *) + +(* Pour générer des transitions propres entre les notes, on a besoin de modifier + nos tons à travers une "enveloppe". La plus classique est l'enveloppe dite + "Attack-Decay-Sustain-Release", cf. Wikipédia. + + https://en.wikipedia.org/wiki/Envelope_(music)#ADSR + + Le noeud ci-dessous produit une telle enveloppe périodiquement, tous les t + instants. L'enveloppe prend la forme d'un gain entre 0 et 1. + + Les paramètres a, d et s doivent-être tels que 0.0 < a + d + s < 1.0. Ils + expriment la fraction de t correspondant à chacune des quatre phases, la + phase d étant la fraction de t définie comme 1 - a - d - s. + + Le paramètre s_level est le niveau de la phase S, entre 0 et 1 donc. + + *) + +node adsr_envelope(t : int; a, d, s : float; s_level : float) + returns (e : float) +var c, a_stop, d_stop, s_stop : int; +let + a_stop = int(float(t) *. a); + d_stop = a_stop + int(float(t) *. d); + s_stop = d_stop + int(float(t) *. s); + c = periodic(t); + automaton + state Attack + do e = float(c) /. float(a_stop); + unless c >= a_stop continue Decay + + state Decay + var f : float; + do e = 1.0 -. (1.0 -. s_level) *. f; + f = float(c - a_stop) /. float(d_stop - a_stop); + unless c >= d_stop continue Sustain + + state Sustain + do e = s_level; + unless c >= s_stop continue Release + + state Release + do e = s_level *. (1.0 -. float(c - s_stop) /. float(t - s_stop)); + until c + 1 >= t continue Attack + end +tel + +node midi_piano<>(keys : int^2^n; time : int^n) returns (o : stereo) +var i, j : int; next : bool; duree_mesure : int; e : float; +let + duree_mesure = 2 * period; (* 1 mesure = 8 noires = 4 sec à 120 BPM. *) + + i = periodic(n); + j = maintain(next, i, 0); + o = stereo_gain(e, stereo_mix<<2>>(map<<2>> tone_of_piano_key(keys[>j<]))); + e = adsr_envelope(duree_mesure / time[> j <], 0.3, 0.1, 0.4, 0.5); + + automaton + state Next + do next = true + until true then Wait + + state Wait + var c : int; + do next = false; + c = 0 fby (c + 1); + until c >= (duree_mesure / time[> j <]) then Next + end +tel + +const num_keys : int = 82 + +node main_pure_4() returns (o : stereo) +var keys : int^2^num_keys; time : int^num_keys; +let + keys = [ + [44, 00], [37, 00], [40, 00], [42, 00], + [44, 00], [37, 00], [40, 00], [42, 00], + [44, 00], [37, 00], [40, 00], [42, 00], + [44, 00], [37, 00], [40, 00], [42, 00], + [44, 00], [37, 00], [41, 00], [42, 00], + [44, 00], [37, 00], [41, 00], [42, 00], + [44, 00], [37, 00], [41, 00], [42, 00], + [44, 00], [37, 00], [41, 00], [42, 00], + + [44, 00], + [37, 00], + [40, 37], [42, 00], [44, 00], + [37, 00], [40, 00], [42, 00], + + [35, 39], [32, 00], [35, 00], [37, 00], + [39, 00], [32, 00], [35, 00], [37, 00], + [35, 39], [32, 00], [35, 00], [37, 00], + [39, 00], [32, 00], [35, 00], + + [42, 00], + [35, 00], + [35, 40], [39, 00], [42, 00], + [35, 00], [40, 00], [39, 00], + [33, 37], [30, 00], [33, 00], [35, 00], + [33, 00], [30, 00], [33, 00], [35, 00], + [33, 37], [30, 00], [33, 00], [35, 00], + [37, 00], [30, 00], [33, 00], + + [0, 0], [0, 0], [0, 0], [0, 0] (* silence *) + ]; + time = [ + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + + 2, + 2, + 8, 8, 2, + 2, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 4, + + 2, + 2, + 8, 8, 2, + 2, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 8, 8, + 4, 4, 4, + + 1, 1, 1, 1 (* silence *) + ]; + o = midi_piano<>(keys, time); +tel + +(* Bonus : la méthode de Karplus-Strong pour la synthèse de son de guitare. + + https://en.wikipedia.org/wiki/Karplus%E2%80%93Strong_string_synthesis + + http://sites.music.columbia.edu/cmc/MusicAndComputers/chapter4/04_09.php +*) + +node flip(i: int) returns (o: float) +let + o = if (i % 2 = 0) then 1.0 else -.1.0 +tel + +node karplus_strong<>() returns (y : float) +var b : float^l; i: int; +let + i = 0 fby ((i+1) % l); + y = 0.5 *. (b[>i<] +. 0.0 fby y); + b = (mapi<> flip ()) fby ([b with [i] = y]); +tel + +node repeat<>(x : stereo) returns (o : stereo) +var last t : stereo^n = silence^n; +let + automaton + state Fill + do o = x; + t = [ last t with [ periodic(n) ] = x ] + until periodic(n) = n - 1 then Repeat + + state Repeat + do o = t[> periodic(n) <] + end +tel + +node saturating_counter(max : int) returns (o : int) +var c : int; +let + c = 0 fby (c + 1); + o = if c < max then c else max; +tel + +node main_kp() returns (o : stereo) +var s : stereo; +let + s = repeat<>({ l = karplus_strong<<115>>(); + r = karplus_strong<<55>>() }); + o = stereo_gain(float(saturating_counter(5 * period)) /. float(5 * period), + s); +tel + +(* Le noeud principal du programme. *) + +(* Vous pouvez choisir un des noeuds principaux main_XXX écrits ci-dessus, ou + bien écrire le votre. *) + +node main() returns (o : stereo) +let + o = main_pure_4(); +tel diff --git a/cours/audio/buffer.c b/cours/audio/buffer.c new file mode 100644 index 0000000..20cfffa --- /dev/null +++ b/cours/audio/buffer.c @@ -0,0 +1,61 @@ +#include "buffer.h" + +#include +#include +#include +#include +#include + +#define max(a, b) ((a) <= (b) ? (a) : (b)) + +void *malloc_checked(size_t size) { + void *result = malloc(size); + if (!result) { + perror("malloc()"); + exit(EXIT_FAILURE); + } + return result; +} + +char *strdup_checked(const char *s) { + char *result = strdup(s); + if (!result) { + perror("strdup()"); + exit(EXIT_FAILURE); + } + return result; +} + +buffer_t *buffer_alloc(size_t initial_size) { + buffer_t *buff = malloc_checked(sizeof *buff); + buff->data = malloc_checked(initial_size * sizeof *buff->data); + buff->size = initial_size; + buff->occupancy = 0; + return buff; +} + +void buffer_free(buffer_t *buffer) { + assert (buffer); + + free(buffer->data); + free(buffer); +} + +void buffer_resize(buffer_t *buff, size_t new_size) { + assert (buff); + assert (new_size >= buff->size); + + unsigned char *new_data = malloc_checked(new_size); + memcpy(new_data, buff->data, buff->occupancy); + free(buff->data); + buff->data = new_data; + buff->size = new_size; +} + +void buffer_write(buffer_t *buff, void *data, size_t data_size) { + assert (buff); + if (buff->occupancy + data_size > buff->size) + buffer_resize(buff, max(buff->size + data_size, 2 * buff->size)); + memcpy(buff->data + buff->occupancy, data, data_size); + buff->occupancy += data_size; +} diff --git a/cours/audio/buffer.h b/cours/audio/buffer.h new file mode 100644 index 0000000..f8a2743 --- /dev/null +++ b/cours/audio/buffer.h @@ -0,0 +1,27 @@ +#ifndef BUFFER_H +#define BUFFER_H + +#include + +/* A simple type of append-only buffers. */ + +typedef struct buffer { + unsigned char *data; + size_t size; + size_t occupancy; +} buffer_t; + +void *malloc_checked(size_t size); +char *strdup_checked(const char *); + +buffer_t *buffer_alloc(size_t initial_size); +void buffer_free(buffer_t *buff); + +void buffer_write(buffer_t *buff, void *data, size_t data_size); + +#define buffer_foreach(ty, var, buffer) \ + for (ty *var = (ty *)buffer->data; \ + var < (ty *)(buffer->data + buffer->occupancy); \ + var++) + +#endif /* BUFFER_H */ diff --git a/cours/audio/hept_ffi.h b/cours/audio/hept_ffi.h new file mode 100644 index 0000000..2d9b81b --- /dev/null +++ b/cours/audio/hept_ffi.h @@ -0,0 +1,52 @@ +#ifndef HEPT_FFI_H +#define HEPT_FFI_H + +#define UNPAREN(...) __VA_ARGS__ + +#define DECLARE_HEPT_FUN(module, name, inputs, outputs) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *) + +#define DECLARE_HEPT_FUN_NULLARY(module, name, outputs) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *) + +#define DEFINE_HEPT_FUN(module, name, inputs) \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *out) + +#define DEFINE_HEPT_FUN_NULLARY(module, name, inputs) \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *out) + +#define DECLARE_HEPT_NODE(module, name, inputs, outputs, state) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + typedef struct { state; } module ## __ ## name ## _mem; \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *, \ + module ## __ ## name ## _mem *); \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *) + +#define DECLARE_HEPT_NODE_NULLARY(module, name, outputs, state) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + typedef struct { state; } module ## __ ## name ## _mem; \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *, \ + module ## __ ## name ## _mem *); \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *) + +#define DEFINE_HEPT_NODE_RESET(module, name) \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *mem) + +#define DEFINE_HEPT_NODE_STEP(module, name, inputs) \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *out, \ + module ## __ ## name ## _mem *mem) + +#define DEFINE_HEPT_NODE_NULLARY_STEP(module, name, inputs) \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *out, \ + module ## __ ## name ## _mem *mem) + +/* FIXME remove when Heptagon's pervasives.h has been fixed. */ +typedef char * string; + +#endif /* HEPT_FFI */ diff --git a/cours/audio/main.c b/cours/audio/main.c new file mode 100644 index 0000000..d4e529a --- /dev/null +++ b/cours/audio/main.c @@ -0,0 +1,153 @@ +#include +#include +#include + +#include +#include + +#include "audio.h" +#include "vcd.h" + +const size_t sample_rate = Audio__period; + +void die(const char *message) { + fprintf(stderr, message); + exit(EXIT_FAILURE); +} + +SNDFILE *file_out = NULL; + +int main(int argc, char** argv) +{ + Audio__main_mem mem; + Audio__main_out res; + SDL_AudioSpec spec; + SDL_AudioDeviceID dev; + int opt = -1; + bool quiet = false; + size_t max_sec = SIZE_MAX; /* largest value of type size_t */ + const char *filename = NULL; + Uint32 buffered; + + while ((opt = getopt(argc, argv, "ho:qm:t:")) != -1) { + switch (opt) { + case 'h': + printf("Usage: %s OPTIONS\n", argv[0]); + printf("Options:\n"); + printf(" -o write samples to \n"); + printf(" -q do not play sound\n"); + printf(" -m play for seconds\n"); + printf(" -t dump traces in \n"); + printf(" -h display this message\n"); + return 0; + case 'q': + quiet = true; + break; + case 'o': + filename = optarg; + break; + case 'm': + max_sec = atoi(optarg); + break; + case 't': + hept_vcd_init(optarg, VCD_TIME_UNIT_US, 20); + break; + default: + fprintf(stderr, "Unknown option '%c'\n", opt); + exit(EXIT_FAILURE); + } + } + + if (SDL_Init(SDL_INIT_AUDIO) < 0) + die("Could not initialize SDL2\n"); + + /* Specification of requested output device. */ + bzero(&spec, sizeof spec); + spec.freq = sample_rate; /* Samples per second */ + spec.format = AUDIO_F32; /* Sample format: IEEE-754 32 bits */ + spec.channels = 2; /* Two channels */ + spec.samples = 4096; /* Buffers sized 4 KiB */ + spec.callback = NULL; + + if (!(dev = SDL_OpenAudioDevice(NULL, 0, &spec, NULL, 0))) + die("Could not open audio device\n"); + + if (filename != NULL) { + /* Specification of requested output file, if any. */ + SF_INFO info_out; + bzero(&info_out, sizeof info_out); + info_out.channels = 2; /* Two channels */ + info_out.samplerate = sample_rate; /* Samples per second */ + info_out.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16; /* File format */ + + if (!(file_out = sf_open(filename, SFM_WRITE, &info_out))) { + fprintf(stderr, "Could not open WAV file %s for writing\n", argv[1]); + SDL_Quit(); + exit(EXIT_FAILURE); + } + } + + Audio__main_reset(&mem); + float *buffer = calloc(spec.samples, sizeof *buffer); + SDL_PauseAudioDevice(dev, 0); + + /* Loop until we've produced the requested amount of samples, that is the + duration in seconds multiplied by the number of samples per second. This + number of samples shall be sent on each of both stereo channels. + + Each iteration sends spec.samples stereo samples to the audio device, + hence we halve it to get the number of generated samples per-channel. */ + for (size_t samples = 0; + samples < max_sec * sample_rate; + samples += spec.samples / 2) { + + /* Print sound progress. */ + printf("\rSent %08zu samples", samples); + if (max_sec != SIZE_MAX) { + printf(" (%2.0f%)", 100. * (double)samples / (max_sec * sample_rate)); + } + fflush(stdout); + + /* Exit immediately if requested, e.g., the user pressed Ctrl-C. */ + if (SDL_QuitRequested()) { + printf("\n"); + return 1; + } + + /* Step the node as much as necessary to fill a buffer. Each step produces + one stereo sample. */ + for (size_t i = 0; i < spec.samples; i += 2) { + Audio__main_step(&res, &mem); + buffer[i+0] = res.o.l; + buffer[i+1] = res.o.r; + } + + /* Send the generated sound to the sound card and/or file. */ + if (!quiet) + SDL_QueueAudio(dev, buffer, spec.samples * sizeof *buffer); + if (file_out) + sf_writef_float(file_out, buffer, spec.samples / 2); + + /* Throttle queued audio, otherwise we will certainly end up consuming all + available memory. */ + buffered = SDL_GetQueuedAudioSize(dev); + while (!quiet && buffered >= 1 << 22) { + SDL_Delay(50); + buffered = SDL_GetQueuedAudioSize(dev); + } + } + printf("\n"); + + /* Wait until the audio buffer is empty. */ + printf("Waiting for queue flush... "); fflush(stdout); + while ((buffered = SDL_GetQueuedAudioSize(dev)) != 0) + SDL_Delay(50); + printf("done.\n"); + + free(buffer); + if (file_out) + sf_close(file_out); + SDL_Quit(); + + return 0; +} diff --git a/cours/audio/mathext.c b/cours/audio/mathext.c new file mode 100644 index 0000000..17b4b42 --- /dev/null +++ b/cours/audio/mathext.c @@ -0,0 +1,55 @@ +#include +#include +#include + +/* Avoid Heptagon's math.h. I don't think there's a single place where to find + math.h on Apple-platforms. */ +#ifndef __APPLE__ +#include +#endif + +#include "mathext.h" + +void Mathext__float_step(int x, Mathext__float_out *o) { + o->o = (float)x; +} + +void Mathext__int_step(float x, Mathext__int_out *o) { + o->o = (int)x; +} + +void Mathext__floor_step(float x, Mathext__floor_out *o) { + o->o = floorf(x); +} + +void Mathext__sin_step(float x, Mathext__sin_out *o) { + o->o = sinf(x); +} + +void Mathext__cos_step(float x, Mathext__cos_out *o) { + o->o = cosf(x); +} + +void Mathext__atan2_step(float y, float x, Mathext__atan2_out *o) { + o->o = atan2f(y, x); +} + +void Mathext__pow_step(float x, float y, Mathext__pow_out *o) { + o->o = powf(x, y); +} + +void Mathext__hypot_step(float x, float y, Mathext__hypot_out *o) { + o->o = hypotf(x, y); +} + +void Mathext__sqrt_step(float x2, Mathext__sqrt_out *o) { + o->o = sqrtf(x2); +} + +void Mathext__modulo_step(int x, int y, Mathext__modulo_out *o) { + o->o = x % y; +} + +void Mathext__piano_freq_of_key_step(int n, Mathext__piano_freq_of_key_out *o) { + o->f = (float)(pow(2, (float)(n - 49) / (float)12) * 440.); +} diff --git a/cours/audio/mathext.epi b/cours/audio/mathext.epi new file mode 100644 index 0000000..daf1fd5 --- /dev/null +++ b/cours/audio/mathext.epi @@ -0,0 +1,16 @@ +external fun float(x : int) returns (o : float) +external fun int(x : float) returns (o : int) +external fun floor(x : float) returns (o : float) + +external fun sin(x : float) returns (o : float) +external fun cos(x : float) returns (o : float) +external fun atan2(y : float; x : float) returns (o : float) +external fun hypot(x : float; y : float) returns (o : float) +external fun sqrt(x2 : float) returns (o : float) +external fun pow(x : float; y : float) returns (o : float) + +external fun modulo(x : int; y : int) returns (o : int) + +external fun piano_freq_of_key(k : int) returns (f : float) + +const pi : float = 3.14115 diff --git a/cours/audio/mathext.h b/cours/audio/mathext.h new file mode 100644 index 0000000..cf2f285 --- /dev/null +++ b/cours/audio/mathext.h @@ -0,0 +1,27 @@ +#ifndef MATHEXT_H +#define MATHEXT_H + +#include "stdbool.h" +#include "assert.h" +#include "pervasives.h" + +#include "hept_ffi.h" + +DECLARE_HEPT_FUN(Mathext, float, (int), float o); +DECLARE_HEPT_FUN(Mathext, int, (float), int o); +DECLARE_HEPT_FUN(Mathext, floor, (float), float o); + +DECLARE_HEPT_FUN(Mathext, sin, (float), float o); +DECLARE_HEPT_FUN(Mathext, cos, (float), float o); +DECLARE_HEPT_FUN(Mathext, atan2, (float, float), float o); +DECLARE_HEPT_FUN(Mathext, hypot, (float, float), float o); +DECLARE_HEPT_FUN(Mathext, sqrt, (float), float o); +DECLARE_HEPT_FUN(Mathext, pow, (float, float), float o); + +DECLARE_HEPT_FUN(Mathext, modulo, (int, int), int o); + +DECLARE_HEPT_FUN(Mathext, piano_freq_of_key, (int), float f); + +static const float Mathext__pi = 3.14115; + +#endif /* MATHEXT_H */ diff --git a/cours/audio/mathext_types.h b/cours/audio/mathext_types.h new file mode 100644 index 0000000..00a5ee7 --- /dev/null +++ b/cours/audio/mathext_types.h @@ -0,0 +1,4 @@ +#ifndef MATHEXT_TYPES_H +#define MATHEXT_TYPES_H + +#endif /* MATHEXT_TYPES_H */ diff --git a/cours/audio/vcd.c b/cours/audio/vcd.c new file mode 100644 index 0000000..0bcbce2 --- /dev/null +++ b/cours/audio/vcd.c @@ -0,0 +1,75 @@ +#include "vcd.h" + +#include +#include +#include + +#include "vcd_lib.h" + +vcd_file_t *vcd = NULL; + +bool enabled = false; + +void hept_vcd_cleanup() { + if (vcd) { + printf("[vcd] saving trace\n"); + vcd_file_write(vcd); + vcd_file_free(vcd); + } +} + +void hept_vcd_init(const char *filename, + size_t time_number, + vcd_time_unit_t unit) { + if (vcd) { + fprintf(stderr, "[vcd] initialization has already been performed\n"); + } else { + vcd = vcd_file_alloc(filename, time_number, unit); + assert (vcd); + if (atexit(hept_vcd_cleanup)) { + perror("[vcd] atexit() failed"); + exit(EXIT_FAILURE); + } + printf("[vcd] will save trace to %s\n", filename); + } +} + +static inline void trace_samples(vcd_signal_t **signal, + const char *name, vcd_signal_type_t type, + void *samples, size_t count) { + if (!vcd) + return; + + if (!*signal) { + *signal = vcd_signal_alloc(name, type, 1 << 17); + if (!vcd_file_add_signal(vcd, *signal)) { + perror("vcd_file_add_signal()\n"); + exit(EXIT_FAILURE); + } + } + vcd_add_samples(*signal, samples, count); +} + +DEFINE_HEPT_NODE_RESET(Vcd, trace_bool) { + mem->signal = NULL; +} + +DEFINE_HEPT_NODE_STEP(Vcd, trace_bool, (string name, int v)) { + trace_samples(&mem->signal, name, VCD_SIGNAL_TYPE_BOOL, &v, 1); +} + +DEFINE_HEPT_NODE_RESET(Vcd, trace_int) { + mem->signal = NULL; +} + +DEFINE_HEPT_NODE_STEP(Vcd, trace_int, (string name, int v)) { + trace_samples(&mem->signal, name, VCD_SIGNAL_TYPE_INT, &v, 1); +} + +DEFINE_HEPT_NODE_RESET(Vcd, trace_float) { + mem->signal = NULL; +} + +DEFINE_HEPT_NODE_STEP(Vcd, trace_float, (string name, float v)) { + trace_samples(&mem->signal, name, VCD_SIGNAL_TYPE_FLOAT, &v, 1); +} diff --git a/cours/audio/vcd.epi b/cours/audio/vcd.epi new file mode 100644 index 0000000..20eb209 --- /dev/null +++ b/cours/audio/vcd.epi @@ -0,0 +1,6 @@ +(* A basic library for producing VCD files from Heptagon code, typically used + for debugging purposes. *) + +external node trace_bool(name : string; v : bool) returns () +external node trace_int(name : string; v : int) returns () +external node trace_float(name : string; v : float) returns () diff --git a/cours/audio/vcd.h b/cours/audio/vcd.h new file mode 100644 index 0000000..f8e8920 --- /dev/null +++ b/cours/audio/vcd.h @@ -0,0 +1,18 @@ +#ifndef VCD +#define VCD + +#include +#include + +#include "hept_ffi.h" +#include "vcd_lib.h" + +void hept_vcd_init(const char *filename, + size_t time_number, + vcd_time_unit_t unit); + +DECLARE_HEPT_NODE(Vcd, trace_bool, (string, int),, vcd_signal_t *signal); +DECLARE_HEPT_NODE(Vcd, trace_int, (string, int),, vcd_signal_t *signal); +DECLARE_HEPT_NODE(Vcd, trace_float, (string, float),, vcd_signal_t *signal); + +#endif /* VCD */ diff --git a/cours/audio/vcd_lib.c b/cours/audio/vcd_lib.c new file mode 100644 index 0000000..d9b329c --- /dev/null +++ b/cours/audio/vcd_lib.c @@ -0,0 +1,198 @@ +#include "vcd_lib.h" + +#include +#include +#include +#include +#include +#include + +#include "buffer.h" + +size_t vcd_sizeof_signal_type(vcd_signal_type_t type) { + switch (type) { + case VCD_SIGNAL_TYPE_BOOL: + return sizeof(int); + case VCD_SIGNAL_TYPE_INT: + return sizeof(int); + case VCD_SIGNAL_TYPE_FLOAT: + return sizeof(float); + } +} + +typedef struct vcd_signal { + char *name; + vcd_signal_type_t type; + buffer_t *samples; +} vcd_signal_t; + +vcd_signal_t *vcd_signal_alloc(const char *name, + vcd_signal_type_t type, + size_t initial_buffer_size) { + vcd_signal_t *res = malloc_checked(sizeof *res); + res->name = strdup_checked(name); + res->type = type; + res->samples = + buffer_alloc(initial_buffer_size * vcd_sizeof_signal_type(type)); + return res; +} + +void vcd_signal_free(vcd_signal_t *signal) { + assert (signal); + + buffer_free(signal->samples); + free(signal->name); + free(signal); +} + +void vcd_add_samples(vcd_signal_t *signal, void *samples, size_t count) { + assert (signal); + assert (samples); + buffer_write(signal->samples, + samples, + count * vcd_sizeof_signal_type(signal->type)); +} + +const char *vcd_time_unit_repr(vcd_time_unit_t u) { + const char *table[] = { "s", "ms", "us", "ns", "ps", "fs" }; + assert (VCD_TIME_UNIT_S <= u && u <= VCD_TIME_UNIT_FS); + return table[u]; +} + +typedef struct vcd_file { + char *filename; + buffer_t *signals; + vcd_time_unit_t time_unit; + size_t time_unit_factor; +} vcd_file_t; + +vcd_file_t *vcd_file_alloc(const char *filename, + vcd_time_unit_t time_unit, + size_t time_unit_factor) { + assert (filename); + + vcd_file_t *vcd = malloc_checked(sizeof *vcd); + vcd->filename = strdup_checked(filename); + vcd->time_unit = time_unit; + vcd->time_unit_factor = time_unit_factor; + vcd->signals = buffer_alloc(10 * sizeof(vcd_signal_t)); + + return vcd; +} + +void vcd_file_free(vcd_file_t *vcd) { + assert (vcd); + + /* Free buffers. */ + buffer_foreach (vcd_signal_t *, psig, vcd->signals) + vcd_signal_free(*psig); + buffer_free(vcd->signals); + + free(vcd->filename); + free(vcd); +} + +vcd_signal_t *vcd_file_lookup_signal(const vcd_file_t *vcd, const char *name) { + assert (vcd); + assert (name); + + vcd_signal_t *res = NULL; + + buffer_foreach (vcd_signal_t *, psig, vcd->signals) { + if (strcmp((*psig)->name, name) == 0) { + res = *psig; + break; + } + } + + return res; +} + +bool vcd_file_add_signal(const vcd_file_t *vcd, vcd_signal_t *signal) { + assert (vcd); + assert (signal); + + if (vcd_file_lookup_signal(vcd, signal->name)) + return false; + + buffer_write(vcd->signals, &signal, sizeof signal); + return true; +} + +bool vcd_file_write(vcd_file_t *vcd) { + FILE *f = fopen(vcd->filename, "w"); + + if (!f) + return false; + + time_t current_time; + time(¤t_time); + + fprintf(f, "$version Generated by vcd.c $end\n"); + fprintf(f, "$date %s $end\n", ctime(¤t_time)); + fprintf(f, "$timescale %zu %s $end\n", + vcd->time_unit_factor, + vcd_time_unit_repr(vcd->time_unit)); + + /* Dump signal declarations. */ + fprintf(f, "$scope module Top $end\n"); + buffer_foreach (vcd_signal_t *, psig, vcd->signals) { + fprintf(f, "$var "); + switch ((*psig)->type) { + case VCD_SIGNAL_TYPE_BOOL: + fprintf(f, "wire 1"); + break; + case VCD_SIGNAL_TYPE_INT: + fprintf(f, "integer %zu", 8 * sizeof(int)); + break; + case VCD_SIGNAL_TYPE_FLOAT: + fprintf(f, "real 32"); + break; + } + fprintf(f, " %p %s $end\n", (*psig), (*psig)->name); + } + fprintf(f, "$upscope $end\n"); + fprintf(f, "$enddefinitions\n"); + + /* Dump samples. */ + fprintf(f, "$dumpvars\n"); + + /* We maintain a pointer to the current sample in each buffer. */ + size_t signal_count = vcd->signals->occupancy / sizeof(vcd_signal_t *); + unsigned char **psamples = calloc(signal_count, sizeof(unsigned char *)); + assert (psamples); + for (size_t i = 0; i < signal_count; i++) + psamples[i] = ((vcd_signal_t **)vcd->signals->data)[i]->samples->data; + + /* We dump */ + bool active = true; + for (size_t step = 0; active; step++) { + fprintf(f, "#%zu\n", step); + active = false; + for (size_t i = 0; i < signal_count; i++) { + vcd_signal_t *sig = ((vcd_signal_t **)vcd->signals->data)[i]; + if (psamples[i] < sig->samples->data + sig->samples->occupancy) { + active = true; + switch (sig->type) { + case VCD_SIGNAL_TYPE_BOOL: + fprintf(f, "%d%p\n", (*(int *)psamples[i] ? 1 : 0), sig); + psamples[i] += sizeof(int); + break; + case VCD_SIGNAL_TYPE_INT: + fprintf(f, "r%d %p\n", *(int *)psamples, sig); + psamples[i] += sizeof(int); + break; + case VCD_SIGNAL_TYPE_FLOAT: + fprintf(f, "r%.16g %p\n", *(float *)psamples, sig); + psamples[i] += sizeof(float); + break; + } + } + } + } + + free(psamples); + + fclose(f); + return true; +} diff --git a/cours/audio/vcd_lib.h b/cours/audio/vcd_lib.h new file mode 100644 index 0000000..cc5189a --- /dev/null +++ b/cours/audio/vcd_lib.h @@ -0,0 +1,47 @@ +#ifndef VCD_LIB_H +#define VCD_LIB_H + +#include +#include + +typedef enum vcd_signal_type { + VCD_SIGNAL_TYPE_FLOAT, + VCD_SIGNAL_TYPE_INT, + VCD_SIGNAL_TYPE_BOOL +} vcd_signal_type_t; + +size_t vcd_sizeof_signal_type(vcd_signal_type_t); + +typedef struct vcd_signal vcd_signal_t; + +vcd_signal_t *vcd_signal_alloc(const char *name, + vcd_signal_type_t type, + size_t initial_buffer_size); +void vcd_signal_free(vcd_signal_t *signal); + +void vcd_add_samples(vcd_signal_t *signal, void *samples, size_t count); + +typedef enum vcd_time_unit { + VCD_TIME_UNIT_S, + VCD_TIME_UNIT_MS, + VCD_TIME_UNIT_US, + VCD_TIME_UNIT_NS, + VCD_TIME_UNIT_PS, + VCD_TIME_UNIT_FS, +} vcd_time_unit_t; + +const char *vcd_time_unit_repr(vcd_time_unit_t); + +typedef struct vcd_file vcd_file_t; + +vcd_file_t *vcd_file_alloc(const char *filename, + vcd_time_unit_t time_unit, + size_t time_unit_factor); +void vcd_file_free(vcd_file_t *); + +vcd_signal_t *vcd_file_lookup_signal(const vcd_file_t *vcd, const char *name); +bool vcd_file_add_signal(const vcd_file_t *vcd, vcd_signal_t *signal); + +bool vcd_file_write(vcd_file_t *); + +#endif /* VCD_LIB_H */ diff --git a/cours/audio/vcd_types.h b/cours/audio/vcd_types.h new file mode 100644 index 0000000..0512dd9 --- /dev/null +++ b/cours/audio/vcd_types.h @@ -0,0 +1,4 @@ +#ifndef VCD_TYPES +#define VCD_TYPES + +#endif /* VCD_TYPES */ diff --git a/cours/inverted-pendulum/Makefile b/cours/inverted-pendulum/Makefile new file mode 100644 index 0000000..60b5ef5 --- /dev/null +++ b/cours/inverted-pendulum/Makefile @@ -0,0 +1,62 @@ +HEPTC=heptc +PYGMENTS=python -m pygments -x -l ../../../notes/heptagon.py:HeptagonLexer + +HEPT_NAME=pendulum + +HEPT_SOURCES=$(HEPT_NAME).ept +HEPT_GENERATED=\ + $(HEPT_NAME)_types.h \ + $(HEPT_NAME).h \ + $(HEPT_NAME)_types.c \ + $(HEPT_NAME).c +C_SOURCES=mathext.c + +SWIG_SOURCE=$(HEPT_NAME).i +SWIG_GENERATED=$(HEPT_NAME).py $(HEPT_NAME)_wrap.c + +PY_SUFFIX=$(shell python -c \ + "import sysconfig as s; print(s.get_config_var('EXT_SUFFIX'))") +TARGET?=_$(HEPT_NAME)$(PY_SUFFIX) + +HTML=pendulum.html + +.PHONY: all clean run runplot repl + +all: $(TARGET) + +clean: + rm -f *.epci *.log *.mls *.obc *.epci *.epo + rm -f $(SWIG_GENERATED) $(HEPT_GENERATED) + rm -rf $(HEPT_NAME)_c build __pycache__ *.so + rm -f $(HTML) + +repl: all + python -i -c 'import $(HEPT_NAME)' + +run: all + ./main.py + +runplot: all + ./main.py --plot + +%.html: %.ept + $(PYGMENTS) -O full -o $@ $^ + +$(TARGET): $(SWIG_GENERATED) $(HEPT_GENERATED) $(C_SOURCES) + python setup.py build_ext --inplace + @echo $@ + +$(SWIG_GENERATED): $(SWIG_SOURCE) $(HEPT_GENERATED) $(C_SOURCES) + swig -python $< + +%_types.h %_types.c %.c %.h %.epci : %.ept + $(HEPTC) -c -target c $< + cp $(foreach ext,c h,$(basename $<)_c/*.$(ext)) . + +%.epci: %.epi + $(HEPTC) $< + +$(HEPT_SOURCES): mathext.epci debug.epci + +%.html: %.ept + $(PYGMENTS) -O full -o $@ $^ diff --git a/cours/inverted-pendulum/README.md b/cours/inverted-pendulum/README.md new file mode 100644 index 0000000..8fd88a7 --- /dev/null +++ b/cours/inverted-pendulum/README.md @@ -0,0 +1,43 @@ +# Pendule inversé et contrôleur PID + +## Présentation + +Ce dossier contient un exemple classique de système physique qu'on peut chercher +à contrôler, le pendule inversé. + +## Dépendances + +Le projet utilise du code Python 3, et le générateur d'interfaces SWIG pour +exposer le code C généré par Heptagon à Python. Pour le lancer, vous devez : + +- avoir Heptagon, à installer via OPAM, + +- avoir SWIG, à installer via le gestionnaire de paquet de votre système + d'exploitation (`apt-get` sous GNU/Linux Debian ou Ubuntu, `pacman` sous Arch + Linux, `brew` sous macOS, etc.), + +- avoir la bibliothèque Python `matplotlib`, à installer `pip install --user + matplotlib`. + +## Utilisation + +Vous pouvez lancer le programme en exécutant `make run`. + +Pour lancer le programme avec l'affichage des courbes, exécutez `make runplot`. +Attention : afficher les courbes engendre un ralentissement assez conséquent. + +L'interface utilisateur fonctione de la façon suivante : + +- **Q** permet de quitter, + +- **P** permet de mettre en pause, + +- **R** permet de réinitialiser la simulation, + +- **M** permet de changer de mode (manuel, contrôleur PID, contrôleur BangBang), + +- **←** et **→** permettent de déplacer le mobile, + +- **clic gauche** permet de déplacer le mobile à la souris. + +Les deux dernières commandes ne fonctionnent qu'en mode manuel. diff --git a/cours/inverted-pendulum/cutils.c b/cours/inverted-pendulum/cutils.c new file mode 100644 index 0000000..d4ab5ca --- /dev/null +++ b/cours/inverted-pendulum/cutils.c @@ -0,0 +1,88 @@ +/* This file is part of SyncContest. + Copyright (C) 2017-2020 Eugene Asarin, Mihaela Sighireanu, Adrien Guatto. */ + +#include "cutils.h" + +#include +#include +#include +#include +#include + +log_verbosity_level level = LOG_INFO; +FILE *f = NULL; +char *filename = NULL; + +void log_set_verbosity_level(log_verbosity_level l) { + level = l; +} + +void log_message_v(log_verbosity_level msg_level, const char *fmt, va_list va) { + va_list vb; + FILE *out = msg_level == LOG_INFO ? stdout : stderr; + + if (msg_level > level) + return; + + va_copy(vb, va); + if (f != NULL) { + vfprintf(f, fmt, va); + } + + vfprintf(out, fmt, vb); + fflush(out); +} + +void log_message(log_verbosity_level msg_level, const char *fmt, ...) { + va_list va; + va_start(va, fmt); + log_message_v(msg_level, fmt, va); + va_end(va); +} + +void log_fatal(const char *fmt, ...) { + va_list va; + va_start(va, fmt); + log_message_v(LOG_FATAL, fmt, va); + va_end(va); + exit(EXIT_FAILURE); +} + +void log_info(const char *fmt, ...) { + va_list va; + va_start(va, fmt); + log_message_v(LOG_INFO, fmt, va); + va_end(va); +} + +void log_debug(const char *fmt, ...) { + va_list va; + va_start(va, fmt); + log_message_v(LOG_DEBUG, fmt, va); + va_end(va); +} + +void log_init(const char *fn) { + if (!fn) { + log_info("[log] initializing, not saving to file\n"); + return; + } + + filename = strdup(fn); + assert (filename); + f = fopen(filename, "w"); + if (!f) + log_fatal("[log] could not open log file %s (fopen)", filename); + + log_info("[log] logging to %s\n", filename); +} + +void log_shutdown() { + if (f) { + log_info("[log] shutting down, closing %s\n", filename); + fclose(f); + free(filename); + } else { + log_info("[log] shutting down\n"); + } +} diff --git a/cours/inverted-pendulum/cutils.h b/cours/inverted-pendulum/cutils.h new file mode 100644 index 0000000..2876046 --- /dev/null +++ b/cours/inverted-pendulum/cutils.h @@ -0,0 +1,30 @@ +/* This file is part of SyncContest. + Copyright (C) 2017-2020 Eugene Asarin, Mihaela Sighireanu, Adrien Guatto. */ + +#ifndef CUTILS_H +#define CUTILS_H + +#include + +typedef enum { + LOG_FATAL = 0, + LOG_INFO = 1, + LOG_DEBUG = 2, +} log_verbosity_level; + +/* Calling `log_init(fn)` initializes the logging subsystem, asking it to save + log messages to the file `fn`. This pointer may be NULL, in which case the + messages are not saved. */ +void log_init(const char *filename); +void log_shutdown(); + +void log_set_verbosity_level(log_verbosity_level level); + +void log_message_v(log_verbosity_level level, const char *fmt, va_list); +void log_message(log_verbosity_level level, const char *fmt, ...); + +void log_fatal(const char *fmt, ...); +void log_info(const char *fmt, ...); +void log_debug(const char *fmt, ...); + +#endif /* CUTILS_H */ diff --git a/cours/inverted-pendulum/debug.c b/cours/inverted-pendulum/debug.c new file mode 100644 index 0000000..6fc6f31 --- /dev/null +++ b/cours/inverted-pendulum/debug.c @@ -0,0 +1,39 @@ +#include "debug.h" + +#include "cutils.h" + +void Debug__dbg_step(char *msg, Debug__dbg_out *o) { + log_info("%s\n", msg); +} + +void Debug__dbg_bool_step(char *msg, bool x, Debug__dbg_bool_out *o) { + log_info("%s %d\n", msg, x); +} + +void Debug__dbg_int_step(char *msg, int x, Debug__dbg_int_out *o) { + log_info("%s %d\n", msg, x); +} + +void Debug__dbg_float_step(char *msg, float x, Debug__dbg_float_out *o) { + log_info("%s %f\n", msg, x); +} + +void Debug__d_init_step(Debug__d_init_out *o) { + /* Empty by design */ +} + +void Debug__d_string_step(Debug__world _w, char *s, Debug__d_string_out *o) { + log_info("%s", s); +} + +void Debug__d_bool_step(Debug__world _w, bool b, Debug__d_bool_out *o) { + log_info("%d", b); +} + +void Debug__d_int_step(Debug__world _w, int i, Debug__d_int_out *o) { + log_info("%d", i); +} + +void Debug__d_float_step(Debug__world _w, float f, Debug__d_float_out *o) { + log_info("%f", f); +} diff --git a/cours/inverted-pendulum/debug.epi b/cours/inverted-pendulum/debug.epi new file mode 100644 index 0000000..56e5014 --- /dev/null +++ b/cours/inverted-pendulum/debug.epi @@ -0,0 +1,11 @@ +external fun dbg(msg : string) returns () +external fun dbg_bool(msg : string; x : bool) returns () +external fun dbg_int(msg : string; x : int) returns () +external fun dbg_float(msg : string; x : float) returns () + +type world +external fun d_init() returns (n : world) +external fun d_string(w : world; string) returns (n : world) +external fun d_bool(w : world; bool) returns (n : world) +external fun d_int(w : world; int) returns (n : world) +external fun d_float(w : world; float) returns (n : world) diff --git a/cours/inverted-pendulum/debug.h b/cours/inverted-pendulum/debug.h new file mode 100644 index 0000000..9a85ac1 --- /dev/null +++ b/cours/inverted-pendulum/debug.h @@ -0,0 +1,23 @@ +#ifndef DEBUG_H +#define DEBUG_H + +#include "stdbool.h" +#include "assert.h" +#include "pervasives.h" + +#include "hept_ffi.h" + +DECLARE_HEPT_FUN(Debug, dbg, (char *),); +DECLARE_HEPT_FUN(Debug, dbg_bool, (char *, bool),); +DECLARE_HEPT_FUN(Debug, dbg_int, (char *, int),); +DECLARE_HEPT_FUN(Debug, dbg_float, (char *, float),); + +typedef struct { } Debug__world; + +DECLARE_HEPT_FUN_NULLARY(Debug, d_init, Debug__world n); +DECLARE_HEPT_FUN(Debug, d_string, (Debug__world, char *), Debug__world n); +DECLARE_HEPT_FUN(Debug, d_bool, (Debug__world, bool), Debug__world n); +DECLARE_HEPT_FUN(Debug, d_int, (Debug__world, int), Debug__world n); +DECLARE_HEPT_FUN(Debug, d_float, (Debug__world, float), Debug__world n); + +#endif /* DEBUG_H */ diff --git a/cours/inverted-pendulum/debug_types.h b/cours/inverted-pendulum/debug_types.h new file mode 100644 index 0000000..d12b30d --- /dev/null +++ b/cours/inverted-pendulum/debug_types.h @@ -0,0 +1,4 @@ +#ifndef DEBUG_TYPES_H +#define DEBUG_TYPES_H + +#endif /* DEBUG_TYPES_H */ diff --git a/cours/inverted-pendulum/hept_ffi.h b/cours/inverted-pendulum/hept_ffi.h new file mode 100644 index 0000000..2d9b81b --- /dev/null +++ b/cours/inverted-pendulum/hept_ffi.h @@ -0,0 +1,52 @@ +#ifndef HEPT_FFI_H +#define HEPT_FFI_H + +#define UNPAREN(...) __VA_ARGS__ + +#define DECLARE_HEPT_FUN(module, name, inputs, outputs) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *) + +#define DECLARE_HEPT_FUN_NULLARY(module, name, outputs) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *) + +#define DEFINE_HEPT_FUN(module, name, inputs) \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *out) + +#define DEFINE_HEPT_FUN_NULLARY(module, name, inputs) \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *out) + +#define DECLARE_HEPT_NODE(module, name, inputs, outputs, state) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + typedef struct { state; } module ## __ ## name ## _mem; \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *, \ + module ## __ ## name ## _mem *); \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *) + +#define DECLARE_HEPT_NODE_NULLARY(module, name, outputs, state) \ + typedef struct { outputs; } module ## __ ## name ## _out; \ + typedef struct { state; } module ## __ ## name ## _mem; \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *, \ + module ## __ ## name ## _mem *); \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *) + +#define DEFINE_HEPT_NODE_RESET(module, name) \ + void module ## __ ## name ##_reset(module ## __ ## name ## _mem *mem) + +#define DEFINE_HEPT_NODE_STEP(module, name, inputs) \ + void module ## __ ## name ##_step(UNPAREN inputs, \ + module ## __ ## name ## _out *out, \ + module ## __ ## name ## _mem *mem) + +#define DEFINE_HEPT_NODE_NULLARY_STEP(module, name, inputs) \ + void module ## __ ## name ##_step(module ## __ ## name ## _out *out, \ + module ## __ ## name ## _mem *mem) + +/* FIXME remove when Heptagon's pervasives.h has been fixed. */ +typedef char * string; + +#endif /* HEPT_FFI */ diff --git a/cours/inverted-pendulum/main.py b/cours/inverted-pendulum/main.py new file mode 100755 index 0000000..2a3e608 --- /dev/null +++ b/cours/inverted-pendulum/main.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python3 + +import pendulum +import math, time, argparse, sys + +from tkinter import * +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg +from matplotlib.figure import Figure +import matplotlib.animation as animation + +parser = argparse.ArgumentParser( + usage = "%(prog)s [OPTION]", + description = "Simulate an inverted pendulum" +) +parser.add_argument('-p', '--plot', dest='plot', action='store_true') +parser.set_defaults(plot=False) +args = parser.parse_args() + +c = pendulum.Cart() +input_x0 = pendulum.max_w / 2 +mouse_diff = False +mode_diff = False +pause = False +last_time = {0: time.time()} + +theta_samples = [0] +error_samples = [0] + +root = Tk() +if sys.platform == 'linux': + root.attributes('-type', 'dialog') +root.title("Inverted pendulum") + +w = Canvas(root, width=pendulum.max_w, height=pendulum.max_h) + +f = Frame(root) +f.pack(side=TOP, fill=X) + +vmode = StringVar() +lmode = Label(f, textvariable=vmode) +lmode.pack(side=LEFT) +vinfo = StringVar() +linfo = Label(f, textvariable=vinfo) +linfo.pack(side=RIGHT) + +def create_line(x1, y1, x2, y2, **kwargs): + w.create_line(x1, pendulum.max_h - y1, x2, pendulum.max_h - y2, **kwargs) + +def create_oval(x1, y1, x2, y2, **kwargs): + w.create_oval(x1, pendulum.max_h - y1, x2, pendulum.max_h - y2, **kwargs) + +def push_sample(data, sample): + data.append(sample) + +def repaint(): + global input_x0, out, pause, error_samples + joint_color = "#4760B2" + mass_color = "#EF1010" + radius = 4 + if not pause: + out = c.step(input_x0, mouse_diff, mode_diff) + push_sample(theta_samples, math.fmod(out.p.theta, 2 * math.pi)) + push_sample(error_samples, out.p.error) + w.delete("all") + create_line(0, out.p.y0, pendulum.max_w, out.p.y0, width=5); + create_line(out.p.x0, out.p.y0, out.p.x, out.p.y, width=3) + create_oval(out.p.x - radius, out.p.y - radius, + out.p.x + radius, out.p.y + radius, fill=mass_color) + create_oval(out.p.x0 - 2 * radius, out.p.y0 - 2 * radius, + out.p.x0 + 2 * radius, out.p.y0 + 2 * radius, fill=joint_color) + input_x0 = out.p.x0 + vmode.set(f"Current mode: {out.p.mode_name}" + + (" (paused)" if pause else "")) + new_time = time.time() + vinfo.set("x: {:07.2f}, theta: {:06.2f}, error: {:06.2f}, fps: {:06.2f}" \ + .format(out.p.x, + out.p.theta, + out.p.error, + 1. / (new_time - last_time[0]))) + last_time.update({0:new_time}) + + w.after(int(1000 * pendulum.dt), repaint) + +def mouse_change_callback(event): + global mouse_diff, input_x0 + input_x0 = event.x + mouse_diff = not mouse_diff + +def key_pressed(event): + global input_x0, mouse_diff, mode_diff, pause + x_step = 3 + if event.keysym == 'q': + quit() + elif event.keysym == 'r': + c.reset() + input_x0 = pendulum.max_w / 2 + theta_samples.clear() + error_samples.clear() + elif event.keysym == 'm': + mode_diff = not mode_diff + elif event.keysym == 'p': + pause = not pause + elif event.keysym == 'Right' and not pause: + input_x0 = input_x0 + x_step + mouse_diff = not mouse_diff + elif event.keysym == 'Left' and not pause: + input_x0 = input_x0 - x_step + mouse_diff = not mouse_diff + +w.pack(expand = YES, fill = BOTH) +root.bind("", mouse_change_callback) +root.bind("", mouse_change_callback) +root.bind("", key_pressed) + +fig = Figure(figsize=(6, 4.5), dpi=100) +ax = fig.add_subplot(111) +text = ax.text(0.97,0.97, "", transform=ax.transAxes, ha="left", va="top") +line1, = ax.plot([]) +line1.set_label("theta") +line2, = ax.plot([]) +line2.set_label("error") +ax.legend(loc = 'upper left') + +def animate(i): + line1.set_data(range(len(theta_samples)), theta_samples) + line2.set_data(range(len(error_samples)), error_samples) + ax.set_xlim([0, len(theta_samples)]) + ax.set_ylim([min(min(theta_samples), min(error_samples)), + max(max(theta_samples), max(error_samples)) + 0.0001]) + + return line1, text + +if args.plot: + canfig = FigureCanvasTkAgg(fig, master=root) + canfig.draw() + canfig.get_tk_widget().pack(side=BOTTOM, fill=X) + ani = animation.FuncAnimation(fig, + animate, + interval=int(1000 * pendulum.dt)) + +root.focus_set() +repaint() +mainloop() diff --git a/cours/inverted-pendulum/mathext.c b/cours/inverted-pendulum/mathext.c new file mode 100644 index 0000000..9cc1169 --- /dev/null +++ b/cours/inverted-pendulum/mathext.c @@ -0,0 +1,55 @@ +#include +#include +#include + +/* Avoid Heptagon's math.h. I don't think there's a single place where to find + math.h on Apple-platforms. */ +#ifndef __APPLE__ +#include +#endif + +#include "mathext.h" + +void Mathext__float_step(int x, Mathext__float_out *o) { + o->o = (float)x; +} + +void Mathext__int_step(float x, Mathext__int_out *o) { + o->o = (int)x; +} + +void Mathext__floor_step(float x, Mathext__floor_out *o) { + o->o = floorf(x); +} + +void Mathext__sin_step(float x, Mathext__sin_out *o) { + o->o = sinf(x); +} + +void Mathext__cos_step(float x, Mathext__cos_out *o) { + o->o = cosf(x); +} + +void Mathext__atan2_step(float y, float x, Mathext__atan2_out *o) { + o->o = atan2f(y, x); +} + +void Mathext__pow_step(float x, float y, Mathext__pow_out *o) { + o->o = powf(x, y); +} + +void Mathext__hypot_step(float x, float y, Mathext__hypot_out *o) { + o->o = hypotf(x, y); +} + +void Mathext__sqrt_step(float x2, Mathext__sqrt_out *o) { + o->o = sqrtf(x2); +} + +void Mathext__fmod_step(float x, float y, Mathext__fmod_out *o) { + o->o = fmodf(x, y); +} + +void Mathext__piano_freq_of_key_step(int n, Mathext__piano_freq_of_key_out *o) { + o->f = (float)(pow(2, (float)(n - 49) / (float)12) * 440.); +} diff --git a/cours/inverted-pendulum/mathext.epi b/cours/inverted-pendulum/mathext.epi new file mode 100644 index 0000000..5fe95cd --- /dev/null +++ b/cours/inverted-pendulum/mathext.epi @@ -0,0 +1,16 @@ +external fun float(x : int) returns (o : float) +external fun int(x : float) returns (o : int) +external fun floor(x : float) returns (o : float) + +external fun sin(x : float) returns (o : float) +external fun cos(x : float) returns (o : float) +external fun atan2(y : float; x : float) returns (o : float) +external fun hypot(x : float; y : float) returns (o : float) +external fun sqrt(x2 : float) returns (o : float) +external fun pow(x : float; y : float) returns (o : float) + +external fun fmod(x : float; y : float) returns (o : float) + +external fun piano_freq_of_key(k : int) returns (f : float) + +const pi : float = 3.14115 diff --git a/cours/inverted-pendulum/mathext.h b/cours/inverted-pendulum/mathext.h new file mode 100644 index 0000000..a2e8624 --- /dev/null +++ b/cours/inverted-pendulum/mathext.h @@ -0,0 +1,27 @@ +#ifndef MATHEXT_H +#define MATHEXT_H + +#include "stdbool.h" +#include "assert.h" +#include "pervasives.h" + +#include "hept_ffi.h" + +DECLARE_HEPT_FUN(Mathext, float, (int), float o); +DECLARE_HEPT_FUN(Mathext, int, (float), int o); +DECLARE_HEPT_FUN(Mathext, floor, (float), float o); + +DECLARE_HEPT_FUN(Mathext, sin, (float), float o); +DECLARE_HEPT_FUN(Mathext, cos, (float), float o); +DECLARE_HEPT_FUN(Mathext, atan2, (float, float), float o); +DECLARE_HEPT_FUN(Mathext, hypot, (float, float), float o); +DECLARE_HEPT_FUN(Mathext, sqrt, (float), float o); +DECLARE_HEPT_FUN(Mathext, pow, (float, float), float o); + +DECLARE_HEPT_FUN(Mathext, fmod, (float, float), float o); + +DECLARE_HEPT_FUN(Mathext, piano_freq_of_key, (int), float f); + +static const float Mathext__pi = 3.14115; + +#endif /* MATHEXT_H */ diff --git a/cours/inverted-pendulum/mathext_types.h b/cours/inverted-pendulum/mathext_types.h new file mode 100644 index 0000000..a046195 --- /dev/null +++ b/cours/inverted-pendulum/mathext_types.h @@ -0,0 +1,7 @@ +#ifndef MATHEXT_TYPES_H +#define MATHEXT_TYPES_H + +/* FIXME remove once Heptagon properly defines its string type. */ +typedef char *string; + +#endif /* MATHEXT_TYPES_H */ diff --git a/cours/inverted-pendulum/node_wrapping.py b/cours/inverted-pendulum/node_wrapping.py new file mode 100644 index 0000000..5e9181f --- /dev/null +++ b/cours/inverted-pendulum/node_wrapping.py @@ -0,0 +1,38 @@ + +from sys import modules +from inspect import getmembers, isfunction, isclass + +# This should probably be done using MetaClasses. + +def create_node_wrapper_class(name): + table = dict([ (n, o) for n, o in getmembers(modules[__name__]) \ + if isfunction(o) or isclass(o) ]) + out_class = table[name + "_out"] + mem_class = table[name + "_mem"] + step_fun = table[name + "_step"] + reset_fun = table[name + "_reset"] + + def wrapper_init(self): + self.mem = mem_class() + reset_fun(self.mem) + + def wrapper_reset(self): + reset_fun(self.mem) + + def wrapper_step(self, *args): + out = out_class() + step_fun(*args, out, self.mem) + return out + + return type(name.capitalize(), (), { + '__init__': wrapper_init, + 'step': wrapper_step, + 'reset': wrapper_reset, + }) + +# For now, we only create wrappers for nodes that have memory. + +for node_mem in [ n for n, o in getmembers(modules[__name__]) \ + if isclass(o) and n.endswith("_mem") ]: + node = node_mem.replace("_mem", "") + globals()[node.capitalize()] = create_node_wrapper_class(node) diff --git a/cours/inverted-pendulum/pendulum.ept b/cours/inverted-pendulum/pendulum.ept new file mode 100644 index 0000000..438367f --- /dev/null +++ b/cours/inverted-pendulum/pendulum.ept @@ -0,0 +1,223 @@ +(******************************************************************************) +(* PENDULE INVERSÉ EN HEPTAGON *) +(******************************************************************************) + +(* Le but de ce fichier est d'illustrer l'automatique par un exemple classique, + le pendule inversé. Il s'agit d'un mobile sur lequel est fixé une tringle + rigide à l'extrémité de laquelle se trouve une sphère d'une certaine + masse. Le but est de déplacer le mobile de sorte à maintenir la sphère à la + verticale. + + https://en.wikipedia.org/wiki/Inverted_pendulum + + Le présent fichier est conçu pour s'interfacer avec l'interface graphique en + Python contenue dans `main.py`. *) + +(* On commence par quelques fonctions et noeuds utilitaires. *) + +fun max(x, y : float) returns (z : float) +let + z = if x <. y then y else x; +tel + +fun min(x, y : float) returns (z : float) +let + z = if x <. y then x else y; +tel + +fun clamp(x, x_min, x_max : float) returns (y : float) +let + y = min(x_max, max(x_min, x)); +tel + +fun abs(x : float) returns (y : float) +let + y = if x >. 0.0 then x else -. x; +tel + +fun ignore_below(x, eps : float) returns (y : float) +let + y = if abs(x) <. eps then 0.0 else x; +tel + +node maintain(ini : float; c : bool; x : float on c) returns (y : float) +let + y = merge c x ((ini fby y) whenot c); +tel + +node edge(b : bool) returns (c : bool) +let + c = b -> (b <> pre b); +tel + +(* Le pendule inversé, en tant que dispositif physique, voit sa dynamique + gouvernée par des équations différentielles. Pour le simuler, nous avons + besoin de résoudre celles-ci, ou du moins d'approcher leur solution par des + moyens numériques. *) + +(* La résolution d'équations différentielles par des moyens numériques constitue + un vaste pan de l'analyse numérique. Nous allons nous contenter de solutions + très simples, dites à pas fixe, où chaque instant synchrone sera supposé + correspondre à une certaine durée de temps physique `dt`. *) + +const dt : float = 0.01 + +(* L'intégration s'effectue par la méthode des rectangles, à la manière de la + définition de l'intégrale de Riemann. Plus `dt` est petit, plus les valeurs + calculées vont être précises. *) + +node integ(x_ini, dx : float) returns (x : float) +let + x = x_ini -> (dt *. dx +. pre x); +tel + +(* La dérivée s'effectue en suivant la formule classique + + x'(t) = lim_{dt → 0} (x(t + dt) - x(t)) / dt + + sauf qu'ici on s'arrête à un `dt` fixé, qu'on espère assez petit. *) + +node deriv(x : float) returns (dx : float) +let + dx = 0.0 -> ((x -. pre x) /. dt); +tel + +(* Pour définir la physique du pendule, on a besoin de certaines constantes : la + constante gravitationnelle sur Terre, la longueur de la tige, et la masse. *) + +const g : float = 9.81 + +const l : float = 100.0 + +const m : float = 10.0 + +(* Une partie de la simulation dépend du référentiel choisi, un rectangle de + `max_w` unités de côté et `max_h` unités de haut. *) + +const max_w : int = 1200 + +const max_h : int = 300 + +(* L'équation différentielle ci-dessous gouverne la physique du pendule inversé. + Ici ℓ est la longueur du pendule, m sa masse, g la gravité à la surface de la + Terre, et θ l'angle par rapport à la verticale. + + ℓ · d²θ/d²t - m · g · sin(θ) = - d²x₀/d²t · cos(θ) + + On peut exprimer d²θ/d²t en fonction de x₀. + + d²θ/dt² = (m · g · sin(θ) + d²x₀/dt² · cos(θ)) / ℓ + + Ensuite, on intègre. + + θ = ∫² (m · g · sin(θ) + d²x₀/dt² · cos(θ)) / ℓ · dt² + + On peut maintenant approcher la solution de l'intégrale et les dérivées + numériquement, à l'aide des noeuds integ() et deriv() définis précédemment, + pour obtenir la valeur de θ en fonction de celle de x₀. + + Attention : comme θ apparaît à gauche et à droite du signe égal, on doit + introduire un délai (`fby`) pour éviter une erreur de causalité. *) + +node pendulum(d2x0 : float) returns (theta : float) +var thetap, d2theta : float; +let + thetap = 0.0 fby theta; + d2theta = (m *. g *. Mathext.sin(thetap) -. d2x0 *. Mathext.cos(thetap)) /. l; + theta = integ(0.0, integ(0.0, d2theta)); +tel + +(* Notre but est maintenant d'écrire des contrôleurs, c'est à dire des noeuds + qui vont chercher à maintenir l'angle θ à 0 en déplaçant le mobile vers la + gauche quand il est positif, et vers la droite quand il est négatif. *) + +(* En pratique, on va écrire nos contrôleurs de manière générique, comme des + noeuds cherchant à calculer une _consigne_ `y` de sorte à ramener leur entrée + d'erreur `e` à 0. *) + +(* Il faut commencer par calculer l'erreur à fournir aux contrôleurs, en + fonction de θ. Pour ce faire, on ramène les angles dans l'intervalle [π, 2π] + dans l'intervalle [-π, 0]. *) + +fun error(theta : float) returns (err : float) +var tm : float; +let + tm = Mathext.fmod(theta, 2.0 *. Mathext.pi); + err = if tm <. Mathext.pi then tm else tm -. 2.0 *. Mathext.pi; +tel + +(* On peut maintenant écrire nos contrôleurs. *) + +(* Le premier contrôleur, le contrôleur bang-bang, est le plus simple. Si + l'erreur est positive, on envoie `act`, si elle est négative on envoie + `-. act`. *) + +node bangbang(e, act : float) returns (y : float) +var eps : float; +let + eps = 0.001; + y = if e >. eps then act + else if e <. -. eps then -. act + else 0.0; +tel + +(* Le contrôleur PID agit de façon proportionnelle à l'erreur, son intégrale et + sa dérivée. Il est paramétré par trois _gains_, c'est à dire les facteurs + qu'on applique respectivement à l'erreur, l'intégrale et la dérivée. *) + +node pid(e : float; kp, ki, kd : float) returns (y : float) +let + y = kp *. e +. ki *. integ(0.0, e) +. kd *. deriv(e); +tel + +(* Le noeud ci-dessous réunit tous les composants du fichier, le pendule et son + mobile, ainsi que la détection d'erreur, et envoie les données nécessaires à + l'interface graphique. Ces données sont représentées par le type `pend`. *) + +type pend = { x0 : float; y0 : float; x : float; y : float; + mode_name : string; theta : float; error : float } + +node cart(mouse_x0 : float; mouse_diff, mode_diff : bool) returns (p : pend) +var x_ini, y_ini, x0user, x0, x, y, theta, error : float; + last d2x0 : float = 0.0; + mouse_click, mode_change : bool; + mode_name : string; +let + mouse_click = edge(mouse_diff); + mode_change = edge(mode_diff); + x_ini = Mathext.float(max_w) /. 2.0; + y_ini = Mathext.float(max_h) /. 2.0; + x0user = maintain(x_ini, mouse_click, mouse_x0 when mouse_click); + + x0 = integ(x_ini, integ(0.0, d2x0)); + theta = pendulum(d2x0); + error = 0.0 fby error(theta); + + x = x0 +. l *. Mathext.sin(theta); + y = y_ini +. l *. Mathext.cos(theta); + + p = { x0 = x0; y0 = y_ini; x = x; y = y; + mode_name = mode_name; + theta = theta; error = error }; + + automaton + state User + do mode_name = "user"; + d2x0 = x0user -. (x_ini fby x0); + unless mode_change then PID + + state PID + var kP, kI, kD : float; + do mode_name = "PID control"; + d2x0 = 10.0 *. pid(error, kP, kI, kD); + kP = 40.0; + kI = 41.0; + kD = 70.0; + unless mode_change then BangBang + + state BangBang + do mode_name = "BangBang control"; + d2x0 = bangbang(error, 50.0) + unless mode_change then User + end; +tel diff --git a/cours/inverted-pendulum/pendulum.i b/cours/inverted-pendulum/pendulum.i new file mode 100644 index 0000000..a19438e --- /dev/null +++ b/cours/inverted-pendulum/pendulum.i @@ -0,0 +1,16 @@ +%module pendulum + +%{ +#define SWIG_FILE_WITH_INIT +#include "pendulum.h" +#include "/usr/include/math.h" // hack to avoid Heptagon's math.h +%} + +// Important: makes SWIG aware of our dumb "string" type. +%include "mathext_types.h" + +%rename("%(strip:[Pendulum__])s") ""; +%include "pendulum_types.h" +%include "pendulum.h" + +%pythoncode "node_wrapping.py" diff --git a/cours/inverted-pendulum/setup.py b/cours/inverted-pendulum/setup.py new file mode 100644 index 0000000..630cf63 --- /dev/null +++ b/cours/inverted-pendulum/setup.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python + +""" +setup.py file for Pendulum example +""" + +from distutils.core import setup, Extension +from subprocess import check_output + +heptagon_headers = \ + check_output("heptc -where", shell=True) \ + .decode("utf-8") \ + .replace("\n", "/c") + +pendulum_module = Extension('_pendulum', + include_dirs = [heptagon_headers, + '.', + 'pendulum_c'], + sources=['mathext.c', + 'debug.c', + 'cutils.c', + 'pendulum_wrap.c', + 'pendulum_c/pendulum_types.c', + 'pendulum_c/pendulum.c'], + ) + +setup (name = 'pendulum', + version = '0.1', + author = "Adrien Guatto", + description = "Wrapper for Heptagon module", + ext_modules = [pendulum_module], + py_modules = ["pendulum"], + ) diff --git a/cours/journal.org b/cours/journal.org new file mode 100644 index 0000000..8e6b235 --- /dev/null +++ b/cours/journal.org @@ -0,0 +1,628 @@ +#+TITLE: Programmation synchrone 2025/2026 -- Journal du cours +#+AUTHOR: Adrien Guatto +#+EMAIL: guatto@irif.org +#+LANGUAGE: fr +#+OPTIONS: ^:nil p:nil +#+LATEX_CLASS: article +#+LATEX_CLASS_OPTIONS: [a4paper,11pt] +#+LATEX_HEADER: \usepackage{a4wide} +#+LATEX_HEADER: \usepackage{microtype} +#+LATEX_HEADER: \hypersetup{hidelinks} +#+LATEX_HEADER: \usepackage[french]{babel} +# (org-latex-export-to-pdf) + +* Cours 01 <2025-09-29> + On lit le syllabus ainsi que le début des notes. + + Le manuel d'Heptagon peut être trouvé sur le dépôt git du projet. On peut + cliquer sur le bouton /download/ pour télécharger une copie PDF. + + https://gitlab.inria.fr/synchrone/heptagon/-/blob/ec26be27b91f3e601b98b8b7e15e8d56d4b9afc7/manual/heptagon-manual.pdf +* Cours 02 <2025-10-06> + On traite des constructions de base du langage. +* Cours 03 <2025-10-13> + On traite des horloges et automates. +* Cours 04 <2025-10-20> + On traite des automates et tableaux. +* Cours 05 <2025-10-27> +** Programmation Heptagon + On traite des tableaux. +** Exemples + On traite la génération d'échantillons audios, disponbile dans son + [dossier](audio/). +* Cours 06 <2025-11-04> +** Projet + On traite du fonctionnement du projet. +* Cours 07 <2025-11-10> +** Exemples + On traite l'exemple du pendule inversé, disponible dans son + [dossier](inverted-pendulum/). +** Implémentation des langages synchrones + On commence à expliquer la compilation d'Heptagon, avec l'explication de + MiniLS vers Obc jusqu'au traitement des horloges *exclu*. +* Cours 08 <2025-11-17> +** Implémentation des langages synchrones + On discute de la compilation des horloges en MiniLS, de la traduction des + automates et des structures de contrôle plus simples (notamment ~switch~) en + Heptagon. On présente l'analyse d'initialisation. +* Cours 09 <2025-11-24> +** Implémentation des langages synchrones + On termine la partie concernant l'implémentation des langages synchrones. +** Vérification automatique + L'outil Kind2 est un vérificateur de modèle (/model checker/) pour un langage + proche d'Heptagon développé à l'université de l'Iowa. + + https://kind.cs.uiowa.edu + + Un tel outil est une des justifications pour l'emploi d'un langage très + limité comme Heptagon, SCADE et leurs variantes. La vérification automatique + est beaucoup plus difficile lorsque le langage d'entrée est un langage + généraliste comme C, Python ou OCaml. +*** Introduction + Le langage accepté par Kind2 est une variante d'Heptagon beaucoup plus + proche de son ancêtre Lustre. Les différences sont assez mineures mais + peuvent rendre l'écriture du code désagréable lorsqu'on est habitué aux + constructions d'Heptagon. + + Par exemple, le langage ne propose pas l'opérateur ~fby~, uniquement ~pre~ + et ~->~. Notre classique compteur doit donc être écrit comme suit. + + #+BEGIN_SRC +node counter() returns (o : int) +let + o = 0 -> (pre o + 1); +tel + #+END_SRC + + Kind2 est capable de traduire ce programme vers du code impératif, jouant + ainsi le même rôle que le compilateur ~heptc~ pour Heptagon. Plutôt que de + générer du code C, Kind2 produit du code Rust. + + #+BEGIN_SRC +$ kind2 --compile true ex0.lus + kind2 v1.1.0-16-g9572d9b + + System counter has no property, skipping verification step. + +-------------------------------------------------------------------------------- +Post-analysis: rust generation + + Compilation to Rust is still a rather experimental feature: + in particular, arrays are not supported. + + Compiling node "counter" to Rust in `ex0.lus.out/counter/implem`. + + Done compiling. + +-------------------------------------------------------------------------------- +$ wc -l ex0.lus.out/counter/implem/src/main.rs +366 ex0.lus.out/counter/implem/src/main.rs + #+END_SRC + + L'intérêt principal de Kind2 est ailleurs : il est capable de /vérifier + formellement/ des propriétés. Par opposition au test, la vérification + formelle permet de vérifier qu'une propriété est vraie pour *toute* entrée, + plutôt que pour une entrée particulière. De plus, dans le cas des langages + synchrones, les entrées sont des suites infinies, dont on ne peut tester + qu'un préfixe fini. La vérification formelle, elle, considère des entrées de + longueur infinie. + + Lorsqu'on découvre un outil de vérification formelle, un des premières + questions à se poser est celle du langage dans lequel les propriétés sont + écrites. Certains outils de vérification offrent un langage spécialisé pour + ce faire, par exemple une logique temporelle pour le model-checker PRISM. + Les outils de vérification de langages synchrones comme Kind2 font un choix + un peu différent. Le langage des propriétés est le même que celui des + programmes, et la propriété à vérifier est exprimée par un flot booléen + désigné par l'utilisateur doit toujours être vrai. Avec cette approche, on + va avoir tendance à écrire des /observateurs synchrones/, c'est-à-dire des + noeuds qui sont utilisés uniquement pour vérifier que le noeud qui va être + déployé sur le système embarqué final va respecter certaines propriétés, + éventuellement sous certaines hypothèses concernant ses entrées. + + #+BEGIN_SRC +node top() returns (ok : bool) +var o : int; +let + o = counter(); + ok = o >= 0; + --%PROPERTY ok; +tel + #+END_SRC + + En lançant Kind2 sur ce fichier, on lui demande de vérifier que le flot ~ok~ + est toujours vrai. Il est capable de vérifier que c'est le cas, même si ~ok~ + contient une infinité de valeurs. + + #+BEGIN_SRC +$ kind2 ex0.lus + kind2 v1.6.0 + +==================================================================================================================== +Analyzing top + with First top: 'top' + subsystems + | concrete: counter + + Property ok is valid by property directed reachability after 0.148s. + +-------------------------------------------------------------------------------------------------------------------- +Summary of properties: +-------------------------------------------------------------------------------------------------------------------- +ok: valid (at 1) +==================================================================================================================== + #+END_SRC + + (Cet exemple montre que Kind2 suppose que la sémantique du type ~int~ de + Lustre est celle des entiers relatifs, sans limite de taille. Pourquoi ? + Remarquons que jusqu'ici on avait supposé que les types des langages + synchrones étaient tous finis, ce qui assurait la finitude de la mémoire. + Comme on le décrira plus bas, ce choix n'est pas sans conséquences + théoriques et pratiques.) + + Cet exemple laisse aussi deviner que la vérification est un processus + coûteux, beaucoup plus coûteux que la compilation. Ce fichier étant écrit + dans le sous-ensemble commun à Heptagon et Lustre, on peut le compiler avec + ~heptc~ et comparer le temps de compilation au temps de vérification. + + #+BEGIN_SRC +$ time heptc -target c ex0.lus +real 0m0,026s +user 0m0,013s +sys 0m0,013s + #+END_SRC + + On peut montrer des propriétés un peu plus intéressantes, par exemple en + écrivant une variante du compteur avec une valeur d'initialisation. + + #+BEGIN_SRC +node counterini(ini : int) returns (o : int) +let + o = ini -> (pre o + 1); +tel + +node top(ini : int) returns (ok : bool) +var o : int; +let + o = counterini(ini); + ok = o >= 0; + check ok; +tel + #+END_SRC + + Cet exemple est plus intéressant que le précédent en ce que le noeud ~top~ + dispose désormais d'une entrée. Kind2 va donc essayer de vérifier que toutes + les valeurs de la suite ~ok~ sont vraies et ce /pour toutes les suites + d'entrée ~ini~/. Va-t-il réussir ? + + #+BEGIN_SRC shell +$ kind2 ex0.lus + kind2 v1.6.0 + +==================================================================================================================== +Analyzing top + with First top: 'top' + subsystems + | concrete: counterini + + Property ok is invalid by bounded model checking for k=0 after 0.041s. + +Counterexample: + Node top () + == Inputs == + ini -1 + == Outputs == + ok false + == Locals == + o -1 + + Node counterini (top[l14c7]) + == Inputs == + ini -1 + == Outputs == + o -1 + + +-------------------------------------------------------------------------------------------------------------------- +Summary of properties: +-------------------------------------------------------------------------------------------------------------------- +ok: invalid after 0 steps +==================================================================================================================== + #+END_SRC + + La vérification a échoué, mais Kind2 nous a fourni un contre-exemple sous la + forme d'une séquence d'entrées qui falsifie la propriété. Il suffit bien sûr + de donner un ~ini~ strictement négatif au premier instant ! + + Notre propriété est donc fausse. Peut-on exprimer que le flot ~o~ est + toujours vrai mais *à condition* que la première valeur de ~ini~ soit + positive ? Oui, en utilisant l'opérateur booléen d'implication pour exprimer + une /précondition/ (hypothèse) sur l'entrée ~ini~. + + (Remarque : l'implication dont il s'agit ici n'est que l'opérateur + combinatoire sur les flots défini par ~x => y = not x or y~, rien de plus. + En particulier, Kind2 n'a aucun traitement spécifique de cet opérateur.) + + Un premier essai pourrait ressembler à ce qui suit. + + #+BEGIN_SRC +node top(ini : int) returns (ok : bool) +var o : int; +let + o = counterini(ini); + ok = ((ini >= 0) -> true) => (o >= 0); + check ok; +tel + #+END_SRC + + La vérification échoue encore, avec le contre-exemple qui suit. + + #+BEGIN_VERBATIM +Counterexample: + Node top () + == Inputs == + ini -2 0 + == Outputs == + ok true false + == Locals == + o -2 -1 + + Node counterini (top[l19c7]) + == Inputs == + ini -2 0 + == Outputs == + o -2 -1 + #+END_VERBATIM + + Comment comprendre ce contre-exemple ? Au premier instant, ~i~ vaut ~-2~ et + la précondition sur ~ini~ est fausse, donc l'implication est vraie. Au + deuxième instant, la précondition sur ini est vraie, donc l'implication + prend la valeur de ~o >= 0~, qui se trouve fausse. + + Quel est le problème ? Ce que nous avions en tête n'est pas ce que nous + avons écrit. On aurait voulu dire "si ~ini < 0~ au premier instant, la + précondition est fausse pour toujours", mais on a écrit "si ~ini < 0~ au + premier instant, la précondition est fausse /au premier instant/". On peut + régler le problème en introduisant un noeud ~always~ dont la sortie devient + fausse pour toujours dès que son entrée est fausse. + + #+BEGIN_SRC +node always(b : bool) returns (o : bool) +let + o = (true -> pre o) and b; +tel + +node top(ini : int) returns (ok : bool) +var o : int; +let + o = counterini(ini); + ok = always((ini >= 0) -> true) => (o >= 0); + check ok; +tel + #+END_SRC + + La vérification réussit. En particulier, le contre-exemple précédent n'est + plus valide. Si la valeur initiale de ~ini~ est ~-2~, la valeur de + ~always((ini >= 0) -> true)~ sera toujours fausse, et donc la vérité de + l'implication dépendra de ~o >= 0~, toujours vraie. + + Le fait d'exprimer une précondition sur les entrées est très fréquent, et + Kind2 propose une syntaxe un peu plus agréable pour ce faire. + + #+BEGIN_SRC +node top(ini : int) returns (ok : bool) +var o : int; +let + o = counterini(ini); + ok = o >= 0; + check ok provided always((ini >= 0) -> true); +tel + #+END_SRC +*** Comparaison de noeuds + On peut utiliser le principe des observateurs synchrones pour comparer le + comportement de deux noeuds, et éventuellement tester leur équivalence. + + Par exemple, comment quelle relation voyez vous entre notre premier noeud + ~counter~ et le noeud ~counterini~ ? Et comment la feriez-vous vérifier par + Kind2 ? + + Le lien est essentiellement qu'ajouter la première valeur de ~ini~ à la + sortie de ~counter~ donne celle de ~counterini~. La tentative donnée + ci-dessous est trop naïve puisqu'il ne faut prendre en compte ~ini~ + uniquement au premier instant. + + #+BEGIN_SRC +node compare(ini : int) returns (ok : bool) +let + ok = counter() + ini = counterini(ini); + check ok; +tel + #+END_SRC + + On peut utiliser un noeud ~constant~ qui maintient la valeur reçue au + premier instant. On obtient le noeud suivant, dont la propriété est + vérifiée. + + #+BEGIN_SRC +node constant(ini : int) returns (o : bool) +let + o = ini -> pre o; +tel + +node compare(ini : int) returns (ok : bool) +let + ok = counter() + constant(ini) = counterini(ini); + check ok; +tel + #+END_SRC +*** Les contrats + Les propriétés qu'on a vérifié jusqu'ici sont des propriétés globales du + système. Il s'agit certainement d'une approche utile, mais elle n'est pas + tout à fait convaincante du point de vue du génie logiciel. On aimerait + plutôt vérifier chaque noeud indépendamment. Kind2 le permet via l'écriture + de /contrats/ qui permettent à chaque noeud de documenter ses attentes + vis-à-vis de ses appelants (mot-clef /assume/) et ses garanties (mot-clef + /guarantee/). + + #+BEGIN_SRC +node counterini(ini : int) returns (o : int) +(*@contract + assume (ini >= 0) -> true; + guarantee o >= 0; + *) +let + o = ini -> (pre o + 1); +tel + #+END_SRC + + Notez que Kind2 assure la bonne sémantique pour l'hypothèse (/assume/), + c'est-à-dire qu'elle doit être valide de toute la séquence d'entrées, et ce + sans avoir à utiliser une noeud comme ~always~. (Pour ceux et celles d'entre + vous qui connaissent la logique temporelle, si $A$ est la formule booléenne + des hypothèses et $G$ est la formule booléenne des garanties, la formule que + l'outil cherche à démontrer est $\square A \Rightarrow \square G$ plutôt que + $\square A \Rightarrow \square G$. + + Je finis cette partie pratique par une difficulté à laquelle tous les + utilisateurs de vérification formelle finissent par se heurter un jour ou + l'autre : des hypothèses trop fortes rendent la garantie triviale. Voir + ci-dessous pour un exemple caricatural ; /garbage in, garbage out/. + + #+BEGIN_SRC +node counterini(ini : int) returns (o : int) +(*@contract + assume false; + guarantee (o = constant(o)); + *) +let + o = ini -> (pre o + 1); +tel + #+END_SRC + + Cependant, Kind2 aide à détecter le problème. L'outil va chercher à vérifier + que tout appelant à ce ~counterini~ satisfait aux hypothèses, ce qui ici est + impossible... Sauf si l'appelant a lui-même fait des hypothèses + contradictoires sur ses entrées, bien entendu ! +* Cours 10 <2025-12-01> +** Vérification automatique, suite +*** Exemples +**** Compteurs + On va s'intéresser à la vérification d'équivalence entre compteurs + périodiques de 0 à 3, formulé avec des entiers, puis en codage binaire, + classique puis de Gray. + + On commence par écrire le compteur sur deux bits avec un entier + + #+BEGIN_SRC +node intcounter (rst : bool; const max : int) returns (out : bool); +var t : int; +let + t = 0 -> if rst or pre t = max then 0 else pre t + 1; + out = t = 3; +tel + +node graycounter (rst : bool) returns (out : bool); +var a, b : bool; +let + a = false -> (not rst and not pre b); + b = false -> (not rst and pre a); + out = not a and b; +tel + +node bincounter (rst : bool) returns (out : bool); +var a, b : bool; +let + a = false -> (not rst and (pre a <> pre b)); + b = false -> (not rst and not pre b); + out = a and b; +tel + +node top (rst : bool) returns (); +var grayc, binc, intc : bool; +let + grayc = graycounter(rst); + intc = intcounter(rst, 3); + binc = bincounter(rst); + check binc = intc; + check grayc = intc; +tel + #+END_SRC +*** Les bases de l'algorithmique de la vérification + Les données du problème sont : + + 1. un programme synchrone à vérifier, + + 2. un /observateur/ pour la propriété à démontrer (~check ok~), + + 3. un /observateur/ pour l'hypothèse supposée (~assert hyp~). + + L'objectif d'un vérificateur de modèle pour un programme synchrone à flots + de données est de démontrer : + + ~always(so_far(hyp) => ok)~. + + Notons qu'on ne traite que les propriétés de sûreté (``quelque-chose de + mauvais n'arrive toujours'') et pas de vivacité (``quelque-chose de bon + finit toujours par arriver''), comme expliqué au cours précédent. + + On va s'intéresser à ce problème dans le cas *purement booléen*, + c'est-à-dire où toutes les variables manipulées sont des suites de + booléens. Notre description est basée sur le rapport de recherche de Pascal + Raymond disponible à + [[https://www.di.ens.fr/~pouzet/cours/mpri/bib/lesar-rapport.pdf][cette + adresse]]. + + Si $X$ est un ensemble, on écrit $\overline{X}$ pour les fonctions de $X$ + dans $\mathbb{B}$, qui sont en bijection avec les parties de $X$. On + identifiera librement les sous-ensembles et leurs fonctions + charactéristiques dans ce qui suit. + + Les données du problème sont les n-uplets ($V$, $S$, $I \subseteq + \overline{S}$, $g : \overline{S} \times \overline{V} \to \overline{S}$, $H + \subseteq \overline{S} \times \overline{V}$, $\Phi \overline{S} \times + \overline{V}$), où : + + - $S$ et $V$ sont des ensembles *finis* dont les éléments sont + respectivement les variables d'état de d'entrée du programme, + + - $I$ est l'ensemble des états initiaux, + + - $H$ est l'hypothèse, + + - $\Phi$ est la propriété à démontrer. + + On écrira $q \to^v q'$ pour $q' = g(q, v)$ avec $q,q' \in \overline{S}$ et + $v \in \overline{V}$. + + On définit les états successeur/prédécesseurs d'un état $q$ ou d'un + ensemble d'états~$X$ par : + + \begin{align*} + post_H(q) & = \{ q' \mid \exists v \in q \to_v q' \wedge H(q, v) \} + \\ + pre_H(q) & = \{ q' \mid \exists v \in q' \to_v q \wedge H(q', v) \} + \\ + POST_H(X) & = \bigcup_{q \in X} post_H(q) + \\ + PRE_H(X) & = \bigcup_{q \in X} pre_H(q) + \end{align*} + + L'ensemble des états d'erreurs immédiate sont ceux qui satisfont + l'hypothèse mais pas la propriété à démontrer. + + \[ + Err = \{ q \mid \exists v, H(q, v) \wedge \neg \Phi(q, v) \} + \] + + On définit maintenant l'ensemble $Acc$ des états accessibles du programme + et l'ensemble $Bad$ des mauvais états par un plus petit point fixe + ensembliste : $\mu X . F(X)$ désigne le plus petit ensemble $Y$ tel que + $F(Y) \subseteq Y$ ; son existence est garantie dès lors que $F$ est + croissante pour l'inclusion. + + \begin{align*} + Acc & = \mu X. (I \cup POST_H(X)) + \\ + Err & = \{ q \mid \exists v, H(q, v) \wedge \neg \Phi(q, v) \} + \\ + Bad & = \mu X . (Err \cup PRE_H(X)) + \end{align*} + + L'objectif est maintenant de démontrer que le système est *sûr*, + c'est-à-dire que $Acc \cap Bad = \emptyset$. Autrement dit, aucun état + accessible du système n'est un état d'erreur. + + Pour ce faire, il existe deux grandes méthodes : + + - l'exploration en *avant* : on démontre $Acc \cap Err = \emptyset$, + + - l'exploration en *arrière* : on démontre $Bad \cap I = \emptyset$. + + L'approche la plus simple pour ce faire est une approche énumérative : on + explore explicitement l'automate induit par le programme, en imitant le + calcul du plus petit point fixe correspondant à $Acc$ ou $Bad$. + + #+BEGIN_SRC +Known := I; Explored := {} +while there is q in Known \ Explored: + for all q' in post(q): + if q' in Err: + return FAILURE + move q' to Known + move q to Explored +} +return SUCCESS + #+END_SRC + + Cet algorithme peut être implémenté et termine puisque tous les ensembles + manipules sont finis, les prédicats calculables, tout comme la fonction de + transition du programme. En particulier, vous pouvez implémenter les + ensembles ~Known~ et ~Explored~ par votre structure d'ensemble fini + favorite, et énumérer les valeurs possibles pour calculer ~post(q)~. Il + s'agit donc en somme d'exécuter le programme sur toutes les entrées + possibles ; on voit que c'est très inefficace. + + On peut aussi implémenter une version en arrière. + + #+BEGIN_SRC +Known := Err; Explored := {} +while there is q in Known \ Explored: + for all q' in pre(q): + if q' in I: + return FAILURE + move q' to Known + move q to Explored +} +return SUCCESS + #+END_SRC + + Les deux algorithmes sont très inefficaces, tellement qu'ils en deviennent + inutilisables en pratique. En effet, les structures de données + représentant les ensembles d'état ou d'entrées deviendront très grosses, + et les calculs attenants (différence ensembliste, ajout, choix d'un + élément...) coûteux à réaliser. + + Notons que le second algorithme est encore moins efficace que le + premier. Voyez-vous pourquoi ? Le calcul de ~pre(q)~ est encore plus + coûteux, puisqu'il faut calculer l'image inverse de la fonction de + transition $g$, c'est-à-dire trouver l'ensemble des $q'$ tel que $g(q', v) + = q$ pour $v$ et $q$ fixés. + + Depuis les années 1990, on utilise plutôt des méthodes dites + *symboliques*. L'idée générale est de représenter les ensembles finis par + des formules logiques plutôt que par des énumérations explicites. Par + exemple, on pourrait représenter l'ensemble $\{ (x, y, z) \mid x + \overline{y} + z = 1 \}$ directement par le *code* de sa fonction + caractéristique $f(x, y, z) = x \overline{y} + z$ plutôt qu'explicitement + par l'ensemble de triplets $\{ (1, 0, 0), (1, 0, 1), (0, 0, 1), (0, 1, 1), + (1, 1, 1) \}$. Les méthodes symboliques sont donc un exemple de compromis + temps-espace : on y gagne une représentation nettement plus succincte des + données que la représentation explicite, en contrepartie de quoi les + opérations (notamment test d'égalité, du vide, d'inclusion, etc.) seront + plus coûteuses. + + Il existe deux grandes méthodes d'implémentation pour exprimer les + opérations et les tests : + + 1. On peut représenter les formules directement par leur code et utiliser + un solveur SAT général pour répondre aux questions intéressantes pour + la vérification. Par exemple, l'ensemble $\{ \vec{x} \mid \phi(\vec{x}) + \}$ est inclus dans l'ensemble $\{ vec{x} \mid \psi(\vec{x}) \}$ si et + seulement si $\forall \vec{x}, \phi(\vec{x}) \implies \psi(\vec{x})$, + question à laquelle un solveur SAT peut répondre, y compris pour de + très grosses formules $\phi$ et $\psi$. D'autres questions, comme par + exemple l'énumération des valeurs appartenant à un ensemble, sont d'un + traitement moins direct mais peuvent être gérée avec un peu de + réflexion et d'effort. + + 2. On peut utiliser une représentation canonique des formules booléennes, + par exemple les /diagrammes de décision binaire/ (ou /binary decision + diagrams/ en anglais, BDD). Ici, toutes les formules équivalentes sont + représentées par la même structure de données. La construction de cette + représentation est donc forcément coûteuse, mais une fois construite on + peut l'utiliser pour répondre efficaement à beaucoup de questions. C'est + cette approche qui est très bien exposée dans le rapport de P. Raymond + cité plus haut. + + Le principe des diagrammes de décision binaire a été brièvement expliqué en + cours mais ne sera pas détaillé dans ces notes. diff --git a/examens/2024/examen-2024.pdf b/examens/2024/examen-2024.pdf new file mode 100644 index 0000000..8bad2c2 Binary files /dev/null and b/examens/2024/examen-2024.pdf differ diff --git a/notes/notes-de-cours-progsync.pdf b/notes/notes-de-cours-progsync.pdf new file mode 100644 index 0000000..e6f8df1 Binary files /dev/null and b/notes/notes-de-cours-progsync.pdf differ diff --git a/src/city.ept b/src/city.ept index b20aa9a..ffd6ef5 100644 --- a/src/city.ept +++ b/src/city.ept @@ -66,13 +66,15 @@ let false); tel +(* The car angle with the road should be in the range -120..120 degree *) node wrong_dir(ph : phase) returns (wrong : bool) -var data : map_data; error : float; ang : float; +var data : map_data; norm, error : float; ang : float; let data = lookup_phase(ph); ang = ph.ph_head *. (pi /. 180.0); + norm = Mathext.sqrt(data.dir_x *. data.dir_x +. data.dir_y *. data.dir_y); error = data.dir_x *. Mathext.cos(ang) +. data.dir_y *. Mathext.sin(ang); - wrong = error <. -. 0.5; + wrong = error <. -. 0.5 *. norm; tel fun aggregate_events(lightRun, speedExcess, exitRoad, collisionEvent, @@ -188,7 +190,7 @@ tel (* Scoring *) node scoringA(e : event; rstatus : status) returns (score : int) -var penalty : int; collision_count : int; +var penalty : int; let score = (10000 fby score) + if Utilities.rising_edge(rstatus = Arrived) then 1000 else 0 @@ -198,9 +200,7 @@ let + (if e.speedExcess then -2 else 0) + (if Utilities.rising_edge(e.exitRoad) then -5000 else 0) + (if Utilities.rising_edge(e.dirEvent) then -2000 else 0) - + (if collision_count = 0 then -500 else 0) - + (if collision_count < 0 then 10 else 0); - collision_count = Utilities.countdown(not e.collisionEvent, 20); + + (if Utilities.rising_edge(e.collisionEvent) then -500 else 0); tel node scoringB(ph : phase) returns (score : int) diff --git a/src/map.c b/src/map.c index b2a08d7..170f1f0 100644 --- a/src/map.c +++ b/src/map.c @@ -683,29 +683,6 @@ bool isOnRoadArc(/* IN */ return true; } -int isOnTLight(int x, int y, int rd, float dir_x, float dir_y) -{ - if (map->tlight_arr == NULL || - map->tlight_sz <= 0) - return -1; - //no traffic light - - for (int i = 0; i < map->tlight_sz; i++) { - tlight_t *tl = &map->tlight_arr[i]; - if (tl->road != rd) - continue; - double d = distance(x, y, tl->tl.ptl_pos.x, tl->tl.ptl_pos.y); - - /* double cosdir = dirCos(tl->tl.ptl_pos.y - y, x - tl->tl.ptl_pos.x, - dir_x, dir_y); */ //MS - double cosdir = dirCos(tl->tl.ptl_pos.x-x, tl->tl.ptl_pos.y-y, - dir_x, dir_y); //EA - if (d < TL_VIEW && cosdir > TL_COSDIR) - return i; - } - return -1; -} - bool isAfterStop(int x, int y, int rid, float dir_x, float dir_y, int* tl) { (*tl) = -1; @@ -822,7 +799,7 @@ isPositionOnPoint(float x, float y, road_t* rd, position_t* p, double pWidth) } } -Globals__color getColorPoint(int rid, float x, float y) { +Globals__color getColorPoint(int rid, float x, float y, int* tl) { // first go through waypoints log_debug("[geometry] looking for waypoints at (%.2f, %.2f) on road %d\n", x, y, rid); @@ -855,6 +832,7 @@ Globals__color getColorPoint(int rid, float x, float y) { // stop points are ruban if (isPositionOnPoint(x, y, rd, &sp->position, RD_SIZE_STOP)) { log_debug("[geometry] (%.2f, %.2f) at stop %d\n", x, y, i); + *tl = sp->sema; return COL_STOP; //one stop by position } @@ -909,28 +887,22 @@ DEFINE_HEPT_FUN(Map, lookup_pos, (Globals__position pos)) { out->data.dir_x = dir_X; out->data.dir_y = dir_Y; out->data.max_speed = map->road_arr[min_rd].max_speed; - /* Update color when a waypoint or stop. */ - col = getColorPoint(rid, x, y); - if (colors_equal(&col, &COL_OUT)) - out->data.color = out->data.color; - else if (colors_equal(&col, &COL_STOP)) { - /* TODO: update red color */ - out->data.color = col; - } - else { - /* TODO: update green color */ - out->data.color = col; - } } } + /* Update color when a waypoint or stop. */ + int tl = -1; + if (min_rd >= 0) { + Globals__color col = getColorPoint(min_rd, x, y, &tl); + if (!colors_equal(&col, &COL_OUT)) + out->data.color = col; + } + /* Compute the return type. */ if (min_rd >= 0) { log_debug("[geometry] (%.2f, %.2f) is on road %d\n", x, y, min_rd); out->data.on_road = true; - int tl = -1; - out->data.tl_number = - isOnTLight(x, y, min_rd, out->data.dir_x, out->data.dir_y); + out->data.tl_number = tl; out->data.tl_required = isAfterStop(x, y, min_rd, out->data.dir_x, out->data.dir_y, &tl); if(out->data.tl_required) { diff --git a/src/simulation_loop.c b/src/simulation_loop.c index 12c3965..2066011 100644 --- a/src/simulation_loop.c +++ b/src/simulation_loop.c @@ -157,7 +157,7 @@ race_result_t simulation_loop(bool show_guide, SDL_WINDOWPOS_UNDEFINED, MAX_X, MAX_Y, - SDL_WINDOW_SHOWN)) == NULL) + SDL_WINDOW_SHOWN|SDL_WINDOW_RESIZABLE)) == NULL) log_fatal("[sdl] could not open window (%s)\n", SDL_GetError()); r = SDL_CreateRenderer(w, -1, SDL_RENDERER_ACCELERATED); @@ -256,6 +256,13 @@ race_result_t simulation_loop(bool show_guide, case SDL_QUIT: quit = true; break; + case SDL_WINDOWEVENT: + if (e.window.event == SDL_WINDOWEVENT_RESIZED) { + float newwidth = e.window.data1; + float newheight = e.window.data2; + SDL_RenderSetScale(r, newwidth / MAX_X, newheight / MAX_Y); + } + break; case SDL_KEYDOWN: switch (e.key.keysym.sym) { case SDLK_q: diff --git a/sujet/sujet-projet.pdf b/sujet/sujet-projet.pdf index cb78c3f..a93d392 100644 Binary files a/sujet/sujet-projet.pdf and b/sujet/sujet-projet.pdf differ diff --git a/tp/TP03.pdf b/tp/TP03.pdf new file mode 100644 index 0000000..e4c2327 Binary files /dev/null and b/tp/TP03.pdf differ