Discussion :

[Matthieu Brucher]

Tout à faut d'accord avec JP, un fit trop précis, et on va droit dans le mur -> over-fitting.
Il peut aussi être intéressant de calculer le meilleur modèle par maximum d'entropie sur le spectre des données.

[User]

ref j.p.mignot,

Oui ton algo m'intéresse, je te mail tout ça !

Merci encore !

Pour le string f je pensais utiliser 1 case (conditions) sur la string:

Je sélectionne f dans une liste déroulante contenant le type de fonction envisagée:

et dans la sub:

Si f="a + bx + cx^2" alors approx par fit polynomiale de degré 2
Si f="a*Exp(-bx)" alors faire routine approprié pour déterminer a et b
Si f="a*Sin(bx)" alors faire routine approprié pour déterminer a et b
Si f="a*sin(bx)*Exp(cx) alors faire routine approprié pour déterminer a, b et c
...

mais l'éventaille des fonctions possibles est tailement grand

Désolé,

@+

[j.p.mignot]

ci-joint un unit en pascal code pour fit polynômial. sous delphi
data_mx pourrair êrtre accru
Less square suivra des que j'ai le temps

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
 
{$n+,g+,a+}
unit math_rx;
interface
const     deg_mx     = 20;               { degre max du fit }
          data_mx    = 16383;             { nombre max de data }
 
type
          a0         = array[0..deg_mx,0..deg_mx] of extended;
          b0         = array[0..deg_mx]      of extended;
          xy         = array[1..data_mx]     of single;
 
function  polynome_fits(
                        px,py : pointer;  { pointe les tableaux de x,y de type single dont
                                            le nombre d'element ne doit pas
                                            exceder data_mx }
                        pf : pointer;     { pointe 1 tableau 0..dg of extended pour retour coef }
                        dg,ndat : byte;   { degre ( doit etre <= dg_mx) et nombre de data (doit etre <= ndata_mx )}
                        var ki2 : single; { pour retour du residu }
                        _weight : pointer) :   { pointeur sur tableau de poids (NIL si non souhait‚)}
                        boolean;
implementation
 
const     epsilon1   = 1e-25; { arbitrairement petit }
 
var
    a_mat     : a0;
    arot,
    piv,
    b_mat     : b0;
    brot      : extended;
    function puiss0(s : extended; v : byte) : extended;
    var se : extended;
       begin
       if v=0 then
          puiss0:=1
       else
       if abs(s) < 1e-40 then
          puiss0:=0
       else
          begin
          se:=exp(v*ln(abs(s)));
          if (se < 0) then if ((v and 1) = 1) then se:=-se;
          puiss0:=se;
          end;
       end;
function polynome_fits( px,py : pointer; pf : pointer ; dg,ndat : byte;
                      var ki2 : single ; _weight : pointer) : boolean;
   procedure make_a_b;
   var i,j,k : integer;
      begin
      if _weight=nil then
         begin
         for i:=0 to dg do
            begin
            b_mat[i]:=0;
            for k:=1 to ndat do
               b_mat[i]:=b_mat[i] + xy(py^)[k] *
                         puiss0( xy(px^)[k],i);
            for j:=0 to dg do
               begin
               a_mat[i,j]:=0;
               for k:=1 to ndat do
                  a_mat[i,j]:=a_mat[i,j] +
                         puiss0( xy(px^)[k],i+j)
               end
            end
         end
      else
         begin
         for i:=0 to dg do
            begin
            b_mat[i]:=0;
            for k:=1 to ndat do
               begin
               b_mat[i]:=b_mat[i] +
                      xy(py^)[k] * puiss0( xy(px^)[k],i) *
                      xy(_weight^)[k];
               end;
            for j:=0 to dg do
               begin
               a_mat[i,j]:=0;
               for k:=1 to ndat do
                  begin
                  a_mat[i,j]:=a_mat[i,j] +
                       puiss0(xy(px^)[k],i+j)  *
                       xy(_weight^)[k];
                  end;
               end
            end
         end
      end;
   var toto : boolean;
       ktest : integer;
   procedure permut (i,ns : integer);
   var j,k : integer;
      begin
      inc(ktest);
      brot:=b_mat[i];
      for j:=i to dg do arot[j]:=a_mat[i,j];
      for j:=i to dg-1 do
         begin
         for k:=i to dg do a_mat[j,k]:=a_mat[j+1,k];
         b_mat[j]:=b_mat[j+1]
         end;
       for k:=i to dg do
          a_mat[dg,k]:=arot[k];
       b_mat[dg]:=brot;
       toto:=false
       end;
   function triangle : boolean;
   var i,j,k : integer;
      begin
      triangle:=false;
      for i:=0 to dg-1 do
         begin
         ktest:=0;
            repeat
            toto:=true;
            if ktest > dg then exit;
            if a_mat[i,i]=0 then permut(i,dg)
            until toto;
         for j:=i+1 to dg do piv[j]:=a_mat[i,j]/a_mat[i,i];
         for j:=i+1 to dg do
            begin
            for k:=i to dg do
               a_mat[j,k]:=a_mat[j,k]-piv[j]*a_mat[i,k];
            b_mat[j]:=b_mat[j]-piv[j]*b_mat[i]
            end
         end;
      triangle:=true
      end;
   procedure solu;
   var u   : extended;
       i,j : integer;
      begin
      for i:=dg downto 0 do
         begin
         u:=b_mat[i];
         if i < dg then
            for j:=i+1 to dg do u:=u-a_mat[i,j] * b0(pf^)[j];
         u:=u/a_mat[i,i];
         b0(pf^)[i]:=u;
         end
      end;
   var s,n : extended;
       i,j : integer;
   begin
   polynome_fits:=false;
   if dg > deg_mx  then exit;
   if ndat <= dg   then exit;
   make_a_b;
   if not triangle then exit;
   solu;
   ki2:=0;
   n:=0;
   for i:=1 to ndat do
      begin
      s:=0;
      for j:=0 to dg do
         s := s + puiss0(xy(px^)[i],j) * b0(pf^)[j];
      n:=n + sqr(xy(py^)[i]);
      ki2:=ki2 + sqr(xy(py^)[i]-s)
      end;
   if n > 1e-100 then
      ki2:=sqrt(ki2/n)*100
   else
      ki2:=1e35;
   polynome_fits:=true
   end;
end.

[User]

Bonjour j.p.,

Je vais essayer de digérer tout ca,

D'ici là,
1 Grand


@ +

[j.p.mignot]

Pour compléter mon message ci-dessus, il convient de dire que cela n'est qu'1 cadre de base. Une application + correcte doit
1- associer à chque paramètre une booléene qui le libère ou non lors du fit. Il convient en effet de fitter en 1er les paramètres significatifs et liberer les autre que lorsque l'on approche de la solution
2- prévoir de pouvoir pondérer les points
3- prévoir de limiter le domaine où les paramètres sont acceptables.

