%% COURTE INTRODUCTION A MATLAB 
% Ceci est une courte introduction a Matlab. 
%
%% Mise en route
% Pour creer un raccourci de l'application Matlab, ajouter dans votre
% fichier |.profile| les lignes suivantes
%
% |# Set PATH variable|
%
% |export PATH=/Applications/MATLAB_R2012b.app/bin:$PATH|
%
% Si le fichier |.profile| n'existe pas vous pouvez le creer. Relancez
% ensuite la fenetre de terminal pour que le raccourci fonctionne.
%
% A l'invite de commande dans votre fenetre de terminal, tapez |matlab|
% pour lancez Matlab en mode graphique. Pour lancez Matlab en ligne de
% commande seulement (sans fenetre), tapez |matlab -nodesktop|. Ce tutoriel suppose que
% vous avez lance Matlab en mode graphique. 

%% Elements de l'interface graphique
% La fenetre Matlab est divisee en plusieurs espaces qui peuvent etre
% configures. L'espace le plus important est la |Command Windows|. Cette
% fenetre permet d'executer des commandes Matlab. L'autre espace important
% est l'editeur de texte, qui vous permet d'editer des fichier qui pourront
% etre executes par Matlab.

%% Philosophie
% Matlab, qui est l'abbreviation de matrix laboratory, est optimise pour
% manipuler des matrices. Pour en tirer partie, il faut vectorialiser les
% expressions. En C, pour initialiser un vecteur |x| qui prendra les
% valeurs de l a 1000000 (un million), on procede de la facon suivante:
%
%   double *x;
%   int i, sizex = 1000000;
%   x = malloc(size*sizeof(double));
%   for ( i = 0; i < sizex; i++)
%   {
%       x[i] = i+1.0;
%   }
%
% Sous Matlab, la meme procedure se fait de la facon suivante:

x = 1:1000000;

%%
% L'initialisation est vectorialisee: le membre de droite (|1:1000000|) est
% le vecteur 1 a 1000000 est l'assignation se fait de vecteur a vecteur.

%%
% Notez qu'il n'est pas necessaire de declarer la variable |x|, elle est
% cree lors de l'initialization. Le type sera automatiquement un array de
% |double|:

class(x)

%%
% Contrairement a C, les expressions n'ont pas a se terminer par un ;. Si
% le ; est omis, le resultat de l'expression sera affiche dans la Command
% Windows. L'expression 

1+1;

%%
% ne produit aucin effet (mais le resultat est stocke dans une variable
% appelee |ans|). L'expression 

1+1

%%
% renvoie le resultat escompte: 2.


%%
% Les valeurs scalaires sont des arrays |1x1|, et ne sont pas different des
% arrays representant des vecteurs ou des matrices.

%% 
% On aurait pu initialiser |x| comme en C, mais cela aurait ete beucoup
% plus long. Les commandes |tic| et |toc| permettent de calculer un
% intervalle de temps ecoule:

clear x; % on efface la variable x
tic
x = 1:1000000;
toc

%%
% Maintenant, avec la boucle |for|

clear x;
tic
for i = 1:1000000
    x(i) = i;
end
toc

%%
% La boucle for est beaucoup plus lente que l'assignation vectorielle!
% Conclusion, si on peut, *on evite les boucles |for| en Matlab*. 

%% Variables, vecteurs, matrices
% La type de variable de base de Matlab est l'array de double. Il existe
% d'autres types (|class| en Matlab) comme les |char|, |cell|, |struct|,
% |int8|, etc, mais les |double| sont utilisees la plupart du temps. Comme
% les variables creees sont de classe |double| par defaut, on suppose qu'on
% travaillera en |double| dans ce tutoriel. Un array de double |A| peut
% comporter plusieurs dimensions. Pour un array de dimension 2, c.-a-d. un
% matrice (le defaut), on accede a l'entree ligne |i|, colonne |j| par
% A(i,j). (Comparer a C: |A[i-1][j-1]|.)
%
% Voici un liste de commandes pour creer des vecteurs (|mx1| pour un
% vecteur colonne, et |1xn| pour un vecteur ligne).

