/*
C code library for Neuron Network training
File : model0Trn.c
Ce code a été crée automatiquement par le module NEURO CODE à partir d'un
modèle neuronal développé avec un outil de la
Suite NEURO ONE
Le présent code source généré par le progiciel NEURO CODE est protégé tant
par les dispositions nationales qu'internationales en matière de droits de
la propriété intellectuelle, dont les droits sont détenus, à titre
exclusif, par la société NETRAL.
L'utilisation et la modification de ce code source est soumise à un contrat
de licence d'utilisation.
La contrefaçon est un délit pénal puni de 2 ans d'emprisonnement et de
150.000 Euros d'amende.
Le modèle du présent code source a fait l'objet d'un dépôt auprès de
l'Agence pour la Protection des Programmes sous le numéro :
IDDN.FR.001.500018.00.S.P.1999.000.20700
NETRAL ne peut en aucun cas être tenu pour responsable des conséquences de
l'utilisation de ce code.
Lisez attentivement le fichier "licfr.txt" joint à ce fichier pour
connaitre vos droits et obligations concernant l'usage de ce code.
Pour obtenir une licence ou tout renseignement complémentaire,
adressez vous à :
NETRAL
14, rue Verdi
9213 Issy-les-Moulineaux
tel : (33) 146 387 512
email: info@netral.com
Date : 25/10/2005
Time : 17:31:15
User : Jean_Luc_PLOIX
Computer : SOPHRONE
Counter : 2
NetworkName : model0
FileName : D:\Program Files\Netral\Data\Exemples\Modeles\Static.NML
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "model0tfr.h"
#include "model0grd.h"
#include "model0lev.h"
#include "model0trn.h"
/* Please check name of included header files */
#define MAXCOST 1e100 /*Coût maximum*/
#define LAMBDAINI 100 /*Initialisation du Lambda*/
#define STEPREDN 0.2 /*Réduction du pas d'apprentissage*/
#define MAXINV 1e24 /*Valeur max pour inversion de matrice*/
#define CTOL 1e-12 /*Tolérance de calcul*/
#define CTOLM 1e12 /*Inverse tolérance de calcul*/
#define RACC 1e-2 /*Seuil d'acceptation d'un rang de matrice.*/
#define LAMBDAMIN 1e-6 /*lambda minimum*/
#define REVERSELIM 100 /*cycle d'inversion*/
#define MAXWEIGHT 100 /*valeur max des poids*/
#define ITMAX 30 /*max des itérations de DVS*/
real* VData; /* Données d'apprentissage*/
real* VDataValid; /* Données de validsation */
real* VHII; /* Table des Hii */
long DataCount; /* Nombre de données d'apprentissage */
long ValidDataCount; /* Nombre de données de validation */
real CurCost; /* Coût d'apprentissage courant */
real GenCost[model0OUTPUTS]; /* Vecteur de coût de généralisation */
real TrainCost[model0OUTPUTS]; /* Coûts d'apprentissage par sorties */
real JacobRank[model0OUTPUTS]; /* Rang du Jacobien par sortie */
long WeightLimited; /* Limitation des poids */
long Initialized = 0;
long InternalData = 0;
long InternalDataValid = 0;
long NDataDim = model0INPUTS + 2*model0OUTPUTS;
long indix = 0;
long CompteAppr = 0;
real mu = 0;
real lambda = LAMBDAINI;
long IsModif = 0;
real* TableResidus;
real* Jacobien;
real model0gethii(long index)
{
return VHII[index];
};
real sign(real a, real b)
{
return ((b >= 0)? fabs(a): -fabs(a));
};
real hypot(real X, real Y)
{
real Xt = fabs(X);
real Yt = fabs(Y);
return (Xt > Yt)? Xt*sqrt(1 + (Yt/Xt)*(Yt/Xt)):((Yt == 0)? 0: Yt*sqrt(1 + (Xt/Yt)*(Xt/Yt)));
};
long DVS(real* Jacobian, real WVct[model0SYNAPSES], real
VMat[model0SYNAPSES][model0SYNAPSES], long* code)
/*Module : model0Trn
Method : DVS
Visibility : Public
Arguments : real* Jacobian -> Le Jacobien à décomposer
real WVct[Q] -> Vecteur des diagonales
real VMat[Q][Q] -> Matrice (Q,Q)
long* code -> code d'erreur
Return : None
Description: Calcule la décomposition en valeurs singulières du Jacobien. Le
Jacobien lui-même est conservé
*/
{
long N = DataCount;
real* WMat;
real VTmp[model0SYNAPSES];
long nm = 0;
long l = 0;
long i, j, jj, k, iterations, mnmin;
real z, y, x, scale, s, h, g, f, c, anorm;
long jp;
WMat = calloc(N*model0SYNAPSES, sizeof(real));
if (!WMat) return 1;
for (i = 0; i < N; i++)
for (j = 0; j < model0SYNAPSES; j++)
WMat[i*model0SYNAPSES + j] = Jacobian[i*model0SYNAPSES + j];
g = 0;
scale = 0;
anorm = 0;
for (i = 0; i < model0SYNAPSES; i++) {
l = i + 1;
VTmp[i] = scale*g;
g = 0;
s = 0;
scale = 0;
if (i < N) {
for (k = i; k < N; k++)
scale += fabs(WMat[k*model0SYNAPSES + i]);
if (scale) {
for (k = i; k < N; k++) {
WMat[k*model0SYNAPSES + i] = WMat[k*model0SYNAPSES + i]/scale;
s = s + WMat[k*model0SYNAPSES + i]*WMat[k*model0SYNAPSES + i];
};
f = WMat[i*model0SYNAPSES + i];
g = -sign(sqrt(s),f);
h = f * g - s;
WMat[i*model0SYNAPSES + i] = f-g;
for (j = l; j < model0SYNAPSES; j++) {
s = 0;
for (k = i; k < N; k++)
s += WMat[k*model0SYNAPSES + i]*WMat[k*model0SYNAPSES + j];
f = s/h;
for (k = i; k < N; k++)
WMat[k*model0SYNAPSES + j] = WMat[k*model0SYNAPSES + j]+ f*WMat[k*model0SYNAPSES + i];
};
for (k = i; k < N; k++)
WMat[k*model0SYNAPSES + i] = scale*WMat[k*model0SYNAPSES + i];
}
};
WVct[i] = scale*g;
g = 0;
s = 0;
scale = 0;
if ((i < N) && (i != model0SYNAPSES - 1)) {
for (k = l; k < model0SYNAPSES; k++)
scale += fabs(WMat[i*model0SYNAPSES + k]);
if (scale) {
for (k = l; k < model0SYNAPSES; k++) {
WMat[i*model0SYNAPSES + k] = WMat[i*model0SYNAPSES + k]/scale;
s += WMat[i*model0SYNAPSES + k] * WMat[i*model0SYNAPSES + k];
};
f = WMat[i*model0SYNAPSES + l];
g = -sign(sqrt(s),f);
h = f*g-s;
WMat[i*model0SYNAPSES + l] = f-g;
for (k = l; k < model0SYNAPSES; k++)
VTmp[k] = WMat[i*model0SYNAPSES + k] / h;
for (j = l; j < N; j++) {
s = 0;
for (k = l; k < model0SYNAPSES; k++)
s += WMat[j*model0SYNAPSES + k] * WMat[i*model0SYNAPSES + k];
for (k = l; k < model0SYNAPSES; k++)
WMat[j*model0SYNAPSES + k] += s * VTmp[k];
};
for (k = l; k < model0SYNAPSES; k++)
WMat[i*model0SYNAPSES + k] *= scale;
}
};
/* anorm = max(anorm, (fabs(WVct[i]) + fabs(VTmp[i])));*/
if (anorm < (fabs(WVct[i]) + fabs(VTmp[i])))
anorm = (fabs(WVct[i]) + fabs(VTmp[i]));
};
for (i = model0SYNAPSES - 1; i >= 0; i--) {
if (i < model0SYNAPSES - 1) {
if (g) {
for (j = l; j < model0SYNAPSES; j++)
VMat[j][i] = (WMat[i*model0SYNAPSES + j]/WMat[i*model0SYNAPSES + l])/g;
for (j = l; j < model0SYNAPSES; j++) {
s = 0;
for (k = l; k < model0SYNAPSES; k++)
s += WMat[i*model0SYNAPSES + k] * VMat[k][j];
for (k = l; k < model0SYNAPSES; k++)
VMat[k][j] += s * VMat[k][i];
}
};
for (j = l; j < model0SYNAPSES; j++) {
VMat[i][j] = 0;
VMat[j][i] = 0;
}
};
VMat[i][i] = 1;
g = VTmp[i];
l = i;
};
if (N < model0SYNAPSES)
mnmin = N - 1;
else
mnmin = model0SYNAPSES - 1;
for (i = mnmin; i >= 0; i--) {
l = i + 1;
g = WVct[i];
for (j = l; j < model0SYNAPSES; j++)
WMat[i*model0SYNAPSES + j] = 0;
if (g) {
g = 1/g;
for (j = l; j < model0SYNAPSES; j++) {
s = 0;
for (k = l; k < N; k++)
s += WMat[k*model0SYNAPSES + i] * WMat[k*model0SYNAPSES + j];
f = (s/WMat[i*model0SYNAPSES + i])*g;
for (k = i; k < N; k++)
WMat[k*model0SYNAPSES + j] += f * WMat[k*model0SYNAPSES + i];
};
for (j = i; j < N; j++)
WMat[j*model0SYNAPSES + i] *= g;
}
else for (j = i; j < N; j++)
WMat[j*model0SYNAPSES + i] = 0;
WMat[i*model0SYNAPSES + i] = WMat[i*model0SYNAPSES + i] + 1;
};
for (k = model0SYNAPSES - 1; k >= 0; k--) {
for (iterations = 0; iterations <= ITMAX; iterations++) {
jp = 0;
for (l = k; l >= 0; l--) {
nm = l-1;
if (fabs(VTmp[l]) + anorm == anorm) {
jp = 1;
break;
};
if ((nm >= 0) && (fabs(WVct[nm]) + anorm == anorm))
break;
};
if (jp == 0) {
c = 0;
s = 1;
for (i = l; i < k+1; i++) {
f = s*VTmp[i];
VTmp[i] *= c;
if (fabs(f) + anorm == anorm)
break;
g = WVct[i];
h = hypot(f,g);
WVct[i] = h;
h = 1/h;
c = (g*h);
s = -(f*h);
for (j = 0; j < N; j++)
{
y = WMat[j*model0SYNAPSES + nm];
z = WMat[j*model0SYNAPSES + i];
WMat[j*model0SYNAPSES + nm] = (y*c)+(z*s);
WMat[j*model0SYNAPSES + i] = -(y*s)+(z*c);
}
};
};
z = WVct[k];
if (l == k) {
if (z < 0) {
WVct[k] = -z;
for (j = 0; j < model0SYNAPSES; j++)
VMat[j][k] = -VMat[j][k];
};
break;
};
if (iterations == ITMAX) {
*code = *code | code_dvserror;
break;
};
x = WVct[l];
nm = k-1;
y = WVct[nm];
g = VTmp[nm];
h = VTmp[k];
f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g = hypot(f,1);
f = ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
c = 1;
s = 1;
for (j = l; j <= nm; j++) {
i = j+1;
g = VTmp[i];
y = WVct[i];
h = s*g;
g = c*g;
z = hypot(f,h);
VTmp[j] = z;
c = f/z;
s = h/z;
f = (x*c)+(g*s);
g = -(x*s)+(g*c);
h = y*s;
y = y*c;
for (jj = 0; jj < model0SYNAPSES; jj++) {
x = VMat[jj][j];
z = VMat[jj][i];
VMat[jj][j] = (x*c)+(z*s);
VMat[jj][i] = -(x*s)+(z*c);
};
z = hypot(f,h);
WVct[j] = z;
if (z) {
z = 1/z;
c = f*z;
s = h*z;
};
f = (c*g)+(s*y);
x = -(s*g)+(c*y);
for (jj = 0; jj < N; jj++) {
y = WMat[jj*model0SYNAPSES + j];
z = WMat[jj*model0SYNAPSES + i];
WMat[jj*model0SYNAPSES + j] = (y*c)+(z*s);
WMat[jj*model0SYNAPSES + i] = -(y*s)+(z*c);
}
};
VTmp[l] = 0;
VTmp[k] = f;
WVct[k] = x;
};
};
free(WMat);
return 0;
};
long tracejac(long OutRnk, real* Source, real WVct[model0SYNAPSES],
real VMat[model0SYNAPSES][model0SYNAPSES], long VRnk[model0SYNAPSES],
long RCnt)
/*Module : model0Trn
Method : tracejac
Visibility : Private
Arguments : long OutRnk -> rang de la sortie du réseau
real* Source -> Jacobien
real WVct[Q] -> vecteur des diagonales
real VMat[Q][Q] -> Matrice (Q,Q)
long VRnk[Q] -> rangs des valeurs singulières
long RCnt -> Nombre de valeurs singulières rangées
Return : long -> 0 si problème, sinon 1
Description: Calcule les Hii et la trace du hessien
*/
{
long i, j, k, m;
real Hii, CumulTrace, Cumul, temp;
long result = 1;
CumulTrace = 0;
for (i = 0; (i < DataCount) && result; i++)
{
Hii = 0;
for (j = 0; (j < RCnt); j++)
{
m = VRnk[j];
Cumul = 0;
for (k = 0; k < model0SYNAPSES; k++)
Cumul += VMat[k][m] * Source[i*model0SYNAPSES + k];
if (WVct[m])
{
temp = Cumul / WVct[m];
Hii += temp*temp;
};
};
VHII[OutRnk*DataCount + i] = Hii;
CumulTrace += Hii;
if ((Hii > 1) && (Hii <= (1 + CTOL))) /* pour les erreurs d'arrondi */
Hii = 1- CTOL;
/* vérifier que Hii doit être inferieur à 1 */
result = result && (Hii <= 1 + CTOL);
};
result = result && (fabs(CumulTrace - RCnt) < RACC);
if (! result)
for (i = 0; i < DataCount; VHII[OutRnk*DataCount + i++] = 0);
return result;
};
real fullcompute(real* ZMat, real* Jacobien, long OutRnk, long* code)
/*Module : model0Trn
Method : fullcompute
Visibility : Private
Arguments : real* ZMat -> matrice ZtZ inverse à calculer
real* Jacobien -> Jacobien
long OutRnk -> rang de la sortie
long* code -> code d'erreur
Return : real -> le coût de généralisation
Description: calcule la décomposition en valeurs singulières du Jacobien, classe
les valeurs singulières, calcule la matrice ZtZ inverse, et le coût de
généralisation
*/
{
long i, j, k, m;
real Cumul;
real cost = 0;
real result = 0;
real WVct[model0SYNAPSES];
real VMat[model0SYNAPSES][model0SYNAPSES];
long VRnk[model0SYNAPSES];
long RCnt;
long Vdty;
real temp;
for (i = 0; i < DataCount; i++)
{
temp = TableResidus[i*model0OUTPUTS + OutRnk];
cost += temp*temp;
};
if (DVS(Jacobien, WVct, VMat, code))
return 0;
VRnk[0] = 0;
RCnt = 1;
for (i = 1; i < model0SYNAPSES; i++)
if (fabs(WVct[i]) > fabs(WVct[VRnk[0]]))
/* VRnk[0] est l'indice de la plus grande des valeurs singulières */
VRnk[0] = i;
for (i = 0; i < model0SYNAPSES; i++)
if ((VRnk[0] != i) && (fabs(WVct[i]/WVct[VRnk[0]]) > CTOL))
{
for(j = 0; (j < RCnt) && (fabs(WVct[i]) < fabs(WVct[VRnk[j]])); j++);
RCnt++;
for (k = RCnt-1; k >= j + 1; k--)
VRnk[k] = VRnk[k - 1];
VRnk[j] = i;
};
for (i = 0; i < model0SYNAPSES; i++)
{
for (j = 0; j <= i; j++)
{
Cumul = 0;
for (k = 0; k < RCnt; k++)
{
m = VRnk[k];
Cumul += VMat[i][m] * VMat[j][m]/WVct[m]/WVct[m];
};
ZMat[j + i*(i+1)/2] = Cumul;
};
};
Vdty = tracejac(OutRnk, Jacobien, WVct, VMat, VRnk, RCnt);
GenCost[OutRnk] = 0;
if (Vdty)
{
for (i = 0; i < DataCount; i++)
{
temp = TableResidus[i*model0OUTPUTS + OutRnk]/(1 - VHII[OutRnk*DataCount + i]);
GenCost[OutRnk] += temp*temp;
};
result = GenCost[OutRnk];
GenCost[OutRnk] = result/DataCount;
};
return result;
};
real JCC(real* weights, long* code)
/*Module : model0Trn
Method : JCC
Visibility : Private
Arguments : real* weights -> les poids courants
long* code -> code d'erreur
Return : real
Description: Calcule le Jacobien et le coût local*/
{
long i, j, k;
long code0 = *code;
real Outputs[model0OUTPUTS];
real* VD;
real* VSD;
real* VG;
real* PCD;
real Temp;
real IniBack[model0OUTPUTS];
real result = 0;
for (i = 0; i < model0OUTPUTS; TrainCost[i++] = 0);
for (i = 0; i < model0OUTPUTS*DataCount; TableResidus[i++] = 0);
for (i = 0; i < model0OUTPUTS*DataCount*model0SYNAPSES; Jacobien[i++] = 0);
for (i = 0; (i < DataCount) && (*code == code0); i++)
{
VD = VData + i*NDataDim;
VSD = VD + model0INPUTS;
PCD = TableResidus + i*model0OUTPUTS;
for (j = 0; j < model0OUTPUTS; j++)
{
VG = Jacobien + (j* DataCount + i)*model0SYNAPSES;
for (k = 0; k < model0OUTPUTS; k++)
IniBack[k] = (j == k);
model0transfergradient(weights, VD, Outputs, VG, IniBack);
PCD[j] = Outputs[j] - VSD[j];
Temp = PCD[j]*PCD[j];
TrainCost[j] = TrainCost[j] + Temp;
result = result + Temp;
if (result > MAXCOST)
{
*code = *code | code_toohighcost;
result = 0;
break;
};
};
};
return result/DataCount;
};
real CalculCoutGradient(real* weights, real Hessien[model0SYNAPSES]
[model0SYNAPSES], real VGradient[model0SYNAPSES], long* code)
/*Module : model0Trn
Method : CalculCoutGradient
Visibility : Private
Arguments : real* weights -> les poids courants.
real Hessien[Q][Q] -> le hessien courant
real VGradient[Q] -> le Gradient
long* code -> code d'erreur
Return : real
Description: Assure une propagation de un pas dans le réseau et effectue
simultanement le calcul du coût et du gradient.*/
{
long i, j, k, l, m;
real Temp;
real* VG;
long code0 = *code;
real result;
for (i = 0; i < model0SYNAPSES; i++)
{
VGradient[i] = 0;
for (j = 0; j < model0SYNAPSES; Hessien[i][j++] = 0);
};
result = JCC(weights, code);
if (*code == code0)
{
for (i = 0; i < DataCount; i++)
{
for (k = model0OUTPUTS - 1; k >= 0; k--)
{
VG = Jacobien + (k*DataCount + i)*model0SYNAPSES;
for (j = 0; j < model0SYNAPSES; j++)
VGradient[j] = VGradient[j] + VG[j]*TableResidus[i*model0OUTPUTS + k];
for (l = 0; l < model0SYNAPSES; l++)
{
Hessien[l][l] = Hessien[l][l] + VG[l] * VG[l];
for (m = 0; m < l; m++)
{
Temp = VG[l] * VG[m];
Hessien[l][m] = Hessien[l][m] + Temp;
Hessien[m][l] = Hessien[m][l] + Temp;
};
};
};
};
};
return result;
};
long CalculAMat(real Hessien[model0SYNAPSES][model0SYNAPSES], real
AMat[model0SYNAPSES][model0SYNAPSES], real lambda)
/*Module : model0Trn
Method : CalculAMat
Visibility : Private
Arguments : real Hessien[Q][Q] -> hessien
real AMat[Q][Q] -> matrice R à calculer
real lambda -> paramètre lambda
Return : long ->
Description: Calcule la matrice A telle que A*Atrans == Hessien + lambda*[I].*/
{
long i, j, k;
long Compte = 0;
long Tourne = 1;
long result = 0;
real fg;
real ua;
real li = 0;
while (Tourne)
{
Tourne = 0;
for (i = 0; i < model0SYNAPSES; i++)
{
fg = Hessien[i][i] + lambda + li;
for (k = 0; k < i; k++)
fg = fg - AMat[i][k] * AMat[i][k];
if (fg >= 0)
{
fg = sqrt(fg);
AMat[i][i] = fg;
for (j = i + 1; j < model0SYNAPSES; j++)
{
ua = Hessien[i][j];
for (k = 0; k < i; k++)
ua = ua - AMat[i][k] * AMat[j][k];
AMat[j][i] = ua / fg;
};
}
else if (Compte < REVERSELIM)
{
Tourne = 0;
result = code_reverseerror;
break;
}
else
{
if (li == 0)
li = LAMBDAMIN;
else
li = sqrt(li);
Tourne = 1;
Compte++;
break;
};
};
};
return result;
};
real norm(real HessienInv[model0SYNAPSES][model0SYNAPSES], real
VGradient[model0SYNAPSES], real VTGradient[model0SYNAPSES], real
*NormeVec)
/*Module : model0Trn
Method : norm
Visibility : Private
Arguments : real HessienInv[Q][Q] -> hessien inverse
real VGradient[Q] -> gradient
real VTGradient[Q] -> gradient de descente à calculer
real* NormeVec -> Norme du gradient
Return : real -> Norme du gradient de descente
Description: Calcule le vecteur de descente, transformé du gradient par le
hessien inverse, ainsi que la norme de ce vecteur.*/
{
real Temp, Elmt;
long i, j;
real result = 0;
*NormeVec = 0;
for (i = 0; i < model0SYNAPSES; i++)
{
Elmt = VGradient[i];
*NormeVec = *NormeVec + Elmt*Elmt;
Temp = 0;
for (j = 0; j < model0SYNAPSES; j++)
Temp += HessienInv[i][j] * VGradient[j];
VTGradient[i] = Temp;
result += Temp * Elmt;
};
return result;
};
long ResetInverse(real HessienInv[model0SYNAPSES][model0SYNAPSES], real
Hessien[model0SYNAPSES][model0SYNAPSES], real lambda)
/*Module : model0Trn
Method : ResetInverse
Visibility : Private
Arguments : real HessienInv[Q][Q] -> Hessien inverse à calculer
real Hessien[Q][Q] -> Hessien en entrée
real lambda -> constante lambda
Return : long
Description: Calcule l'inverse de [Hessien] + lambda*[I], et met le résultat
dans HessienInv.*/
{
long i, j, k;
long test;
real temp;
real AMat[model0SYNAPSES][model0SYNAPSES];
real BMat[model0SYNAPSES][model0SYNAPSES];
real DML[model0SYNAPSES];
long result = 0;
for (i = 0; i < model0SYNAPSES; i++)
for (j = 0; j < model0SYNAPSES; j++)
{
AMat[i][j] = 0;
BMat[i][j] = 0;
HessienInv[i][j] = 0;
};
test = lambda > 0;
i = 0;
while (test && (i < model0SYNAPSES))
{
test = lambda > fabs(Hessien[i][i]) * MAXINV;
i++;
};
if (test)
for (i = 0; i < model0SYNAPSES; i++)
{
HessienInv[i][i] = (1 - Hessien[i][i]/lambda)/lambda;
for (j = 1; j < model0SYNAPSES; j++)
HessienInv[i][j] = - Hessien[i][j];
}
else {
result = CalculAMat(Hessien, AMat, lambda);
if (result == 0)
{
for (i = 0; i < model0SYNAPSES; i++) DML[i] = AMat[i][i] ? 1.0/AMat[i][i]:0;
for (j = 0; j < model0SYNAPSES; j++)
{
for (i = 0; i < model0SYNAPSES; i++)
{
if (i == j)
temp = 1;
else
temp = 0;
for (k = 0; k < i; k++)
temp -= AMat[i][k]*BMat[k][j];
BMat[i][j] = temp*DML[i];
};
};
for (j = 0; j < model0SYNAPSES; j++)
{
for (i = model0SYNAPSES - 1; i >= 0; i--)
{
temp = BMat[i][j];
for (k = i + 1; k < model0SYNAPSES; k++)
temp -= AMat[k][i] * HessienInv[k][j];
HessienInv[i][j] = temp*DML[i];
};
};
};
};
return result;
};
real model0getcost(long* code)
{
return model0getcostw(model0weights, code);
};
real model0getcostw(real* weights, long* code)
/*Module : model0Trn
Method : model0getcostw
Visibility : Public
Arguments : real* weights -> poids courants
long* code -> code d'erreur
Return : real
Description: calcul du coût sur l'ensemble d'apprentissage, avec les poids
fournis. Les données d'apprentissage doivent avoir été enregistrées*/
{
long i, j;
long code0 = *code;
real* VD;
real* VS;
real* VSD;
real* VC;
real temp;
real result;
for (i = 0; i < model0OUTPUTS; TrainCost[i++] = 0);
for (i = 0; i < model0OUTPUTS*DataCount; TableResidus[i++] = 0);
result = 0;
for(i = 0; (i < DataCount) && (*code == code0); i++)
{
VD = VData + i*NDataDim;
VSD = VD + model0INPUTS;
VS = VD + model0INPUTS + model0OUTPUTS;
VC = TableResidus + i*model0OUTPUTS;
model0transferw(weights, VD, VS);
for (j = 0; j < model0OUTPUTS; j++)
{
VC[j] = VS[j] - VSD[j];
temp = VC[j]*VC[j];
TrainCost[j] += temp;
result = result + temp;
if (result > MAXCOST)
{
*code = *code | code_toohighcost;
result = 0;
break;
};
};
};
return result/DataCount;
};
long model0getdata(long Index, real* DataLoc)
/*Module : model0Trn
Method : model0GetData
Visibility : Public
Arguments : long Index -> indice de la donnée dans le tableau
real* DataLoc -> pointeur sur tableau de real. Le tableau
pointé doit porter les entrées puis les sorties.
Return : long -> 0 si OK, sinon code_dataaccessdenied
Description: Copie les données désignées dans le tableau de données*/
{
long i;
if (InternalData && (Index >= 0) && (Index < DataCount))
{
for (i = 0; i < model0INPUTS + model0OUTPUTS; i++)
VData[Index*NDataDim + i] = DataLoc[i];
return 0;
}
else
return code_dataaccessdenied;
};
long model0getdatavalid(long Index, real* DataLoc)
/*Module : model0Trn
Method : model0getdatavalid
Visibility : Public
Arguments : long Index -> indice de la donnée dans le tableau
real* DataLoc -> pointeur sur tableau de real. Le tableau
pointé doit porter les entrées puis les sorties.
Return : long -> 0 si OK, sinon code_dataaccessdenied
Description: Copie les données désignées dans le tableau de données de
validation*/
{
long i;
if ((!VDataValid) && (Index >= 0) && (Index < ValidDataCount))
{
for (i = 0; i < model0INPUTS + model0OUTPUTS; i++)
VDataValid[Index*NDataDim + i] = DataLoc[i];
return 0;
}
else
return code_dataaccessdenied;
};
real model0getgencost(real* Homogen, long* code)
{
return model0getgencostw(model0weights, model0matrix, Homogen, code);
};
real model0getgencostw(real* weights, real* ZMat, real* Homogen, long* code)
/*Module : model0Trn
Method : model0GetGenCost
Visibility : Public
Arguments : real* weights -> poids courants
real* ZMat -> matrice ZtZ inverse courante
real* Homogen -> homogénéité des Hii
long* code -> code d'erreur
Return : real -> le cout de généralisation.
Description: Calcule le Coût de généralisation et l'homogénéité. */
{
long i, j;
long code0 = *code;
real outputs[model0OUTPUTS];
real leverage[model0OUTPUTS];
real* VD;
real* VSD;
real temp;
real result = 0;
*Homogen = 0;
if (DataCount > model0SYNAPSES)
{
for(i = 0; (i < DataCount) && (*code == code0); i++)
{
VD = VData + i*NDataDim;
VSD = VD + model0INPUTS;
model0transferlevW(weights, ZMat, VD, outputs, leverage);
for (j = 0; j < model0OUTPUTS; j++)
{
*Homogen = *Homogen + sqrt(leverage[j]);
temp = outputs[j] - VSD[j];
result += temp*temp/(1 - leverage[j]);
if (result > MAXCOST)
{
*code = *code | code_toohighcost;
result = 0;
*Homogen = 0;
break;
};
};
};
*Homogen = *Homogen/DataCount;
result = result/(DataCount - model0SYNAPSES);
};
return result;
};
real model0getvalidcost(long* code)
{
return model0getvalidcostw(model0weights, code);
};
real model0getvalidcostw(real* weights, long* code)
/*Module : model0Trn
Method : model0getvalidcost
Visibility : Public
Arguments : real* weights -> poids courants
long* code -> code d'erreur
Return : real
Description: calcul du coût sur l'ensemble de validation*/
{
long i, j;
real outputs[model0OUTPUTS];
real* VD;
real* VSD;
long code0 = *code;
real temp;
real result = 0;
if (!VDataValid)
for(i = 0; (i < ValidDataCount) && (*code == code0); i++)
{
VD = VDataValid + i*NDataDim;
VSD = VD + model0INPUTS;
model0transferw(weights, VD, outputs);
for (j = 0; j <model0OUTPUTS; j++)
{
temp = outputs[j] - VSD[j];
result += temp*temp;
if (result > MAXCOST)
{
*code = *code | code_toohighcost;
result = 0;
break;
};
};
};
return result/ValidDataCount;
};
void model0free()
{
if (Initialized)
{
free(VHII);
free(Jacobien);
free(TableResidus);
if (InternalData)
free(VData);
if (InternalDataValid)
free(VDataValid);
}
Initialized = 0;
}
long model0init(long ADataCount, long AValidDataCount, real* AVData,
real* AVDataValid)
/*Module : model0Trn
Method : model0Init
Visibility : Public
Arguments : long ADataCount -> Nombre d'exemples de l'ensemble d'apprentissage
long AValidDataCount -> Nombre d'exemples de l'ensemble de validation
real* AVData -> tableau des données d'apprentissage
real* AVDataValid -> tableau des données de validation
Return : None
Description: procédure d'initialisation. Le nombre de données d'apprentissage
est fixé par ADataCount, et le nombre de donné"es de validation par
AValidDataCount.
AVData est un pointeur sur le tableau de données d'apprentisage.
AVDataValid est un pointeur sur le tableau de données de validation.
Si AVData est non nul, le programme pointe VData sur le tableau fourni
crée un tableau interne de données
Cette procédure doit être appellée avant toute autre utilisation de ce code */
{
if (Initialized)
{
free(VHII);
free(Jacobien);
free(TableResidus);
if (InternalData)
free(VData);
if (InternalDataValid)
free(VDataValid);
};
Initialized = 1;
WeightLimited = 0;
DataCount = ADataCount;
ValidDataCount = AValidDataCount;
TableResidus = calloc(DataCount*model0OUTPUTS, sizeof(real));
if (!TableResidus) return 1;
Jacobien = calloc(DataCount*model0OUTPUTS*model0SYNAPSES, sizeof(real));
if (!Jacobien) return 1;
VHII = calloc(DataCount*model0OUTPUTS, sizeof(real));
if (!VHII) return 1;
InternalData = (!AVData);
InternalDataValid = (!AVDataValid);
if (InternalData)
{
VData = calloc(DataCount* NDataDim, sizeof(real));
if (!VData) return 1;
}
else
VData = AVData;
if ((InternalDataValid) && (ValidDataCount))
{
VDataValid = calloc(ValidDataCount* NDataDim, sizeof(real));
if (!VDataValid) return 1;
}
else
VDataValid = AVDataValid;
return 0;
};
void model0iniweights(real env)
{
model0iniweightsw(model0weights, env);
};
void model0iniweightsw(real* weights, real env)
/*Module : model0Trn
Method : model0iniweights
Visibility : Public
Arguments : real* weights -> les pids courants à initialiser
env: real -> le paramètre d'initialisation
Return : None
Description: Initialise les poids par un valeur aléatoire de densité de
probabilité uniforme entre -env et env.
note : Cette méthode d'initialisation utilise le générateur de nombre aléatoire
de l'implémentation utilisée. */
{
long i;
for (i = 0; i < model0SYNAPSES; weights[i++] = env*2*((real)rand()/RAND_MAX - 0.5));
model0setwinit(1);
};
real model0setleverages(long* code)
{
return model0setleveragesw(model0weights, model0matrix, code);
};
real model0setleveragesw(real* weights, real* ZMat, long* code)
/*Module : model0Trn
Method : model0setleverages
Visibility : Public
Arguments : real* weights -> les poids courants
real* ZMat -> la matrice ZtZ inverse à calculer
long* code -> code d'erreur
Return : real -> le coût de généralisation
Description: Calcule le Jacobien, la matrice (ZtZ)Inv, le rang du Jacobien,
puis les leviers. */
{
long i, j;
real result = 0;
for (i = 0; i < model0OUTPUTS; JacobRank[i++] = 0);
if (DataCount > model0SYNAPSES) {
JCC(weights, code);
for (i = model0OUTPUTS - 1; i >= 0; i--) {
if (*code == 0)
result += fullcompute(ZMat, (Jacobien + i*DataCount*model0SYNAPSES), i,
code);
JacobRank[i] = 0;
for (j = 0; j < DataCount; j++)
JacobRank[i] += VHII[i*DataCount + j];
};
};
return result;
};
long model0trainingstep(real* weights, long indixmax, long* fini, long* TrainEnd,
long* code)
/*Module : model0Trn
Method : model0trainingstep
Visibility : Public
Arguments : real* weights -> les poids courants
long* fini -> vrai: apprentissage terminé
long* code -> code d'erreur
Return : long -> la valeur de *fini
*/
{
long i, j, k;
long Compteur;
long code0 = *code;
long OKAccuracy = 1;
long OKweights = 1;
long FirstTry = 1;
long AccPoint = 0;
long AccuracyCount = 0;
real LocalCurrentCost;
real ProjGrad;
real NormeVg = 0;
real stp;
real newval;
real VMWeight[model0SYNAPSES];
real VGradient[model0SYNAPSES];
real VTGradient[model0SYNAPSES];
real MTrainCost[model0OUTPUTS];
real Hessien[model0SYNAPSES][model0SYNAPSES];
real HessienInv[model0SYNAPSES][model0SYNAPSES];
if (*fini == 0) {
for (i = 0; i < model0SYNAPSES; i++) {
VMWeight[i] = 0;
VTGradient[i] = 0;
};
CurCost = 0;
for (j = 0; j < model0OUTPUTS; CurCost += TrainCost[j++]);
CurCost = CurCost/DataCount;
for(i = 0; WeightLimited && OKweights && (i < model0SYNAPSES); i++)
OKweights = (fabs(weights[i]) < MAXWEIGHT);
if (OKweights == 0)
*TrainEnd = teHighWeight;
*fini = (OKweights == 0) || (indix >= indixmax);
};
if (*fini == 0) {
CompteAppr++;
mu = 1;
for (i = 0; i < model0SYNAPSES; i++)
VMWeight[i] = weights[i];
LocalCurrentCost = CalculCoutGradient(weights, Hessien, VGradient, code);
if (*code == code0) {
for (j = 0; j < model0OUTPUTS; j++)
MTrainCost[j] = TrainCost[j];
while (OKAccuracy && (AccPoint == 0) && (*code == code0)) {
*code = code0 | ResetInverse(HessienInv, Hessien, lambda);
if (*code != code0) break;
/*calcul du vecteur de descente VTGradient, de la projection
du gradient sur la direction de descente, et du module carré du Gradient*/
ProjGrad = norm(HessienInv, VGradient, VTGradient, &NormeVg);
/*la projection doit être positive pour correspondre à une direction de descente,
et le facteur multiplicatif du module du gradient transformé est limité*/
AccPoint = (ProjGrad >= 0) && (ProjGrad <= CTOLM * NormeVg);
if (AccPoint) {
/*correction du vecteur de poids par les composantes du vecteur de descente
pondérées par l'incrément MU. Les corrections sont comparées à la tolérance.
Si toutes les corrections sont plus petites que la tolérance, on met à faux
le booléen OKAccuracy*/
Compteur = 0;
for (i = 0; i < model0SYNAPSES; i++) {
stp = mu * VTGradient[i];
newval = VMWeight[i] - stp;
weights[i] = newval;
if (fabs(stp) < CTOL) Compteur++;
};
OKAccuracy = Compteur != model0SYNAPSES;
AccPoint = (model0getcostw(weights, code) < LocalCurrentCost) &&
(*code == code0);
};
OKAccuracy = OKAccuracy && (lambda <= CTOLM) && (lambda != 0);
if (AccPoint == 0) {
if (OKAccuracy)
lambda = lambda/STEPREDN;
else if (FirstTry && (lambda > 0)) {
lambda = LAMBDAINI;
FirstTry = 0;
OKAccuracy = 1;
};
};
};
if (AccPoint) {
if (lambda > LAMBDAMIN)
lambda = lambda * STEPREDN;
IsModif = 1;
indix++;
}
else {
/*sinon, on revient à l'ancien vecteur de coefficients*/
for (k = 0; k < model0SYNAPSES; k++)
weights[k] = VMWeight[k];
for (k = 0; k < model0OUTPUTS; k++)
TrainCost[k] = MTrainCost[k];
};
IsModif = IsModif || OKAccuracy;
if (OKAccuracy)
AccuracyCount = 0;
else {
AccuracyCount++;
if (lambda == 0)
AccuracyCount++;
*TrainEnd = teAccuracy;
};
*fini = *fini || (AccuracyCount >= 2) || (*code != code0);
if (*fini)
AccuracyCount = 0;
}
else if (lambda != 0)
lambda = LAMBDAINI;
};
if (*TrainEnd == teNone)
{
if (*code)
*TrainEnd = teError;
else if (indix == indixmax)
*TrainEnd = teEpoch;
};
*fini = *fini || (*TrainEnd != teNone);
return *fini;
};
long model0train(long* Nb, long* TrainEnd, real* cost)
{
return model0trainw(model0weights, Nb, TrainEnd, cost);
};
long model0trainw(real* weights, long* Nb, long* TrainEnd, real* cost)
/*Module : model0Trn
Method : model0train
Visibility : Public
Arguments : real* weights -> les poids courants
long* Nb -> nombre d'époques remandées/ réalisées
long* TrainEnd -> Motif de l'arrêt d'apprentisage.
real* Cost -> valeur du coût en fin d'apprentisage.
Return : 0 si Ok, sinon, code d'erreur.
Description: Lancement d'un apprentissage. Le nombre d'époques demandées est Nb.
en retour, Nb porte le nombre d'époques effectivement réalisées. Le motif de
l'arret est porté par la variable TrainEnd.
*/
{
long code = 0;
long fini = 0;
long indixmax = *Nb;
*TrainEnd = teNone;
indix = 0;
*cost = model0getcostw(weights, &code);
while ((indix < *Nb) && (!fini) && (!code))
model0trainingstep(weights, indixmax, &fini, TrainEnd, &code);
*Nb = indix;
*cost = CurCost;
return code;
};
NETRAL - Neuro Code 6