* Bugs 5824, 7732, 15344 fixed: datafit() upgraded + page overhauled 61/19661/18
Samuel GOUGEON [Wed, 10 Jan 2018 23:52:32 +0000 (00:52 +0100)]
  http://bugzilla.scilab.org/5824
  http://bugzilla.scilab.org/7732
  http://bugzilla.scilab.org/15344

  datafit.pdf: to be removed before merging

Change-Id: Ic994529d79b3f3192887c9c01db17db7caf992a5

20 files changed:
scilab/CHANGES.md
scilab/modules/helptools/data/pages/homepage-en_US.html
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/datafit_1.png
scilab/modules/helptools/images/datafit_2.png
scilab/modules/helptools/images/datafit_3.png
scilab/modules/helptools/images/datafit_4.png [new file with mode: 0644]
scilab/modules/optimization/demos/datafit/datafit.dem.gateway.sce
scilab/modules/optimization/demos/datafit/datafit.dem.sce
scilab/modules/optimization/demos/datafit/demo_datafit.sci [deleted file]
scilab/modules/optimization/demos/datafit/ellipse.dem.sce [new file with mode: 0644]
scilab/modules/optimization/demos/datafit/slopingGaussian.dem.sce [new file with mode: 0644]
scilab/modules/optimization/demos/optimization.dem.gateway.sce
scilab/modules/optimization/help/en_US/nonlinearleastsquares/datafit.xml
scilab/modules/optimization/macros/datafit.sci
scilab/modules/optimization/tests/nonreg_tests/bug_15117.tst
scilab/modules/optimization/tests/nonreg_tests/bug_2330.dia.ref [deleted file]
scilab/modules/optimization/tests/nonreg_tests/bug_2330.tst
scilab/modules/optimization/tests/nonreg_tests/bug_244.dia.ref [deleted file]
scilab/modules/optimization/tests/nonreg_tests/bug_244.tst

index 8f55c22..541ee1c 100644 (file)
@@ -171,17 +171,18 @@ Feature changes and additions
   - The probing direction can now be numeric: 1, 2, ..
   - Returned indices can now be formatted with the new option `indType`.
   - There were no unit tests. More than 100 tests are added.
-
+* `datafit` is now able to fit weighted data. It now supports any gap function vectorized for Data points, and so is much faster. It now accepts any matrix of parameters (not necessarily a colum). It now returns the average Mode-to-data distance, and the termination status for the quasi-Newton algo.
 
 Help pages:
 -----------
 
-* overhauled / rewritten: `bitget`, `edit`, `factorial`, `vectorfind`
+* overhauled / rewritten: `bitget`, `edit`, `factorial`, `vectorfind`, `datafit`
 * fixed / improved:  `bench_run` `M_SWITCH`, `comet`, `comet3d`
 * Rewritten: `weekday`
 * Translations added:
   - (ru): `weekday`
 
+
 User Interface improvements:
 ----------------------------
 