m = 10; n = 10;
xcol  = ones(m,1);      % cree un vecteur colonne de dimension 10 avec UNs partout
xline = ones(1,n);      % cree un vecteur ligne de dimension 10 avec UNs partout
ycol = (1:m)';          % cree un vecteur colonne des entiers de 1 a m. Le ' est l'adjoint
yline = 1:n;            % cree un vecteur ligne des entiers de 1 a n.
z1 = linspace(0,1,n);   % cree le vecteur ligne des n doubles espaces lineairement entre 0 et 1.
h = 1/(n-1);            % definit la distance entre deux elements de z1
z2 = 0:h:1;             % autre facon de definir le meme vecteur, mais attention aux erreurs d'arrondis!
isequal(z1,z2)          % on teste pour savoir si |z1 == z2|, le resultat est FAUX.


%% 
% pour connaitre les dimensions des variables |xcol| et |xline|, on peut
% utiliser les commandes |length| et |size|. La commande |size| renvoie la
% taille de chaque dimension

size(xcol)
size(xline)

%%
% La commande |length| renvoie la plus grande taille seulement, ce qui est
% pratique pour des expressions du type
%
%   mcols = length(xcol);
%
% L'autre facon de connaitre l'etat des variables est la commande |whos|

whos

%%
% Cette commande affiche toutes les variables du workspace courant. Le workspace
% contient les variables definies par l'utilisateur.

%% 
% De la meme facon que les vecteurs, on peut creer des matrices |mxn|.

A  = ones(m,n);         % cree une matrice 10x10 avec UNs partout
Ar = rand(m,n);         % cree une matrice 10x10 avec nombres aleatoires entre 0 et 1
A  = vander(z1);        % cree une matrice carree de Vandermonde de taille length(z1) 

%% Aide 
% Pour de l'aide sur une commande, tapez
%
%   help nom_de_la_commande
%
% Par exemple, pour afficher l'aide sur la commande |ones|

help ones

%%
% La ligne |See also| est utile pour trouver des commandes similaire. Ici
% on trouve les commandes |zeros| et |eye| qui permettent d'initialiser des
% matrices (ou vecteurs) de zeros et la matrice identite.

%% Principales commandes et syntaxe sur les vecteurs et matrices
%
%   a = 2.0         % assigne la valeur 2 a la varialble a
%   (a == 2)        % teste l'egalite entre 2 et a
%   find(x == 200)  % trouve les entrees i qui verifient l'expression logique x(i) == 200
%   find(x < 10)    % trouve les entrees i qui verifient l'expression logique x(i) < 10
%   A*xcol          % multiplication matrice*vecteur
%   xline*A         % multiplication vecteur*matrice
%   A\x             % division a gauche (equivalent a A^(-1)*x)
%   A(:,1:2:end)    % colonnes impaires de A
%   A([1 2 5 6],:)  % lignes 1, 2, 4, 5, 6