Tous ses ajouts sont très faciles à implémenter sur le unit Less_2 et peuvent faire toute la différence entre un soft qui fonctionne + ou - et un soft robuste ( en regard de divergences, d'over/unser flow, ...).


HELP! pour mieux encapsuler ce unit, je prefairerais ne pas avoir à programmer " Ma_fonction_a_fitter( a : one_col_param; x: single ) : single; " dans le unit mais la laisser en externe (mode $f+)et définir dans l'implémentation du unit un pointer sur la fonction afin de l'executer depuis l'interne lors du fit.

N'etant pas informatitien, j'ai un petit problème de syntaxe en regard du passage de paramètre. Si qq sait comment faire cela en pascal ou en C cela serait le bien venu.

L'idée serait d'encapsuler définitivement de unit dans une classe et d'y acceder que par des pointers de fonction et de de tableaux.

Merci!

[User]

Bonsoir,

J'ai pu tester le moindre carré 'standard',
en copiant tes unit dans 1 projet Delphi 7

Le résultat est impressionnant !



J'ai comparé les paramètres de départ (série1) avec ceux obtenus pour la série2 et visionner les courbes dans le Graph ca colle toujours bien:

Pour Une_fonction_a_fitter du style: polynômes, Exp, Log, Sin*Expo, Cos*Exp etc.. le résultat est parfait.

Il y a juste un bémol pour les fonctions trigo. sin, cos ou la ca colle moins bien mais tu dois certainement avoir une explication.

Mon problème maintenant va être d'incorporer ca sous VB (Excel, Access ...)

J'ai commencé à traduire ton 1er algo fits-Polynomial..
j'ai juste 1 problème au niveau de la traduction des pointeurs ..

Ceci dit encore bravo !!!

Quel bol de tomber sur 1 expert comme toi !

L'idée serait d'encapsuler définitivement de unit dans une classe et d'y acceder que par des pointers de fonction et de de tableaux.

Bonne idée, je me suis trop éloigner de delphi, il faudrait que je reprenne mon ancien code sous Delphi ou j'avais crée des classes,
peut-être avec un 1 truc du genre:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
interface
 
Uses Windows, Classes,SysUtils,ExtCtrls, StdCtrls, Dialogs;
 
Type
 
  TFit=Class(TObject)
  Private
 
  px, py        : pointer; 
  FCoef_a_fitter : one_col_param ;  
  N_Cy         : integer;       
  Residu       : single;        
 
  procedure SetCoef_a_fitter (); { méthode write }  
 
 
  public                                           { déclaration de la propriété }
  property Coef_a_fitter: one_col_param read FCoef_a_fitter write SetCoef_a_fitter;
 
  function Ma_fonction_a_fitter( a : one_col_param; x: single ) : single; 
 
  procedure FIT_ma_fonction( n_param_a_fitter, 
                           n_data           : integer; 
                           var Matrice_covariance : matrix_ma_ma); 
 
  constructor Create(AOwner: TComponent);
 
 
implementation 
 
function TFit.Ma_fonction_a_fitter( a : one_col_param; x: single ) : single;
 
....

et tu fais un truc du style:

f:=new TFit(Form1);


@ +

[j.p.mignot]

1- Il n'y a aucune raison que cos ne fonctionne pas! peut-être faut-il liberer la phase et fitter quelquechose du type A1 * cos( A2*x + A3) + A4.
2- Si vous translatez ces units en VB je vous serais recommaissant si vous m'en adressiez une copie ( jpmignot@freesurf.ch )
3- Je n'ai pas mis dans cette toute petite demo un affichage des erreurs sur les résultats: en general cela est toujours utile! de plus le nombre de cycles pour atteindre la convergence devrait aussi être analysé. En particulier si il vaut sa valeur max ( defini par une constante de l'interface ) alors on a pas convergé et si il vaut -1 alors le fit n'a pas eu lieu car il n'a pas pu être correctement initialisé
4- Pointers de fonction: je ne me suis pas correctement fais comprendre.
Actuellement la fonction à fitter est dans Less_2
je souhaite qu'elle soit écrite par l'utilisateur HORS de Less_2 avec le prototype correct ( (A ,x) ), que Less_2, la supprime de l'interface et mette à disposition un pointeurs, dison P_Fct, l'utilisateur devra charger dans P_Fct l'adresse de sa routine ( comme on le fait déjà pour l'adresse des tableaux de points ). Dans Less_2, je ne derai plus alors executer une fonction interne lors du fit, mais une fonction dont j'aurrai l'adresse via P_Fct..

donc sil'utilisateur écrit

function Ma_gaussienne(a : one_col_param; x: single ) : single;
begin
Ma_gaussienne:=a[1] * exp(-sqr(x/a[2]));
end;

puis dans ses initialisations

P_Fct := @Ma_Gaussienne;

comment dans Less_2 executer Ma_Gaussienne connaissant uniquement P_Fct?

Interet : => a- encapsulation Less_2
b- Si on a plusieurs types de fit, ne rien changer sur Less_2, mais uniquement changer la fonction en allant pointer sur une autre

[User]

Rebonjour,

J'ai traduit le 1er Algo,

Ca marche bien pour des fonctions polynômes dans xy mais sinon ca ne marche pas trop bien
j'utilise peut être mal les arguments weight, et le résidu Ki2 comment le prendre en compte dans le calcul ?

Voici ma fonction test:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
Public Function Test_polynome_fits()
Dim i As Integer
Dim Result As Boolean
Dim u As Double
Dim x As Double
Dim ki2 As Single
 
For i = 1 To 600
xy(1, i) = -1 + 0.035 * i   ' j'ai du supprimer les pointeurs
xy(2, i) = Sin(xy(1, i))  ' je teste avec sin 
Next i
 
Result = polynome_fits(20, 600, ki2, 0) ' (dg,ndat,Ki2,weight)
 
Debug.Print ("---------------------------------")
 
x = 4  ' argument de la fonction
 
u = 0
 
For i = 0 To 20
u = u + b0(i) * (x ^ i)   ' évaluation de la valeur du polynome de degré dg avec les coef b0 (me suis-je trompé ?)
Next i
 
Debug.Print ("Valeur de la fonction:")
Debug.Print (Sin(x))  ' affiche la valeur de la fonction
 
Debug.Print ("Valeur du polynôme:")
Debug.Print (u) ' affiche la valeur du polynôme
 
End Function

J'ai fait une traduction simple: ca bug peut-être ???
comment utiliser le résidu dans la fonction ???:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
Function polynome_fits(dg As Byte, ndat As Integer, ki2 As Single, weight As Integer) As Boolean
Dim toto As Boolean
Dim ktest As Integer
Dim s As Double, n As Double
Dim i As Integer, j As Byte
 
polynome_fits = False
 
   If dg > deg_mx Then Exit Function
   If ndat <= dg Then Exit Function
   make_a_b dg, ndat, 0
   If Not triangle(dg) Then Exit Function
   solu (dg)
   ki2 = 0
   n = 0
   For i = 1 To ndat
   s = 0
      For j = 0 To dg
      s = s + puiss0(xy(1, i), j) * b0(j)
      Next j
   n = n + (xy(2, i) * xy(2, i))
   ki2 = ki2 + (xy(2, i) - s) * (xy(2, i) - s)
   Next i
 
   If n > 1E-100 Then
      ki2 = Sqr(ki2 / n) * 100
   Else
      ki2 = 1E+35
   polynome_fits = True
   End If
 
End Function

Désolé pour mon amateurisme


Merci pour vos réponses

@+

[j.p.mignot]

code C mi à jour + loin

[User]

J'ai pu tester !
Ce qui m'impressionne dans ce code c'est la rapidité alliée à la justesse ,

[User]

Bonsoir,

Commentaire sur le dernier code de j.p.mignot:

Comme indiqué par j.p.mignot:

En fait cela montre une des difficultés de moinde carré : Il faut partir de valeurs relativement proche de la solution pour eviter de diverger. Il ne faut pas les fixer 'au pifomètrre' mais réellement les estimer à partir de féquences moyeene, de pente et ou valeurs moyenne, de constantes de temps sur une enveloppe, ...

Partant de la:
Tout le problème est donc de faire avant tout une bonne évaluations des coefficients de départ, que la proposition de départ soit la plus proche de la courbe recherchée.

