function DemoRegression
%% Resolution de systeme d'equation simple avec optimisation
%
% Jean-Nicolas Ouellet Automne 2007, revisé 2009


%% DONNEES 
Pts = [[1:10]' [1:10]' ones(10,1)]; %genere des points
Pts(:,1:2) = Pts(:,1:2)+rand(10,2)*2; %bruit au coordonnees x,y des points


%% SOLUTION INITIALE
%un point (x,y) passant par une droite L=(a,b,c)' doit respecter l'equation:
%ax+by+c = 0
Eqn = Pts; %voir l'introduction a matlab sur le site web du cours pour plus de details
[U,S,V] = svd(Eqn);
L = V(:,end);
%on normalise L pour avoir un vecteur normal unitaire.
L = L./sqrt(L(1)^2+L(2)^2);


%% CHANGEMENT DE REPRESENTATION
% on ne peut pas minimiser les parametres d'un systeme d'equation homogenes
% avec lsqnonlin de matlab, on doit avoir une representation minimale
% Pour une droite, la representation rho,theta est minimale 
% ax+by+c=0  , 3 parametres (a,b,c) est equavalent a
% sin(theta)*x + cos(theta)*y + rho = 0,   2 parametres (rho,theta)
% lorsque le vecteur (a,b) est de norme unitaire. L'angle Theta est
% simplement atan2(b,a).  Si vous calculez theta avec atan, faites
% attention au signe de rho.

% passage en representation rho,theta
theta = atan2(L(2),L(1));
rho = L(3);

initialLineCoefficients = [theta rho];


Err = [Pts*L];
ErrRMS = abs(Err);
disp('--------- Erreur avant optimisation ---------')
disp(mean(ErrRMS))


%% OPTIMISATION
optset = optimset('Display','off',...
                  'Jacobian','off',...
                  'LevenbergMarquardt','off',...
                  'Diagnostics','off',...
                  'TolX',1e-4,...
                  'TolFun', 1e-4,...
                  'MaxFunEvals',1500);
optimalLineCoefficients = lsqnonlin(@FonctionErreurSurDroite,initialLineCoefficients,[],[],optset,Pts);

% passage en representation cartesienne
theta = optimalLineCoefficients(1);
rho = optimalLineCoefficients(2);
L = [cos(theta) sin(theta) rho]';


Err = [Pts*L];
ErrRMS = abs(Err);
disp('--------- Erreur apres optimisation ---------')
disp(mean(ErrRMS))








function Err = FonctionErreurSurDroite(L, Pts)
%calcul l'erreur entre la doite L et les points Pts

% passage en representation cartesienne
theta = L(1);
rho = L(2);
L = [cos(theta) sin(theta) rho]';


Err = [Pts*L];