%% Fonctions et scripts
% Les fonctions Matlab sont le plus souvent ecrites dans un fichier texte
% separe. Le nom du fichier est le nom de la fonction. La fonction du
% fichier |mafonction.m| sera donc appelee par la commande |mafonction|.
% Construisons un fonction qui prend en argument deux vecteurs de meme
% taille, x et y, et qui renvoie les coefficients du polynome de degre |n|
% qui interpole au sens des moindres carres y aux points x: y ~ p(x). 
% 
%   edit monpolyfit.m
%
% La premiere ligne doit etre
%
%   function p = monpolyfit(x,y,n)
%
% Le mot cle |function| indique que le fichier est une fonction. |p| est la
% la variable retournee par la fonction et |(x,y,n)| sont les arguments de
% la fonction. La notation |function [p1, p2] = mafonction(x,y)| est utilisee s'il y a
% plusieurs variables de sortie et |function mafunction(x,y)| est utilisee
% si il n'y a pas de variable de sortie (pour une fonction qui genererait
% un graphique ou un fichier, par exemple).
%
% La deuxieme ligne devrait etre un ligne descriptive de la fonction, et
% est en commentaire:
%
%   % MONPOLYFIT interpolation polynomiale par moindres carres
%
% Cette ligne est utilisee pour afficher les descriptions courtes des
% fonction lorque, par exemple la commande |lookfor| est utilisee. Les
% lignes suivantes, en commentaire, decrivent l'usage de la fonction. C'est
% le moment de mettre le |See also|, pour faire un lien vers d'autres
% fonctions Matlab. Ici on pourrait faire un lien vers la fonction
% |polyfit| de Matlab.
%
%   % See also POLYFIT
%
% (Les majuscules ne se verront pas avec la commande |help|).
%
% Ensuite viens le corps de la fonction. Tout d'abord on verifie que les
% deux vecteurs |x| et |y| sont des vecteurs de meme taille, et que |n+1| est
% inferieur ou egal a la taille de |y|. 
%
%   sx = size(x);
%   sy = size(y);
%   if min(sx)>1 % alors pas un vecteur 
%       error('monpolyfit:tailledonnees','X n''est pas un vecteur');
%   end
%   if min(sy)>1 % alors pas un vecteur 
%       error('monpolyfit:tailledonnees','Y n''est pas un vecteur');
%   end
%   if sx(2) ~= 1 % alors pas un vecteur colonne (~= : pas egal)
%       x = x.'; % on transpose (pas conjugue)
%       sx = size(x); % on recalcule la taille de x
%   end
%   if sy(2) ~= 1 % alors pas un vecteur colonne (~= : pas egal)
%       y = y.'; % on transpose (pas conjugue)
%       sy = size(y); % on recalcule la taille de y
%   end
%   if (sx(1)~=sy(1)) % verifions les tailles de X et Y
%       error('monpolyfit:tailledonnees','Les tailles de X et Y ne correspondent pas');
%   end
%   if (n+1)>length(x)
%       error('monpolyfit:tailledonnees','Le degre n du polynome a ajuster est trop eleve, il doit etre inferieur a la longueur de X');
%   end
%
% Maintenant on etablit notre equation normale pour les moindres carres. On
% genere un matrice de Vandermonde |A| a l'aide de |x|:
%   
%   A = vander(x);
%
% |A| est une matrice de taille |length(x)|. On ne retient de la matrice que les |n+1| dernieres lignes
% 
%   m = sx(1);
%   A = A(:,(m-n):m);
%
% On cherche les coefficients du polynome |p| qui approximent le mieux
% |A*p=y|. L'equation normale (en notation Matlab) est A'*A*p = A'*y. Si on
% fait une decomposition QR de la matrice |A|, on obtient l'equation R*p =
% Q'*y. On resoud ces equations a l'aide des fonctions |householder|,
% |Qadj| et |backsubs|.
%
%   [W,R]=householder(A);
%   Qy = Qadjb(W,y);
%   p = backsubs(R(1:(n+1),1:(n+1)),Qy(1:(n+1)));
%
% Le code est termine, pas besoin de terminer le fichier avec un |return|
% et un |end|. Le vecteur |p| est a ete defini et sera retourne par la
% fonction. Toutes les variables creees dans la fonction sont locales a
% cette fonction et ne peuvent etre vues a l'exterieur de celle-ci. Les
% variables sont detruites lors du retour de la fonction, a moins de 
% declarer ces variables |persistent|. Cela pourrait etre utile, par
% exemple si on voulait faire plusieurs interpolations de meme degre avec
% le meme vecteur |x|. Il serait alors plus efficace de ne pas recalculer
% |A| est sa decomposition a chaque appel de la fonction.
%
% Pour resoudre un probleme, il est utile d'ecrire dans un fichier la liste
% des commandes a faire. Ce fichier n'est qu'une liste d'expressions a
% evaluer par Matlab. Tout fichier avec l'extension |.m| dans le chemin de
% Matlab peut etre appele de l'invite de commande. Un fichier |.m| est un
% fichier script si il ne commence pas par |function|. Les fichiers script
% sont executes ligne par ligne et les variables sont dans le workspace.
%
% On cree un fichier script appele |fit_polynome.m|
%
%   edit fit_polynome.m
%
% La premiere ligne devrait etre une description du fichier en commentaire:
%
%   %% fit_polynome: script pour ajuster un polynome a des donnees
%
% Ensuite on definit les vecteur |x| et |y|.
%
%   x = 10*rand(34,1);      % 34 points aleatoires dans l'intervalle [0, 10];
%   y = 10*rand(34,1) + x;  % 34 points aleatoires correles avec x
%   n = 3;                  % interpolation avec un polynome de degre 3
%   p = monpolyfit(x,y,n);  % on trouve le polynome qui s'ajuste le mieux
%
% En executant le script |fit_polynome|, on obtient le polynome |p| demande. Toutes les variables definies dans le
% script se retrouvent dans le workspace.