@@ -244,10 +245,12 @@ Bug Fixes
 
 ### Bugs fixed in 6.1.0:
 * [#2694](https://bugzilla.scilab.org/2694): `bitget` did not accept positive integers of types int8, int16 or int32.
+* [#5824](https://bugzilla.scilab.org/5824): The `datafit` algorithm was not documented.
 * [#6070](https://bugzilla.scilab.org/6070): How to make multiscaled plots was not documented.
 * [#7562](https://bugzilla.scilab.org/7562): `factorial` could use a huge memory amount even for a scalar argument.
 * [#7657](https://bugzilla.scilab.org/7657): `lstsize` was a duplicate of `size` and should be removed.
 * [#7724](https://bugzilla.scilab.org/7724): When a figure is created in .auto_resize="on" mode, its .axes_size sets its .figure_size accordingly, not the reverse. But this was not documented.
+* [#7732](https://bugzilla.scilab.org/7732): The `datafit` help page needed to be fixed and overhauled.
 * [#7765](https://bugzilla.scilab.org/7765): `champ1` is useless. `champ().colored` is available for a long time.
 * [#7948](https://bugzilla.scilab.org/7948): `gsort` could not perform multilevel sorting, and could not sort complex numbers completely.
 * [#7967](https://bugzilla.scilab.org/7967): The tricky size `[ny,nx]` of `meshgrid(x,y)` results and usages with graphics was not enough documented.
index ff4e55f..271a26a 100644 (file)
                     <li>There were no unit tests. More than 100 tests are added.</li>
                 </ul>
             </li>
+            <li><code>datafit</code> is now able to fit weighted data. It now supports any gap function vectorized for Data points, and so is much faster. It now accepts any matrix of parameters (not necessarily a colum). It now returns the average Mode-to-data distance, and the termination status for the quasi-Newton algo.</li>
         </ul>
         <h2 class="title">Help pages</h2>
         <ul>
index fbc93dc..9bc17db 100644 (file)
@@ -67,6 +67,7 @@ Graphics_29.png=87ec396c51a7cc6aafb148e1aa5f9256
 Graphics_3.png=9d4c072a9f27f8f0e07f7d098ee533ec
 Graphics_30.png=67f6ae9c0d36df93eb50e514746aac13
 Graphics_31.png=26826b399f567100ed0118a54c896b5a
+Graphics_32.png=26826b399f567100ed0118a54c896b5a
 Graphics_36.png=76950c3a78d543018ca046b7e5709d88
 Graphics_4.png=ca17ee46fae89a343407d090111051b5
 Graphics_5.png=6e915789427344f2111be90b463014b3
@@ -625,9 +626,10 @@ damp_fr_FR_1.png=e81873a345f84e42ce2c4f4779d3421d
 damp_ja_JP_1.png=e81873a345f84e42ce2c4f4779d3421d
 damp_pt_BR_1.png=e81873a345f84e42ce2c4f4779d3421d
 damp_ru_RU_1.png=e81873a345f84e42ce2c4f4779d3421d
-datafit_1.png=cd52087997c788d3b0942584d38f9e8a
-datafit_2.png=0df0ac6cce3ab8e656471675e248b36d
-datafit_3.png=7b4a378cf730d4e181aa7a6f76448863
+datafit_1.png=5c2eaec5b6e7fc82764cf6ee163d1db7
+datafit_2.png=0cfb5856f4b4be19a2e4ea36b5ac97e9
+datafit_3.png=1cddd11ca3b1dbaaf4d77b95d4af88ca
+datafit_4.png=56e06b267e0f9035f6ade4f350fd3b8f
 datatipCreate_1.png=8ec8bfcd7f9dfd2a04b65f08309ee8fe
 datatipGetEntities_1.png=71ddeec8bd4accf61d6e6aac4d089514
 datatipGetStruct_1.png=b2a89937e8fee4f020fc6aa88de8d065
@@ -854,7 +856,7 @@ nanreglin_2.png=5a3c75b736030778ba4a1717051810cb
 ndgrid_1.png=881b8c6e4e1ccf21fce8004766f2df8d
 ndgrid_2.png=02c9c5b309ed4b25ee57c8b35e6a1fe1
 ndgrid_ru_RU_2.png=4ea2b3c6a3035633dbbdf9ffc2c7a676
-nicholschart_2.png=0f5b59c05f28aad6beae03208bc3fa07
+nicholschart_2.png=d24a24247f03f61ff1a2c6647064221e
 nicholschart_en_US_1.png=364f970d6bb93b0f839013928c7c006a
 nicholschart_en_US_2.png=0f5b59c05f28aad6beae03208bc3fa07
 nicholschart_fr_FR_1.png=364f970d6bb93b0f839013928c7c006a
index a638f7f..80c9cb0 100644 (file)
Binary files a/scilab/modules/helptools/images/datafit_1.png and b/scilab/modules/helptools/images/datafit_1.png differ
index cad4822..e612c35 100644 (file)
Binary files a/scilab/modules/helptools/images/datafit_2.png and b/scilab/modules/helptools/images/datafit_2.png differ
index 71691fb..03479a7 100644 (file)
Binary files a/scilab/modules/helptools/images/datafit_3.png and b/scilab/modules/helptools/images/datafit_3.png differ
diff --git a/scilab/modules/helptools/images/datafit_4.png b/scilab/modules/helptools/images/datafit_4.png
new file mode 100644 (file)
index 0000000..22474c1
Binary files /dev/null and b/scilab/modules/helptools/images/datafit_4.png differ
index b9fc8db..1f4b1c3 100644 (file)
@@ -5,5 +5,9 @@
 // This file is released under the 3-clause BSD license. See COPYING-BSD.
 
 
-subdemolist = [_("Non linear data fitting"), "datafit.dem.sce"];
+subdemolist = [
+    _("Parabolic model (3 params)"), "datafit.dem.sce"
+    _("Sloping gaussian model (5 params)"), "slopingGaussian.dem.sce"
+    _("Tilted ellipse (5 params)"), "ellipse.dem.sce"
+    ];
 subdemolist(:,2) = SCI + "/modules/optimization/demos/datafit/" + subdemolist(:,2);
index 18f5969..76441d7 100644 (file)
@@ -4,6 +4,51 @@
 //
 // This file is released under the 3-clause BSD license. See COPYING-BSD.
 
-exec("SCI/modules/optimization/demos/datafit/demo_datafit.sci");
+function demo_datafit()
+
+    function y = FF(x)
+        // parametric function model
+        y = a*(x-b)+c*x.*x;
+    endfunction
+
+    function e = G(p, z)
+        // datafit external computes the error
+        a = p(1),
+        b = p(2),
+        c = p(3),
+        y = z(1),
+        x = z(2),
+        e = y - FF(x)
+    endfunction
+
+    // create the experimental data
+    X = [];
+    Y = [];
+    a = 34;
+    b = 12;
+    c = 14;
+    for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
+    Z = [Y; X];
+    //show the data points
+    my_handle = scf(100001);
+    clf(my_handle, "reset");
+    demo_viewCode(SCI + "/modules/optimization/demos/datafit/datafit.dem.sce");
+    plot(X, Y, "+");
+    f = gcf();
+    l=legend(_("Experimental data"),2);
+    sleep(500);
+    // solve the non linear data fitting
+    [p,err] = datafit(G,Z,[3;5;10]);
+    if is_handle_valid(f) then // If the window is still open after the sleep
+        // show the fitting curve
+        drawlater()
+        plot(X,FF(X), "r")
+        delete(l);
+        l = legend([_("Experimental data"); _("Fitting function")],2);
+        drawnow()
+    end
+
+endfunction
+
 demo_datafit();
 clear demo_datafit;
diff --git a/scilab/modules/optimization/demos/datafit/demo_datafit.sci b/scilab/modules/optimization/demos/datafit/demo_datafit.sci
deleted file mode 100644 (file)
index f705206..0000000
+++ /dev/null
@@ -1,54 +0,0 @@
-//
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) ????-2008 - INRIA
-// Copyright (C) 2010 - DIGITEO - Allan CORNET
-//
-// This file is distributed under the same license as the Scilab package.
-//
-
-function demo_datafit()
-
-    function y = FF(x)
-        // parametric function model
-        y = a*(x-b)+c*x.*x;
-    endfunction
-
-    function e = G(p, z)
-        // datafit external computes the error
-        a = p(1),
-        b = p(2),
-        c = p(3),
-        y = z(1),
-        x = z(2),
-        e = y - FF(x)
-    endfunction
-
-    // create the experimental data
-    X = [];
-    Y = [];
-    a = 34;
-    b = 12;
-    c = 14;
-    for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
-    Z = [Y; X];
-    //show the data points
-    my_handle = scf(100001);
-    clf(my_handle, "reset");
-    demo_viewCode(SCI + "/modules/optimization/demos/datafit/demo_datafit.sci");
-    plot(X, Y, "+");
-    f = gcf();
-    l=legend(_("Experimental data"),2);
-    sleep(500);
-    // solve the non linear data fitting
-    [p,err] = datafit(G,Z,[3;5;10]);
-    if is_handle_valid(f) then // If the window is still open after the sleep
-        // show the fitting curve
-        drawlater()
-        plot(X,FF(X), "r")
-        delete(l);
-        l = legend([_("Experimental data"); _("Fitting function")],2);
-        drawnow()
-    end
-
-endfunction
-
diff --git a/scilab/modules/optimization/demos/datafit/ellipse.dem.sce b/scilab/modules/optimization/demos/datafit/ellipse.dem.sce
new file mode 100644 (file)
index 0000000..6d8537b
--- /dev/null
@@ -0,0 +1,107 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2018 - Samuel GOUGEON
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+function demo_datafit()
+    function g = Gap(p, Data)
+        // p: a, b, center xc, yc, tilt °
+        C = cosd(p(5));
+        S = -sind(p(5));
+        x = Data(1,:) - p(3);
+        y = Data(2,:) - p(4);
+        g = ((x*C + y*S )/p(1)).^2 + ((-x*S + y*C)/p(2)).^2 - 1;
+    endfunction
+
+    // Generating the data
+    // -------------------
+    // Actual parameters :
+    //    a   b   xc  yc  tilt
+    a = grand(1,1,"unf",1,10);
+    b = grand(1,1,"unf",a/6, a);
+    C = grand(1,2,"unf",-10,10);
+    tilt = grand(1,1,"unf", -85,85);
+    Pa = [a b C(1) C(2) tilt];
+    Np = 30;            // Number of data points
+    Theta_min = grand(1,1,"unf",-180,-90)
+    Theta = grand(1,Np, "unf",Theta_min, Theta_min+270);
+    // untilted centered noisy arc
+    x = Pa(1)*(cosd(Theta) + grand(Theta, "unf",-0.07, 0.07));
+    y = Pa(2)*(sind(Theta) + grand(Theta, "unf",-0.07, 0.07));
+    // Tilting and off-centering the arc
+    A = Pa(5);
+    C = cosd(A); S = sind(A);
+    xe =  C*x + y*S + Pa(3);
+    ye = -S*x + y*C + Pa(4);
+    Data = [xe ; ye];
+
+    // Retrieving parameters from noisy data
+    // -------------------------------------
+    // Initial estimation
+    ab0 = (strange(xe)+strange(ye))/2;
+    xc0 = mean(xe);
+    yc0 = mean(ye);
+    A0 = -atand(reglin(xe,ye));
+    P0 = [ab0 ab0/2 xc0 yc0 A0];
+    // Parameters bounds
+    Pinf = [ 1     0.2     xc0-ab0/2, yc0-ab0/2 -89];
+    Psup = [1.2*ab0 1.2*ab0  xc0+ab0/2, yc0+ab0/2  89];// Fitting
+
+    // FITTING
+    [P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
+
+    // P(1) is the longest axis: postprocessing the orientation angle
+    if P(1)<P(2) then
+        P([1 2]) = P([2 1]);
+        if P(5)>0
+            P(5) = 90 - P(5);
+        else
+            P(5) = 90 + P(5);
+        end
+    end
+    if (P(5)*Pa(5))<0 & abs(-P(5)/Pa(5)-1)<0.3 & P(5)>5 then
+        P(5) = -P(5)
+    end
+
+    // COMPUTING THE FITTING CURVE
+    Theta = 0:2:360;
+    x = P(1) * cosd(Theta);
+    y = P(2) * sind(Theta);
+    A = P(5);
+    [x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
+
+    // COMPUTING THE TRUE ELLIPSE CURVE
+    xt = Pa(1) * cosd(Theta);
+    yt = Pa(2) * sind(Theta);
+    At = Pa(5);
+    [xt, yt] = (xt*cosd(At)+yt*sind(At)+Pa(3), -xt*sind(At)+yt*cosd(At)+Pa(4));
+
+    // DISPLAY
+    // -------
+    my_handle = scf(100001);
+    drawlater
+    clf(my_handle, "reset");
+    gcf().figure_name = _("datafit() random demo : ellipse fitting (unweighted)");
+    demo_viewCode(get_absolute_file_path("ellipse.dem.sce")+"ellipse.dem.sce");
+
+    plot(xt, yt, "g")   // True ellipse
+    plot(xe, ye, "+")   // Data
+    plot(x, y, "r")     // Fitting ellipse
+    isoview
+    legend([_("True elliptical law")
+            _("Data = true + noise dx|y = +/-7%")
+            _("Best ellipse fitting Data")], "in_lower_right");
+    // Parameters values:
+    Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
+    Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
+    xstring(P(3), P(4), ["" "    a           b        xc       yc     A°"
+    "True:" Patxt
+    "Fit:" Ptxt
+    "" "    dmin = "+ msprintf("%5.3f", dmin)]);
+    set(gce(),"text_box",[0 0], "text_box_mode", "centered");
+    title(_("datafit: Fitting on 30 random unweighted data points"),"fontsize",3)
+    drawnow
+endfunction
+
+demo_datafit();
+clear demo_datafit;
diff --git a/scilab/modules/optimization/demos/datafit/slopingGaussian.dem.sce b/scilab/modules/optimization/demos/datafit/slopingGaussian.dem.sce
new file mode 100644 (file)
index 0000000..fb1e035
--- /dev/null
@@ -0,0 +1,75 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2018 - Samuel GOUGEON
+//
+// This file is released under the 3-clause BSD license. See COPYING-BSD.
+
+function demo_datafit()
+    function Y = model(P,X)
+        Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
+    endfunction
+    // -------------------------------
+    // Gap function
+    function r = gap(P,XY)
+        r = XY(2,:) - model(P,XY(1,:))
+    endfunction
+    // -------------------------------
+    // True parameters
+    Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
+
+    // Generating data
+    // ---------------
+    X = 0:0.2:10;
+    Np = length(X);
+    Y0 = model(Pg, X);              // True law
+    noise = (rand(Y0)-.5)*Pg(3)/4;
+    Y = Y0 + noise                  // Noised data
+    Z = [X ; Y];
+
+    // Let's try to find the original parameters from the (noisy) data:
+    // ---------------------------------------------------------------
+    // Performing the fitting = parameters optimization
+    [v, k] = max(Y0);
+    P0 = [X(k) 1 v 1 1];
+    [P, fmin, s]  = datafit(gap, Z, P0);
+    Wz =  1./abs(noise);
+    [Pw, fminw, s] = datafit(gap, Z, Wz, P0);
+
+    // Computing fitting curves from found parameters
+    fY = model(P, X);
+    fYw = model(Pw,X);
+
+    // Display
+    // -------
+    my_handle = scf(100001);
+    clf(my_handle, "reset");
+    gcf().figure_name = "datafit() demo : Sloping gaussian";
+    demo_viewCode(get_absolute_file_path("slopingGaussian.dem.sce")+"slopingGaussian.dem.sce");
+    xsetech([0 0 1 0.75])
+
+    plot2d(X, Y0, 24)  // True underlaying law
+    plot2d(X, Y,  -1)  // Noisy data
+    plot2d(X, fY, 15)  // best unweighted fit
+    plot2d(X, fYw,18)  // best weighted fit
+    legend(["True law" "Noisy Data" ..
+           msprintf("Unweighted fit (%5.2f /%5.2f)",mean(abs(fY-Y0)), sqrt(fmin/Np)) ..
+           msprintf("Weighted fit (%5.2f /%5.2f)", mean(abs(fYw-Y0)), sqrt(fminw/sum(Wz)))],1)
+           // Average vertical Model-Data distance: True one / according to the residual cost
+    // Params: row#1: true params row#2: from unweighted fit row#3: from weighted fit
+    [xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
+    xnumb(xs, ys, [Pg ; P ; Pw])
+    xstring(2,1.5,"$\Large p_3\,.\,exp\left({-(x-p_1)^2}\over{2\,p_2^2} \right)+p_4.x+p_5$")
+
+    // Plotting residus:
+    xsetech([0 .75 1 0.25])
+    plot2d(X, model(P ,X)-Y0, 15)  // remaining gap with best unweighted fit
+    plot2d(X, model(Pw,X)-Y0, 18)  // remaining gap best weighted fit
+    ylabel("Fit - True law")
+    ax = gca();
+    ax.x_location = "top";
+    tmp = ax.x_ticks.labels;
+    gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+    xgrid(color("grey70"))
+endfunction
+
+demo_datafit();
+clear demo_datafit;
index 16ff9c1..8271c62 100644 (file)
@@ -10,10 +10,12 @@ function subdemolist = demo_gateway()
     _("Optimization and Simulation");  // lets gettext() harvesting it
     add_demo("Optimization and Simulation", demopath+"optimization.dem.gateway.sce");
 
-    subdemolist = [_("Non linear data fitting"), "datafit/datafit.dem.gateway.sce"; ..
-    _("Optimisation"),            "optim/optim.dem.gateway.sce"; ..
-    _("fminsearch"),              "neldermead/neldermead.dem.gateway.sce"; ..
-    _("ICSE"),              "icse/icse.dem.gateway.sce"];
+    subdemolist = [
+    _("datafit: Nonlinear fitting"), "datafit/datafit.dem.gateway.sce"
+    _("Optimisation"),  "optim/optim.dem.gateway.sce"
+    _("fminsearch"),    "neldermead/neldermead.dem.gateway.sce"
+    _("ICSE"),          "icse/icse.dem.gateway.sce"
+    ];
 
     subdemolist = [subdemolist; ..
     _("Genetic algorithms"), "genetic/genetic_algorithms.dem.gateway.sce"];
index 5f79d12..064a580 100644 (file)
@@ -3,6 +3,7 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
  * Copyright (C) 2011 - DIGITEO - Michael Baudin
+ * Copyright (C) 2018 - Samuel GOUGEON
  *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * along with this program.
  *
  -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="datafit" xml:lang="en">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+        xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml"
+        xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
+        xmlns:scilab="http://www.scilab.org" xml:id="datafit" xml:lang="en">
     <refnamediv>
         <refname>datafit</refname>
-        <refpurpose>Parameter identification based on measured data</refpurpose>
+        <refpurpose>Non linear (constrained) parametric fit of measured (weighted) data</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
-        <synopsis>[p,err]=datafit([iprint,] G [,DG],Z [,W],[contr],p0,[algo],[df0,[mem]],
-            [work],[stop],['in'])
+        <synopsis>[p, fmin, status] = datafit([iprint,] G [,DG], Data [,Wd] [,Wg][,'b',p_inf,p_sup], p0 [,algo][,stop])
         </synopsis>
     </refsynopsisdiv>
     <refsection>
             <varlistentry>
                 <term>iprint</term>
                 <listitem>
-                    <para>scalar argument used to set the trace mode.
-                        <literal>iprint=0</literal> nothing (except errors) is reported,
-                        <literal>iprint=1</literal> initial and final reports,
-                        <literal>iprint=2</literal> adds a report per iteration,
-                        <literal>iprint&gt;2</literal> add reports on linear search.
-                    </para>
-                    <para>
+                    <para>scalar used to set the trace mode:
+                        <table>
+                            <tr>
+                                <th>Value</th><th>Action</th>
+                            </tr>
+                            <tr>
+                                <td>0</td><td>nothing (except errors) is reported</td>
+                            </tr>
+                            <tr>
+                                <td>1</td><td>initial and final reports</td>
+                            </tr>
+                            <tr>
+                                <td>2</td><td>adds a report per iteration</td>
+                            </tr>
+                            <tr>
+                                <td>>2</td><td>adds reports on linear search</td>
+                            </tr>
+                        </table>
                         <warning>
-                            Warning: most of these reports are written on the Scilab standard output.
+                            Most of the reports are displayed in the Scilab console.
                         </warning>
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>G</term>
+                <term>p0</term>
                 <listitem>
-                    <para>function descriptor (e=G(p,z), e: ne x 1, p: np x 1, z: nz x
-                        1)
-                    </para>
+                    Matrix of initial guess of the model parameters.
+                  <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>DG</term>
+                <term>"b"</term>
                 <listitem>
-                    <para>partial of G wrt p function descriptor (optional; S=DG(p,z),
-                        S: ne x np)
-                    </para>
+                    keyword heading the <literal>p_inf, p_sup</literal> sequence, where
+                    <literal>p_inf</literal> and <literal>p_sup</literal> are real matrices
+                    with the <literal>p0</literal> shape and sizes.
+                    <literal>p_inf</literal> is the vector of lower bounds of respective
+                    <literal>p</literal> components. <literal>p_sup</literal> are upper bounds.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>Z</term>
+                <term>p</term>
                 <listitem>
-                    <para>matrix [z_1,z_2,...z_n] where z_i (nz x 1) is the ith
-                        measurement
-                    </para>
+                    Matrix of the best fitting model parameters, with the <varname>p0</varname> sizes.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>W</term>
+                <term>Data</term>
                 <listitem>
-                    <para>weighting matrix of size ne x ne (optional; default no
-                        ponderation)
-                    </para>
+                    (ndc x nd) matrix of experimental data points, where <literal>Data(:,i)</literal>
+                    are the #nzc coordinates of the ith measurement.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>contr</term>
+                <term>Wd</term>
                 <listitem>
-                    <para>
-                        <literal>'b',binf,bsup</literal> with
-                        <literal>binf</literal> and <literal>bsup</literal> real vectors
-                        with same dimension as <literal>p0</literal>.
-                        <literal>binf</literal> and <literal>bsup</literal> are lower and
-                        upper bounds on <literal>p</literal>.
-                    </para>
+                    Data weights (optional). Row vector with size(Data,2) elements.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>p0</term>
+                <term>G</term>
                 <listitem>
-                    <para>initial guess (size np x 1)</para>
+                    <para>
+                        Gap function: handle of a function in Scilab language or of a compiled
+                        external function. Computed gap(s) between the fitting model and data
+                        may be of the first order, not quadratic.
+                    </para>
+                    <para>
+                        <emphasis role="bold">
+                            If <literal>G()</literal>
+                        </emphasis>
+                        does not need
+                        any input parameters (like constants), only arguments, its syntax must be
+                        <literal>g = G(p, Data)</literal>.
+                    </para>
+                    <para>
+                        <emphasis role="bold">
+                            If <literal>G()</literal>
+                        </emphasis>
+                        needs
+                        some input parameters in addition to input arguments, its syntax must be
+                        <literal>g = G(p, Data, param1, param2, .., paramN)</literal>.
+                        Then, <literal>list(G, param1, param2, .., paramN)</literal> must be provided
+                        to <literal>datafit()</literal> instead of only <literal>G</literal>.
+                    </para>
+                    <para>
+                        The following definitions will be used:
+                        <itemizedlist>
+                            <listitem>
+                                <literal>ng = size(g,1)</literal>: number of fitting criteria (= of
+                                gap definitions).
+                            </listitem>
+                            <listitem>
+                                <literal>np = size(p,"*")</literal>: total number of fitting parameters.
+                            </listitem>
+                            <listitem>
+                                <literal>nd = size(Data,2)</literal>: number of data points.
+                            </listitem>
+                            <listitem>
+                                <literal>ndc = size(Data,1)</literal>: number of coordinates of each
+                                data point.
+                            </listitem>
+                        </itemizedlist>
+                    </para>
+                    <para>
+                        Then, the sizes of G arguments are the following ones:
+                        <table>
+                            <tr valign="top">
+                                <th>g</th>
+                                <td>: (ng x nd) matrix
+                                    <itemizedlist>
+                                        <listitem>
+                                            Each column is related to a given data point and retuns its
+                                            gaps to the fitting model.
+                                        </listitem>
+                                        <listitem>
+                                            Each row is about a gap definition. <literal>G</literal> may
+                                            compute several types of gaps. For instance, the vertical
+                                            gaps may be returned in <literal>g(1,:)</literal> ; the
+                                            horizontal gaps in <literal>g(2,:)</literal> ; the nearest
+                                            distances in <literal>g(3,:)</literal> ; etc.
+                                        </listitem>
+                                    </itemizedlist>
+                                </td>
+                            </tr>
+                            <tr valign="top">
+                                <th>p</th>
+                                <td>: fitting parameters to be optimized.
+                                    <para/>
+                                </td>
+                            </tr>
+                            <tr valign="top">
+                                <th>Data</th>
+                                <td>: (ndc x nd) matrix: Each column provides the
+                                    <literal>#ndc</literal> coordinates of a given data point.
+                                    If <literal>G</literal> computes gaps only for a single
+                                    data point, <literal>datafit()</literal> will automatically
+                                    loop over the set of <literal>Data</literal> columns.
+                                    This makes <literal>datafit()</literal> slower.
+                                </td>
+                            </tr>
+                        </table>
+                    </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>algo</term>
+                <term>DG</term>
                 <listitem>
-                    <para>
-                        <literal>'qn'</literal> or <literal>'gc'</literal> or
-                        <literal>'nd'</literal> . This string stands for quasi-Newton
-                        (default), conjugate gradient or non-differentiable respectively.
-                        Note that <literal>'nd'</literal> does not accept bounds on
-                        <literal>x</literal> ).
-                    </para>
+                    (optional) Handle of a function in Scilab language or of a compiled external
+                    function, computing the partial of <varname>G</varname> wrt <varname>p</varname>.
+                    Its syntax must be like <literal>S = DG(p, Data)</literal>. The expected sizes
+                    of the result are (ng x np x 1) (Data-unvectorized) or (ng x np x nd)
+                    (Data-vectorized).
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>df0</term>
+                <term>Wg</term>
                 <listitem>
-                    <para>
-                        real scalar. Guessed decreasing of <literal>f</literal> at
-                        first iteration. (<literal>df0=1</literal> is the default
-                        value).
-                    </para>
+                    Optional (ng x ng) matrix of weights of gaps definitions. Default = eye().
+                    <varname>G()</varname> defines as many types of gaps as required.
+                    Most often, a unique gap will be computed -- most often the vertical one --,
+                    but other ones can be computed (horizontal, nearest distance, etc).
+                    This matrix allows to weight the different types.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>mem :</term>
+                <term>algo</term>
                 <listitem>
-                    <para>integer, number of variables used to approximate the Hessian,
-                        (<literal>algo='gc' or 'nd'</literal>). Default value is around
-                        6.
-                    </para>
+                    Word <literal>'qn'</literal> (quasi-Newton, default), or
+                    <literal>'gc'</literal> (conjugate gradient), or
+                    <literal>'nd'</literal> (non-differentiable. Then <varname>p_inf</varname>
+                    and <varname>p_sup</varname> parameters bounds are not accepted).
+                    Selects the algorithm used to perform the fitting.
+                    See <function>optim()</function> for more details.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>stop</term>
                 <listitem>
-                    <para>sequence of optional parameters controlling the convergence of
-                        the algorithm. <literal> stop= 'ar',nap, [iter [,epsg [,epsf
-                            [,epsx]]]]
-                        </literal>
+                    <para>sequence of optional arguments controlling the convergence of
+                        the fitting algorithm:
                     </para>
-                    <variablelist>
-                        <varlistentry>
-                            <term>"ar"</term>
-                            <listitem>
-                                <para>reserved keyword for stopping rule selection defined as
-                                    follows:
-                                </para>
-                            </listitem>
-                        </varlistentry>
-                        <varlistentry>
-                            <term>nap</term>
-                            <listitem>
-                                <para>
-                                    maximum number of calls to <literal>fun</literal>
-                                    allowed.
-                                </para>
-                            </listitem>
-                        </varlistentry>
-                        <varlistentry>
-                            <term>iter</term>
-                            <listitem>
-                                <para>maximum number of iterations allowed.</para>
-                            </listitem>
-                        </varlistentry>
-                        <varlistentry>
-                            <term>epsg</term>
-                            <listitem>
-                                <para>threshold on gradient norm.</para>
-                            </listitem>
-                        </varlistentry>
-                        <varlistentry>
-                            <term>epsf</term>
-                            <listitem>
-                                <para>threshold controlling decreasing of
-                                    <literal>f</literal>
-                                </para>
-                            </listitem>
-                        </varlistentry>
-                        <varlistentry>
-                            <term>epsx</term>
-                            <listitem>
-                                <para>
-                                    threshold controlling variation of <literal>x</literal>.
-                                    This vector (possibly matrix) of same size as
-                                    <literal>x0</literal> can be used to scale
-                                    <literal>x</literal>.
-                                </para>
-                            </listitem>
-                        </varlistentry>
-                    </variablelist>
+                    <para>
+                        <literal>'ar',nap, [iter [,epsg [,epsf [,epsx]]]]</literal>, with
+                    </para>
+                    <table>
+                        <tr>
+                            <th>"ar"</th>
+                            <td>
+                                reserved keyword heading the list of stopping criteria as defined below.
+                            </td>
+                        </tr>
+                        <tr>
+                            <th>nap</th>
+                            <td>
+                                maximum number of calls to the cost f function allowed (default: 100).
+                            </td>
+                        </tr>
+                        <tr>
+                            <th>iter</th>
+                            <td>maximum number of iterations allowed (default: 100).</td>
+                        </tr>
+                        <tr>
+                            <th>epsg</th>
+                            <td>threshold on gradient norm (default: %eps).</td>
+                        </tr>
+                        <tr>
+                            <th>epsf</th>
+                            <td>
+                                threshold controlling the decreasing of <literal>f</literal>
+                                (defined in the description section). Default: 0.
+                            </td>
+                        </tr>
+                        <tr>
+                            <th>epsx</th>
+                            <td>
+                                threshold controlling the variation of <varname>p</varname>.
+                                This vector (possibly matrix) of same size as <varname>p0</varname>
+                                can be used to scale <varname>p</varname>. Default: 0.
+                            </td>
+                        </tr>
+                    </table>
                 </listitem>
             </varlistentry>
-            <varlistentry>
-                <term>"in"</term>
+            <!--            <varlistentry>
+                <term>fmin</term>
                 <listitem>
-                    <para>reserved keyword for initialization of parameters used when
-                        <literal>fun</literal> in given as a Fortran routine (see
-                        below).
+                    <para>
+                        <literal>f</literal> quadratic residu (minimal <literal>f</literal> value reached).
+                        <literal>sqrt(fmin/ng/nz)</literal> or
+                        <literal>sqrt(fmin/ng/sum(Wd))</literal> estimates the Model-to-Data average
+                        distance.
                     </para>
                 </listitem>
             </varlistentry>
+-->
             <varlistentry>
-                <term>p</term>
+                <term>dmin</term>
                 <listitem>
-                    <para>Column vector, optimal solution found</para>
+                    Average Model-to-data residual distance. It is equal to
+                    <literal>sqrt(fmin/ng/nd)</literal> or
+                    <literal>sqrt(fmin/ng/sum(Wd))</literal>, where fmin is the quadratic residue
+                    of the cost function <literal>f</literal> (see the description).
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>err</term>
+                <term>status</term>
                 <listitem>
-                    <para>scalar, least square error.</para>
+                  return flag = termination status (only for the 'qn' default algorithm.
+                  Set to %nan otherwise): <literal>9</literal> is successful. Possible values are:
+                  <table>
+                    <tr valign="top">
+                      <th>1</th><td>: Norm of projected gradient lower than...</td>
+                      <td>....</td>
+                      <th>6</th><td>: Optim stops: too small variations in gradient direction.</td>
+                    </tr>
+                    <tr valign="top">
+                      <th>2</th><td>: At last iteration, f decreases by less than...</td>
+                      <td></td>
+                      <th>7</th><td>: Stop during calculation of descent direction.</td>
+                    </tr>
+                    <tr valign="top">
+                      <th>3</th><td>: Optim stops: Too small variations for p.</td>
+                      <td></td>
+                      <th>8</th><td>: Stop during calculation of estimated hessian.</td>
+                    </tr>
+                    <tr valign="top">
+                      <th>4</th><td>: Optim stops: maximum number of calls to f is reached.</td>
+                      <td></td>
+                      <th>9</th><td>: Successful completion.</td>
+                    </tr>
+                    <tr valign="top">
+                      <th>5</th><td>: Optim stops: maximum number of iterations is reached.</td>
+                      <td></td>
+                      <th>10</th><td>: Linear search fails.</td>
+                    </tr>
+                  </table>
+                  <para/>
                 </listitem>
             </varlistentry>
         </variablelist>
+        <para>
+            Other <function>optim()</function> input arguments such that <varname>df0</varname>,
+            <varname>mem</varname>, or <varname>'in'</varname> may be used with
+            <literal>datafit()</literal>. They will be passed as is to <function>optim()</function>.
+        </para>
     </refsection>
     <refsection>
         <title>Description</title>
         <para>
-            <literal>datafit</literal> is used for fitting data to a model. For
-            a given function <literal>G(p,z)</literal>, this function finds the best
-            vector of parameters <literal>p</literal> for approximating
-            <literal>G(p,z_i)=0</literal> for a set of measurement vectors
-            <literal>z_i</literal>. Vector <literal>p</literal> is found by minimizing
-            <literal>G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)</literal>
+            <literal>datafit</literal> is used to fit a parametrized model to given data.
+        </para>
+        <para>
+            A function <literal>G(p, Data)</literal> must be defined to compute the gaps between the
+            <varname>Data</varname> points and a model whose parameters to be tuned are provided
+            through the matrix <varname>p</varname>.
+
+        </para>
+        <para>
+            <literal>datafit()</literal> checks whether <literal>G()</literal> is vectorized or not
+            over <varname>Data</varname>. If <literal>G()</literal> works only with a single
+            <literal>data</literal> point as a single column of coordinates, then <literal>datafit()</literal>
+            loops over <literal>data=Data(:,i)</literal> points. Otherwise, <literal>G(p, Data)</literal>
+            is called for the whole <varname>Data</varname> array of data coordinates, which is
+            way faster if <varname>G()</varname>'s code is truly vectorized.
+        </para>
+        <para>
+            Starting from initial values <varname>p0</varname> of the model parameters,
+            <literal>datafit()</literal> varies them in order to minimize the whole Model-to-Data
+            distance defined as
+            <screen>
+f = (G(p, Data(:,1))' * Wg * G(p, Data(:,1)))* Wd(1) + ..
+    (G(p, Data(:,2))' * Wg * G(p, Data(:,2)))* Wd(2) + ..
+    ...
+    (G(p, Data(:,nz))' * Wg * G(p,Data(:,nz)))* Wd(nz)
+</screen>
+        </para>
+        <para>
+            When <literal>Wg = eye()</literal> (default), this definition is equivalent to
+<screen>
+f = sum(G(p, Data(:,1)).^2) * Wd(1) + ..
+    sum(G(p, Data(:,2)).^2) * Wd(2) + ..
+    ...
+    sum(G(p, Data(:,nz)).^2) * Wd(nz)
+</screen>
+        </para>
+        <para>
+            If in addition <varname>G</varname> returns only one gap type, this even simplifies to
+            <screen>
+f = G(p, Data(:,1))^2 * Wd(1) + ..
+    G(p, Data(:,2))^2 * Wd(2) + ..
+    ...
+    G(p, Data(:,nz))^2 * Wd(nz)
+</screen>
+        </para>
+        <para>
+            <literal>datafit()</literal> calls <function>optim()</function> to minimize f.
+        </para>
+        <para>
+            Choosing appropriately <varname>p0</varname> may improve and faster
+            <literal>datafit()</literal>'s convergence to the best fit.
         </para>
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[
-//generate the data
-function y=FF(x,p)
-  y=p(1)*(x-p(2))+p(3)*x.*x
+        <table>
+            <tr>
+                <td>
+                    <para>
+                        <emphasis role="bold">Simple example: Polynomial fitting (parabola = 3 parameters)</emphasis>.
+                    </para>
+                    <programlisting role="example"><![CDATA[
+// Model (parabola)
+function y = FF(x,p)
+    y = (p(1)*x + p(2)).*x + p(3)
 endfunction
 
-X=[];
-Y=[];
-pg=[34;12;14] //parameter used to generate data
-for x=0:.1:3
-  Y=[Y,FF(x,pg)+100*(rand()-.5)];
-  X=[X,x];
-end
-Z=[Y;X];
+// Data creation
+pa = [2;-4;14] // parameters of the actual law, to generate data
+X = 0:0.1:3.5;
+Y0 = FF(X, pa);
+noise = 10*(rand(Y0)-.5);   // We add some noise
+Y = Y0 + noise
+Data = [X ; Y];
 
-//The criterion function
-function e=G(p,z),
-  y=z(1),x=z(2);
-  e=y-FF(x,p),
+// Gap function definition
+function e = G(p, Data)
+    e = Data(2,:) - FF(Data(1,:), p);  // vertical distance
 endfunction
 
-//Solve the problem
-p0=[3;5;10]
-[p,err]=datafit(G,Z,p0);
-
-scf(0);clf()
-plot2d(X,FF(X,pg),5) //the curve without noise
-plot2d(X,Y,-1)  // the noisy data
-plot2d(X,FF(X,p),12) //the solution
- ]]></programlisting>
-        <scilab:image>
-            function y=FF(x,p)
-            y=p(1)*(x-p(2))+p(3)*x.*x
-            endfunction
-
-            X=[];Y=[];
-            pg=[34;12;14]
-            for x=0:.1:3
-            Y=[Y,FF(x,pg)+100*(rand()-.5)];
-            X=[X,x];
-            end
-            Z=[Y;X];
-
-            function e=G(p,z),
-            y=z(1),x=z(2);
-            e=y-FF(x,p),
-            endfunction
-
-            p0=[3;5;10]
-            [p,err]=datafit(G,Z,p0);
-
-            scf(0);clf()
-            plot2d(X,FF(X,pg),5)
-            plot2d(X,Y,-1)
-            plot2d(X,FF(X,p),12)
-        </scilab:image>
-        <programlisting role="example"><![CDATA[
-//generate the data
-function y=FF(x,p)
-  y=p(1)*(x-p(2))+p(3)*x.*x
-endfunction
+// Performing the fitting
+// ----------------------
+// Without weighting data:
+p0 = [1;-1;10]
+[p, dmin] = datafit(G, Data, p0);
+// The uncertainty is identified to the noise value.
+// The weight is set inversely proportional to the uncertainty
+[pw, dmin] = datafit(G, Data, 1./abs(noise), p0);
 
-//the gradient of the criterion function
-function s=DG(p,z),
-  a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),
-  s=-[x-b,-a,x*x]
-endfunction
+// Display
+// -------
+scf(0);
+clf()
+xsetech([0 0 1 0.83])
+plot2d(X, Y0, 24)       // True underlaying law
+plot2d(X, Y, -1)        // Noisy data
+plot2d(X, FF(X, p), 15)  // best unweighted fit
+plot2d(X, FF(X, pw), 18) // best weighted fit
+legend(["True law" "Noisy data" "Unweighted fit" "Weighted fit"],2)
 
-function e=G(p,z),
-  y=z(1),x=z(2);
-  e=y-FF(x,p),
+xsetech([0 .75 1 0.25])
+plot2d(X, FF(X,p)-Y0, 15)  // remaining gap <= best unweighted fit
+plot2d(X, FF(X,pw)-Y0, 18) // remaining gap <= best weighted fit
+ylabel("Data fit - True law")
+ax = gca(); ax.x_location = "top";
+tmp = ax.x_ticks.labels;
+gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+xgrid(color("grey70"))
+         ]]></programlisting>
+                </td>
+                <td>
+                    <scilab:image>
+                        // Model (parabola)
+                        function y = FF(x,p)
+                        y = p(1)*x.*x + p(2)*x + p(3)
+                        endfunction
+
+                        // Data creation
+                        pa = [2;-4;14] // parameters of the actual law, used to generate data
+                        X = 0:0.1:3.5;
+                        Y0 = FF(X, pa);
+                        noise = 10*(rand(Y0)-.5);   // We add some noise
+                        Y = Y0 + noise
+                        Data = [X ; Y];
+
+                        // Gap function definition
+                        function e = G(p, Data)
+                        e = Data(2,:) - FF(Data(1,:), p);  // vertical distance
+                        endfunction
+
+                        // Performing the fitting
+                        // ----------------------
+                        // Without weighting data:
+                        p0 = [1;-1;10]
+                        [p, dmin] = datafit(G, Data, p0);
+
+                        // The uncertainty is identified to the noise value.
+                        // The weight is set inversely proportional to the uncertainty
+                        [pw, dmin] = datafit(G, Data, 1 ./abs(noise), p0);
+
+                        // Display
+                        // -------
+                        scf(0);
+                        clf()
+                        xsetech([0 0 1 0.83])
+                        plot2d(X, Y0, 24)       // True underlaying law
+                        plot2d(X, Y, -1)        // Noisy data
+                        plot2d(X, FF(X,p), 15)  // best unweighted fit
+                        plot2d(X, FF(X,pw), 18) // best weighted fit
+                        legend(["True law" "Noisy data" "Unweighted fit" "Weighted fit"],2)
+
+                        xsetech([0 .75 1 0.25])
+                        plot2d(X, FF(X,p)-Y0, 15)  // remaining gap with best unweighted fit
+                        plot2d(X, FF(X,pw)-Y0, 18) // remaining gap best weighted fit
+                        ylabel("Data fit - True law")
+                        ax = gca(); ax.x_location = "top";
+                        tmp = ax.x_ticks.labels;
+                        gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+                        xgrid(color("grey70"))
+                        // ===============
+                        gcf().axes_size = [500 500];
+                        gcf().children(1).margins(1) = 0.07;
+                        gcf().children(2).margins(1) = 0.07;
+                    </scilab:image>
+                </td>
+            </tr>
+            <tr>
+                <td>
+                    <para/>
+                    <para>
+                        <emphasis role="bold">Fitting a gaussian curve + sloping base (5 parameters)</emphasis>:
+                    </para>
+                    <programlisting role="example"><![CDATA[
+function Y = model(P, X)
+    Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
+endfunction
+// -------------------------------
+// Gap function
+function r = gap(P,XY)
+    r = XY(2,:) - model(P,XY(1,:))
 endfunction
+// -------------------------------
+// True parameters
+Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
 
-X=[];Y=[];
-pg=[34;12;14]
-for x=0:.1:3
-  Y=[Y,FF(x,pg)+100*(rand()-.5)];
-  X=[X,x];
-end
-Z=[Y;X];
+// Generating data
+// ---------------
+X = 0:0.2:10;
+Np = length(X);
+Y0 = model(Pg, X);              // True law
+noise = (rand(Y0)-.5)*Pg(3)/4;
+Y = Y0 + noise                  // Noised data
+Data = [X ; Y];
+
+// Trying to recover original parameters from the (noisy) data:
+// -----------------------------------------------------------
+// Performing the fitting = parameters optimization
+[v, k] = max(Y0);
+P0 = [X(k) 1 v 1 1];
+[P, dmin, s]  = datafit(gap, Data, P0);
+Wd =  1./abs(noise);
+[Pw, dminw, s] = datafit(gap, Data, Wd, P0);
+
+// Computing fitting curves from found parameters
+fY = model(P, X);
+fYw = model(Pw, X);
 
-p0=[3;5;10]
-[p,err]=datafit(G,DG,Z,p0);
-scf(1);
+// Display
+// -------
+scf(0);
 clf()
-plot2d(X,FF(X,pg),5) //the curve without noise
-plot2d(X,Y,-1)  // the noisy data
-plot2d(X,FF(X,p),12) //the solution
- ]]></programlisting>
-        <scilab:image>
-            function y=FF(x,p)
-            y=p(1)*(x-p(2))+p(3)*x.*x
-            endfunction
-
-            function s=DG(p,z),
-            a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),
-            s=-[x-b,-a,x*x]
-            endfunction
-
-            function e=G(p,z),
-            y=z(1),x=z(2);
-            e=y-FF(x,p),
-            endfunction
-
-            X=[];Y=[];
-            pg=[34;12;14]
-            for x=0:.1:3
-            Y=[Y,FF(x,pg)+100*(rand()-.5)];
-            X=[X,x];
-            end
-            Z=[Y;X];
-
-            p0=[3;5;10]
-            [p,err]=datafit(G,DG,Z,p0);
-            scf(1);
-            clf()
-            plot2d(X,FF(X,pg),5)
-            plot2d(X,Y,-1)
-            plot2d(X,FF(X,p),12)
-        </scilab:image>
-        <programlisting role="example"><![CDATA[
-//generate the data
-function y=FF(x,p)
-  y=p(1)*(x-p(2))+p(3)*x.*x
+xsetech([0 0 1 0.8])
+plot2d(X, Y0, 24)  // True underlaying law
+plot2d(X, Y,  -1)  // Noisy data
+plot2d(X, fY, 15)  // best unweighted fit
+plot2d(X, fYw,18)  // best weighted fit
+
+// Legends: Average vertical Model-Data distance:
+// True one v/s according to the residual cost
+L_1 = "Unweighted fit (%5.2f /%5.2f)";
+L_1 = msprintf(L1, mean(abs(fY-Y0)), sqrt(dmin/Np));
+L_2 = "Weighted fit (%5.2f /%5.2f)";
+L_2 = msprintf(L1, mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)));
+legend(["True law", "Noisy Data", L_1, L_2],1)
+
+// Params: row#1: true params row#2:
+//  from unweighted fit row#3: from weighted fit
+[xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
+xnumb(xs, ys, [Pg ; P ; Pw]);
+LawTxt = "$\Large p_3\,.\,"+..
+         "exp\left({-(x-p_1)^2}\over{2\,p_2} \right)+p_4.x+p_5$";
+xstring(2, 1.3, LawTxt)
+
+// Plotting residus:
+xsetech([0 .75 1 0.25])
+plot2d(X, model(P ,X)-Y0, 15) // remaining gap with best unweighted fit
+plot2d(X, model(Pw,X)-Y0, 18) // remaining gap best weighted fit
+ylabel("Data fit - True law")
+ax = gca();
+ax.x_location = "top";
+tmp = ax.x_ticks.labels;
+gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+xgrid(color("grey70"))
+         ]]></programlisting>
+                </td>
+                <td>
+                    <scilab:image>
+                        function Y = model(P,X)
+                        Y = P(3) * exp(-(X-P(1)).^2 / (2*P(2)*P(2))) + P(4)*X + P(5);
+                        endfunction
+                        // -------------------------------
+                        // Gap function
+                        function r = gap(P,XY)
+                        r = XY(2,:) - model(P,XY(1,:))
+                        endfunction
+                        // -------------------------------
+                        // True parameters
+                        Pg = [3.5 1.5 7 0.4 -0.5]; // mean stDev ymax a(.x) b
+
+                        // Generating data
+                        // ---------------
+                        X = 0:0.2:10;
+                        Np = length(X);
+                        Y0 = model(Pg, X);              // True law
+                        noise = (rand(Y0)-.5)*Pg(3)/4;
+                        Y = Y0 + noise                  // Noised data
+                        Data = [X ; Y];
+
+                        // Let's try to find the original parameters from the (noisy) data:
+                        // ---------------------------------------------------------------
+                        // Performing the fitting = parameters optimization
+                        [v, k] = max(Y0);
+                        P0 = [X(k) 1 v 1 1];
+                        [P, dmin, s]  = datafit(gap, Data, P0);
+                        Wd =  1./abs(noise);
+                        [Pw, dminw, s] = datafit(gap, Data, Wd, P0);
+
+                        // Computing fitting curves from found parameters
+                        fY = model(P, X);
+                        fYw = model(Pw,X);
+
+                        // Display
+                        // -------
+                        scf(0);
+                        clf()
+                        xsetech([0 0 1 0.82])
+                        plot2d(X, Y0, 24)  // True underlaying law
+                        plot2d(X, Y,  -1)  // Noisy data
+                        plot2d(X, fY, 15)  // best unweighted fit
+                        plot2d(X, fYw,18)  // best weighted fit
+                        legend(["True law" "Noisy Data" ..
+                        msprintf("Unweighted fit (%5.2f /%5.2f)",mean(abs(fY-Y0)), sqrt(dmin/Np)) ..
+                        msprintf("Weighted fit (%5.2f /%5.2f)", mean(abs(fYw-Y0)), sqrt(dminw/sum(Wd)))],1)
+                        // Average vertical Model-Data distance: True one / according to the residual cost
+                        // Params: row#1: true params row#2: from unweighted fit row#3: from weighted fit
+                        [xs, ys] = meshgrid(linspace(2,6,5), linspace(0.5,-0.5,3));
+                        xnumb(xs, ys, [Pg ; P ; Pw])
+                        xstring(2,1.3,"$\Large p_3\,.\,exp\left({-(x-p_1)^2}\over{2\,p_2} \right)+p_4.x+p_5$")
+
+                        // Plotting residus:
+                        xsetech([0 .75 1 0.25])
+                        plot2d(X, model(P ,X)-Y0, 15)  // remaining gap with best unweighted fit
+                        plot2d(X, model(Pw,X)-Y0, 18)  // remaining gap best weighted fit
+                        ylabel("Data fit - True law")
+                        ax = gca();
+                        ax.x_location = "top";
+                        tmp = ax.x_ticks.labels;
+                        gca().x_ticks.labels = emptystr(size(tmp,"*"),1);
+                        xgrid(color("grey70"))
+                        // ===============
+                        gcf().axes_size = [500 600];
+                    </scilab:image>
+                </td>
+            </tr>
+            <tr>
+                <td>
+                    <para/>
+                    <para>
+                        <emphasis role="bold">Extraction of 3 overlapping gaussian curves (9 parameters)</emphasis>:
+                        example with constrained bounds and a matrix of parameters.
+                    </para>
+                    <programlisting role="example"><![CDATA[
+// Basic gaussian component:
+function Y = gauss(X, average, stDev, ymax)
+    Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
 endfunction
 
-//the gradient of the criterion function
-function s=DG(p,z),
-  a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),
-  s=-[x-b,-a,x*x]
+// Model: cumulated gaussian laws
+function r = model(P, X)
+    // P(1,:): averages  P(2,:): standard deviation  P(3,:): ymax
+    r = 0;
+    for i = 1:size(P,2)
+        r = r + gauss(X, P(1,i), P(2,i), P(3,i));
+    end
 endfunction
 
-function e=G(p,z),
-  y=z(1),x=z(2);
-  e=y-FF(x,p),
+// Gap function:
+function r = gap(P,XY)
+    r = XY(2,:) - model(P,XY(1,:))
 endfunction
 
-X=[];Y=[];
-pg=[34;12;14]
-for x=0:.1:3
-  Y=[Y,FF(x,pg)+100*(rand()-.5)];
-  X=[X,x];
-end
-Z=[Y;X];
+// Actual parameters :
+Pa = [ 1.1  4    7    // averages
+        1  0.8  1.5   // standard deviations
+        3   2   5.5   // ymax
+      ];
+Nc = size(Pa,2);      // Number of summed curves
+Np = 200;             // Number of data points
+
+// Data generation (without noise)
+xmin = 0;
+xmax = 10;
+X = linspace(xmin-2, xmax+2, Np);
+Y = model(Pa, X);
+Data = [X ; Y];          // Curve to analyze / fit
 
-p0=[3;5;10]
+// FITTING
+// -------
+// Fitting parameters:
+// Initial values & bounds (amplitudes must be > 0)
+Po = [ linspace(xmin,xmax,Nc)
+       ones(1,Nc)*0.5
+       ones(1,Nc)*max(Y)/2 ];
+Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
+Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3   ; ones(1,Nc)*max(Y)*1.2];
+// Fitting
+[P, dmin] = datafit(gap,Data,'b',Pinf,Psup,Po,'qn','ar',200,200);
 
-// Add some bounds on the estimate of the parameters
-// We want positive estimation (the result will not change)
-[p,err]=datafit(G,DG,Z,'b',[0;0;0],[%inf;%inf;%inf],p0,algo='gc');
-scf(1);
+// Display
+// -------
 clf()
-plot2d(X,FF(X,pg),5) //the curve without noise
-plot2d(X,Y,-1)  // the noisy data
-plot2d(X,FF(X,p),12) //the solution
- ]]></programlisting>
-        <scilab:image>
-            function y=FF(x,p)
-            y=p(1)*(x-p(2))+p(3)*x.*x
-            endfunction
-
-            function s=DG(p,z),
-            a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),
-            s=-[x-b,-a,x*x]
-            endfunction
-
-            function e=G(p,z),
-            y=z(1),x=z(2);
-            e=y-FF(x,p),
-            endfunction
-
-            X=[];Y=[];
-            pg=[34;12;14]
-            for x=0:.1:3
-            Y=[Y,FF(x,pg)+100*(rand()-.5)];
-            X=[X,x];
-            end
-            Z=[Y;X];
-
-            p0=[3;5;10]
-
-            [p,err]=datafit(G,DG,Z,'b',[0;0;0],[%inf;%inf;%inf],p0,algo='gc');
-            scf(1);
-            clf()
-            plot2d(X,FF(X,pg),5)
-            plot2d(X,Y,-1)
-            plot2d(X,FF(X,p),12)
-        </scilab:image>
+xsetech([0 0 1 0.8]);
+plot(X, Y, "-b", X, model(P,X), "-r");
+gca().children.children(1:2).thickness = 2;
+for i = 1:Nc
+    // True gaussian components
+    plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c");
+    // Extracted gaussian components
+    plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g");
+end
+legend("Data", "Parametric Model", ..
+       "True components", "Extracted components",  2);
+
+xsetech([0 0.75 1 0.25])
+plot(X, model(P,X)-Y);
+gca().x_location = "top";
+gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
+ylabel("Model - Data")
+xgrid(color("grey70"))
+         ]]></programlisting>
+                </td>
+                <td>
+                    <scilab:image>
+                        // Basic gaussian component:
+                        function Y = gauss(X, average, stDev, ymax)
+                        Y = ymax * exp(-(X-average).^2 / (2*stDev*stDev))
+                        endfunction
+
+                        // Model: cumulated gaussian laws
+                        function r = model(P, X)
+                        // P(1,:): averages  P(2,:): standard deviation  P(3,:): ymax
+                        r = 0;
+                        for i = 1:size(P,2)
+                        r = r + gauss(X, P(1,i), P(2,i), P(3,i));
+                        end
+                        endfunction
+
+                        // Gap function:
+                        function r = gap(P,XY)
+                        r = XY(2,:) - model(P,XY(1,:))
+                        endfunction
+
+                        // True parameters :
+                        Pa = [ 1.1  4    7    // averages
+                        1  0.8  1.5   // standard deviations
+                        3   2   5.5   // ymax
+                        ];
+                        Nc = size(Pa,2);       // Number of summed curves
+                        Np = 200;             // Number of data points
+
+                        // Data generation (without noise)
+                        xmin = 0;
+                        xmax = 10;
+                        X = linspace(xmin-2, xmax+2, Np);
+                        Y = model(Pa, X);
+                        Data = [X ; Y];          // Curve to analyze / fit
+
+                        // FITTING
+                        // -------
+                        // Fitting parameters: initial values and bounds (amplitudes must be > 0)
+                        Po = [ linspace(xmin,xmax,Nc) ; ones(1,Nc)*0.5 ; ones(1,Nc)*max(Y)/2];
+                        Pinf = [ones(1,Nc)*-2 ; ones(1,Nc)*0.1 ; ones(1,Nc)*1];
+                        Psup = [ones(1,Nc)*12 ; ones(1,Nc)*3   ; ones(1,Nc)*max(Y)*1.2];
+                        // Fitting
+                        [P, dmin] = datafit(gap, Data, 'b',Pinf,Psup, Po, 'qn','ar',200,200);
+
+                        // Display
+                        // -------
+                        clf()
+                        xsetech([0 0 1 0.8]);
+                        plot(X, Y, "-b", X, model(P,X), "-r");
+                        gca().children.children(1:2).thickness = 2;
+                        for i = 1:Nc
+                        plot(X, gauss(X, Pa(1,i), Pa(2,i), Pa(3,i)), "-c"); // True gaussian components
+                        plot(X, gauss(X, P(1,i), P(2,i), P(3,i)), "-g");    // Extracted gaussian components
+                        end
+                        legend("Data", "Parametric Model","True components", "Extracted components",  2);
+
+                        xsetech([0 0.75 1 0.25])
+                        plot(X, model(P,X)-Y);
+                        gca().x_location = "top";
+                        gca().x_ticks.labels = emptystr(size(gca().x_ticks.labels,"*"),1);
+                        ylabel("Model - Data")
+                        xgrid(color("grey70"))
+                        // ===============
+                        gcf().axes_size = [500 500];
+                    </scilab:image>
+                </td>
+            </tr>
+            <tr>
+                <td>
+                    <para>
+                        <emphasis role="bold">With a parametric curve: Fitting a off-centered tilted 2D elliptical arc (5 parameters)</emphasis>.
+                    </para>
+                    <programlisting role="example"><![CDATA[
+function g = Gap(p, Data)
+    // p: a, b, center xc, yc, alpha tilt °
+    C = cosd(p(5)); S = -sind(p(5));
+    x = Data(1,:) - p(3);
+    y = Data(2,:) - p(4);
+    g = ((x*C - y*S )/p(1)).^2 + ((x*S + y*C)/p(2)).^2 - 1;
+endfunction
+
+// Generating the data
+// -------------------
+// Actual parameters :
+Pa = [3, 1.3, -2, 1.5, 30];
+Np = 30;            // Number of data points
+Theta = grand(1,Np, "unf",-100, 170);
+// untilted centered noisy arc
+x = Pa(1)*(cosd(Theta) + grand(Theta, "unf", -0.07, 0.07));
+y = Pa(2)*(sind(Theta) + grand(Theta, "unf", -0.07, 0.07));
+// Tilting and off-centering the arc
+A = Pa(5);
+C = cosd(A); S = sind(A);
+xe =  C*x + y*S + Pa(3);
+ye = -S*x + y*C + Pa(4);
+Data = [xe ; ye];
+
+// Retrieving parameters from noisy data
+// -------------------------------------
+// Initial estimation
+ab0 = (strange(xe)+strange(ye))/2;
+xc0 = mean(xe);
+yc0 = mean(ye);
+A0 = -atand(reglin(xe,ye));
+P0 = [ab0 ab0/2 xc0 yc0 A0];
+// Parameters bounds
+Pinf = [ 0     0     xc0-2*ab0, yc0-2*ab0 -90];
+Psup = [3*ab0 3*ab0  xc0+2*ab0, yc0+2*ab0  90];// Fitting
+
+// FITTING
+[P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
+
+// P(1) is the longest axis:
+if P(1)<P(2) then
+    P([1 2]) = P([2 1]);
+    P(5) = 90 - P(5);
+end
+
+// DISPLAY
+// -------
+disp([Pa;P], dmin);
+// A successful fit should lead to dmin ~ noise amplitude
+
+clf
+plot(xe, ye, "+")   // Data
+// Model
+Theta = 0:2:360;
+x = P(1) * cosd(Theta);
+y = P(2) * sind(Theta);
+A = P(5);
+[x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
+plot(x, y, "r")
+isoview
+// Parameters values:
+Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
+Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
+xstring(-3.7, 1.2, ["" "    a           b        xc       yc     A°"
+                    "Actual:" Patxt
+                    "Fit:" Ptxt ])
+xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
+         ]]></programlisting>
+                </td>
+                <td>
+                    <scilab:image><![CDATA[
+function g = Gap(p, Data)
+    // p: a, b, center xc, yc, alpha tilt °
+    C = cosd(p(5)); S = -sind(p(5));
+    x = Data(1,:) - p(3);
+    y = Data(2,:) - p(4);
+    g = ((x*C - y*S )/p(1)).^2 + ((x*S + y*C)/p(2)).^2 - 1;
+endfunction
+
+// Generating the data
+// -------------------
+// Actual parameters :
+Pa = [3, 1.3, -2, 1.5, 30];
+Np = 30;            // Number of data points
+Theta = grand(1,Np, "unf",-100, 170);
+// untilted centered noisy arc
+x = Pa(1)*(cosd(Theta) + grand(Theta, "unf",-0.07, 0.07));
+y = Pa(2)*(sind(Theta) + grand(Theta, "unf",-0.07, 0.07));
+// Tilting and off-centering the arc
+A = Pa(5);
+C = cosd(A); S = sind(A);
+xe =  C*x + y*S + Pa(3);
+ye = -S*x + y*C + Pa(4);
+Data = [xe ; ye];
+
+// Retrieving parameters from noisy data
+// -------------------------------------
+// Initial estimation
+ab0 = (strange(xe)+strange(ye))/2;
+xc0 = mean(xe);
+yc0 = mean(ye);
+A0 = -atand(reglin(xe,ye));
+P0 = [ab0 ab0/2 xc0 yc0 A0];
+// Parameters bounds
+Pinf = [ 0     0     xc0-2*ab0, yc0-2*ab0 -90];
+Psup = [3*ab0 3*ab0  xc0+2*ab0, yc0+2*ab0  90];// Fitting
+
+// FITTING
+[P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
+
+// P(1) is the longest axis:
+if P(1)<P(2) then
+    P([1 2]) = P([2 1]);
+    P(5) = 90 - P(5);
+end
+
+// DISPLAY
+// -------
+clf
+plot(xe, ye, "+")   // Data
+// Model
+Theta = 0:2:360;
+x = P(1) * cosd(Theta);
+y = P(2) * sind(Theta);
+A = P(5);
+[x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
+plot(x, y, "r")
+isoview
+// Parameters values:
+Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
+Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
+xstring(-3.7, 1.2, ["" "    a           b        xc       yc     A°"
+                    "Actual:" Patxt
+                    "Fit:" Ptxt ]);
+xstring(-2.5, 0.9, "dmin = " + msprintf("%5.3f", dmin))
+title("datafit() example #4", "fontsize", 3)
+// ===============
+gcf().axes_size = [500 500];
+gcf().children(1).margins(1) = 0.05;
+]]>             </scilab:image>
+                </td>
+            </tr>
+        </table>
     </refsection>
     <refsection role="see also">
         <title>See also</title>
@@ -424,4 +956,31 @@ plot2d(X,FF(X,p),12) //the solution
             </member>
         </simplelist>
     </refsection>
+    <refsection role="history">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1</revnumber>
+                <revdescription>
+                    <para>
+                        New option <varname>Wd</varname> to weight measured data <varname>Data</varname>.
+                    </para>
+                    <para>
+                        The gap function G() may now be vectorized over <literal>Data</literal> points.
+                    </para>
+                    <para>
+                        Initial parameters <varname>p0</varname> and their bounds may now be
+                        provided as a matrix instead of a mandatory column.
+                    </para>
+                    <para>
+                        The <literal>err</literal> output argument is replaced with
+                        <literal>dmin</literal> = average Model-to-data distance.
+                    </para>
+                    <para>
+                        New output argument <literal>status</literal> added (for the "qn" algorithm).
+                    </para>
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index 9bd398f..3122ad3 100644 (file)
@@ -1,7 +1,7 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA
-//
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2018 - Samuel GOUGEON
 //
 // This file is hereby licensed under the terms of the GNU GPL v2.0,
 // pursuant to article 5.3.4 of the CeCILL v.2.1.
 // and continues to be available under such terms.
 // For more information, see the COPYING file which you should have received
 // along with this program.
-//
-function [p,err]=datafit(iprint,G,varargin)
-    //
-    //         [p,err]=datafit([iprint,] G [,DG],Z [,W],...)
-    //
-    //         Function used for fitting data to a model.
-    // For a given function G(p,z), this function finds the best vector
-    // of parameters p for approximating G(p,z_i)=0 for a set of measurement
-    // vectors z_i. Vector p is found by minimizing
-    //    G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)
-    //
-    //      G: Scilab function (e=G(p,z), e: nex1, p: npx1, z: nzx1)
-    //     p0: initial guess (size npx1)
-    //      Z: matrix [z_1,z_2,...z_n] where z_i (nzx1) is the ith measurement
-    //      W: weighting matrix of size nexne (optional)
-    //     DG: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
-    //
-    //                     Examples
+
+function [p, dmin, status] = datafit(iprint, G, varargin)
     //
-    //deff('y=FF(x)','y=a*(x-b)+c*x.*x')
-    //X=[];Y=[];
-    //a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
-    //Z=[Y;X];
-    //deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
-    //[p,err]=datafit(G,Z,[3;5;10])
-    //scf(0)
-    //clf();
-    //plot2d(X',Y',-1)
-    //plot2d(X',FF(X)',5,'002')
-    //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
+    //         [p,err]=datafit([iprint,] G [,DG], Data [,Wd] [,Wg],...)
     //
-    //a=34;b=12;c=14;
-    //deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
-    //[p,err]=datafit(G,DG,Z,[3;5;10])
-    //scf(1)
-    //clf();
-    //plot2d(X',Y',-1)
-    //plot2d(X',FF(X)',5,'002')
-    //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
+    // Data: matrix where Data(:,i) are the nzc coordinates of the ith measurement point
+    //       nz  = size(Data,"c") : number of measured available points.
+    //       nzc = size(Data,"r") : number of coordinates for each Data point.
+    //    G: Scilab function (e = G(p,Data) (e: ne x nz)
+    //   p0: initial guess of parameters values (vector or matrix)
+    //       np = length(p0) = length(p) : number of fitting parameters
+    //   Wd: Data weights (optionnal). Row of length nz.
+    //   Wg: weighting matrix of gaps definitions, of size ng x ng (optional)
+    //   DG: partial of G wrt p (optional; S=DG(p,Data), S: ne x np x nz)
     //
+    //      p: vector of best parameters values (minimizing dmin)
+    // status: integer as returned by optim() when the default "qn" algo is used.
+    //   dmin: scalar = average distance of data points to model.
 
-    [lhs,rhs]=argn(0)
+    [lhs,rhs] = argn(0)
 
-    if type(iprint)<>1 then
-        wflag=warning("query") // Disable warnings as the following lines may produce some
+    // Processing the 2 first argin imp and G:
+    if type(iprint)<>1 then      // 1st argin iprint is actually G (no iprint directive)
+        wflag = warning("query") // Disable warnings as the following lines may produce some
         warning("off")
-        varargin(0)=G
-        G=iprint
-        iprint=0
+        varargin(0) = G     // so the argin G is the first arg after the true G
+        G = iprint          // G gets its true value
+        iprint = 0          // iprint gets its default value
         warning(wflag)
     end
 
-    if type(G)==15 then
-        Gparams=G;Gparams(1)=null();
+    if type(G)==15 then     // G = list(G, params)
+        Gparams = G;
+        Gparams(1) = null();
         G=G(1)
     else
-        Gparams=list()
+        Gparams = list()    // no explicit parameter
     end
 
-
-    DG=varargin(1)
-    if type(DG)==10|type(DG)==13 then
-        GR=%t  //Jacobian provided
-        varargin(1)=null()
+    // Next argins: optional dG (Jacobian), or mandatory points Data
+    DG = varargin(1)
+    if type(DG)==10 | type(DG)==13 then   // is a function
+        GR = %t      // Jacobian provided
+        varargin(1) = null()
     elseif type(DG)==15 then
-        error(msprintf(gettext("%s: Jacobian cannot be a list, parameters must be set in G."),"datafit"));
+        msg = gettext("%s: Jacobian cannot be a list, parameters must be set in G.")
+        error(msprintf(msg, "datafit"));
     else
-        GR=%f
+        GR = %f  // Jacobian not provided
     end
 
-    Z=varargin(1);
-    varargin(1)=null()
-    if type(Z)<>1 then
-        error(msprintf(gettext("%s: Wrong measurement matrix."),"datafit"));
+    Data = varargin(1);          // Data points to fit
+    varargin(1) = null()
+    if type(Data) <>1 | ~isreal(Data,0) then
+        msg = gettext("%s: Data points must be real numbers.");
+        error(msprintf(msg, "datafit"));
     end
+    [mz,nz] = size(Data);    // nz: Number of points
 
-    nv=size(varargin)
-    if nv>=1 then
-        if size(varargin(1),2)==1 then // p0 ou 'b'
-            W=1
-        else
-            W=varargin(1);varargin(1)=null()
-            if size(W,1)~=size(W,2) then
-                if size(W,1)==1 then
-                    error(msprintf(gettext("%s: Initial guess must be a column vector."),"datafit"));
-                else
-                    error(msprintf(gettext("%s: Weighting matrix must be square."),"datafit"));
-                end
+    // Data weights: Must be a line vector of length nz (Wg can _neither_ be so long)
+    if size(varargin(1),2)==nz
+        Wd = varargin(1)
+        varargin(1) = null()
+    else
+        Wd = ones(1,nz);
+    end
+    // Wg assignment / initialization  (from [,Wg][,'b', pinf, psup], p0 [,algo])
+    nv = size(varargin);
+    if nv==0 then
+        msg =  gettext("%s: Initial guess p0 is missing.");
+        error(msprintf(msg, "datafit"));
+    elseif nv==1 then   // p0 expected
+        Wg = 1;
+    else
+        v1 = varargin(1);
+        v2 = varargin(2);
+        if v1=="b"| (nv>=5 & or(varargin(5)==["qn" "gc" "nd" "ar" "in"]))
+            Wg = 1
+        elseif (nv>1 & v2=="b") | ..
+               (nv>=6 & or(varargin(6)==["qn" "gc" "nd" "ar" "in"])) | ..
+               (type(v1)==1 & issquare(v1) & nv>1 & type(v2)==1 & (size(v2,"*")>1 | nv==2))
+            Wg = varargin(1)
+            if size(Wg,1) ~= size(Wg,2) then
+                msg = gettext("%s: Weighting matrix must be square.");
+                error(msprintf(msg, "datafit"));
             end
+            varargin(1) = null()
+            nv = nv-1
+        else
+            Wg = 1;
         end
     end
-    if type(varargin(1))==1 then // p0
-        p0=varargin(1)
+    // next: [,'b', pinf, psup], p0 [,algo])
+    if varargin(1)=="b" then
+        if nv<4
+            msg = gettext("%s: error with parameters bounds or initial guess\n");
+            error(msprintf(msg, "datafit"));
+        end
+        p0 = varargin(4);
     else
-        p0=varargin(4)
+        p0 = varargin(1);
     end
+    sp0 = size(p0);
+    np = size(p0,"*");    // Number of parameters
+    _format = format()
+    format("v",23)        // for string()
 
-    [mz,nz]=size(Z);np=size(p0,"*");
-
-    if type(G)==10 then  //form function to call hard coded external
+    if type(G)==10 then   // form function to call hard coded external
         if size(Gparams)==0 then
-            error(msprintf(gettext("%s: With hard coded function, user must give output size of G."),"datafit"));
+            format(_format(2),_format(1))
+            msg = gettext("%s: With hard coded function, user must give output size of G.")
+            error(msprintf(msg, "datafit"));
         end
-        m=Gparams(1);Gparams(1)=null()
+        m = Gparams(1);
+        Gparams(1) = null();
 
-        // foo(m,np,p,mz,nz,Z,pars,f)
-        deff("f=G(p,Z)","f=call(''"+G+"'',"+..
-        "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
-        "pars,7,''out'',["+string(m)+",1],8,''d'')")
+        // foo(m,np,p,mz,nz,Data[,Wd],pars,f)
+        tmp = "f=call(''" + G + "''," + ..
+              "m,1,''i'', np,2,''i'', p,3,''d'', mz,4,''i'', nz,5,''i''," + ..
+              "Data,6,''d'', pars,7,''out'',[" + string(m) + ",1],8,''d'')"
+        deff("f = G(p,Data)", tmp)
 
-        pars=[];
-        for k=1:size(Gparams)
-            p=Gparams(k)
-            pars=[pars;p(:)]
+        pars = [];
+        for k = 1:size(Gparams)
+            p = Gparams(k)
+            pars = [pars;p(:)]
         end
-        Gparams=list()
+        Gparams = list()
     end
 
-    if type(DG)==10 then //form function to call hard coded external
-        // dfoo(m,np,p,mz,nz,Z,pars,f)
-        deff("f=DG(p,Z)","f=call(''"+DG+"'',"+..
-        "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
-        "pars,7,''out'',["+string(m)+","+string(np)+"],8,''d'')")
+    if type(DG)==10 then // form function to call hard coded external
+        // dfoo(m,np,p,mz,nz,Data[,Wd],pars,f)
+        tmp = "f = call(''" + DG + "''," + ..
+          "m,1,''i'', np,2,''i'', p,3,''d'', mz,4,''i'', nz,5,''i'',Data,6,''d''," + ..
+          "pars,7,''out'',[" + string(m) + "," + string(np) + "],8,''d'')"
+        deff("f = DG(p,Data)", tmp)
     end
 
+    // Is the G function vectorized?
+    try
+        if Gparams==list() then
+            e = G(p0,[Data(:,1) Data(:,1)]);
+        else
+            e = G(p0,[Data(:,1) Data(:,1)], Gparams(:));
+        end
+        isGvec = size(e,2)==2
+    catch
+        isGvec = %f;
+    end
 
     // form square cost gradient function DGG
-
     if Gparams==list() then
-        GP   = "G(p,Z(:,i))"
-        GPPV = "G(p+v,Z(:,i))"
-        DGP  = "DG(p,Z(:,i))"
+        GP   = "G(p,  Data(:,i))"
+        GPPV = "G(p+v,Data(:,i))"
+        DGP  = "DG(p, Data(:,i))"
     else
-        GP   = "G(p,Z(:,i),Gparams(:))"
-        GPPV = "G(p+v,Z(:,i),Gparams(:))"
-        DGP  = "DG(p,Z(:,i),Gparams(:))"
+        GP   = "G(p,  Data(:,i), Gparams(:))"
+        GPPV = "G(p+v,Data(:,i), Gparams(:))"
+        DGP  = "DG(p, Data(:,i), Gparams(:))"
     end
 
     if ~GR then // finite difference
-        DGG=["g=0*p";
-        "pa=sqrt(%eps)*(1+1d-3*abs(p))"
-        "f=0;"
-        "for i=1:"+string(nz)
-        "  g1="+GP
-        "  f=f+g1''*W*g1"
-        "end"
-        "for j=1:"+string(np)
-        "  v=0*p;v(j)=pa(j),"
-        "  e=0;"
-        "  for i=1:"+string(nz)
-        "    g1="+GPPV
-        "    e=e+g1''*W*g1"
-        "  end"
-        "  g(j)=e-f;"
-        "end;"
-        "g=g./pa;"]
-    else // using Jacobian of G
-        DGG="g=0;for i=1:nz,g=g+2*"+DGP+"''*W*"+GP+",end"
+        // Version vectorized over Data:
+        if isGvec then
+            DGG = [
+            "g  = 0*p"
+            "pa = sqrt(%eps)*(1+1d-3*abs(p))"
+            "i = : "
+            "g1 = " + GP
+            "f = sum(sum((g1''*Wg) .* g1'',""c"").*Wd'');"
+            "for j = 1:" + msprintf("%d\n",np)
+            "    v = 0*p;"
+            "    v(j) = pa(j);"
+            "    g1 = " + GPPV
+            "    e = sum(sum((g1''*Wg) .* g1'',""c"").*Wd'');"
+            "    g(j) = e - f;"
+            "end"
+            "g = g ./ pa;"
+            ];
+        else
+            DGG = [
+            "g  = 0*p"
+            "pa = sqrt(%eps)*(1+1d-3*abs(p))"
+            "f = 0;"
+            "for i = 1:" + msprintf("%d\n", nz)
+            "    g1 = " + GP
+            "    f = f + g1''*Wg*g1 * Wd(i)"
+            "end"
+            "for j = 1:" + msprintf("%d\n", np)
+            "    v = 0*p;"
+            "    v(j) = pa(j),"
+            "    e = 0;"
+            "    for i = 1:" + msprintf("%d\n", nz)
+            "        g1 = " + GPPV
+            "        e = e + g1''*Wg*g1 * Wd(i)"
+            "    end"
+            "    g(j) = e - f;"
+            "end"
+            "g = g ./ pa;"
+            ];
+        end
+    else // using Jacobian of G  (UNVECTORIZED)
+        DGG = "g=0; for i=1:nz, g = g + 2*"+DGP+"''*Wg*"+GP+", end"
     end
 
-    // form cost function for optim
-    deff("[f,g,ind]=costf(p,ind)",[
-    "if ind==2|ind==4 then "
-    "  f=0;"
-    "   for i=1:"+string(nz)
-    "     g1="+GP
-    "     f=f+g1''*W*g1"
-    "   end"
-    "else "
-    "  f=0;"
-    "end";
-    "if ind==3|ind==4 then"
-    DGG
-    "else"
-    " g=0*p;"
-    "end"])
-
-    [err,p]=optim(costf,varargin(:),iprint=iprint)
+    // Defining the costf() function for optim:
+    if isGvec then
+        tmp = [
+            "    i = : "
+            "    g1 = " + GP
+            "    f = sum(sum((g1''*Wg) .* g1'',""c"") .* Wd'');"
+            ]
+    else
+        tmp = [
+            "    f = 0;"
+            "    for i = 1:" + msprintf("%d\n", nz)
+            "        g1 = " + GP
+            "        f = f + g1''*Wg*g1 * Wd(i)"
+            "    end"
+            ]
+    end
+    deff("[f,g,ind] = costf(p,ind)",[
+        "if ind==2 | ind==4 then "
+        tmp
+        "else"
+        "    f = 0;"
+        "end";
+        "if ind==3 | ind==4 then"
+        DGG
+        "else"
+        "    g = 0*p;"
+        "end"
+        ])
+
+    format(_format(2), _format(1))
+    qn = %t;
+    for v = varargin
+        if or(v==['gc' 'nd'])
+            qn = %f;
+            break
+        end
+    end
+    if qn then
+        [fmin, p, t,t,t,t, status] = optim(costf, varargin(:), iprint=iprint)
+    else
+        [fmin, p] = optim(costf, varargin(:), iprint=iprint)
+        status = %nan;
+    end
 
+    p = matrix(p, sp0);
 
+    // fmin => dmin: average Model-to-data distance
+    ng = size(G(p, Data(:,1)),1);
+    dmin = sqrt(fmin/ng/sum(Wd));
 endfunction
index 86ef079..fe03375 100644 (file)
@@ -70,4 +70,4 @@ assert_checkequal(v, zeros(6, 1));
 // datafit
 p0 = x0;
 [p, err] = datafit(G, Z, p0);
-assert_checkalmostequal([p err], [4 0], [], 1e-9);
+assert_checkalmostequal([p err], [4 0], [], 2e-9);
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_2330.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_2330.dia.ref
deleted file mode 100644 (file)
index b74edf0..0000000
+++ /dev/null
@@ -1,33 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2007-2008 - INRIA - Serge STEER <serge.steer@inria.fr>
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// <-- Non-regression test for bug 2330 -->
-//
-// <-- Bugzilla URL -->
-// http://bugzilla.scilab.org/show_bug.cgi?id=2330
-//
-// <-- Short Description -->
-//    datafit does not work equivalently on WinXP and Linux. With a given dataset and
-//    same routines centered on datafit function, it works perfectly on Linux and
-//    partially (some data are fitted some others not) on WinXP.
-//build the data to fit
-//---------------------
-function Xcalc=biexp(p,t)
-  Xcalc=p(1).*exp(-p(2).*t)+p(3).*exp(-p(4).*t)+p(5);
-endfunction;
-t=(0:100:36000)';
-p=[0.1;0.0001;0.2;0.0002;0.3];
-X=biexp(p,t);
-//try to fit the data
-//-------------------
-//the error function
-function e=myerf(p,X,t),e=X-biexp(p,t),endfunction
-// the initial point
-p0=[0.01;0.001;0.01;0.001;0.1];
-//call datafit
-[pr,err]=datafit(list(myerf,t),X,p0);
-if err>=5d-4 then bugmes();quit;end
index fffd5be..2032a95 100644 (file)
@@ -4,9 +4,10 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-
+//
 // <-- CLI SHELL MODE -->
-
+// <-- NO CHECK REF -->
+//
 // <-- Non-regression test for bug 2330 -->
 //
 // <-- Bugzilla URL -->
 //build the data to fit
 //---------------------
 function Xcalc=biexp(p,t)
-  Xcalc=p(1).*exp(-p(2).*t)+p(3).*exp(-p(4).*t)+p(5);
+  Xcalc = p(1).*exp(-p(2).*t) + p(3).*exp(-p(4).*t) + p(5);
 endfunction;
-t=(0:100:36000)'; 
-p=[0.1;0.0001;0.2;0.0002;0.3];
-X=biexp(p,t);
+t = (0:100:36000);
+p = [0.1 ; 0.0001 ; 0.2 ; 0.0002 ; 0.3];
+X = biexp(p,t);
 
 //try to fit the data
 //-------------------
 //the error function
-function e=myerf(p,X,t),e=X-biexp(p,t),endfunction
+function e=myerf(p,X)
+    t = X(1,:)
+    e = sum(abs(X(2,:) - biexp(p,t)))
+endfunction
 // the initial point
-p0=[0.01;0.001;0.01;0.001;0.1];
+p0 = [0.01 ; 0.001 ; 0.01 ; 0.001 ; 0.1];
 //call datafit
-[pr,err]=datafit(list(myerf,t),X,p0);
-if err>=5d-4 then pause,end
+[pr,err] = datafit(myerf,[t ; X], p0);
+if err >= 6d-4 then pause,end
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_244.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_244.dia.ref
deleted file mode 100644 (file)
index 6462829..0000000
+++ /dev/null
@@ -1,160 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Vincent COUVERT
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- TEST WITH GRAPHIC -->
-// <-- Non-regression test for bug 244 -->
-//
-// <-- Bugzilla URL -->
-// http://bugzilla.scilab.org/show_bug.cgi?id=244
-//
-// <-- Short Description -->
-//   recursion problems with fsolve
-// Titration of a dibase with HCl (at the beginning, a little amount of NaOH was added)
-// parameters
-//C0=2.39e-4
-C0=2e-4;
-C02=0.02;
-na0=2e-7;
-v0=0.01;
-pK1=6.9;
-pK2=6.9;
-// function to get the calculated pH
-deff('[x]=f0(pH)','x=(10^(-pH-pK1)+2*10^(-2*pH))/(10^(-pK1-pK2)+10^(-pH-pK1)+10^(-2*pH))+((v0+v)*(10^(-pH)-10^(pH-14))+na0-C02*v)/(C0*v0)');
-deff('[pH]=fpH(p)','v=p;pH=fsolve(7,f0)');
-// data
-X=[];Y=[];
-X=[
-//0    
-0.00001        
-0.00002        
-0.00003        
-0.00004        
-0.00005        
-0.000055       
-0.00006        
-0.000065       
-0.00007        
-0.000075       
-0.00008        
-0.000085       
-0.00009        
-0.000095       
-0.0001 
-0.000105       
-0.00011        
-0.000115       
-0.00012        
-0.000125       
-0.00013        
-0.000135       
-0.00014        
-0.000145       
-0.00015        
-0.000155       
-0.00016        
-0.000165       
-0.00017        
-0.000175       
-0.00018        
-0.000185       
-0.00019        
-0.000195       
-0.0002 
-0.000205       
-0.00021        
-0.000215       
-0.00022        
-0.000225       
-0.00023        
-0.000235       
-0.00024        
-0.000245       
-0.00025        
-0.000255       
-0.00026        
-0.000265       
-0.00027        
-0.000275       
-0.00028        
-0.000285       
-0.000295       
-0.000305       
-0.00032        
-0.000345       
-0.000395       
-0.000445       
-];
-Y=[
-//8.21 
-7.82   
-7.64   
-7.48   
-7.37   
-7.22   
-7.20   
-7.17   
-7.13   
-7.12   
-7.12   
-7.10   
-7.07   
-7.04   
-7.07   
-7.04   
-7.01   
-6.98   
-6.94   
-6.91   
-6.87   
-6.84   
-6.80   
-6.76   
-6.73   
-6.68   
-6.61   
-6.57   
-6.51   
-6.45   
-6.36   
-6.27   
-6.14   
-6.02   
-5.84   
-5.64   
-5.34   
-5.00   
-4.80   
-4.59   
-4.46   
-4.40   
-4.33   
-4.24   
-4.18   
-4.11   
-4.07   
-4.02   
-3.99   
-3.95   
-3.91   
-3.88   
-3.86   
-3.80   
-3.76   
-3.69   
-3.60   
-3.47   
-3.36   
-];
-// fitting
-Z=[Y;X];
-deff('e=G(p,z)','pK1=p(1),pK2=p(2),v=z(2),pHexp=z(1),e=pHexp-fpH(v)');
-[p,err]=datafit(G,Z,[6;7]);
-// graphic part
-clf()
-//v=[0:1e-5:4.5e-4]
-v=X;
-fplot2d(v,fpH);
-plot2d(X,Y,[-2]);
index 786c3ec..5653117 100644 (file)
@@ -6,6 +6,7 @@
 // =============================================================================
 
 // <-- TEST WITH GRAPHIC -->
+// <-- NO CHECK REF -->
 
 // <-- Non-regression test for bug 244 -->
 //
@@ -37,131 +38,131 @@ deff('[pH]=fpH(p)','v=p;pH=fsolve(7,f0)');
 X=[];Y=[];
 
 X=[
-//0    
-0.00001        
-0.00002        
-0.00003        
-0.00004        
-0.00005        
-0.000055       
-0.00006        
-0.000065       
-0.00007        
-0.000075       
-0.00008        
-0.000085       
-0.00009        
-0.000095       
-0.0001 
-0.000105       
-0.00011        
-0.000115       
-0.00012        
-0.000125       
-0.00013        
-0.000135       
-0.00014        
-0.000145       
-0.00015        
-0.000155       
-0.00016        
-0.000165       
-0.00017        
-0.000175       
-0.00018        
-0.000185       
-0.00019        
-0.000195       
-0.0002 
-0.000205       
-0.00021        
-0.000215       
-0.00022        
-0.000225       
-0.00023        
-0.000235       
-0.00024        
-0.000245       
-0.00025        
-0.000255       
-0.00026        
-0.000265       
-0.00027        
-0.000275       
-0.00028        
-0.000285       
-0.000295       
-0.000305       
-0.00032        
-0.000345       
-0.000395       
-0.000445       
+//0
+0.00001
+0.00002
+0.00003
+0.00004
+0.00005
+0.000055
+0.00006
+0.000065
+0.00007
+0.000075
+0.00008
+0.000085
+0.00009
+0.000095
+0.0001
+0.000105
+0.00011
+0.000115
+0.00012
+0.000125
+0.00013
+0.000135
+0.00014
+0.000145
+0.00015
+0.000155
+0.00016
+0.000165
+0.00017
+0.000175
+0.00018
+0.000185
+0.00019
+0.000195
+0.0002
+0.000205
+0.00021
+0.000215
+0.00022
+0.000225
+0.00023
+0.000235
+0.00024
+0.000245
+0.00025
+0.000255
+0.00026
+0.000265
+0.00027
+0.000275
+0.00028
+0.000285
+0.000295
+0.000305
+0.00032
+0.000345
+0.000395
+0.000445
 ];
 
 Y=[
-//8.21 
-7.82   
-7.64   
-7.48   
-7.37   
-7.22   
-7.20   
-7.17   
-7.13   
-7.12   
-7.12   
-7.10   
-7.07   
-7.04   
-7.07   
-7.04   
-7.01   
-6.98   
-6.94   
-6.91   
-6.87   
-6.84   
-6.80   
-6.76   
-6.73   
-6.68   
-6.61   
-6.57   
-6.51   
-6.45   
-6.36   
-6.27   
-6.14   
-6.02   
-5.84   
-5.64   
-5.34   
-5.00   
-4.80   
-4.59   
-4.46   
-4.40   
-4.33   
-4.24   
-4.18   
-4.11   
-4.07   
-4.02   
-3.99   
-3.95   
-3.91   
-3.88   
-3.86   
-3.80   
-3.76   
-3.69   
-3.60   
-3.47   
-3.36   
+//8.21
+7.82
+7.64
+7.48
+7.37
+7.22
+7.20
+7.17
+7.13
+7.12
+7.12
+7.10
+7.07
+7.04
+7.07
+7.04
+7.01
+6.98
+6.94
+6.91
+6.87
+6.84
+6.80
+6.76
+6.73
+6.68
+6.61
+6.57
+6.51
+6.45
+6.36
+6.27
+6.14
+6.02
+5.84
+5.64
+5.34
+5.00
+4.80
+4.59
+4.46
+4.40
+4.33
+4.24
+4.18
+4.11
+4.07
+4.02
+3.99
+3.95
+3.91
+3.88
+3.86
+3.80
+3.76
+3.69
+3.60
+3.47
+3.36
 ];
 
 // fitting
-Z=[Y;X];
+Z=[Y,X];
 deff('e=G(p,z)','pK1=p(1),pK2=p(2),v=z(2),pHexp=z(1),e=pHexp-fpH(v)');
 [p,err]=datafit(G,Z,[6;7]);