Donc il ne suffit pas de donner à l'algorythme l'expression général de la fonction (ex:f=a1*(a2*Exp(x)+a3), il faut aussi lui donner des coefficients de départ proche des coefficients à trouver ce qui n'est pas une mince affaire.

Je propose donc de complèter le code de j.p.mignot de la facon suivante:

1. saisie des mesures (x1,yi) dans 1 bd ou dans une plage Excel

2. L'utilisateurs visionne ces points de mesures seuls (courbe en rouge)

3. L'utilisateurs cherche une fonction proche:
3.1. il estime la forme générale de la fonction (ex:f=a1*(a2*Exp(x)+a3) + a4,)
3.2. puis il essaye successivement différentes valeurs de coefficient
et à chaque essaie il peux comparer la courbe rouge (points de mesures) à la courbe bleu de la fonction essai, jusqu'a obtenir une courbe assez proche.

4. L'utilisateur applique en dernier l'algorythme de j.p.m. pour affiner définitivement la courbe (coube en vert)


après ces 4 points la courbe finale en vert devrait passer définitivement par les points de mesures !

Vos critiques ou remarque à 1 telle méthode sont les bien venus !


@+

[j.p.mignot]

Il est bien evidant que le succes d'1 moindre carrés réside en
1- fournir 1 point de départ suffisemment proche de la solution pour ne pas se trapper dans des minimums locaux parfois fort loin de la réalité
2- Estimer en 1er les paramètres significatifs et libérer les autres plus tard.

En fonction du modèle il peut y avoit plusieurs moyens de chercher ses points de départ. Par contre J'ECARTERAI une méthode par visualisation. Un soft DOIT se suffire à lui même.

Prenons l'exemple qui vous interesse
Y = a * exp(b*x+c) + d

si on lisse cette courbe ( FFT, fit poly, ...) on peut alors peut-être estimer

Y' = ab * exp(b*x+c)

soit log(y') = ln(ab) + b*x + c. si Y' <0 , <=> a*b <0, imultiplier les data par -1 au préalable. si Y' change de signe => le modèle est faux( l'exp est monotome) ou la méthode d'estimarion de Y' l'est)
Un fit linéaire de ces point (par la 1ere méthode que je vous ai fourni => Fit polynômoal parmis n points, choisir ici degré = 1 ) fixe alors

les valeures de départ pour
b et ln(|ab|) + c

En estimant les rayons de courbure ( (1+y'^2) / y'' ) on en

déduit une approximation de c

puiqu'il est minimum pour

X = -( ln(b^2) + 2c ) / ( 2*b )

dans le cadre de ce modèle et que l'on a déjà une estimation de b

On connait à présent ln(|ab|)+c, b et c

on en déduit un point de départ pour a

( son signe devra être ajusté en fonction d'une éventuelle multipllication par -1 du tableau Y' avant ln)

et finalement un

point de départ pour d

en prenant la moyenne des Yi - a*exp(bXi+c)

Votre remarque "apres cela la courbe devrait passer definitivement par les points experimentaux"

FAUX : elle passe probablement par AUCUN point! par contre elle ajute globalement au mieux le jeux de points

[User]

Votre remarque "apres cela la courbe devrait passer definitivement par les points experimentaux"

FAUX : elle passe probablement par AUCUN point! par contre elle ajute globalement au mieux le jeux de points

Il va sans dire, il s'agit d'1 ajustement de points !!!

en tous cas je vais étudier vos commentaires !

[User]

Bonjour,

J'ai pu tester le dernier code sous Delphi de j.p.mignot :
Un seul mot:
Bravo pour les progrès apportés depuis la 1ère version !!!

Ca marche pour pratiquement toutes les fonctions sin,cos, exp, Polynôme etc... l'ajustement est excellent.

La convergence est atteinte autour des 23 cycles

Sauf biensur pour le log si l'argument est négatif (ce qui peut arriver dans certain cas avec les coef)

il m'indique "OPERATION EN VIRGULE FLOTTANTE INCORRECTE !"
ce qui est normal !

Par contre j'ai testé pour une fonction du style:
Ma_fonction_a_fitter1:=a[1]*exp(a[2]*x)*Sin(a[3]*x) + a[4];

et par moment il m'indique "DEBORDEMENT EN VIRGULE FLOTTANTE !"

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Sur Covar[jb,kb] := Covar[jb,kb] + wt * derivee_partielle_y_da[kb];

Pour cette même fonction:
assez rarement j'ai aussi le message "ERROR GAUSS_JORDAN=0 !"
et aussi rarement "CONVERGENCE NON ATTEINTE APRES N CYCLES !"

EN TOUS CAS CETTE VERSION SE RAPPROCHE DU BUT FIXE !!!

@+

[j.p.mignot]

pour le log il n'y a pas de probleme: il suffit dans la definition de la fonction à fitter d'ecire if ( x < 1e-50 ) then FCT:=ln(1e-130) else ...

Vous avez pu voir que le programme ouvre les matrice[1..$ff] et [1..$ff,1..$ff] pour les coef, les poids, covar. or il est assez rare d'avoir plus que quelques paramètres ( a l'exeption de quelques problèmes de physique que l'affinement de structures cristallographique )

J'ai donc une version qui ne travaille que par pointeur et alloue la mémoire nécessaire et rien d'autre pour toutes ses matrices.

un probleme de underflow/overflow peut parfois se resoudtr en changeant spread. En Pascal il est facile d'utiliser le type extended dans gauss_jordan pour resoudre se probleme. Ce type est extrement precis et couvre des floating point sur une etendue nettement + grande que le type real. Il a le defaut d'être peu exportable, mais tant qu'on l'utilise en interne pas de problème. Il ralentit aussi les opérations: Il faut donc l'utiliser que dans les routines qui en ont réellement besoin. Je pense essayer de l'implémenter dans la prochaine mouture que je metterai à disposition. Son translatage en VB ne sera peut-être pas si aisé...

Je suis toujours intéressé à votre transcription en VB. - ne serai-ce que pour ma culture -


Pour des fonctions du type
Ma_fonction_a_fitter1:=a[1]*exp(a[2]*x)*Sin(a[3]*x) + a[4];
il se peut aussi que l'exponemtielle soit la raison d'1 under/over flow

il seriat plutôt prferable d'ecire

function Ma_fonction_a_fitter1(a,x) : real;
var arg : real;
begin
arg:=a[2]*x;
if ( arg < -25 ) then Ma_fonction_a_fitter1:=a[4] else
if arg > 50 ) then Ma_fonction_a_fitter1:=a[1]*exp(50)*Sin(a[3]*x) + a[4]
else
Ma_fonction_a_fitter1:=a[1]*exp(a[2]*x)*Sin(a[3]*x) + a[4];
end;

seuils de l'expo a controler !!

cela stabilise Ma_fonction_a_fitter1 mais peut induire des problèmes sur les derivées car alors la fonction est continue mais plus dérivable - surtout au niveau de arg=50- on peut alors imaginer pour eviter une divergence du systeme, écrire une fonction qui va arrondir cet angle.

[User]

Bonsoir,

Voici quelques copies d'écran des démos d'ajustements de points de mesures par la méthode des moindre carrés, en partant d'une proposition de départ:

Dans ce cas on part d'1 modèle du type a*cos(bx+c)+d
puis on choisit les coefs de départ (proposition..)
puis l'algo trouve les coefficients a, b, c et d qui ajuste au mieux les points de mesures:

Sous Delphi:



Sous Excel:



Série1oints de mesures
Série2: Proposition de départ
Série3: courbe d'ajustement des points de mesures

Vos commentaires ou propositions sont les bien venus !

[j.p.mignot]

FIT MOINDRE CARRES
A toutes fins utiles, voici une version C du module. J'y ai de + implémenté un le calcul d'un facteur de confiance.


Header LESS2_C.H

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
 
 
 
 
//---------------------------------------------------------------------------
 
#ifndef LESS2_CH
#define LESS2_CH
// les codes de sortie possibles
#define Error_Code_pas_assez_de_data            0xff
#define Error_Code_Wrong_Nbr_of_Cycle           0xfe
#define Error_Code_Incertitudes_non_definie     0xfa
#define Error_Code_Fct_non_definie              0xf0
#define Error_Code_somme_des_poids_nulle        80
#define Error_Code_Poids_negatif_trouve         81
#define Error_Code_pas_de_convergence           1
#define Error_Code_GJ_Nulle                     2
#define Error_Code_NONE                         0
// constantes                  
#define max_allowed_cycles                      0xff  /* limité à 1 byte (U8) */
#define min_allowed_cycles                      2
#define max_param                               0xff   /* nombre max de paramètres du modele, doit etre <= $ff */
#define max_points                              0xffff /* nombre max de couples (x,y) <= $ffff ( word)     */
#define SQRT10                                  3.1622776601683793319988935444327 /* sqrt(10)  */
#define Shift_Ini                               1e-4
// types
#define U16                                     unsigned short
#define U8                                      unsigned char
#define Real                                    double   /* mettre long double si necessaire */
// macros
#define sqr(f)                                  (   (f) * (f)   )
#define Maxs(a, b)                              (((a) > (b)) ? (a) : (b))
#define Mins(a, b)                              (((a) < (b)) ? (a) : (b))
class LESS2_Classe
   {
   private:
   //Variable internes :
        U8                      n_par;          // nombre de paramètre(s) du modèle <=  max_param
        int                     SNR2;           // sizeof(Real) * sqr(n_param) peut etre > que 0xFFFF
        U16                     SNR,            // sizeof(Real) * (long)n_param < 0xFFFF
                                W_dat;          // nombre de couples (x,y) <= max_points, doit aussi être >= n_param+2
        Real                    MAG,            // 'gain' pour normaliser y entre 1 et 2 lors du calcul
                                Ki_Square,      // valeur courrante du résidu
                                YM,             // Offset pour normaliser y entre 1 et 2 lors du calcul
                                Norme,          // valeur courrante de l'expention sur les dérivées
                                sm_weigh,       // facteur de pondération
                                Shift,          // décalage
                                PCT_OK_IN,      // critère de convergenge sur l'évolution du résidu ( en % ) si 0 on laisse le moindre carrés finir. initialisé par FIT_ma_fonction
                                Norme_ini,      // valeur initiale de l'expention sur les dérivées ( 10^-3 à 10^-5 semble OK,initialisé par FIT_ma_fonction)
                                _p_,            // _p_,_q_, paramètre de la droite qui ajuste les écarts dY ( x) normalement on doit truver _p_ #0 et _q_ #0
                                _q_,            // "    "
                                _r_;            // répartition des écarts pour essayer d'estimer si 1 fit est acceptable ou non
        bool                    bottom_reached, // valeur min de spred atteinte
                                _2nd_deriv,     // on utilise [ ou non ] -dy* @2R@a@aj dans la construction de Hessian
                                converged;      // le fit a convergé
   // Fonctions internes
        bool    Solve_Systeme
                (
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres d'être ajustés. Les autres sont fixes, maintenus à leurs valeurs initiales pas necessairement nulle
                Real Try_delta_a[],             // tableau de n_param Real définissant les (petites) variations sur les paramètres libres du modèle
                Real Covar[],                   // tableau n_param X n_param de réels définissant la matrice de covariance @R/@ai*@R/@aj. Elle est symétrique
                bool Use_snd_deriv              // si true alors covar n'est + symetrique => scan complet necessaire
                );
        void Covar_Builder
                (
                Real Try_delta_a[],             // tableau de n_param Real définissant les (petites) variations sur les paramètres libres du modèle
                Real Tab_W[],                   // tableau de n_param Real définissant les poids de chaque points ( 1 partout en génértal ) (initialisé par FIT_ma_fonction)
                Real Essai_New_a[],             // tableau de n_param Real définissant le jeu de paramètres actuellement estimés
                Real Tab_X[],                   // Tableau des coordonnées X
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres
                Real derivee_partielle_y_da[],  // tableau de n_param Real définissant les dérivées partielles
                Real Tab_Y[],                   // tableau des ordonnés Y
                Real p_er0[],                   // Tableau de n_data Real pour mémoriser les écarts modèle - observé
                Real Covar[],                   // tableau n_param X n_param de réels définissant la matrice de covariance @R/@ai*@R/@aj
                Real True_Covar[],              // Covariance sans pondération
                bool Use_snd_deriv   ,          // permettre d'utiliser @2R/@ai@aj dans matrice hessian
                Real d2y_daidaj[]               // derivées secondes : Tableau n_paramXn_param de réels. Il n'est pas necessairement symétrique
                );
        void make_deriv
                (
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres
                Real X[],                       // Tableau des coordonnées X
                Real Essai_New_a[],             // tableau de n_param Real définissant le jeu de paramètres actuellement estimés
                Real derivee_partielle_y_da[],  // tableau de n_param Real définissant les dérivées partielles
                U16  w1                         // N° du point
                );
        void make_2nd_deriv
                (
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres
                Real Tab_X[],                   // Tableau des coordonnées X
                Real Essai_New_a[],             // tableau de n_param Real définissant le jeu de paramètres actuellement estimés
                Real d2y_daidaj[],              //  tableau de n_param X n_param de réels contenant en sortie @2R/@ai1@ai2
                U16  w1                         // N° du point
                );
        bool Essai_de_param_modifies
                (
                Real Essai_new_a[],             // tableau n_param Real des coefs à ajuster que l'on vient d'easser
                Real Coef_a_fitter[],           // tableau n_param Real des coefs à ajuster au meilleur fit actuel
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres
                Real Try_delta_a[],             // tableau de n_param Real définissant les (petites) variations sur les paramètres libres du modèle
                Real Tab_W[],                   // tableau de n_param Real définissant les poids de chaque points ( 1 partout en génértal ) (initialisé par FIT_ma_fonction)
                Real Essai_New_a[],             // tableau de n_param Real définissant le jeu de paramètres actuellement estimés
                Real Tab_X[],                   // Tableau des coordonnées X
                Real Tab_Y[],                   // tableau des ordonnés Y
                Real derivee_partielle_y_da[],  // tableau de n_param Real définissant les dérivées partielles
                Real p_er0[],                   // Tableau de n_data Real pour mémoriser les écarts actuels modèle - observé
                Real p_er1[],                   // Tableau de n_data Real pour mémoriser les écarts modèle - observé dans le meilleur fit
                Real Covar[],                   // tableau n_param X n_param de réels définissant la matrice de covariance @R/@ai*@R/@aj
                Real True_Covar[],              // matrice de covariance uniquement pondérée par les poids ( sans Shift ni [sigma Y^2]^-1, ... )
                Real Best_Covar[],              // matrice de covariance au meilleur fit
                U8   *N_Cy_b,                   // itération
                Real *Residu_du_fit,            // Résidu obtenu au meilleur fit
                Real d2y_daidaj[]               // dérivées 2nd,
                );
        U8      MLi                             // si true alor Covar contient des @2R@ai@aj
                (                               // non necessairement symétrique => la scanner
                U8 Li,                          // en entier
                bool b
                );
        bool    distribution_des_ecart
                (
                Real Tab_X[],                   // Tableau des coordonnées X
                Real Tab_Y[],                   // tableau des ordonnés Y
                Real Tab_W[],                   // tableau de n_param Real définissant les poids de chaque points ( 1 partout en génértal ) (initialisé par FIT_ma_fonction)
                Real p_er0[],                   // Tableau de n_data Real pour mémoriser les écarts actuels modèle - observé
                Real p_er1[]                    // Tableau de n_data Real pour mémoriser les écarts modèle - observé au meilleur fit
                );
        short   init_au
                (
                U8   np,                        // nombre de paramètres du modèle
                U16  ndata,                     // nombre de données pour lr fit
                Real Tab_Y[],                   // tableau des ordonnés Y
                Real Tab_W[],                   // tableau de n_param Real définissant les poids de chaque points ( 1 partout en génértal ) (initialisé par FIT_ma_fonction)
                Real Essai_new_a[],             // tableau de n_param Real définissant le jeu de paramètres actuellement estimés
                Real Coef_a_fitter[],           // tableau de n_param Real des valeurs initiales des param à fitter
                Real Try_delta_a[],             // tableau de n_param Real définissant les (petites) variations sur les paramètres libres du modèle
                Real Tab_X[],                   // Tableau des coordonnées X
                bool Free[],                    // tableau de n_param boolean definissant les paramètreslibres
                Real derivee_partielle_y_da[],  // tableau de n_param Real définissant les dérivées partielles
                Real p_er0[],                   // tableau de n_data ecarts modèle-obs. Il est initialisé avec le modèle basésur les param proposé d'origine.
                Real Covar[],                   // tableau n_param X n_param de réels définissant la matrice de covariance @R/@ai*@R/@aj
                Real True_Covar[],              // matrice de covariance uniquement pondérée par les poids ( sans Shift ni [sigma Y^2]^-1, ... )
                U8   *N_Cy_b,                   // nombre de cycles necessaire pour converger initialisé à 0
                Real *Residu_du_fit,            // Résidu du fit initialisé à 1 reel arbitrairement grand
                Real d2y_daidaj[]               // derivées 2nd
                );
        Real    Calibre                         // normalise entre 1 et 2 les donnés Y proposées
                (
                Real y                          // pas fait "en bloc" comme pour UnCalibre car on a besoin
                );                              // de l'acceder point par point durant le fit
        void    UnCalibre
                (
                Real Tab_y[]                    // restitution des données en fin de calcul
                );
        bool    Compute_Incertitudes
                (
                Real    p_er1[],                // tableau de Real, contient dY
                Real    Covar[],                // Buffer pour triangulariser Best_Covar
                Real    Best_Covar[],           // Covariance non pondérée par Y ( mais par W tout de même ) au meilleur fit
                Real    Incertitude[],          // incertitudes en retour
                Real    Tab_w[],                // poids de chaque point
                bool    Free[]                  // parametres libre ( => où Covar est defini )
                );
   public:
   // variables de l'interface
        bool    DONE;                           // passe à TRUE quand le calcul est fini, se met à false lors de l'appel de FIT_ma_fonction
   // fonctions de l'interface
        short   FIT_ma_fonction                 // pour fitter le modèle
               (
                Real    Norme_ini_caller,       // 10^3 à 10^-5
                Real    PCT_OK,                 // en général 0 ( laisser finir le fit sans interruption si dR/R < x% )
                U8      n_params_a_fitter,      // nombre de paramètre(s) du modèle
                U8      max_cycle,              // limitation sur le nombre de cycle ( <= FF )
                U16     n_data_a_fitter,        // nombre de couples de points experimentaux ( >= n_param+2, <= FFFF )
                Real    Tab_x[],                // Tableau des coordonnées X
                Real    Tab_y[],                // tableau des ordonnés Y
                Real    Tab_w[],                // tableau de n_param Real définissant les poids de chaque points ( 1 partout en génértal )
                bool    Param_Free[],           // tableau de n_param boolean definissant les paramètreslibres du modèle
                Real    Coef_a_fitter[],        // valeurs initiales SENSEES de TOUS les paramètres - y compris ceux que l'on ne libère pas
                                                // au retour, valeurs affinées
                Real    Incertitudes[],         // au retour, Tableau de n_paran Real qui, en retour, contient les incertitudes
                Real    Best_Covar[],           // au retour, matrice de covariance au meilleur fit
                U8      *N_Cy_b,                // nombre de cycles effectivement utilisés
                Real    *Residu_du_fit,         // Résidu du fit
                Real    *Index_de_confience,    // index de confiance 0 à 1
                bool    allow_2nd_deriv         // permet d'utiliser la derivé 2nd dans la matrice Hessian
                );
   //---------------------------------------
       LESS2_Classe();                          //      Creation
       ~LESS2_Classe();                         //      destroy
   //---------------------------------------
 
   };
 
void            Init_LESS2         ( void );    // Initialisation de la classe
void            Kill_LESS2         ( void );    // Destruction de la classe
 
extern          LESS2_Classe      * LESS2_Device; // pour acceder à la classe
 
//------------ Less2_FCT : modèle à fitter----------------------------------------
// doit être fixé par l'appelant en definissant une fonction Real Fonc1(Real[],x)
// et en assignant Less2_FCT = Fonc1
// prototype des fonctions à fitter   Real MY_fit( Real a[], Real x)
extern Real    (*Less2_FCT )      ( Real a[], Real x );
 
//---------------------------------------------------------------------------
#endif

Code LESS2_C.CPP

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
 
 
//---------------------------------------------------------------------------
#pragma hdrstop
 
#include        "LESS2_C.h"
#include        "mem.h"
#include        <stdlib.h>
#include        <math.h>
 
//---------------------------------------------------------------------------
 
#pragma package(smart_init)
 
static  bool    NO_FUNC = true, LESS2_ON = false;
LESS2_Classe      * LESS2_Device;
 
Real    (*Less2_FCT )      ( Real a[] , Real x );   // prototype fonction à fitter
 
static Real Ma_fonction_non_definie( Real a[] , Real x )
   {
   NO_FUNC  =  true;  // pour dectecter si l'utilisateur a oublié de définir sa fonction
   return      x;     // sans interet si non eviter un warning du compilateur
   }
// création de la classe
void Init_LESS2 (void)
   {
   if ( ! LESS2_ON )
      {
      LESS2_ON          =         true;
      LESS2_Device      =         new (LESS2_Classe);
      if ( LESS2_Device == NULL ) exit(1);
      Less2_FCT         =         Ma_fonction_non_definie;  // pas de fonction à fitter
      }
   }
// suppression de la classe
void  Kill_LESS2 (void)
   {
   if ( LESS2_ON )
      {
      LESS2_ON  =       false;
      LESS2_Device->~LESS2_Classe();
      }
   }
LESS2_Classe::LESS2_Classe()    //      Creation
   {
   }
LESS2_Classe::~LESS2_Classe()   //      destroy
   {
   }
static void swap( Real a[], U16 i1, U16 i2, int n ) // pas U16 pour n car *sizeof(Real) peut etre > 0xFFFF
   { //  echange n Real entres les points de départ a[i1] et a[i2]
   Real *Buf = new Real[n];
   n    *=    sizeof(Real);
   memmove(Buf,&a[i1],n);
   memmove(&a[i1],&a[i2],n);
   memmove(&a[i2],Buf,n);
   delete [] Buf;
   }
//--- dérivées du modéle au point X[nw1] par rapport aux différents paramères libre
void LESS2_Classe::make_deriv(  bool Free[],
                                Real X[],
                                Real Essai_New_a[],
                                Real derivee_partielle_y_da[],
                                U16  w1 )    // N° du point
   {
   #define      epsilon         1e-8
   #define      Un_div_2epsi    5e7  /* 1 / 2 / epsilon   */
   for ( U8 i1 =0 ; i1 < n_par; i1++ ) if ( Free[i1] )
      {
      Real z            =       Essai_New_a[i1];
      Essai_New_a[i1]   =       z - epsilon;
      Real u1           =       Less2_FCT( Essai_New_a,X[w1] );
      Essai_New_a[i1]   =       z + epsilon;
      u1                =       (Less2_FCT(Essai_New_a,X[w1]) - u1) * Un_div_2epsi * MAG;
      derivee_partielle_y_da[i1]        =       u1;
      Essai_New_a[i1]   =       z;
      }
   }
void LESS2_Classe::make_2nd_deriv(
                                bool Free[],
                                Real Tab_X[],
                                Real Essai_New_a[],
                                Real d2y_daidaj[],
                                U16  w1)
   {
   for ( U8 i1 =0 ; i1 < n_par; i1++ ) if ( Free[i1] ) // ici matrice non symetrique donc
   for ( U8 i2 =0 ; i2 < n_par; i2++ ) if ( Free[i2] ) // scanner la matrice en entier
      {
      Real z1              =       Essai_New_a[i1];
      if ( i1 == i2 )
         {
         Essai_New_a[i1]   =       z1 - 2*epsilon;
         Real u1           =       Less2_FCT( Essai_New_a,Tab_X[w1] );
         Essai_New_a[i1]   =       z1 ;
         u1                =       (Less2_FCT(Essai_New_a,Tab_X[w1]) - u1) * Un_div_2epsi * MAG;  // @y/@a1 pour a1-epsilon
         Essai_New_a[i1]   =       z1 + 2*epsilon;
         Real u2           =       Less2_FCT( Essai_New_a,Tab_X[w1] );
         Essai_New_a[i1]   =       z1 ;
         u2                =       (Less2_FCT(Essai_New_a,Tab_X[w1]) - u2) * Un_div_2epsi * MAG;  // @y/@a1 pour a1+epsilon
         d2y_daidaj[(U16)i1*(U16)n_par+i2] = (u2-u1)* Un_div_2epsi ;    // @2y/@ai1^2
         }
      else
         {
         Real z2           =       Essai_New_a[i2];
         Essai_New_a[i2]   =       z2 - epsilon;
         Essai_New_a[i1]   =       z1 - epsilon;
         Real u1           =       Less2_FCT( Essai_New_a,Tab_X[w1] );
         Essai_New_a[i1]   =       z1 + epsilon;
         u1                =       (Less2_FCT(Essai_New_a,Tab_X[w1]) - u1) * Un_div_2epsi * MAG;  // @y/@ai1 pour ai2-epsilon
         Essai_New_a[i2]   =       z2 + epsilon;
         Essai_New_a[i1]   =       z1 - epsilon;
         Real u2           =       Less2_FCT( Essai_New_a,Tab_X[w1] );
         Essai_New_a[i1]   =       z1 + epsilon;
         u2                =       (Less2_FCT(Essai_New_a,Tab_X[w1]) - u2) * Un_div_2epsi * MAG;  // @y/@ai1 pour ai2+epsilon
         Essai_New_a[i1]   =       z1;
         Essai_New_a[i2]   =       z2;
         d2y_daidaj[(U16)i1*(U16)n_par+i2] = (u2-u1)* Un_div_2epsi ;    // @2y/@ai1@ai2
         }
      }
   }
#define Cov(i,j)        Covar[(U16)n_par*(U16)i+j]
#define Diag(i)         Covar[(U16)i*((U16)n_par+1)]
// calcul des modifications sur les paramètres pour essayer d'amelliorer le fit
U8 LESS2_Classe::MLi ( U8 Li, bool b)
   { // la matrice est symétrique ssi on néglige @2f/@aiaj
   if (b) return n_par-1;
   return Li;
   }
bool LESS2_Classe::Solve_Systeme(
        bool Free[],
        Real Try_delta_a[],
        Real Covar[],
        bool Use_snd_deriv
        )
   {
   bool Full_Matrix = Use_snd_deriv && _2nd_deriv;
   U16  li,co;  // pas U8 car * n_param pei monter à 0xFFFF
   for ( U8 ib1 = 0 ; ib1 < n_par ; ib1 ++ ) if ( Free[ib1] )
      {
      Real      Sg1       =  0;
      for ( U8 Li   =  0; Li  <   n_par  ; Li++  ) if ( Free[Li]  )
      for ( U8 Col  =  0; Col <=  MLi(Li,Full_Matrix) ; Col++ ) if ( Free[Col] )
      if ( fabs( Cov(Li,Col)) >= Sg1 )
         {   // la matrice est symétrique => chercher le max sur 1/2 matrice [ si !Full_Matrix]
         Sg1    =       fabs(Cov(Li,Col));
         li     =       Li;
         co     =       Col;
         }
      if ( Sg1 < 1e-305 )    // Real 5.0 x 10^-324 .. 1.7 x 10^308
         // toutes les derivee_partielle_y_da = 0 lors de la construction de covar dans Covar_Builder
         return false;
      if ( li != co ) // inverse les 2 lignes pour ramener sur diagonale
         {
         swap(Covar,(U16)li*(U16)n_par,(U16)co*(U16)n_par,n_par);
         swap(Try_delta_a,li,co,1);
         }
      // maitenant le + grand (en abs ) est sur la diagonale [co,co]
      Sg1               =       1 / Diag(co);
      Diag(co)          =       1;
      // on divise toute la ligne Covar + Try_delta_a par le max
      for ( U8 Col = 0; Col < n_par; Col++) if ( Free[Col] )
         Cov(co,Col)    *=      Sg1; // l'element (co,co) est l'inverse de son original
      Try_delta_a[co]   *=      Sg1;
      // maintenant la colonne a l'exeption de l'element de la ligne deja traitée
      for ( U8 Li = 0 ; Li < n_par; Li++) if ((Li != co ) && Free[Li] )
         { // scanne les ligne pour la colonne co
         Sg1            =       Cov(Li,co);     // sauve l'element
         Cov(Li,co)     =       0;       // puis l'annule
         for ( U8 Col = 0 ; Col < n_par; Col++) if ( Free[Col])
            Cov(Li,Col) -=      Cov(co,Col) * Sg1;
         Try_delta_a[Li]-=      Try_delta_a[co] * Sg1;
         }
      }
   return true;  // attention ici Covar a été détruit!  ( inversé )
   }
Real  LESS2_Classe::Calibre  ( Real y )
   {
   return  ( y - YM ) * MAG    +    1;    // signal normé entre 1 et 2
   }
void LESS2_Classe::UnCalibre     ( Real Tab_y[] )
   {
   for ( U16 w= 0 ; w < W_dat; w++) Tab_y[w] = ( Tab_y[w] - 1 ) / MAG + YM;    // restoration du signal normé entre 1 et 2
   }
// création de la matrice de covariance @R/@ai*@R/@aj pour les paramètres libres
void LESS2_Classe::Covar_Builder (
                              Real Try_delta_a[], Real Tab_W[],
                              Real Essai_New_a[], Real Tab_X[], bool Free[],
                              Real derivee_partielle_y_da[],
                              Real Tab_Y[],       Real p_er0[],
                              Real Covar[],       Real True_Covar[],
                              bool Use_snd_deriv,
                              Real d2y_daidaj[]
                                 )
   {
   bool Full_Matrix = Use_snd_deriv && _2nd_deriv;
   Real y_selon_model,wt,sig2i,dy,WG;
   memset(Covar,0,SNR2);
   memset(True_Covar,0,SNR2);
   memset(Try_delta_a,0,SNR);
   Ki_Square = 0;
   for ( U16 iw1 = 0 ; iw1 <  W_dat ; iw1++)
      {
      y_selon_model   =  Less2_FCT(Essai_New_a,Tab_X[iw1]);
      y_selon_model   =  Calibre(y_selon_model);
      make_deriv(Free,Tab_X,Essai_New_a,derivee_partielle_y_da,iw1);
      if ( Use_snd_deriv && _2nd_deriv  )
         make_2nd_deriv(Free,Tab_X,Essai_New_a,d2y_daidaj,iw1);
      else if ( Shift < 0 ) Shift = Shift_Ini;
      WG  = Tab_W[iw1];
      sig2i = WG / sqr( Tab_Y[iw1] * Norme ); // ici mes_data dans [1,2]
      dy = Tab_Y[iw1] - y_selon_model;
      for ( U8 li = 0 ; li < n_par; li++) if ( Free[li] )
         {
         wt = derivee_partielle_y_da[li] * sig2i;
         for ( U8 col = 0 ; col <= MLi(li,Full_Matrix); col++) // cacul 1/2 matrice car elle est symétrique
            {
            Cov(li,col) += ( wt * derivee_partielle_y_da[col]) ;
            if ( Full_Matrix )
               Cov(li,col) -= dy *  sig2i * d2y_daidaj[(U16)li*(U16)n_par + col];
            True_Covar[(U16)li*(U16)n_par + col] += WG*derivee_partielle_y_da[li]*derivee_partielle_y_da[col];
            }
         Try_delta_a[li] += wt * dy;
         }
      // alors conserver les écarts qui seront utiles pour l'estimateur et pour incertitudes
      p_er0[iw1] =  dy / MAG;   // div MAG pour avoir les ecarts après UnCalibre
      Ki_Square += sqr( dy ) * WG;
      }
   Ki_Square *= sm_weigh; //1/[sigma Pi*Yi^2] * Ki
   if ( ! Full_Matrix )
   for ( U8 li = 1; li < n_par; li++) for (U8 col = 0 ; col < li; col ++ )
      {
      Cov(col,li)      = Cov(li,col);
      True_Covar[(U16)col*(U16)n_par + li] = True_Covar[(U16)li*(U16)n_par + col];
      }
   }
// voyons si la nouvelle tentative sur les paramètres amene une amelliorartion du Ki^2
bool LESS2_Classe::Essai_de_param_modifies
                             (
                                Real Essai_new_a[],
                                Real Coef_a_fitter[],
                                bool Free[],
                                Real Try_delta_a[],
                                Real Tab_W[],
                                Real Essai_New_a[],
                                Real Tab_X[],
                                Real Tab_Y[],
                                Real derivee_partielle_y_da[],
                                Real p_er0[],
                                Real p_er1[],
                                Real Covar[],
                                Real True_Covar[],
                                Real Best_Covar[],
                                U8   *N_Cy_b,
                                Real *Residu_du_fit,
                                Real d2y_daidaj[]
                             )
   {
   bool OK;
   if ( (*N_Cy_b) == 0 )
      memmove(Essai_new_a,Coef_a_fitter,SNR );
   if ( Shift > 0 ) for ( U8 bj1 = 0 ; bj1 <  n_par ; bj1++) // initialise -1, continue avec -SQRT10^p tant
      Diag(bj1) *= (1 + Shift);                                 // que @2f/@ai@aj est permis
   if ( Solve_Systeme(Free,Try_delta_a,Covar,*N_Cy_b <=1) )
      {
      for ( U8 bj1 = 0 ; bj1 < n_par ; bj1++)
         Essai_new_a[bj1] +=  Try_delta_a[bj1];
      Covar_Builder( Try_delta_a, Tab_W, Essai_New_a,Tab_X,Free,
                derivee_partielle_y_da,Tab_Y, p_er0,Covar,True_Covar,*N_Cy_b <=1 ,d2y_daidaj);
      if ( Ki_Square < (*Residu_du_fit) )
         {
         if ( Shift < 1e10 )    // la multiplication n'a pas de rôle tant que Shift < 0
            Shift       *=  SQRT10;     // c'est peut-être 1 peu trop ???
         (*Residu_du_fit)  =  Ki_Square;
         bottom_reached =  false;
         memcpy(Coef_a_fitter,Essai_new_a,SNR);
         memcpy(Best_Covar,True_Covar,SNR2);
         converged = (*Residu_du_fit) < PCT_OK_IN;
         memmove(p_er1,p_er0,sizeof(Real)* (long)W_dat );
         }
      else
         {
         OK =  ( Norme > Norme_ini / 10 );
         if ( bottom_reached )
            {
            if ( OK ) // accroitre Norme de increment^2 car sera divive + loin par increment
               Norme /=  sqr( SQRT10 );
            else
               converged = true;  // si OK alors on ne peut + rien faire désolé !!!
            }
         }
      if ( Norme < 1 )
         Norme *= SQRT10; // 10^-3, 3 10^-3, 10^-2, .. , 1
      else
         bottom_reached =  true;
      return true;
      }
   else
      return false; //  false, matrice nulle, tous les derivee_partielle_y_da = 0
   }
bool LESS2_Classe::distribution_des_ecart
        (
        Real Tab_X[],
        Real Tab_Y[],
        Real Tab_W[],
        Real p_er0[],
        Real p_er1[]
        )
   {
   Real WG,s1,sx,sx2,sy,sxy,D;
   s1   =  0;
   sx   =  0;
   sx2  =  0;
   sy   =  0;
   sxy  =  0;
   for ( U16 iw2 =0 ; iw2 <  W_dat ; iw2++)
      {
      WG   =  Tab_W[iw2];
      s1  +=  WG;
      sx  +=  WG  * Tab_X[iw2];
      sx2 +=  WG  * sqr(Tab_X[iw2]);
      sy  +=  WG  * p_er1[iw2];
      sxy +=  WG  * Tab_X[iw2] * p_er1[iw2];
      }
   D  =  s1 * sx2 - sqr (sx);
   if ( fabs(D) > 1e-250 )    // idealement on doit trouver la droite y=0 => p=q=0
      {
      D =  1 / D;
      _p_ =  D * (sxy*s1 - sx*sy);
      _q_ =  D * (sx2*sy - sx*sxy);
      if ( W_dat > 20 ) // si non la distribution des ecarts n'est pas exploitable
         {
         for ( U16 iw2=0 ; iw2 < W_dat ; iw2++) // on met dans p_er[0] les ecart entre la droite et p_er[1]
            p_er0[iw2] = fabs( p_er1[iw2] - _p_ * Tab_X[iw2]  - _q_ );
         // p=q=0 n'est pas suffisant il faut encore que cela soit 'bien' réparti
         int i10  =  W_dat / 200;
         sxy = 0;
         for ( long iw2 =i10 ; iw2 < i10 + (long)W_dat ; iw2++)
            {
            s1 = 0;
            for ( long iw1 = iw2-i10; iw1 <= iw2+i10 ; iw1++)
               s1 += p_er0[iw1 % W_dat];
            p_er1[iw2-i10]= s1 / ( 2 * i10 + 1 );
            sxy += p_er1[iw2-i10];
            }
         sx = p_er1[0];
         sx2=sx;
         sy=Tab_Y[0];
         sxy=sy;
         for ( U16 iw2=1 ; iw2 < W_dat; iw2++)
            {
            sx  =       Maxs( sx , p_er1[iw2]);
            sy  =       Maxs( sy , Tab_Y[iw2]);
            sx2 =       Mins( sx2, p_er1[iw2]);
            sxy =       Mins( sxy, Tab_Y[iw2]);
            }
         _r_ = (sx-sx2)/Maxs(sy-sxy,1e-200)*100;
         }
      return true;
      }
   else
      return false;
   }
short LESS2_Classe::init_au
        (
        U8   np,
        U16  ndata,
        Real Tab_Y[],
        Real Tab_W[],
        Real Essai_new_a[],
        Real Coef_a_fitter[],
        Real Try_delta_a[],
        Real Tab_X[],
        bool Free[],
        Real derivee_partielle_y_da[],
        Real p_er0[],
        Real Covar[],
        Real True_Covar[],
        U8   *N_Cy_b,
        Real *Residu_du_fit,
        Real d2y_daidaj[]
        )
   {
   Real R;
   n_par             =       np;
   YM                   =       Tab_Y[0];
   MAG                  =       YM;
   for ( U16 w1= 1; w1 < ndata; w1++)
      {
      YM                =       Mins  (YM  , Tab_Y[w1]);
      MAG               =       Maxs  (MAG , Tab_Y[w1]);
      }
   MAG                  -=      YM;   // expension du signal
   MAG                  =       1.0 / Maxs( MAG, 1e-200 ); // limitation de principe
   sm_weigh             =       0;
   for ( U16 w1= 0 ; w1 < ndata; w1++)  // il faut calibrer le signal pour eviter des zero divide
      // de plus travailler avec des data de memes ordre de grandeur statbilise des
      // variations dues au codage
      {
      Tab_Y[w1]         =       Calibre(Tab_Y[w1]);
      if ( Tab_W[w1] >= 0 )
         sm_weigh       +=    Tab_W[w1] * sqr(Tab_Y[w1]);
      else
         return Error_Code_Poids_negatif_trouve;
      }
   if ( sm_weigh < 1e-300 )
      return Error_Code_somme_des_poids_nulle;
   sm_weigh             =       1 / sm_weigh;
   W_dat                =       ndata;
   converged            =       false;
   bottom_reached       =       false;
   Norme                =       Norme_ini;
   Shift                =       -1;    // passe à 10^-3 des que on arrete l'utilisation de @2R/@ai@aj
   (*Residu_du_fit)     =       1.7e308;  // Real = Double max - cf help -
   (*N_Cy_b)            =       0;
   memmove(Essai_new_a,Coef_a_fitter,SNR);
   Covar_Builder(Try_delta_a,Tab_W,Essai_new_a,Tab_X,Free,derivee_partielle_y_da,
            Tab_Y,p_er0,Covar,True_Covar,true,d2y_daidaj); // initialise Covar
   return Error_Code_NONE;
   }
bool LESS2_Classe::Compute_Incertitudes (
                                                Real    p_er1[],
                                                Real    Covar[],
                                                Real    Best_Covar[],
                                                Real    Incertitude[],
                                                Real    Tab_w[],
                                                bool    Free[]
                                        )
   {
   bool ok;
   Real R2;
   Real *Bj     =       new     Real    [n_par];
   bool *Fr     =       new     bool    [n_par];
   memcpy(Covar,Best_Covar,SNR2); // contient True_Covar  ( sigma p=0..W_data-1   Wp * @f/@ai *  @f/@aj ( x= xp) )
   memcpy(Fr,Free,sizeof(Fr));   // if faut preserver Free car permut peut le changer
   R2           =       0;
   for ( U16 i=0; i < W_dat; i++ )
      R2        +=      p_er1[i];
   for ( U8 i=0; i < n_par; i++ )
      Bj[i] = R2 * sqrt( Maxs(0,Diag(i))) ;  // |@f/@ai| * sigma(y). Maxs: on n'est jamais trop prudent avec les arrondis...
   ok = Solve_Systeme( Free,Bj, Covar, false ); // pour ce calcul pas de @2R@ai@aj
   memset(Incertitude,0,SNR);
   if ( ok )  // le det du systeme est non nul => on peut calculer la solution
      {
      for (U8 i = 0; i< n_par ;i++ ) if ( Free[i] )
         Incertitude[i] = Bj[i];
      R2        =       0;
      for ( U16 i = 0 ; i < W_dat ; i++ ) R2       +=      Tab_w[i];
      for ( U8  i = 0 ; i < n_par; i++ ) Incertitude[i] = fabs(Incertitude[i]) / R2;
      }
   delete [] Fr;
   delete [] Bj;
   return ok;
   }
short LESS2_Classe::FIT_ma_fonction
   (
   Real    Norme_ini_caller,
   Real    PCT_OK,
   U8      n_param_a_fitter,
   U8      max_cycle,
   U16     n_data_a_fitter,
   Real    Tab_x[],
   Real    Tab_y[],
   Real    Tab_w[],
   bool    Param_Free[],
   Real    Coef_a_fitter[],
   Real    Incertitudes[],
   Real    Best_Covar[],
   U8      *N_Cy_b,
   Real    *Residu_du_fit,
   Real    *Index_de_confience,
   bool    allow_2nd_deriv
   )
   {
   short LESS2_Error_Code;
   DONE            =  false;
   SNR             =  sizeof(Real) * (long)n_param_a_fitter;
   SNR2            =  (long)SNR    * (long)n_param_a_fitter;
   PCT_OK_IN       =  sqr(PCT_OK * 1e-2);
   _2nd_deriv      =  allow_2nd_deriv;
   Norme_ini       =  Norme_ini_caller;
   if  ( max_cycle < min_allowed_cycles )           return Error_Code_Wrong_Nbr_of_Cycle;
   if  ( n_data_a_fitter   <   n_param_a_fitter+2 ) return Error_Code_pas_assez_de_data;
   Real *Essai_new_a            =       new Real [n_param_a_fitter];
   Real *Try_delta_a            =       new Real [n_param_a_fitter];
   Real *derivee_partielle_y_da =       new Real [n_param_a_fitter];
   Real *p_er0                  =       new Real [n_data_a_fitter];
   Real *p_er1                  =       new Real [n_data_a_fitter];
   Real *Covar                  =       new Real [sqr(n_param_a_fitter)];
   Real *True_Covar             =       new Real [sqr(n_param_a_fitter)];
   Real *d2y_daidaj             =       new Real[sqr(n_param_a_fitter)];
   LESS2_Error_Code             =
      init_au
      (
      n_param_a_fitter,         n_data_a_fitter,        Tab_y,
      Tab_w,                    Essai_new_a,            Coef_a_fitter,
      Try_delta_a,              Tab_x,                  Param_Free,
      derivee_partielle_y_da,   p_er0,                  Covar,
      True_Covar,               N_Cy_b,                 Residu_du_fit,
      d2y_daidaj
      );
   if ( LESS2_Error_Code == Error_Code_NONE )
      {
      NO_FUNC = false;
      Less2_FCT(Coef_a_fitter,Tab_x[0]); // no_func repasse a true si
      if ( NO_FUNC ) // Less2_FCT n'est pas redefini
         LESS2_Error_Code =  Error_Code_Fct_non_definie;  // pas return car doit liberer la memoire
      else
         {
         #define        Loop_GoOn       -1
         LESS2_Error_Code  =  Loop_GoOn;
         while ( LESS2_Error_Code == Loop_GoOn )
            {
            if (
               Essai_de_param_modifies
                  (     Essai_new_a,    Coef_a_fitter,          Param_Free,
                        Try_delta_a,    Tab_w,Essai_new_a,      Tab_x,
                        Tab_y,          derivee_partielle_y_da, p_er0,
                        p_er1,          Covar,                  True_Covar,
                        Best_Covar,     N_Cy_b,                 Residu_du_fit,
                        d2y_daidaj
                  )
               )
               {
               if ( ! converged ) // on ameliore le Ki 2 toujours vrai au 1er passage
                  {
                  if  ( (*N_Cy_b) < max_cycle )
                     (*N_Cy_b)++;  // 1 au 1er passage
                  else
                     LESS2_Error_Code  =  Error_Code_pas_de_convergence;
                  }
               else // a ce cycle, Ki2 s'est mis a croitre => stopper
                  LESS2_Error_Code     =  Error_Code_NONE ;
               }
            else // probleme dans l'inversion de GAUSS car toutes des derivées partielles sont nulles
               LESS2_Error_Code        =  Error_Code_GJ_Nulle;
            }
         UnCalibre(Tab_y);  // il a falu calibrer le signal pour eviter des zero divide
         MAG    =       1;
         YM     =       0;
         if (!Compute_Incertitudes ( p_er1, Covar,Best_Covar, Incertitudes,Tab_w, Param_Free))
            LESS2_Error_Code    =       Error_Code_Incertitudes_non_definie;
         (*Residu_du_fit)       =       100 * sqrt(Maxs(0,(*Residu_du_fit))); // Maxs au cas ou...
         if (
            distribution_des_ecart ( Tab_x,Tab_y,Tab_w,p_er0,p_er1)
            )
            {
            _p_ =  abs(_p_);
            _q_ =  abs(_q_);
            if ( _p_ < 0.2 )
               (*Index_de_confience) = 1;
            else
               (*Index_de_confience) = exp( - sqr(( _p_-0.2)/(0.3-0.2)));
            if (_q_ > 0.6 )
               (*Index_de_confience) *= exp( - sqr(( _q_-0.6)/(2-0.6)));
            if ( _r_ > 8)
               (*Index_de_confience) *= exp( - sqr(( _r_-8)/(20-8)));
            }
         else
               (*Index_de_confience)  =  0;
         }
      }
   delete [] d2y_daidaj;
   delete [] True_Covar;
   delete [] Covar;
   delete [] p_er1;
   delete [] p_er0;
   delete [] derivee_partielle_y_da;
   delete [] Try_delta_a;
   delete [] Essai_new_a;
   DONE                 =  true;
   return          LESS2_Error_Code;
   }

Pour a 1voir démo de ce code, me contacter. J'ai 1 soft de demo mais on m'a fait remarquer à juste titre que cela devenait trop long à tout introduire ici.

[Matthieu Brucher]

C'est quand même bien long
En gros, ce que tu as, c'est juste un optimiseur avec une fonction qui est une somme de carrés, on peut donc utiliser toutes les astuces d'optimisation de ce type de fonction de coût, ce qui permettrait de faire un code plus réutilisable ?

[j.p.mignot]

1- oui c'est long et il aurrait été préférable de mettre un zip attaché option qui ne semble pas être à disposition.
2- non ce n'est pas un "optimisateur avec somme de carrés" il s'agit d'une méthode de moindre carrés basé sur du calcul variationnel et l'inversion de la matrice de gauss-jordan @2R/@ai@aj.
Cette technique est connue mais n'est pas forcement maitrisée par tous. Ce unit a déjà rendu un certain nombre de services.
3- la classe du unit a un interface limité et simple. Elle peut et a déjà été réutilisé dans maintes applications telle qu'elle - en tout cas dans sa version pascal ( delphi) -. J'en profite pour encore féliciter l'initiateur de ce topic qui en a fait une macro exel qui fonctionne étonnament ( compte tenu du cadre exel ! )
4- Le programme ne sert en fait à rien si non illustrer les résultats et donner un exemple d'implémentation.

[Matthieu Brucher]

Tu essaies tout de même de minimiser une somme de carrés, non ?
Tu as un lien sur l'optimisation que tu utilises ? Je suis en train de me renseigner sur ces bestiaux pour mon sujet de thèse.