%% Graphiques
% Matlab possede plusieurs commandes pour visualiser des donnees. Pour
% tracer des donnes en 2D, on utilise la commande |plot|. Rajoutons dans le
% fichier script fit_polynome les commandes pour visualiser les nuages de
% points |(x,y)| et le graphe du polynome |p|. 
%
%   %% Visualisation 2D des donnees
%   figure(1);              % cree ou active la figure 1
%   clf;                    % clf efface la figure
%   plot(x,y,'s')           % trace les points avec des carres ('s' pour square)
%   hold on                 % les prochains appels a plot n'effaceront pas les autres
%   r = max(x)-min(x);      % l'etendue de x
%   a = min(x) - 0.1*r;     % limite inferieure sur x pour tracer le polyome
%   b = max(x) + 0.1*r;     % limite superieure sur x pour tracer le polyome
%   xp = linspace(a,b,100); % 100 points equirepartis sur [a,b]
%   plot(xp,polyval(p,xp))  % trace le graphe de p aux valeurs xp
%   xlabel('x')             % etiquette abscisse
%   ylabel('y')             % etiquette ordonnees
%   legend('donnees','polynome',2) % legende, les etiquettes sont dans
%                                  % l'ordre d'aparition des plots
%
% Pour voir le resultat, on execute le script,

fit_polynome

%% Resoudre un systeme d'equations differentielles ordinaires
% Matlab possede plusieurs solveurs EDO, dont une suite |edo*|. Le solveur
% multi-usage a utiliser par default est |ode45| (prononce ode quatre cinq)
% un solveur Runge-Kutta a pas de temps adaptatif d'ordre 4 et 5. Les
% autres solveurs utiles sont |ode23|, un solveur d'ordre  |ode15s|. Tous
% les solveurs ont la meme syntaxe:
%
%   [t,x] = ode45(@monedo,[t0 tfinal],CI);
%   sol   = ode45(@monedo,[t0 tfinal],CI);
%
% La premiere syntaxe renvoie la solution |x| au temps |t|. La deuxieme
% syntaxe renvoie une structure |sol| qui contient la solution ainsi que
% des informations sur le solveur utilise. La fonction |monedo| est une
% fonction Matlab qui definit le systeme d'EDO. Le systeme est integre sur
% l'intervalle |[t0 tfinal]|, avec conditions initiales |CI|. L'equation
% differentielle est de la forme
% 
% $$\frac{dx}{dt}(t) = F(t,x).$$
%
% avec $x(t)$ un vecteur de taille $n$. La fonction EDO doit suivre une
% structure bien definie, avec une premiere ligne
%
%   function dxdt = monedo(t,x)
%
% La fonction prend en argument le temps |t| et l'etat du systeme |x|. La
% variable |x| est un vecteur colonne de taille |n|. La fonction retourne
% la derivee de |x| au temps |t| dans |dxdt| un vecteur colonne de taille
% |n|. Les conditions initiales sont un vecteur de la meme dimension que
% |x|. Comme premier exemple, on prend le systeme de Lotka-Volterra. La
% fonction <lotkavolterra.m> specifie le systeme suivant:
%
% $$\frac{dx}{dt} =  x (a - b y),$$ 
% $$\frac{dy}{dt} = -y (c - d x).$$
%
% Le systeme est resolu de la facon suivante:

