function [y] = gauss_dissym(A, x)
    // generates a dissymetrical gaussian peak
    // inputs: A(1): peak position
    //     A(2): height of the peak
    //     A(3): width on the left side
    //     A(4): width on the right side
    //     x: vector of numbers
    // oututs: y: vector of numbers
    index = (x < A(1)); // vector of %t at the left, %f at the right
    y = zeros(x); // initialisation
    y(index) = A(2)*exp(-(x(index) - A(1)).^2/A(3)); // left profile
    y(~index) = A(2)*exp(-(x(~index) - A(1)).^2/A(4)); // right profile
endfunction

function [B] = approximative_linearisation(f, A, x, deltaA)
    // computes an approximation of the partial derivatives of f
    // with respect to its parameters A(i)
    // f: function of x which parameters are described by the vector A
    // as following: f = (A,x)
    // A: vector of parameters (numbers)
    // x: vector of numbers
    sizex = size(x, '*');
    sizeA = size(A, '*');
    B = zeros(sizex, sizeA);
    for i = 1:sizeA // only the parameter A(i) is modified
        Aleft = [A(1:(i-1)), A(i) - deltaA, A(i+1:$)];
        Aright = [A(1:(i-1)), A(i) + deltaA, A(i+1:$)];
        B(:,i) = (f(Aright, x) - f(Aleft, x))/2/deltaA; // df/dA(i)
    end
endfunction

// **********
// Constants and initialisation
// **********

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

// parameters of the noisy curve
paramgauss(1) = 0; // peak position
paramgauss(2) = 10; // height of the peak
paramgauss(3) = 1; // width on the left side of the peak
paramgauss(4) = 0.3; // width on the right side of the peak
var = 0.5; // variance of the normal law (noise)
nbpts = 200 // number of points
halfwidth = 3*max(paramgauss(3), paramgauss(4)); // to determine the x range
step = 2*halfwidth/nbpts;

// **********
// Functions
// **********

exec('common_functions.sce', -1)

// **********
// Main program
// **********

// Data generation

for i = 1:nbpts
    x = step*i - halfwidth;
    Xinit(i) = x;
    noise(i) = var*rand(1, 'normal');
end

Ybasis = gauss_dissym(paramgauss, Xinit);
Yinit = Ybasis + noise;

// Data saving

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

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

// **********
// Constants and initialisation
// **********

clear;
clf;

chdir('mypath/')

// Newton-Raphson parameters
precision = 1e-7; // stop condition
itermax = 30; // idem

// Precision of approximated derivation
epsilon = 1e-6;

// **********
// Functions
// **********

exec('common_functions.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: set of parameters of the model optimised by regression (vector)
    // R: list of the quality factors for the regression (vector)
    // X: explanatory variable (vector)
    // Yexp : response variable, measured values (vector)
    // A0 : paramètres d'initialisation du modèle (vector)
    // epsilon : stop value (scalar)
    k = 1; // initial damping factor, <=1,
    // avoids divergence 
    n = size(X, '*');
    e0 = sqrt(sum(Yexp.^2)); // normalisation of the quality factor
    Ycal = f(A0, X); // nitial model
    R(1) = res(Yexp, Ycal)/e0; // initial quality factor
    disp('i = 1 ; R = '+string(R(1))) // display of initial param
    i = 1;
    B = A0;
    flag = %t
    while (i < imax) & flag // tests the global convergence
        i = i+1;
        deltay = Yexp - Ycal;
        J = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        tJ = J'; // transposed
        deltap0 = inv((tJ*J))*tJ*deltay;
        flag2 = %t // for the 1st execution
        while flag2 & (k>0.1)
            deltap = k*deltap0;
            Bnew = B + deltap';
            Ycal = f(Bnew, X);
            R(i) = res(Yexp, Ycal)/e0;
            flag2 = (R(i) >= R(i-1)) // true if it diverges
            if flag2 then k = k*0.75; // increase the damping on divergence
            else k0 = k; // to display the value
                k = (1 + k)/2; // reduces the damping on convergence
            end
        end
        B = Bnew;
        flag = abs(R(i-1) - R(i)) > epsilon;
        disp('i = '+string(i)+' ; R = '+string(R(i)))
    end
    A = B;
endfunction

// **********
// Main program
// **********

// Data reading
data = read('noisy_asym_gauss_peak.txt',-1,2);

// Data characteristics
Xdef = data(:,1);
Ydef = data(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];

// Regression
tic();
[Aopt, Rnr] =...
    gaussnewton(gauss_dissym, Xdef, Ydef,...
    Ainit, itermax, precision)
t = toc();

// Calculated curve

Yopt = gauss_dissym(Aopt, Xdef);

// Display

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 = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        tJ = J'; // transposed
        deltap0 = inv((tJ*J))*tJ*deltay;

        J = approximative_linearisation(f, B, X, epsilon); // jacobian matrix
        deltap0 = pinv(J)*deltay;

// **********
// Constants and initialisation
// **********

clear;
clf;

chdir('mypath/')

// **********
// Functions
// **********

exec('common_functions.sce', -1)

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

// **********
// Main program
// **********

// Data reading
data = read('noisy_asym_gauss_peak.txt',-1,2);

// Data characteristics
Xdef = data(:,1);
Ydef = data(:,2);
Ainit = [-0.03, 9.7, 8*((0.84 - 0.03)/2.35)^2, 8*((0.45 + 0.03)/2.35)^2];

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

// Calculated curve

Yopt = gauss_dissym(Aopt, Xdef);

// Display

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

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