function [y] = gauss_dissym(A, x)
    // génère un pic gaussien dissymétrique
    // entrées : A(1) : position du pic
    //     A(2) : hauteur de la courbe
    //     A(3) : largeur de la courbe à gauche
    //     A(4) : largeur de la courbe à droite
    //     x : vecteur de nombres
    // sorties : y : vecteur de nombres
    indice = (x < A(1)); // vecteur de %t à gauche, %f à droite
    y = zeros(x); // initialisation
    y(indice) = A(2)*exp(-(x(indice) - A(1)).^2/A(3)); // profil gauche
    y(~indice) = A(2)*exp(-(x(~indice) - A(1)).^2/A(4)); // profil droit
endfunction

function [B] = linearisation_approchee(f, A, x, deltaA)
    // calcule une valeur approchée de sdérivées partielles de f
    // par rapport à ses paramètres A(i)
    // f : fonction de x acceptant un vecteur de paramètres A
    // sous la forme f = (A,x)
    // A : vecteur de paramètres (nombres)
    // x : vecteur de nombres
    taillex = size(x,'*');
    tailleA = size(A, '*');
    B = zeros(taillex, tailleA);
    for i = 1:tailleA // on ne modifie que le paramètre A(i)
        Agauche = [A(1:(i-1)), A(i) - deltaA, A(i+1:$)];
        Adroite = [A(1:(i-1)), A(i) + deltaA, A(i+1:$)];
        B(:,i) = (f(Adroite, x) - f(Agauche, x))/2/deltaA; // df/dA(i)
    end
endfunction

// **********
// Constantes et initialisation
// **********

clear;
clf;
chdir('monchemin/');

// paramètres de la courbe bruitée
parmgauss(1) = 0; // position du pic
paramgauss(2) = 10; // hauteur de la courbe gaussienne
paramgauss(3) = 1; // largeur de la courbe gaussienne à gauche
paramgauss(4) = 0.3; // largeur de la courbe gaussienne à droite
var = 0.5; // variance de la loi normale
nbpts = 200 // nombre de points
demielargeur = 3*max(paramgauss(3), paramgauss(4)); // pour intervalle x
pas = 2*demielargeur/nbpts;

// **********
// Fonctions
// **********

exec('fonctions_communes.sce', -1)

// **********
// Programme principal
// **********

// Génération des données

for i = 1:nbpts
    x = pas*i - demielargeur;
    Xinit(i) = x;
    bruit(i) = var*rand(1, 'normal');
end

Ybase = gauss_dissym(paramgauss, Xinit);

Yinit = Ybase + bruit;

// Enregistrement des données

// plot(Xinit, Ybase, "r")
// plot(Xinit, Yinit, "b")

write ('pic_gauss_dissym_bruite.txt', [Xinit, Yinit])

// **********
// Constantes et initialisation
// **********

clear;
clf;

chdir('monchemin/')

// Paramètres de Newton-Raphson
precision = 1e-7; // condition d'arrêt
itermax = 30; // idem
 
// Précision de la linéarisation approchée
epsilon = 1e-6;
 
// **********
// Fonctions
// **********
 
exec('fonctions_communes.sce', -1)
 
function [e] = res(Yexp, Ycal)
    e = sqrt(sum((Yexp-Ycal).^2));
endfunction
 
function [A, R] = gaussnewton(f, X, Yexp, A0, imax, epsilon)
    // A : jeu de paramètres optimisé par régression (vecteur)
    // R : liste des facteurs de qualité de la régression
    // pour chaque étape (vecteur)
    // X : variable explicative (vecteur)
    // Yexp : variable expliquée, valeurs mesurées (vecteur)
    // A0 : paramètres d'initialisation du modèle (vecteur)
    // epsilon : valeur d'arrêt (scalaire)
    k = 1; // facteur d'amortissement initial, <=1,
    // évite la divergence 
    n = size(X,'*');
    e0 = sqrt(sum(Yexp.^2)); // normalisation du facteur de qualité
    Ycal = f(A0, X); // modèle initial
    R(1) = res(Yexp, Ycal)/e0; // facteur de qualité initial
    disp('i = 1 ; k = 1 ; R = '+string(R(1))) // affichage param initiaux
    i = 1;
    B = A0;
    drapeau = %t;
    while (i < imax) & drapeau // teste la convergence globale
        i = i+1;
        deltay = Yexp - Ycal;
        J = linearisation_approchee(f, B, X, epsilon); // matrice jacobienne
        tJ = J'; // transposée
        drapeau2 = %t // pour une 1re exécution
        deltap0 = inv((tJ*J))*tJ*deltay;
        while drapeau2 & (k>0.1) // teste la divergence sur 1 étape
            deltap = k*deltap0;
            Bnouveau = B + deltap';
            Ycal = f(Bnouveau, X);
            R(i) = res(Yexp, Ycal)/e0;
            drapeau2 = (R(i) >= R(i-1)) // vrai si diverge
            if drapeau2 then k = k*0.75; // atténue si diverge
            else k0 = k; // pour affichage de la valeur
                k = (1 + k)/2; // réduit l'atténuation si converge
            end
        end
        B = Bnouveau;
        drapeau = abs(R(i-1) - R(i)) > epsilon
        disp('i = '+string(i)+' ; k = '+string(k0)+' ; R = '+string(R(i)))
    end
    A = B;
endfunction
 
// **********
// Programme principal
// **********
 
// lecture des données
donnees = read('pic_gauss_dissym_bruite.txt',-1,2);
 
// carcatéristiques des données
Xdef = donnees(:,1);
Ydef = donnees(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];
 
// Régression
tic();
[Aopt, Rnr] =...
    gaussnewton(gauss_dissym, Xdef, Ydef,...
    Ainit, itermax, precision)
t = toc();

// Courbe calculée
 
Yopt = gauss_dissym(Aopt, Xdef);
 
// Affichage
 
print(%io(2),Ainit)
print(%io(2),Aopt)
print(%io(2),t)
 
clf
 
subplot(2,1,1)
plot(Xdef, Ydef, "-b")
plot(Xdef, Yopt, "-r")
 
subplot(2,1,2)
plot(Rnr)

        J = linearisation_approchee(f, B, X, epsilon); // matrice jacobienne
        tJ = J'; // transposée
        deltap0 = inv((tJ*J))*tJ*deltay;

        J = linearisation_approchee(f, B, X, epsilon); // matrice jacobienne
        deltap0 = pinv(J)*deltay;

// **********
// Constantes et initialisation
// **********

clear;
clf;

chdir('monchemin/')

// **********
// Fonctions
// **********

exec('fonctions_communes.sce', -1)

function [e] = res(A, X, Yexp)
    Ycal = gauss_dissym(A, X);
    e = sqrt(sum((Yexp-Ycal).^2));
endfunction

// **********
// Programme principal
// **********

// lecture des données
donnees = read('pic_gauss_dissym_bruite.txt',-1,2);

// carcatéristiques des données
Xdef = donnees(:,1);
Ydef = donnees(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];

// Régression
tic();
[Rnr, Aopt] =...
    leastsq(list(res, Xdef, Ydef), Ainit);
t = toc();

// Courbe calculée

Yopt = gauss_dissym(Aopt, Xdef);

// Affichage

print(%io(2),Ainit)
print(%io(2),Aopt)
print(%io(2),t)

clf

plot(Xdef, Ydef, "-b")
plot(Xdef, Yopt, "-r")