CI = [30; 2];
tspan = [0 20];
sol = ode23(@lotkavolterra,tspan,CI);

%%
% On trace la solution a l'aide de |deval|, qui permet d'evaluer la
% solution |sol| a n'importe quel temps. 

xint = linspace(0,20,200); 
yint = deval(sol,xint);
figure(2); clf;
plot(xint,yint)
xlabel('temps'); ylabel('population'); 

%%
% Dans cet exemple, il n'est pas possible de facilement changer les
% parametres du systeme, qui sont code dans le fichier <lotkavolterra.m>.
% De plus, pour chaque resolution du systeme, il faut definir les
% conditions initiales et l'intervalle d'integration. Un fichier script
% peut aider, mais ne permet pas d'acceder aux parametres, qui sont
% variables locales a |lotkavolterra|. Un moyen de s'en sortir est
% d'utiliser une fonction imbriquee. Une fonction imbriquee est une
% sous-fonction, mais qui partage certaines variables avec la fonction
% principale dans laquelle elle se trouve.
%
% A lieu de definir un script pour resoudre le systeme d'EDO, on definit
% une fonction run_lotka.m 
%
%   function sol = run_lotka()
%   % RUN_LOTKA simulation du systeme de Lotka-Volterra
% 
%   % parametres d'integration
%   CI = [30; 2];
%   tspan = [0 20];
% 
%   % parametres du systeme
%   pars = [0.1 1 2 ];      % Les variables a, b, c, d sont partagees
%   b = 0.01;   % avec la fonction imbriquee lotka
%   c = 1;
%   d = 0.02;
% 
%   % resolution du systeme pour 3 jeux de parametres differents
%   for i = 1:3
%       a = pars(i); % on met a jour la valeur de a
%       sol{i} = ode23(@lotka,tspan,CI);
%   end
% 
%   % plot des 3 solutions
%   nx = 200;
%   xint = linspace(tspan(1),tspan(2),nx);
%   figure(2); clf;
%   newplot;            % met en place un nouveau plot
%   hold on             % qui ne s'effacera pas a chaque appel de plot
%   cols = 'kbr';       % couleurs k: noir, b: bleu, r:rouge
%   for i = 1:3
%       yint = deval(sol{i},xint);
%       plot(xint,yint,'Color',cols(i))
%   end
%   xlabel('temps'); ylabel('population');
% 
%       function dxdt = lotka(~,xx)  % fonction imbriquee
%       % LOTKA systeme d'EDO
%         
%           x = xx(1); y = xx(2);
%         
%           dxdt = [ x * (a - b * y)
%               -y * (c - d * x)];
%       end                         % on termine d'abord la fonction imbriquee
% 
%   end                             % on termine ensuite la fonction principale 
%
% De cette facon, on peut resoudre le systeme pour n'importe quelle
% condition initiale, intervalle de temps et parametres du systeme. Le
% fichier <run_lotka.m> peut etre lance:

run_lotka;

%% Codes
% Les fichiers matlab peuvent etre telecharges ici:
% <householder.m householder.m>,
% <Qadjb.m Qadjb.m>,
% <backsubs.m backsubs.m>,
% <lotkavolterra.m lotkavolterra.m>,
% <run_lotka.m run_lotka.m>