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

clear;
clf;

chdir('monchemin/')

// Paramètres de Newton-Raphson
precision = 1e-7; // condition d'arrêt
itermax = 60; // 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;
        subplot(2,1,1)
        plot2d(X, Yexp, rect=[-3, -2, 3, 12])
        plot(X, Ycal, "-r")
        xstring(-2.8, -1.5, string(B))
        subplot(2,1,2)
        plot2d(R, rect=[1, 0, 10, 1])
        xstring(1.2, 0.1, 'α = '+string(k)+' ; R = '+string(R(i)))
        nom = 'picassym'+string(i)+'.gif';
        xs2gif(0,nom)
    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
        deltap0 = inv((tJ*J))*tJ*deltay;
        drapeau2 = %t // pour une 1re exécution
        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
        clf;
        subplot(2,1,1)
        plot2d(X, Yexp, rect=[-3, -2, 3, 12])
        plot(X, Ycal, "-r")
        xstring(-2.8, -1.5, string(B))
        subplot(2,1,2)
        plot2d(R, rect=[1, 0, 10, 1])
        xstring(1.2, 0.1, 'α = '+string(k0)+' ; R = '+string(R(i)))
        nom = 'picassym'+string(i)+'.gif';
        xs2gif(0,nom)
//        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];
Ainit = [1, 1, 1, 1];

// 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)