CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
mlp_gen.cc
Go to the documentation of this file.
1 #include <math.h>
2 #include <string.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 
6 #include "mlp_gen.h"
7 #include "mlp_sigmoide.h"
8 
9 #ifdef __cplusplus
10 extern "C" {
11 #endif
12 
13 #define NPMAX 100000
14 #define NLMAX 1000
15 
17 struct learn_ learn_ MLP_HIDDEN;
18 struct pat_ pat_ MLP_HIDDEN;
20 struct stat_ stat_ MLP_HIDDEN;
21 
22 int MessLang = 0;
23 int OutputWeights = 100;
25 int WeightsMemory = 0;
26 int PatMemory[2] = {0,0};
27 int BFGSMemory = 0;
29 int LearnMemory = 0;
30 int NetMemory = 0;
31 float MLPfitVersion = (float) 1.40;
34 
35 dbl ***dir;
42 
43 /* The following lines are needed to use the dgels routine from the LAPACK
44  library in Reslin() */
45 
46 #include "mlp_lapack.h"
47 /* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer *
48  nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
49  doublereal *work, integer *lwork, integer *info);
50 
51 /***********************************************************/
52 /* MLP_Out */
53 /* */
54 /* computes the output of the Neural Network */
55 /* inputs: double *rrin = inputs of the MLP */
56 /* outputs: double *rrout = outputs of the MLP */
57 /* */
58 /* Author: J.Schwindling 29-Mar-99 */
59 /* Modified: J.Schwindling 16-Jul-99 unrolled loops */
60 /* Modified: J.Schwindling 20-Jul-99 fast sigmoid */
61 /***********************************************************/
62 
63 /* extern "C"Dllexport */void MLP_Out(type_pat *rrin, dbl *rrout)
64 {
65 // static int i, il, in, j, ilm1, m, mp1;
66  static int i, il, in, j, m, mp1;
67  dbl **deriv1;
68 
69 /* input layer */
70 
71  deriv1 = NET.Deriv1;
72  m = NET.Nneur[0]%4;
73  if(m==0) goto L10;
74  for(j=0;j<m;j++) NET.Outn[0][j] = rrin[j];
75 L10:
76  mp1 = m+1;
77  for(i=mp1; i<=NET.Nneur[0]; i+=4)
78  {
79  NET.Outn[0][i-1] = rrin[i-1];
80  NET.Outn[0][i] = rrin[i];
81  NET.Outn[0][i+1] = rrin[i+1];
82  NET.Outn[0][i+2] = rrin[i+2];
83  }
84 
85 /* hidden and output layers */
86 
87  MLP_MatrixVectorBias(NET.vWeights[1],NET.Outn[0],
88  NET.Outn[1],NET.Nneur[1],NET.Nneur[0]);
89 
90  for(il=2; il<NET.Nlayer; il++)
91  {
92  MLP_vSigmoideDeriv(NET.Outn[il-1],
93  deriv1[il-1],NET.Nneur[il-1]);
94  MLP_MatrixVectorBias(NET.vWeights[il],NET.Outn[il-1],
95  NET.Outn[il],NET.Nneur[il],
96  NET.Nneur[il-1]);
97  }
98  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
99  {
100  deriv1[NET.Nlayer-1][in] = 1;
101  }
102 }
103 
104 
105 /***********************************************************/
106 /* MLP_Out_T */
107 /* */
108 /* computes the output of the Neural Network when called */
109 /* from MLP_Test */
110 /* inputs: double *rrin = inputs of the MLP */
111 /* */
112 /* Author: J.Schwindling 23-Jul-1999 */
113 /***********************************************************/
114 
115 /* extern "C"Dllexport */void MLP_Out_T(type_pat *rrin)
116 {
117  static int i, il, in, j, ilm1, m, mp1;
118  register dbl a;
119 
120 /* input layer */
121 
122  m = NET.Nneur[0]%4;
123  if(m==0) goto L10;
124  for(j=0;j<m;j++) NET.Outn[0][j] = rrin[j];
125 L10:
126  mp1 = m+1;
127  for(i=mp1; i<=NET.Nneur[0]; i+=4)
128  {
129  NET.Outn[0][i-1] = rrin[i-1];
130  NET.Outn[0][i] = rrin[i];
131  NET.Outn[0][i+1] = rrin[i+1];
132  NET.Outn[0][i+2] = rrin[i+2];
133  }
134 
135 /* hidden and output layers */
136 
137 /* for(in=0;in<NET.Nneur[0]; in++) printf("%e %e\n",
138  NET.Outn[0][in],NET.Weights[1][0][in]);
139  printf("\n"); */
140  for(il=1; il<NET.Nlayer; il++)
141  {
142  ilm1 = il-1;
143  m = NET.Nneur[ilm1]%4;
144  for(in=0; in<NET.Nneur[il]; in++)
145  {
146  a = NET.Weights[il][in][0];
147  if(m==0) goto L20;
148  for(j=1;j<=m;j++) a +=
149  NET.Weights[il][in][j]*NET.Outn[ilm1][j-1];
150 L20:
151  mp1 = m+1;
152  for(j=mp1; j<=NET.Nneur[ilm1]; j+=4)
153  {
154  a +=
155  NET.Weights[il][in][j+3]*NET.Outn[ilm1][j+2]+
156  NET.Weights[il][in][j+2]*NET.Outn[ilm1][j+1]+
157  NET.Weights[il][in][j+1]*NET.Outn[ilm1][j]+
158  NET.Weights[il][in][j]*NET.Outn[ilm1][j-1];
159  }
160  switch(NET.T_func[il][in])
161  {
162  case 2: NET.Outn[il][in] = MLP_Sigmoide(a);
163  break;
164  case 1: NET.Outn[il][in] = a;
165  break;
166  case 0: NET.Outn[il][in] = 0;
167  break;
168  }
169  }
170  }
171 }
172 
173 
174 /***********************************************************/
175 /* MLP_Out2 */
176 /* */
177 /* computes the output of the Neural Network */
178 /* inputs: double *rrin = inputs of the MLP */
179 /* outputs: double *rrout = outputs of the MLP */
180 /* */
181 /* Author: J.Schwindling 29-Mar-99 */
182 /* Modified: J.Schwindling 16-Jul-99 unrolled loops */
183 /* Modified: J.Schwindling 20-Jul-99 fast sigmoid */
184 /***********************************************************/
185 
186 /* extern "C"Dllexport */void MLP_Out2(type_pat *rrin)
187 {
188 // static int il, in, m, mp1, i0, ilm1;
189  static int il, in, m, mp1;
190  register int i;
191  dbl **rrout, **deriv1;
192  register dbl *prrout;
193  type_pat *prrin;
194  int nhid = NET.Nneur[1];
195  int nin = NET.Nneur[0];
196 
197  rrout = NET.Outn;
198  deriv1 = NET.Deriv1;
199 
200  m = NET.Nneur[0]%4;
201  if(m==0) goto L10;
202  if(m==1)
203  {
204  rrout[0][0] = rrin[1];
205  goto L10;
206  }
207  else if(m==2)
208  {
209  rrout[0][0] = rrin[1];
210  rrout[0][1] = rrin[2];
211  goto L10;
212  }
213  else if(m==3)
214  {
215  rrout[0][0] = rrin[1];
216  rrout[0][1] = rrin[2];
217  rrout[0][2] = rrin[3];
218  goto L10;
219  }
220 L10:
221  mp1 = m+1;
222  prrout = &(rrout[0][mp1]);
223  prrin = &(rrin[mp1+1]);
224  for(i=mp1; i<=NET.Nneur[0]; i+=4, prrout+=4, prrin+=4)
225  {
226  *(prrout-1) = *(prrin-1);
227  *prrout = *prrin;
228  *(prrout+1)= *(prrin+1);
229  *(prrout+2) = *(prrin+2);
230  }
231 
232 /* input layer */
233 
234  MLP_MatrixVectorBias(NET.vWeights[1],NET.Outn[0],
235  NET.Outn[1],nhid,nin);
236 
237 
238 /* hidden and output layers */
239 
240  for(il=2; il<NET.Nlayer; il++)
241  {
242  MLP_vSigmoideDeriv(NET.Outn[il-1],deriv1[il-1],NET.Nneur[il-1]);
243  MLP_MatrixVectorBias(NET.vWeights[il],NET.Outn[il-1],
244  NET.Outn[il],NET.Nneur[il],NET.Nneur[il-1]);
245  }
246  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
247  deriv1[NET.Nlayer-1][in] = 1;
248 }
249 
250 
251 /***********************************************************/
252 /* MLP_Test_MM */
253 /* */
254 /* computes the MLP error function using matrix-matrix */
255 /* multiplications */
256 /* inputs: int ifile = file number: 0=learn, 1=test */
257 /* dbl *tmp = a pointer to an array of size */
258 /* 2 x number of neurons in first */
259 /* hidden layer */
260 /* */
261 /* return value (dbl) = error value */
262 /* */
263 /* Author: J.Schwindling 25-Jan-2000 */
264 /***********************************************************/
265 
267 {
268  int ipat;
269  int npat = PAT.Npat[ifile];
270  int nhid = NET.Nneur[1];
271  int nin = NET.Nneur[0];
272  int jpat, j, il, ilm1, m, in, mp1;
273  dbl err, a, rrans;
274  dbl *pweights, *ptmp;
275 
276  err = 0;
277  for(ipat=0; ipat<npat-1; ipat+=2)
278  {
279  MLP_MM2rows(tmp, &(PAT.vRin[ifile][ipat*(nin+1)]),
280  NET.vWeights[1], 2, nhid, nin+1,
281  nin+1, nin+1);
282 
283  switch(NET.T_func[1][0])
284  {
285  case 2:
286  ptmp = &(tmp[0]);
287  MLP_vSigmoide(ptmp,2*nhid);
288  break;
289 
290  case 1:
291  break;
292 
293  case 0:
294  for(jpat=0; jpat<2; jpat++)
295  {
296  for(j=0; j<nhid; j++)
297  {
298  tmp[j+jpat*nhid] = 0;
299  }
300  }
301  break;
302  }
303 
304  for(jpat=0; jpat<2; jpat++)
305  {
306  for(in=0; in<nhid; in++)
307  {
308  NET.Outn[1][in] = tmp[jpat*nhid+in];
309  }
310  for(il=2; il<NET.Nlayer; il++)
311  {
312  ilm1 = il-1;
313  m = NET.Nneur[ilm1]%4;
314  for(in=0; in<NET.Nneur[il]; in++)
315  {
316  pweights = &(NET.Weights[il][in][0]);
317  a = *pweights;
318  pweights++;
319  if(m==0) goto L20;
320  for(j=1;j<=m;j++,pweights++) a +=
321  (*pweights)*NET.Outn[ilm1][j-1];
322 L20:
323  mp1 = m+1;
324  for(j=mp1; j<=NET.Nneur[ilm1];
325  j+=4, pweights+=4)
326  {
327  a +=
328  *(pweights+3)*NET.Outn[ilm1][j+2]+
329  *(pweights+2)*NET.Outn[ilm1][j+1]+
330  *(pweights+1)*NET.Outn[ilm1][j]+
331  *(pweights )*NET.Outn[ilm1][j-1];
332  }
333  switch(NET.T_func[il][in])
334  {
335  case 2: NET.Outn[il][in] = MLP_Sigmoide(a);
336  break;
337  case 1: NET.Outn[il][in] = a;
338  break;
339  case 0: NET.Outn[il][in] = 0;
340  break;
341  }
342  }
343  if(il == NET.Nlayer-1)
344  {
345  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
346  {
347  rrans = (dbl) PAT.Rans[ifile][ipat+jpat][in];
348  err += (rrans-NET.Outn[NET.Nlayer-1][in])*
349  (rrans-NET.Outn[NET.Nlayer-1][in])*
350  PAT.Pond[ifile][ipat+jpat];
351  }
352  }
353  }
354  }
355  }
356 
357 /* cas npat impair */
358  for(ipat=ipat; ipat<npat; ipat++)
359  {
360  MLP_MatrixVector(NET.vWeights[1],
361  &(PAT.vRin[ifile][ipat*(nin+1)]),tmp,
362  nhid,nin+1);
363 
364  switch(NET.T_func[1][0])
365  {
366  case 2:
367  ptmp = &(tmp[0]);
368  MLP_vSigmoide(ptmp,2*nhid);
369  break;
370 
371  case 1:
372  break;
373 
374  case 0:
375  for(j=0; j<nhid; j++)
376  {
377  tmp[j] = 0;
378  }
379  break;
380  }
381 
382  for(in=0; in<nhid; in++)
383  {
384  NET.Outn[1][in] = tmp[in];
385  }
386  for(il=2; il<NET.Nlayer; il++)
387  {
388  ilm1 = il-1;
389  m = NET.Nneur[ilm1]%4;
390  for(in=0; in<NET.Nneur[il]; in++)
391  {
392  pweights = &(NET.Weights[il][in][0]);
393  a = *pweights;
394  pweights++;
395  if(m==0) goto L25;
396  for(j=1;j<=m;j++,pweights++) a +=
397  (*pweights)*NET.Outn[ilm1][j-1];
398 L25:
399  mp1 = m+1;
400  for(j=mp1; j<=NET.Nneur[ilm1];
401  j+=4, pweights+=4)
402  {
403  a +=
404  *(pweights+3)*NET.Outn[ilm1][j+2]+
405  *(pweights+2)*NET.Outn[ilm1][j+1]+
406  *(pweights+1)*NET.Outn[ilm1][j]+
407  *(pweights )*NET.Outn[ilm1][j-1];
408  }
409  switch(NET.T_func[il][in])
410  {
411  case 2: NET.Outn[il][in] = MLP_Sigmoide(a);
412  break;
413  case 1: NET.Outn[il][in] = a;
414  break;
415  case 0: NET.Outn[il][in] = 0;
416  break;
417  }
418  }
419  if(il == NET.Nlayer-1)
420  {
421  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
422  {
423  rrans = (dbl) PAT.Rans[ifile][ipat][in];
424  err += (rrans-NET.Outn[NET.Nlayer-1][in])*
425  (rrans-NET.Outn[NET.Nlayer-1][in])*
426  PAT.Pond[ifile][ipat];
427  }
428  }
429  }
430  }
431  return(err);
432 }
433 
434 
435 /***********************************************************/
436 /* MLP_Test */
437 /* */
438 /* computes the MLP error function */
439 /* inputs: int ifile = file number: 0=learn, 1=test */
440 /* int regul = 0: no regularisation term */
441 /* 1: regularisation term */
442 /* (for hybrid learning method) */
443 /* */
444 /* return value (dbl) = error value */
445 /* */
446 /* Author: J.Schwindling 31-Mar-99 */
447 /***********************************************************/
448 
449 /* extern "C"Dllexport */dbl MLP_Test(int ifile,int regul)
450 {
451  dbl err, rrans;
452  int in,jn,ipat,ipati;
453 
454  dbl *tmp;
455 
456  tmp = (dbl *) malloc(2 * NET.Nneur[1] * sizeof(dbl));
457  if(tmp == 0) /* not enough memory */
458  {
459  printf("not enough memory in MLP_Test\n");
460  err = 0;
461  for(ipat=0; ipat<PAT.Npat[ifile]; ipat++)
462  {
463  if(ifile==0)
464  {
465  ipati = ExamplesIndex[ipat];
466  }
467  else
468  {
469  ipati = ipat;
470  }
471  MLP_Out_T(PAT.Rin[ifile][ipati]);
472  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
473  {
474  rrans = (dbl) PAT.Rans[ifile][ipati][in];
475  err += (rrans-NET.Outn[NET.Nlayer-1][in])*
476  (rrans-NET.Outn[NET.Nlayer-1][in])*
477  PAT.Pond[ifile][ipati];
478  }
479  }
480 
481  if(regul>=1)
482  {
483  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
484  for(jn=0; jn<=NET.Nneur[NET.Nlayer-2]; jn++)
485  {
486  err += LEARN.Alambda*NET.Weights[NET.Nlayer-1][in][jn]*
487  NET.Weights[NET.Nlayer-1][in][jn];
488  }
489  }
490  free(tmp);
491  return(err);
492  }
493  else /* computation using matrix - matrix multiply */
494  {
495  err = MLP_Test_MM(ifile, tmp);
496  if(regul>=1)
497  {
498  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
499  for(jn=0; jn<=NET.Nneur[NET.Nlayer-2]; jn++)
500  {
501  err += LEARN.Alambda*NET.Weights[NET.Nlayer-1][in][jn]*
502  NET.Weights[NET.Nlayer-1][in][jn];
503  }
504  }
505  free(tmp);
506  return(err);
507  }
508 }
509 
510 
511 /***********************************************************/
512 /* MLP_Stochastic */
513 /* */
514 /* one epoch of MLP stochastic training */
515 /* (optimized for speed) */
516 /* */
517 /* Author: J.Schwindling 08-Jul-99 */
518 /* Modified: J.Schwindling 20-Jul-99 remove unused cases */
519 /***********************************************************/
520 
521 /* extern "C"Dllexport */dbl MLP_Stochastic()
522 {
523  int ipat, ii, inm1;
524  dbl err = 0;
525  int il, in1, in, itest2;
526  dbl deriv, deriv1, deriv2, deriv3, deriv4, pond;
527  dbl eta, eps, epseta;
528  register dbl a, b, dd, a1, a2, a3, a4;
529  dbl *pout, *pdelta, *pw1, *pw2, *pw3, *pw4;
530  dbl ***weights;
531 
532  if(NET.Debug>=5) printf(" Entry MLP_Stochastic\n");
533  weights = NET.Weights;
534 /* shuffle patterns */
535  ShuffleExamples(PAT.Npat[0],ExamplesIndex);
536 
537 /* reduce learning parameter */
538  if(LEARN.Decay<1) EtaDecay();
539 
540  eta = -LEARN.eta;
541  eps = LEARN.epsilon;
542  epseta = eps/eta;
543 
544 /* loop on the examples */
545  for(ipat=0;ipat<PAT.Npat[0];ipat++)
546  {
547  ii = ExamplesIndex[ipat];
548  pond = PAT.Pond[0][ii];
549 
550  MLP_Out2(&(PAT.vRin[0][ii*(NET.Nneur[0]+1)]));
551 
552 /* next lines are equivalent to DeDwSum */
553  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
554  {
555  deriv = NET.Deriv1[NET.Nlayer-1][in];
556  a = (dbl) PAT.Rans[0][ii][in];
557  b = NET.Outn[NET.Nlayer-1][in]-a;
558  err += b*b*pond;
559  NET.Delta[NET.Nlayer-1][in] = b*deriv*pond*eta;
560  }
561 
562  for(il=NET.Nlayer-2; il>0; il--)
563  {
564  dd = NET.Delta[il+1][0];
565  for(in=0; in<NET.Nneur[il]-3; in+=4)
566  {
567  deriv1 = NET.Deriv1[il][in];
568  deriv2 = NET.Deriv1[il][in+1];
569  deriv3 = NET.Deriv1[il][in+2];
570  deriv4 = NET.Deriv1[il][in+3];
571  itest2 = (NET.Nneur[il+1]==1);
572  a1 = dd*weights[il+1][0][in+1];
573  a2 = dd*weights[il+1][0][in+2];
574  a3 = dd*weights[il+1][0][in+3];
575  a4 = dd*weights[il+1][0][in+4];
576  if(itest2) goto L1;
577  pdelta = &(NET.Delta[il+1][1]);
578  for(in1=1; in1<NET.Nneur[il+1];
579  in1++, pdelta++)
580  {
581  a1 += *pdelta * weights[il+1][in1][in+1];
582  a2 += *pdelta * weights[il+1][in1][in+2];
583  a3 += *pdelta * weights[il+1][in1][in+3];
584  a4 += *pdelta * weights[il+1][in1][in+4];
585  }
586 L1: NET.Delta[il][in] = a1*deriv1;
587  NET.Delta[il][in+1] = a2*deriv2;
588  NET.Delta[il][in+2] = a3*deriv3;
589  NET.Delta[il][in+3] = a4*deriv4;
590  }
591  for(in=in; in<NET.Nneur[il]; in++)
592  {
593  deriv = NET.Deriv1[il][in];
594  itest2 = (NET.Nneur[il+1]==1);
595  a = dd*weights[il+1][0][in+1];
596  if(itest2) goto L2;
597  pdelta = &(NET.Delta[il+1][1]);
598  for(in1=1; in1<NET.Nneur[il+1];
599  in1++, pdelta++)
600  {
601  a += *pdelta *
602  weights[il+1][in1][in+1];
603  }
604 L2: NET.Delta[il][in] = a*deriv;
605  }
606 
607  } /* end of loop on layers */
608 
609 
610 /* update the weights */
611  if(eps==0)
612  {
613  for(il=1; il<NET.Nlayer; il++)
614  {
615  inm1 = NET.Nneur[il-1];
616  for(in=0; in<NET.Nneur[il]-3; in+=4)
617  {
618  a1 = NET.Delta[il][in];
619  a2 = NET.Delta[il][in+1];
620  a3 = NET.Delta[il][in+2];
621  a4 = NET.Delta[il][in+3];
622  pout = &(NET.Outn[il-1][0]);
623  weights[il][in][0] += a1;
624  weights[il][in+1][0] += a2;
625  weights[il][in+2][0] += a3;
626  weights[il][in+3][0] += a4;
627  weights[il][in][1] += a1* (*pout);
628  weights[il][in+1][1] += a2* (*pout);
629  weights[il][in+2][1] += a3* (*pout);
630  weights[il][in+3][1] += a4* (*pout);
631  pout++;
632  pw1 = &(weights[il][in][2]);
633  pw2 = &(weights[il][in+1][2]);
634  pw3 = &(weights[il][in+2][2]);
635  pw4 = &(weights[il][in+3][2]);
636  for(in1=2; in1<=inm1;
637  ++in1, ++pout, ++pw1, ++pw2,
638  ++pw3, ++pw4)
639  {
640  *pw1 += a1 * *pout;
641  *pw2 += a2 * *pout;
642  *pw3 += a3 * *pout;
643  *pw4 += a4 * *pout;
644  }
645  }
646  for(in=in; in<NET.Nneur[il]; in++)
647  {
648  a1 = NET.Delta[il][in];
649  pout = &(NET.Outn[il-1][0]);
650  weights[il][in][0] += a1;
651  weights[il][in][1] += a1* (*pout);
652  pout++;
653  pw1 = &(weights[il][in][2]);
654  for(in1=2; in1<=inm1;
655  ++in1, ++pout, ++pw1)
656  {
657  *pw1 += a1 * *pout;
658  }
659  }
660  }
661  }
662  else
663  {
664  for(il=1; il<NET.Nlayer; il++)
665  {
666  for(in=0; in<NET.Nneur[il]; in++)
667  {
668 
669  a = NET.Delta[il][in];
670  LEARN.Odw[il][in][0] = a + eps * LEARN.Odw[il][in][0];
671  NET.Weights[il][in][0] += LEARN.Odw[il][in][0];
672 
673  b = a*NET.Outn[il-1][0];
674  LEARN.Odw[il][in][1] = b + eps*LEARN.Odw[il][in][1];
675  NET.Weights[il][in][1] += LEARN.Odw[il][in][1];
676 
677  for(in1=2; in1<=NET.Nneur[il-1]; in1++)
678  {
679  b = a*NET.Outn[il-1][in1-1];
680  LEARN.Odw[il][in][in1] = b + eps*LEARN.Odw[il][in][in1];
681  NET.Weights[il][in][in1] += LEARN.Odw[il][in][in1];
682  }
683  }
684  }
685  }
686 
687  } /* end of loop on examples */
688  return(err);
689 }
690 
691 
692 /***********************************************************/
693 /* MLP_Epoch */
694 /* */
695 /* one epoch of MLP training */
696 /* inputs: int iepoch = epoch number */
697 /* dbl *alpmin = optimal alpha from Line Search*/
698 /* */
699 /* return value (dbl) = error value on learning sample */
700 /* BEFORE changing the weights */
701 /* */
702 /* Author: J.Schwindling 31-Mar-99 */
703 /* Modified: J.Schwindling 09-Apr-99 */
704 /* re-organize routine */
705 /* J.Schwindling 13-Apr-99 */
706 /* remove Quickprop algorithm */
707 /* implement Ribiere-Polak conjugate grad. */
708 /***********************************************************/
709 
710 /* extern "C"Dllexport */dbl MLP_Epoch(int iepoch, dbl *alpmin, int *Ntest)
711 {
712  dbl err, Norm, ONorm, beta, prod, ddir;
713 /* int *index;*/
714  int Nweights, Nlinear, ipat, ierr;
715  int nn;
716 
717  err = 0;
718  Norm = 0;
719  *alpmin = 0.;
720 
721  Nweights = NET.Nweights;
722  Nlinear = NET.Nneur[NET.Nlayer-2] + 1;
723 
724  if(NET.Debug>=5) printf(" Entry MLP_Epoch\n");
725 /* stochastic minimization */
726  if(LEARN.Meth==1)
727  {
728 
729  err = MLP_Stochastic();
730 
731  }
732  else
733  {
734  if(iepoch==1 && LEARN.Meth==7)
735  {
736  SetLambda(10000);
737  MLP_ResLin();
738  if(NET.Debug>=2) PrintWeights();
739  }
740 
741 /* save previous gradient and reset current one */
742  DeDwSaveZero();
743  if(LEARN.Meth==16)
744  {
745  ShuffleExamples(PAT.Npat[0],ExamplesIndex);
746  nn = PAT.Npat[0];
747  PAT.Npat[0] = nn/10;
748  for(ipat=0;ipat<nn;ipat++)
749  {
750  ierr = MLP_Train(&ExamplesIndex[ipat],&err);
751  if(ierr!=0) printf("Epoch: ierr= %d\n",ierr);
752  }
753  }
754  else
755  {
756  for(ipat=0;ipat<PAT.Npat[0];ipat++)
757  {
758  ierr = MLP_Train(&ipat,&err);
759  if(ierr!=0) printf("Epoch: ierr= %d\n",ierr);
760  }
761  }
762  DeDwScale(PAT.Npat[0]);
763  if(LEARN.Meth==2) StochStep();
764  if(LEARN.Meth==3)
765  {
766  SteepestDir();
767  if(LineSearch(alpmin,Ntest,err)==1) StochStep();
768  }
769 
770 /* Conjugate Gradients Ribiere - Polak */
771  if(LEARN.Meth==4)
772  {
773  if((iepoch-1)%LEARN.Nreset==0)
774  {
775  LEARN.Norm = DeDwNorm(); /* for next epoch */
776  SteepestDir();
777  }
778  else
779  {
780  ONorm = LEARN.Norm;
781  LEARN.Norm = DeDwNorm();
782  prod = DeDwProd();
783  beta = (LEARN.Norm-prod)/ONorm;
784  CGDir(beta);
785  }
786  if(LineSearch(alpmin,Ntest,err)==1) StochStep();
787  }
788 
789 /* Conjugate Gradients Fletcher - Reeves */
790  if(LEARN.Meth==5)
791  {
792  if((iepoch-1)%LEARN.Nreset==0)
793  {
794  LEARN.Norm = DeDwNorm(); /* for next epoch */
795  SteepestDir();
796  }
797  else
798  {
799  ONorm = LEARN.Norm;
800  LEARN.Norm = DeDwNorm();
801  beta = LEARN.Norm/ONorm;
802  CGDir(beta);
803  }
804  if(LineSearch(alpmin,Ntest,err)==1) StochStep();
805  }
806  if(LEARN.Meth==6)
807  {
808  if((iepoch-1)%LEARN.Nreset==0)
809  {
810  SteepestDir();
811  InitBFGSH(Nweights);
812  }
813  else
814  {
815  GetGammaDelta();
816  ierr = GetBFGSH(Nweights);
817  if(ierr)
818  {
819  SteepestDir();
820  InitBFGSH(Nweights);
821  }
822  else
823  {
824  BFGSdir(Nweights);
825  }
826  }
827  ddir = DerivDir();
828  if(ddir>0)
829  {
830  SteepestDir();
831  InitBFGSH(Nweights);
832  ddir = DerivDir();
833  }
834  if(LineSearch(alpmin,Ntest,err)==1)
835  {
836  InitBFGSH(Nweights);
837  SteepestDir();
838  if(LineSearch(alpmin,Ntest,err)==1)
839  {
840  printf("Line search fail \n");
841  }
842  }
843  }
844  if(LEARN.Meth==7)
845  {
846  if((iepoch-1)%LEARN.Nreset==0)
847  {
848  SteepestDir();
849  InitBFGSH(Nweights-Nlinear);
850  }
851  else
852  {
853  if(NET.Debug>=5) printf("Before GetGammaDelta \n");
854  GetGammaDelta();
855  if(NET.Debug>=5) printf("After GetGammaDelta \n");
856  ierr = GetBFGSH(Nweights-Nlinear);
857  if(NET.Debug>=5) printf("After GetBFGSH \n");
858  if(ierr)
859  {
860  SteepestDir();
861  InitBFGSH(Nweights-Nlinear);
862  }
863  else
864  {
865  BFGSdir(Nweights-Nlinear);
866  }
867  if(NET.Debug>=5) printf("After BFGSdir \n");
868  }
869  SetLambda(10000);
870  if(LineSearchHyb(alpmin,Ntest)==1)
871  {
872  InitBFGSH(Nweights-Nlinear);
873  SteepestDir();
874  if(LineSearchHyb(alpmin,Ntest)==1)
875  {
876  printf("Line search fail \n");
877  }
878  }
879  }
880  }
881 
882  if(NET.Debug>=5) printf(" End MLP_Epoch\n");
883  return(err);
884 }
885 
886 
887 /***********************************************************/
888 /* MLP_Train */
889 /* */
890 /* Train Network: compute output, update DeDw */
891 /* inputs: int *ipat = pointer to pattern number */
892 /* input/output: dbl *err = current error */
893 /* */
894 /* return value (int) = error value: 1 wrong pattern number*/
895 /* 2 *ipat < 0 */
896 /* */
897 /* Author: J.Schwindling 09-Apr-99 */
898 /***********************************************************/
899 
900 /* extern "C"Dllexport */int MLP_Train(int *ipat, dbl *err)
901 {
902  int in;
903 
904 /* if(*ipat>=PAT.Npat[0]) return(1);*/
905  if(*ipat<0) return(2);
906 
907 /* MLP_Out(PAT.Rin[0][*ipat],NET.Outn[NET.Nlayer-1]);*/
908  MLP_Out2(&(PAT.vRin[0][*ipat*(NET.Nneur[0]+1)]));
909  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
910  {
911  *err += ((dbl) PAT.Rans[0][*ipat][in]-NET.Outn[NET.Nlayer-1][in])
912  *((dbl) PAT.Rans[0][*ipat][in]-NET.Outn[NET.Nlayer-1][in])*
913  PAT.Pond[0][*ipat];
914  }
915  DeDwSum(PAT.Rans[0][*ipat],NET.Outn[NET.Nlayer-1],*ipat);
916  return(0);
917 }
918 
919 
920 /***********************************************************/
921 /* StochStepHyb */
922 /* */
923 /* Update the weights according to stochastic minimization */
924 /* formula (for hybrid methods) */
925 /* */
926 /* return value (int) = 0 */
927 /* */
928 /* Author: J.Schwindling 09-Apr-99 */
929 /***********************************************************/
930 
931 /* extern "C"Dllexport */int StochStepHyb()
932 {
933  int il, in1, in;
934  dbl eta, eps;
935 
936  eta = LEARN.eta;
937  eps = LEARN.epsilon;
938  for(il=NET.Nlayer-2; il>0; il--) {
939 
940  for(in=0; in<NET.Nneur[il]; in++) {
941 
942  /* compute delta weights */
943  for(in1=0; in1<=NET.Nneur[il-1]; in1++) {
944  LEARN.Odw[il][in][in1] = -eta * LEARN.DeDw[il][in][in1]
945  + eps * LEARN.Odw[il][in][in1];
946  }
947 
948  /* update weights */
949  for(in1=0; in1<=NET.Nneur[il-1]; in1++) {
950  NET.Weights[il][in][in1] += LEARN.Odw[il][in][in1];
951  }
952  }
953  }
954  MLP_ResLin();
955  return(0);
956 }
957 
958 
959 /***********************************************************/
960 /* StochStep */
961 /* */
962 /* Update the weights according to stochastic minimization */
963 /* formula */
964 /* */
965 /* return value (int) = 0 */
966 /* */
967 /* Author: J.Schwindling 09-Apr-99 */
968 /***********************************************************/
969 
970 /* extern "C"Dllexport */int StochStep()
971 {
972  int il, in1, in;
973  dbl eta, eps, epseta;
974 
975  eta = -LEARN.eta;
976  eps = LEARN.epsilon;
977  epseta = eps/eta;
978  for(il=NET.Nlayer-1; il>0; il--) {
979  for(in1=0; in1<=NET.Nneur[il-1]; in1++) {
980 
981  /* compute delta weights */
982  for(in=0; in<NET.Nneur[il]; in++) {
983  LEARN.Odw[il][in][in1] = eta * (LEARN.DeDw[il][in][in1]
984  + epseta * LEARN.Odw[il][in][in1]);
985  NET.Weights[il][in][in1] += LEARN.Odw[il][in][in1];
986  }
987 
988  }
989  }
990 
991  return(0);
992 }
993 
994 
995 /***********************************************************/
996 /* DeDwNorm */
997 /* */
998 /* computes the norm of the gradient */
999 /* */
1000 /* Author: J.Schwindling 31-Mar-99 */
1001 /***********************************************************/
1002 
1003 /* extern "C"Dllexport */dbl DeDwNorm()
1004 {
1005  int il,in,jn;
1006  dbl dd=0;
1007  for(il=1; il<NET.Nlayer; il++)
1008  for(in=0; in<NET.Nneur[il]; in++)
1009  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1010  dd += LEARN.DeDw[il][in][jn]*
1011  LEARN.DeDw[il][in][jn];
1012  return(dd);
1013 }
1014 
1015 
1016 /***********************************************************/
1017 /* DeDwProd */
1018 /* */
1019 /* scalar product between new gradient and previous one */
1020 /* (used in Polak-Ribiere Conjugate Gradient method */
1021 /* */
1022 /* Author: J.Schwindling 26-Mar-99 */
1023 /***********************************************************/
1024 
1025 /* extern "C"Dllexport */dbl DeDwProd()
1026 {
1027  int il,in,jn;
1028  dbl dd=0;
1029  for(il=1; il<NET.Nlayer; il++)
1030  for(in=0; in<NET.Nneur[il]; in++)
1031  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1032  dd += LEARN.DeDw[il][in][jn]*
1033  LEARN.ODeDw[il][in][jn];
1034  return(dd);
1035 }
1036 
1037 
1038 /***********************************************************/
1039 /* DeDwZero */
1040 /* */
1041 /* resets to 0 the gradient (should be done before DeDwSum)*/
1042 /* */
1043 /* Author: J.Schwindling 31-Mar-99 */
1044 /***********************************************************/
1045 
1046 /* extern "C"Dllexport */void DeDwZero()
1047 {
1048  int il, in, jn;
1049  for(il=1; il<NET.Nlayer; il++)
1050  for(in=0; in<NET.Nneur[il]; in++)
1051  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1052  LEARN.DeDw[il][in][jn] = 0;
1053 }
1054 
1055 
1056 /***********************************************************/
1057 /* DeDwScale */
1058 /* */
1059 /* divides the gradient by the number of examples */
1060 /* inputs: int Nexamples = number of examples */
1061 /* */
1062 /* Author: J.Schwindling 31-Mar-99 */
1063 /***********************************************************/
1064 
1065 /* extern "C"Dllexport */void DeDwScale(int Nexamples)
1066 {
1067  int il, in, jn;
1068  for(il=1; il<NET.Nlayer; il++)
1069  for(in=0; in<NET.Nneur[il]; in++)
1070  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1071  LEARN.DeDw[il][in][jn] /= (dbl) Nexamples;
1072 }
1073 
1074 
1075 /***********************************************************/
1076 /* DeDwSave */
1077 /* */
1078 /* copies the gradient DeDw into ODeDw */
1079 /* */
1080 /* Author: J.Schwindling 31-Mar-99 */
1081 /***********************************************************/
1082 
1083 /* extern "C"Dllexport */void DeDwSave()
1084 {
1085  int il, in, jn;
1086  for(il=1; il<NET.Nlayer; il++)
1087  for(in=0; in<NET.Nneur[il]; in++)
1088  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1089  LEARN.ODeDw[il][in][jn] = LEARN.DeDw[il][in][jn];
1090 }
1091 
1092 
1093 /***********************************************************/
1094 /* DeDwSaveZero */
1095 /* */
1096 /* copies the gradient DeDw into ODeDw */
1097 /* resets DeDw to 0 */
1098 /* */
1099 /* Author: J.Schwindling 23-Apr-99 */
1100 /***********************************************************/
1101 
1102 /* extern "C"Dllexport */void DeDwSaveZero()
1103 {
1104  int il, in, jn;
1105  for(il=1; il<NET.Nlayer; il++)
1106  for(in=0; in<NET.Nneur[il]; in++)
1107  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1108  {
1109  LEARN.ODeDw[il][in][jn] = LEARN.DeDw[il][in][jn];
1110  LEARN.DeDw[il][in][jn] = 0;
1111  }
1112 }
1113 
1114 
1115 /***********************************************************/
1116 /* DeDwSum */
1117 /* */
1118 /* adds to the total gradient the gradient in the current */
1119 /* example */
1120 /* inputs: int Nexamples = number of examples */
1121 /* */
1122 /* Author: J.Schwindling 31-Mar-99 */
1123 /* Modified: B.Mansoulie 23-Apr-99 */
1124 /* faster and still correct way to compute the */
1125 /* sigmoid derivative */
1126 /***********************************************************/
1127 
1128 /* extern "C"Dllexport */int DeDwSum(type_pat *ans, dbl *out, int ipat)
1129 {
1130  int il, in1, in, ii;
1131 /* dbl err[NMAX][4]; */
1132  dbl deriv;
1133  dbl *pout, *pdedw, *pdelta;
1134  register dbl a, b;
1135 /* char buf[50];*/
1136 
1137 /* output layer */
1138  b = (dbl) PAT.Pond[0][ipat];
1139  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
1140  {
1141  deriv = NET.Deriv1[NET.Nlayer-1][in];
1142  NET.Delta[NET.Nlayer-1][in] =
1143  (out[in] - (dbl) ans[in])*deriv*b;
1144  }
1145 
1146  for(il=NET.Nlayer-2; il>0; il--)
1147  {
1148 
1149  for(in=0; in<NET.Nneur[il]; in++)
1150  {
1151  deriv = NET.Deriv1[il][in];
1152  a = NET.Delta[il+1][0] * NET.Weights[il+1][0][in+1];
1153  pdelta = &(NET.Delta[il+1][1]);
1154  for(in1=1; in1<NET.Nneur[il+1]; in1++, pdelta++)
1155  {
1156  a += *pdelta * NET.Weights[il+1][in1][in+1];
1157  }
1158  NET.Delta[il][in] = a * deriv;
1159  }
1160  }
1161 
1162  for(il=1; il<NET.Nlayer; il++)
1163  {
1164  ii = NET.Nneur[il-1];
1165  for(in=0; in<NET.Nneur[il]; in++)
1166  {
1167  a = NET.Delta[il][in];
1168  LEARN.DeDw[il][in][0] += a;
1169  LEARN.DeDw[il][in][1] += a * NET.Outn[il-1][0];
1170  pout = &(NET.Outn[il-1][1]);
1171  pdedw = &(LEARN.DeDw[il][in][2]);
1172  for(in1=1; in1<ii; ++in1, ++pout, ++pdedw)
1173  {
1174  (*pdedw) += a * (*pout);
1175  }
1176  }
1177  }
1178 
1179  return(0);
1180 }
1181 
1182 
1183 /***********************************************************/
1184 /* SetTransFunc */
1185 /* */
1186 /* to set the transfert function of a neuron */
1187 /* inputs: int layer = layer number (1 -> Nlayer) */
1188 /* int neuron = neuron number (1 -> Nneur) */
1189 /* int func = transfert function: */
1190 /* 0: y=0 */
1191 /* 1: y=x */
1192 /* 2: y=1/(1+exp(-x)) */
1193 /* 3: y=tanh(x) */
1194 /* 4: y=delta*x+1/(1+exp(-x)) */
1195 /* 5: y=exp(-x**2) */
1196 /* */
1197 /* return value (int) = error value: */
1198 /* 0: no error */
1199 /* 1: layer > 4 */
1200 /* 2: neuron > NMAX */
1201 /* */
1202 /* Author: J.Schwindling 02-Apr-99 */
1203 /***********************************************************/
1204 
1205 /* extern "C"Dllexport */int SetTransFunc(int layer, int neuron,
1206  int func)
1207 {
1208  if(layer>NLMAX) return(1);
1209 /* if(neuron>NMAX) return(2);*/
1210 
1211  NET.T_func[layer-1][neuron-1] = func;
1212 
1213  return(0);
1214 }
1215 
1216 
1217 /***********************************************************/
1218 /* SetDefaultFuncs */
1219 /* */
1220 /* - sets the default transfer functions to sigmoid for */
1221 /* hidden layers and linear for output layer */
1222 /* - sets temperatures to 1 */
1223 /* */
1224 /* */
1225 /* Author: J.Schwindling 08-Apr-99 */
1226 /***********************************************************/
1227 
1229 {
1230  int il,in;
1231  for(il=0; il<NET.Nlayer; il++) {
1232  for(in=0; in<NET.Nneur[il]; in++) {
1233  NET.T_func[il][in] = 2;
1234  if(il==NET.Nlayer-1) NET.T_func[il][in] = 1;
1235  }
1236  }
1237 
1238 }
1239 
1240 
1241 /***********************************************************/
1242 /* SteepestDir */
1243 /* */
1244 /* sets the search direction to steepest descent */
1245 /* */
1246 /* Author: J.Schwindling 02-Apr-99 */
1247 /***********************************************************/
1248 
1249 /* extern "C"Dllexport */void SteepestDir()
1250 {
1251  int il,in,jn;
1252  for(il=1; il<NET.Nlayer; il++)
1253  for(in=0; in<NET.Nneur[il]; in++)
1254  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1255  dir[il][in][jn] = -LEARN.DeDw[il][in][jn];
1256 }
1257 
1258 
1259 /***********************************************************/
1260 /* CGDir */
1261 /* */
1262 /* sets the search direction to conjugate gradient dir */
1263 /* inputs: dbl beta : d(t+1) = -g(t+1) + beta d(t) */
1264 /* beta should be: */
1265 /* ||g(t+1)||^2 / ||g(t)||^2 (Fletcher-Reeves) */
1266 /* g(t+1) (g(t+1)-g(t)) / ||g(t)||^2 (Polak-Ribiere) */
1267 /* */
1268 /* Author: J.Schwindling 02-Apr-99 */
1269 /***********************************************************/
1270 
1271 /* extern "C"Dllexport */void CGDir(dbl beta)
1272 {
1273  int il,in,jn;
1274  for(il=1; il<NET.Nlayer; il++)
1275  for(in=0; in<NET.Nneur[il]; in++)
1276  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1277  {
1278  dir[il][in][jn] = -LEARN.DeDw[il][in][jn]+
1279  beta*dir[il][in][jn];
1280  }
1281 }
1282 
1283 
1284 /***********************************************************/
1285 /* DerivDir */
1286 /* */
1287 /* scalar product between gradient and direction */
1288 /* = derivative along direction */
1289 /* */
1290 /* Author: J.Schwindling 01-Jul-99 */
1291 /***********************************************************/
1292 
1293 /* extern "C"Dllexport */dbl DerivDir()
1294 {
1295  int il,in,jn;
1296  dbl ddir = 0;
1297 
1298  for(il=1; il<NET.Nlayer; il++)
1299  for(in=0; in<NET.Nneur[il]; in++)
1300  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1301  {
1302  ddir += LEARN.DeDw[il][in][jn]*dir[il][in][jn];
1303  }
1304  return(ddir);
1305 }
1306 
1307 
1308 /***********************************************************/
1309 /* GetGammaDelta */
1310 /* */
1311 /* sets the vectors Gamma (=g(t+1)-g(t)) */
1312 /* and delta (=w(t+1)-w(t)) */
1313 /* (for BFGS and Hybrid learning methods) */
1314 /* */
1315 /* Author: J.Schwindling 02-Apr-99 */
1316 /***********************************************************/
1317 
1318 /* extern "C"Dllexport */void GetGammaDelta()
1319 {
1320  int i=0;
1321  int il,in,jn;
1322  for(il=1; il<NET.Nlayer; il++)
1323  for(in=0; in<NET.Nneur[il]; in++)
1324  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1325  {
1326  Gamma[i] = LEARN.DeDw[il][in][jn]-
1327  LEARN.ODeDw[il][in][jn];
1328  delta[i] = LEARN.Odw[il][in][jn];
1329  i++;
1330  }
1331 }
1332 
1333 
1334 /***********************************************************/
1335 /* BFGSDir */
1336 /* */
1337 /* sets the search direction to BFGS direction from the */
1338 /* BFGS matrix */
1339 /* */
1340 /* inputs: int Nweights = number of weights */
1341 /* */
1342 /* Author: J.Schwindling 02-Apr-99 */
1343 /***********************************************************/
1344 
1345 /* extern "C"Dllexport */void BFGSdir(int Nweights)
1346 {
1347  dbl *g, *s;
1348  int kk=0;
1349  int il,i,j,in,jn;
1350 
1351  g = (dbl*) malloc(NET.Nweights*sizeof(dbl));
1352  s = (dbl*) malloc(Nweights*sizeof(dbl));
1353 
1354  for(il=1; kk<Nweights; il++)
1355  for(in=0; in<NET.Nneur[il]; in++)
1356  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1357  {
1358  g[kk] = LEARN.DeDw[il][in][jn];
1359  kk++;
1360  }
1361  for(i=0; i<Nweights; i++)
1362  {
1363  s[i] = 0;
1364  for(j=0; j<Nweights; j++)
1365  {
1366  s[i] += BFGSH[i][j] * g[j];
1367  }
1368  }
1369 
1370  kk = 0;
1371  for(il=1; kk<Nweights; il++)
1372  for(in=0; in<NET.Nneur[il]; in++)
1373  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1374  {
1375  dir[il][in][jn] = -s[kk];
1376  kk++;
1377  }
1378  free(g);
1379  free(s);
1380 }
1381 
1382 
1383 /***********************************************************/
1384 /* InitBFGS */
1385 /* */
1386 /* initializes BFGS matrix to identity */
1387 /* */
1388 /* inputs: int Nweights = number of weights */
1389 /* */
1390 /* Author: J.Schwindling 02-Apr-99 */
1391 /***********************************************************/
1392 
1393 /* extern "C"Dllexport */void InitBFGSH(int Nweights)
1394 {
1395  int i,j;
1396  for(i=0; i<Nweights; i++)
1397  for(j=0; j<Nweights; j++)
1398  {
1399  BFGSH[i][j] = 0;
1400  if(i==j) BFGSH[i][j] = 1;
1401  }
1402 }
1403 
1404 
1405 /***********************************************************/
1406 /* GetBFGSH */
1407 /* */
1408 /* sets the BFGS matrix */
1409 /* */
1410 /* inputs: int Nweights = number of weights */
1411 /* */
1412 /* return value (int) = 0 if no problem */
1413 /* 1 is deltaTgamma = 0 -> switch to */
1414 /* steepest dir */
1415 /* */
1416 /* Author: J.Schwindling 02-Apr-99 */
1417 /* Modified: J.Schwindling 04-May-99 */
1418 /* computations as Nw^2 , matrices removed */
1419 /* Modified: J.Schwindling 11-Jun-99 */
1420 /* return value */
1421 /***********************************************************/
1422 
1423 /* extern "C"Dllexport */int GetBFGSH(int Nweights)
1424 {
1425  typedef double dble;
1426  dble deltaTgamma=0;
1427  dble factor=0;
1428  dble *Hgamma;
1429  dble *tmp;
1430  register dble a, b;
1431  int i,j;
1432 
1433  Hgamma = (dble *) malloc(Nweights*sizeof(dble));
1434  tmp = (dble *) malloc(Nweights*sizeof(dble));
1435 
1436  for(i=0; i<Nweights; i++)
1437  {
1438  deltaTgamma += (dble) delta[i] * (dble) Gamma[i];
1439  a = 0;
1440  b = 0;
1441  for(j=0; j<Nweights; j++)
1442  {
1443  a += (dble) BFGSH[i][j] * (dble) Gamma[j];
1444  b += (dble) Gamma[j] * (dble) BFGSH[j][i];
1445  }
1446  Hgamma[i] = a;
1447  tmp[i] = b;
1448  factor += (dble) Gamma[i]*Hgamma[i];
1449  }
1450  if(deltaTgamma == 0) return 1;
1451  a = 1 / deltaTgamma;
1452  factor = 1 + factor*a;
1453 
1454  for(i=0; i<Nweights; i++)
1455  {
1456  b = (dble) delta[i];
1457  for(j=0; j<Nweights; j++)
1458  BFGSH[i][j] += (dbl) (factor*b* (dble)
1459  delta[j]-(tmp[j]*b+Hgamma[i]*(dble)delta[j]))*a;
1460  }
1461  free(Hgamma);
1462  free(tmp);
1463  return 0;
1464 }
1465 
1466 
1467 /***********************************************************/
1468 /* LineSearch */
1469 /* */
1470 /* search along the line defined by dir */
1471 /* */
1472 /* outputs: dbl *alpmin = optimal step length */
1473 /* */
1474 /* Author: B.Mansoulie 01-Jul-98 */
1475 /***********************************************************/
1476 
1477 /* extern "C"Dllexport */
1478 int LineSearch(dbl *alpmin, int *Ntest, dbl Err0)
1479 {
1480  dbl ***w0;
1481  dbl alpha1, alpha2, alpha3;
1482  dbl err1, err2, err3;
1483  dbl tau;
1484  int icount, il, in, jn;
1485 
1486  tau=LEARN.Tau;
1487 
1488 /* store weights before line search */
1489 
1490  *Ntest = 0;
1491  w0 = (dbl ***) malloc(NET.Nlayer*sizeof(dbl**));
1492  for(il=1; il<NET.Nlayer; il++)
1493  {
1494  w0[il] = (dbl **) malloc(NET.Nneur[il]*sizeof(dbl*));
1495  for(in=0; in<NET.Nneur[il]; in++)
1496  {
1497  w0[il][in] = (dbl *) malloc((NET.Nneur[il-1]+1)*
1498  sizeof(dbl));
1499  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1500  {
1501  w0[il][in][jn] = NET.Weights[il][in][jn];
1502  }
1503  }
1504  }
1505 
1506 /* compute error(w0) */
1507 
1508 /* err1 = MLP_Test(0,0);
1509  (*Ntest) ++;*/
1510  err1 = Err0;
1511 
1512  if(NET.Debug>=4) printf("err depart= %f\n",err1);
1513 
1514  *alpmin = 0;
1515  alpha1 = 0;
1516 /* alpha2 = 0.05;
1517  if(LastAlpha != 0) alpha2 = LastAlpha;*/
1518  alpha2 = LastAlpha;
1519  if(alpha2 < 0.01) alpha2 = 0.01;
1520  if(alpha2 > 2.0) alpha2 = 2.0;
1521  MLP_Line(w0,alpha2);
1522  err2 = MLP_Test(0,0);
1523  (*Ntest) ++;
1524  if(NET.Debug>=4) printf("alpha, err= %e %e\n",alpha2,err2);
1525 
1526  alpha3 = alpha2;
1527  err3 = err2;
1528 
1529 /* try to find a triplet (alpha1, alpha2, alpha3) such that
1530  Error(alpha1)>Error(alpha2)<Error(alpha3) */
1531 
1532  if(err1>err2)
1533  {
1534  for(icount=1;icount<=100;icount++)
1535  {
1536  alpha3 = alpha3*tau;
1537  MLP_Line(w0,alpha3);
1538  err3 =MLP_Test(0,0);
1539  if(NET.Debug>=4) printf("alpha, err= %e %e\n",alpha3,err3);
1540  (*Ntest) ++;
1541  if(err3>err2) break;
1542  alpha1 = alpha2;
1543  err1 = err2;
1544  alpha2 = alpha3;
1545  err2 = err3;
1546  }
1547  if(icount>=100) /* line search fails */
1548  {
1549  MLP_Line(w0,0); /* reset weights */
1550  free(w0);
1551  return(1);
1552  }
1553  }
1554  else
1555  {
1556  for(icount=1;icount<=100;icount++)
1557  {
1558  alpha2 = alpha2/tau;
1559  MLP_Line(w0,alpha2);
1560  err2 = MLP_Test(0,0);
1561  if(NET.Debug>=4) printf("alpha, err= %e %e\n",alpha2,err2);
1562  (*Ntest) ++;
1563  if(err1>err2) break;
1564  alpha3 = alpha2;
1565  err3 = err2;
1566  }
1567  if(icount>=100) /* line search fails */
1568  {
1569  MLP_Line(w0,0); /* reset weights */
1570  free(w0);
1571  LastAlpha = 0.05; /* try to be safe */
1572  return(1);
1573  }
1574  }
1575 
1576 /* find bottom of parabola */
1577 
1578  *alpmin = 0.5*(alpha1+alpha3-(err3-err1)/((err3-err2)/(alpha3-alpha2)
1579  -(err2-err1)/(alpha2-alpha1)));
1580  if(*alpmin>10000) *alpmin=10000;
1581 
1582 /* set the weights */
1583  MLP_Line(w0,*alpmin);
1584  LastAlpha = *alpmin;
1585 
1586 /* store weight changes */
1587  for(il=1; il<NET.Nlayer; il++)
1588  for(in=0; in<NET.Nneur[il]; in++)
1589  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1590  LEARN.Odw[il][in][jn] = NET.Weights[il][in][jn]
1591  - w0[il][in][jn];
1592 
1593  for(il=1; il<NET.Nlayer; il++)
1594  for(in=0; in<NET.Nneur[il]; in++)
1595  free(w0[il][in]);
1596  for(il=1; il<NET.Nlayer; il++)
1597  free(w0[il]);
1598  free(w0);
1599 
1600  return(0);
1601 }
1602 
1603 
1604 /***********************************************************/
1605 /* DecreaseSearch */
1606 /* */
1607 /* search along the line defined by dir for a point where */
1608 /* error is decreased (faster than full line search) */
1609 /* */
1610 /* outputs: dbl *alpmin = step length */
1611 /* */
1612 /* Author: J.Schwindling 11-May-99 */
1613 /***********************************************************/
1614 
1615 /* extern "C"Dllexport */
1616 int DecreaseSearch(dbl *alpmin, int *Ntest, dbl Err0)
1617 {
1618  dbl ***w0;
1619  dbl alpha1, alpha2, alpha3;
1620  dbl err1, err2, err3;
1621  dbl tau;
1622  int icount, il, in, jn;
1623 
1624  tau=LEARN.Tau;
1625 
1626 /* store weights before line search */
1627 
1628  *Ntest = 0;
1629  w0 = (dbl ***) malloc(NET.Nlayer*sizeof(dbl**));
1630  for(il=1; il<NET.Nlayer; il++)
1631  {
1632  w0[il] = (dbl **) malloc(NET.Nneur[il]*sizeof(dbl*));
1633  for(in=0; in<NET.Nneur[il]; in++)
1634  {
1635  w0[il][in] = (dbl *) malloc((NET.Nneur[il-1]+1)*
1636  sizeof(dbl));
1637  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1638  {
1639  w0[il][in][jn] = NET.Weights[il][in][jn];
1640  }
1641  }
1642  }
1643 
1644 /* compute error(w0) */
1645 
1646 /* err1 = MLP_Test(0,0);
1647  (*Ntest) ++;*/
1648  err1 = Err0;
1649 
1650  if(NET.Debug>=4) printf("err depart= %f\n",err1);
1651 
1652  *alpmin = 0;
1653  alpha1 = 0;
1654  alpha2 = 0.05;
1655  MLP_Line(w0,alpha2);
1656  err2 = MLP_Test(0,0);
1657  (*Ntest) ++;
1658 
1659  if(err2<err1)
1660  {
1661  *alpmin = alpha2;
1662  }
1663  else
1664  {
1665  alpha3 = alpha2;
1666  err3 = err2;
1667 
1668 
1669  for(icount=1;icount<=100;icount++)
1670  {
1671  alpha2 = alpha2/tau;
1672  MLP_Line(w0,alpha2);
1673  err2 = MLP_Test(0,0);
1674  (*Ntest) ++;
1675  if(err1>err2) break;
1676  alpha3 = alpha2;
1677  err3 = err2;
1678  }
1679  if(icount>=100) /* line search fails */
1680  {
1681  MLP_Line(w0,0); /* reset weights */
1682  free(w0);
1683  return(1);
1684  }
1685  *alpmin = alpha2;
1686  }
1687 
1688 /* set the weights */
1689  MLP_Line(w0,*alpmin);
1690 
1691 /* store weight changes */
1692  for(il=1; il<NET.Nlayer; il++)
1693  for(in=0; in<NET.Nneur[il]; in++)
1694  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1695  LEARN.Odw[il][in][jn] = NET.Weights[il][in][jn]
1696  - w0[il][in][jn];
1697 
1698  for(il=1; il<NET.Nlayer; il++)
1699  for(in=0; in<NET.Nneur[il]; in++)
1700  free(w0[il][in]);
1701  for(il=1; il<NET.Nlayer; il++)
1702  free(w0[il]);
1703  free(w0);
1704 
1705  return(0);
1706 }
1707 
1708 
1709 /* extern "C"Dllexport */int FixedStep(dbl alpha)
1710 {
1711  dbl ***w0;
1712  int il, in, jn;
1713 
1714  w0 = (dbl ***) malloc(NET.Nlayer*sizeof(dbl**));
1715  for(il=1; il<NET.Nlayer; il++)
1716  {
1717  w0[il] = (dbl **) malloc(NET.Nneur[il]*sizeof(dbl*));
1718  for(in=0; in<NET.Nneur[il]; in++)
1719  {
1720  w0[il][in] = (dbl *) malloc((NET.Nneur[il-1]+1)*
1721  sizeof(dbl));
1722  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1723  {
1724  w0[il][in][jn] = NET.Weights[il][in][jn];
1725  }
1726  }
1727  }
1728 
1729 
1730 /* set the weights */
1731  MLP_Line(w0,alpha);
1732 
1733 /* store weight changes */
1734  for(il=1; il<NET.Nlayer; il++)
1735  for(in=0; in<NET.Nneur[il]; in++)
1736  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1737  LEARN.Odw[il][in][jn] = NET.Weights[il][in][jn]
1738  - w0[il][in][jn];
1739 
1740  for(il=1; il<NET.Nlayer; il++)
1741  for(in=0; in<NET.Nneur[il]; in++)
1742  free(w0[il][in]);
1743  for(il=1; il<NET.Nlayer; il++)
1744  free(w0[il]);
1745  free(w0);
1746 
1747  return(0);
1748 }
1749 
1750 /***********************************************************/
1751 /* MLP_Line */
1752 /* */
1753 /* sets the weights to a point along a line */
1754 /* */
1755 /* inputs: dbl ***w0 = initial point */
1756 /* dbl alpha = step length */
1757 /* */
1758 /* Author: B.Mansoulie 01-Jul-98 */
1759 /***********************************************************/
1760 
1761 /* extern "C"Dllexport */
1762 void MLP_Line(dbl ***w0, dbl alpha)
1763 {
1764  register int il,in,jn;
1765 
1766  for(il=1; il<NET.Nlayer; il++)
1767  for(in=0; in<NET.Nneur[il]; in++)
1768  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1769  NET.Weights[il][in][jn] = w0[il][in][jn]+
1770  alpha*dir[il][in][jn];
1771 
1772 }
1773 
1774 
1775 /***********************************************************/
1776 /* LineSearchHyb */
1777 /* */
1778 /* search along the line defined by dir */
1779 /* */
1780 /* outputs: dbl *alpmin = optimal step length */
1781 /* */
1782 /* Author: B.Mansoulie 01-Jul-98 */
1783 /***********************************************************/
1784 
1785 /* extern "C"Dllexport */
1786 int LineSearchHyb(dbl *alpmin, int *Ntest)
1787 {
1788  dbl ***w0;
1789  dbl alpha1, alpha2, alpha3;
1790  dbl err1, err2, err3;
1791  dbl tau;
1792  int icount, il, in, jn;
1793 
1794 /* char buf [50];
1795  sprintf (buf,"entree linesearchhyb\n");
1796  MessageBoxA (0,buf,"dans FreePatterns",MB_OK);*/
1797 
1798  if(NET.Debug>=4){
1799  printf(" entry LineSearchHyb \n");
1800  }
1801  tau=LEARN.Tau;
1802 
1803 /* store weights before line search */
1804 
1805  *Ntest = 0;
1806  w0 = (dbl ***) malloc((NET.Nlayer-1)*sizeof(dbl**));
1807  for(il=1; il<NET.Nlayer-1; il++)
1808  {
1809  w0[il] = (dbl **) malloc(NET.Nneur[il]*sizeof(dbl*));
1810  for(in=0; in<NET.Nneur[il]; in++)
1811  {
1812  w0[il][in] = (dbl *) malloc((NET.Nneur[il-1]+1)*
1813  sizeof(dbl));
1814  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1815  {
1816  w0[il][in][jn] = NET.Weights[il][in][jn];
1817  }
1818  }
1819  }
1820 
1821 /* compute error(w0) */
1822  err1 = MLP_Test(0,1);
1823  (*Ntest) ++;
1824  if(NET.Debug>=4) printf("LinesearchHyb err depart= %f\n",err1);
1825 
1826  *alpmin = 0;
1827  alpha1 = 0;
1828 /* alpha2 = 0.05;
1829  if(LastAlpha != 0) alpha2 = LastAlpha;*/
1830  alpha2 = LastAlpha;
1831  if(alpha2 < 0.01) alpha2 = 0.01;
1832  if(alpha2 > 2.0) alpha2 = 2.0;
1833  MLP_LineHyb(w0,alpha2);
1834  err2 = MLP_Test(0,1);
1835  (*Ntest) ++;
1836 
1837  alpha3 = alpha2;
1838  err3 = err2;
1839 
1840 /* try to find a triplet (alpha1, alpha2, alpha3) such that
1841  Error(alpha1)>Error(alpha2)<Error(alpha3) */
1842 
1843  if(err1>err2)
1844  {
1845  for(icount=1;icount<=100;icount++)
1846  {
1847  alpha3 = alpha3*tau;
1848  MLP_LineHyb(w0,alpha3);
1849  err3 = MLP_Test(0,1);
1850  (*Ntest) ++;
1851  if(err3>err2) break;
1852  alpha1 = alpha2;
1853  err1 = err2;
1854  alpha2 = alpha3;
1855  err2 = err3;
1856  }
1857  if(icount>=100) /* line search fails */
1858  {
1859  MLP_LineHyb(w0,0); /* reset weights */
1860  free(w0);
1861  return(1);
1862  }
1863  }
1864  else
1865  {
1866  for(icount=1;icount<=100;icount++)
1867  {
1868  alpha2 = alpha2/tau;
1869  MLP_LineHyb(w0,alpha2);
1870  err2 = MLP_Test(0,1);
1871  (*Ntest) ++;
1872  if(err1>err2) break;
1873  alpha3 = alpha2;
1874  err3 = err2;
1875  }
1876  if(icount>=100) /* line search fails */
1877  {
1878  MLP_LineHyb(w0,0); /* reset weights */
1879  free(w0);
1880  return(1);
1881  }
1882  }
1883 
1884 /* find bottom of parabola */
1885 
1886  *alpmin = 0.5*(alpha1+alpha3-(err3-err1)/((err3-err2)/(alpha3-alpha2)
1887  -(err2-err1)/(alpha2-alpha1)));
1888  if(*alpmin>10000) *alpmin=10000;
1889 
1890 /* set the weights */
1891  MLP_LineHyb(w0,*alpmin);
1892  LastAlpha = *alpmin;
1893 
1894 /* store weight changes */
1895  for(il=1; il<NET.Nlayer-1; il++)
1896  for(in=0; in<NET.Nneur[il]; in++)
1897  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1898  LEARN.Odw[il][in][jn] = NET.Weights[il][in][jn]
1899  - w0[il][in][jn];
1900 
1901  for(il=1; il<NET.Nlayer-1; il++)
1902  for(in=0; in<NET.Nneur[il]; in++)
1903  free(w0[il][in]);
1904  for(il=1; il<NET.Nlayer-1; il++)
1905  free(w0[il]);
1906  free(w0);
1907  if(NET.Debug>=4){
1908  printf(" exit LineSearchHyb \n");
1909  }
1910 
1911  return(0);
1912 }
1913 
1914 
1915 /***********************************************************/
1916 /* MLP_LineHyb */
1917 /* */
1918 /* sets the weights to a point along a line */
1919 /* (for hybrid methods) */
1920 /* */
1921 /* inputs: dbl ***w0 = initial point */
1922 /* dbl alpha = step length */
1923 /* */
1924 /* Author: B.Mansoulie 01-Jul-98 */
1925 /***********************************************************/
1926 
1927 /* extern "C"Dllexport */
1928 void MLP_LineHyb(dbl ***w0, dbl alpha)
1929 {
1930  int il,in,jn;
1931  for(il=1; il<NET.Nlayer-1; il++)
1932  for(in=0; in<NET.Nneur[il]; in++)
1933  for(jn=0; jn<=NET.Nneur[il-1]; jn++)
1934  {
1935  NET.Weights[il][in][jn] = w0[il][in][jn]+
1936  alpha*dir[il][in][jn];
1937  }
1938  MLP_ResLin();
1939 }
1940 
1941 
1942 /***********************************************************/
1943 /* SetLambda */
1944 /* */
1945 /* sets the coefficient of the regularisation term for */
1946 /* hybrid learning method */
1947 /* */
1948 /* input: double Wmax = typical maximum weight */
1949 /* */
1950 /* Author: J.Schwindling 13-Apr-99 */
1951 /***********************************************************/
1952 
1953 /* extern "C"Dllexport */
1954 void SetLambda(double Wmax)
1955 {
1956  dbl err;
1957  err = MLP_Test(0,0);
1958  LEARN.Alambda =
1959  LEARN.Lambda*err/(Wmax*Wmax*(dbl)(NET.Nneur[NET.Nlayer-2]+1));
1960 }
1961 
1962 
1963 /***********************************************************/
1964 /* MLP_ResLin */
1965 /* */
1966 /* resolution of linear system of equations for hybrid */
1967 /* training method */
1968 /* */
1969 /* */
1970 /* Author: B.Mansoulie end-98 */
1971 /* Modified: J.Schwindling 29-APR-99 */
1972 /* use dgels from LAPACK */
1973 /***********************************************************/
1974 
1975 /* extern "C"Dllexport */
1977 {
1978 /* dbl rrans[NMAX], rrout[NMAX];*/
1979 /* type_pat rrin[NMAX];*/
1980  doublereal *HR,*dpat; //,*wlin,*SV;
1981  double err,lambda,lambda2;
1982  integer Nl,M,Nhr,khr,nrhs,iret,ierr;
1983  int il, in, inl, ipat;
1984  register dbl a;
1985  char Trans = 'N';
1986 
1987 
1988 /* integer rank; */
1989 // doublereal rcond = -1; /* use machine precision */
1990 
1991  lambda2 = LEARN.Alambda;
1992 
1993 /* Nl = number of linear weights
1994  M = number of terms in linear system = number of examples + regularisation*/
1995  Nl = NET.Nneur[NET.Nlayer-2] + 1;
1996  M = PAT.Npat[0]+Nl;
1997 
1998  integer Lwork = 5 * M;
1999  doublereal *Work = (doublereal*) malloc((int) Lwork*sizeof(doublereal));
2000 
2001 /* memory allocation */
2002  dpat = (doublereal*) malloc((int) M*sizeof(doublereal));
2003 // wlin = (doublereal*) malloc((int) Nl*sizeof(doublereal));
2004 // SV = (doublereal*) malloc((int) Nl*sizeof(doublereal));
2005 
2006  Nhr = M * Nl;
2007  HR = (doublereal*) malloc((int) Nhr*sizeof(doublereal));
2008  err = 0.;
2009  for(ipat=0;ipat<PAT.Npat[0];ipat++)
2010  {
2011 /* *** Filling dpat and HR *** */
2012 /* for(in=0; in<NET.Nneur[0]; in++)
2013  {
2014  rrin[in] = PAT.Rin[0][ipat][in];
2015  }*/
2016 
2017  MLP_Out(PAT.Rin[0][ipat],NET.Outn[NET.Nlayer-1]);
2018 /* MLP_Out(rrin,rrout);*/
2019  for(in=0; in<NET.Nneur[NET.Nlayer-1]; in++)
2020  {
2021  a = (dbl) PAT.Rans[0][ipat][in];
2022  }
2023  il = NET.Nlayer-2;
2024  dpat[ipat] = (dbl) PAT.Rans[0][ipat][0]*sqrt(PAT.Pond[0][ipat]);
2025  khr = ipat;
2026  HR[khr] = (dbl) sqrt(PAT.Pond[0][ipat]);
2027  for(in=0;in<NET.Nneur[il];in++)
2028  {
2029  khr = M *(in+1) + ipat;
2030  HR[khr] = NET.Outn[il][in]*
2031  (dbl) sqrt(PAT.Pond[0][ipat]);
2032  }
2033  }
2034  il = NET.Nlayer-2;
2035  lambda = sqrt(lambda2);
2036  for(ipat=0;ipat<=NET.Nneur[il];ipat++)
2037  {
2038  dpat[ipat+PAT.Npat[0]] = 0;
2039  for(in=0;in<=NET.Nneur[il];in++)
2040  {
2041  khr = M *in + ipat + PAT.Npat[0];
2042  HR[khr] = 0;
2043  if(in==ipat) HR[khr]=lambda;
2044  }
2045  }
2046  if(NET.Debug>=4)
2047  {
2048  err = MLP_Test(0,0);
2049  printf("entry ResLin, err=MLP_Test(0,0), err= %f\n",err);
2050  }
2051 /* */
2052 /* Trouve les poids lineaires par resolution lineaire */
2053 /* */
2054  nrhs = 1;
2055  ierr = dgels_(&Trans,&M,&Nl,&nrhs,HR,&M,dpat,&M,Work,
2056  &Lwork,&iret);
2057  if(iret != 0) printf("Warning from dgels: iret = %d\n",(int)iret);
2058  if(ierr != 0) printf("Warning from dgels: ierr = %d\n",(int)ierr);
2059 
2060 /* ierr = dgelss_(&M,&Nl,&nrhs,HR,&M,dpat,&M,SV,&rcond,&rank,Work,&Lwork,
2061  &iret);
2062  if(iret != 0) printf("Warning from dgelss: iret = %d\n",iret);
2063  if(ierr != 0) printf("Warning from dgelss: ierr = %d\n",ierr);*/
2064 
2065  il = NET.Nlayer-1;
2066  for (inl=0; inl<=NET.Nneur[il-1];inl++)
2067  {
2068  NET.Weights[il][0][inl] = dpat[inl];
2069  }
2070  if(NET.Debug>=4)
2071  {
2072  err = MLP_Test(0,0);
2073  printf("ResLin, apres tlsfor, err= %f\n",err);
2074  }
2075  free(Work);
2076  free(dpat);
2077 // free(wlin);
2078  free(HR);
2079 // free(SV);
2080 }
2081 
2082 /***********************************************************/
2083 /* EtaDecay */
2084 /* */
2085 /* decreases the learning parameter eta by the factor */
2086 /* LEARN.Decay */
2087 /* */
2088 /* Author: J.Schwindling 14-Apr-99 */
2089 /***********************************************************/
2090 
2091 void EtaDecay()
2092 {
2093  LEARN.eta *= LEARN.Decay;
2094 }
2095 
2096 
2097 /***********************************************************/
2098 /* ShuffleExamples */
2099 /* */
2100 /* Shuffles the learning examples (for stochastic method) */
2101 /* */
2102 /* Author: J.Schwindling 14-Apr-1999 */
2103 /* Modified: J.Schwindling 21-Jul-1999 inline MLP_Rand */
2104 /***********************************************************/
2105 
2106 /* extern "C"Dllexport */
2107 int ShuffleExamples(int n, int *index)
2108 {
2109  int i,ii,itmp;
2110  dbl a = (dbl) (n-1);
2111 
2112  for(i=0;i<n;i++)
2113  {
2114  ii = (int) MLP_Rand(0.,a);
2115  itmp = index[ii];
2116  index[ii] = index[i];
2117  index[i] = itmp;
2118  }
2119  return 0;
2120 }
2121 
2122 
2123 /***********************************************************/
2124 /* MLP_Rand */
2125 /* */
2126 /* Random Numbers (to initialize weights) */
2127 /* */
2128 /* inputs : dbl min, dbl max = random number between */
2129 /* min and max */
2130 /* return value: (double) random number */
2131 /* */
2132 /* Author: J.Schwindling 14-Apr-99 */
2133 /***********************************************************/
2134 
2135 /* extern "C"Dllexport */
2136 double MLP_Rand(dbl mini, dbl maxi)
2137 {
2138 return mini+(maxi-mini)*random()/RAND_MAX;
2139 }
2140 
2141 
2142 /***********************************************************/
2143 /* InitWeights */
2144 /* */
2145 /* initializes the weights to random numbers between */
2146 /* -0.5 : 0.5 */
2147 /* */
2148 /* Author: J.Schwindling 14-Apr-99 */
2149 /***********************************************************/
2150 
2151 /* extern "C"Dllexport */
2153 {
2154  int ilayer,ineur,i;
2155 
2156  for(ilayer=1;ilayer<NET.Nlayer;ilayer++)
2157  for(ineur=0;ineur<NET.Nneur[ilayer];ineur++)
2158  for(i=0;i<=NET.Nneur[ilayer-1];i++)
2159  NET.Weights[ilayer][ineur][i]=
2160  (dbl) MLP_Rand(-0.5, 0.5);
2161 }
2162 
2163 
2164 /***********************************************************/
2165 /* PrintWeights */
2166 /* */
2167 /* Print weights on the screen */
2168 /* */
2169 /* Author: J.Schwindling 14-Apr-99 */
2170 /***********************************************************/
2171 
2172 /* extern "C"Dllexport */
2174 {
2175  int ilayer,ineur,i;
2176 
2177  for(ilayer=1; ilayer<NET.Nlayer; ilayer++)
2178  {
2179  if(MessLang==1)
2180  {
2181  printf("Couche %d\n",ilayer);
2182  }
2183  else
2184  {
2185  printf("Layer %d\n",ilayer);
2186  }
2187  for(ineur=0; ineur<NET.Nneur[ilayer]; ineur++)
2188  {
2189  if(MessLang==1)
2190  {
2191  printf("Neurone %d",ineur);
2192  }
2193  else
2194  {
2195  printf("Neuron %d",ineur);
2196  }
2197  for(i=0; i<=NET.Nneur[ilayer-1]; i++)
2198  {
2199  printf(" %f",
2200  (double) NET.Weights[ilayer][ineur][i]);
2201  }
2202  printf("\n");
2203  }
2204  printf("\n");
2205  }
2206 }
2207 
2208 
2209 /***********************************************************/
2210 /* ReadPatterns */
2211 /* */
2212 /* parser for learn.pat or test.pat files */
2213 /* */
2214 /* inputs: char *filename = name of the file to read */
2215 /* int ifile = 0: learning examples */
2216 /* 1: test examples */
2217 /* */
2218 /* outputs: int *inet = 1 if a network is defined */
2219 /* int *ilearn = 1 if a learning method is defined*/
2220 /* int *iexamples = 1 if examples are defined */
2221 /* */
2222 /* return value (int) = 0: no error */
2223 /* = -1: file could not be opened */
2224 /* */
2225 /* Author: J.Schwindling 20-Apr-99 */
2226 /* Modified: J.Schwindling 01-Jun-99 return inet, ilearn */
2227 /* iexamples */
2228 /* J.Schwindling 21-Sep-99 return value = -1 */
2229 /***********************************************************/
2230 
2231 #define CLEN 1024
2232 
2233 /* extern "C"Dllexport */
2234 int ReadPatterns(char *filename, int ifile,
2235  int *inet, int *ilearn, int *iexamples)
2236 {
2237 char s[CLEN], s2[CLEN], cc[6], cc2[6];
2238 char otherfile[CLEN];
2239 double p;
2240 //int line,i,j;
2241 int line,i;
2242 //int l,ll,ipat,nmax,il,in,tf;
2243 int l,ll,ipat,nmax;
2244 int np=0; /* nombre d'exemples */
2245 int nin=0; /* nombre d'entrees */
2246 int nout=0; /* nombre de sorties */
2247 int npon=0;
2248 int ntot, ierr;
2249 //char **ss;
2250 char **ss=0;
2251 FILE *LVQpat;
2252 int nlayer, nneur[NLMAX];
2253 
2254 printf("\nLoading file %s\n",filename);
2255 LVQpat=fopen(filename,"r");
2256 if(LVQpat == 0) return -1;
2257 
2258 line=0;
2259 
2260 while(fgets(s,CLEN,LVQpat))
2261  {
2262  if(*s=='N')
2263  {
2264  if(*(s+1)=='N') /* NNEU */
2265  {
2266  printf("Number of neurons %s",s);
2267  *inet = 1;
2268  sscanf(s,"%s %s",cc,s2);
2269  ierr = GetNetStructure(s2,&nlayer,nneur);
2270  if(ierr != 0) return ierr;
2271  ierr = MLP_SetNet(&nlayer,nneur);
2272  if(ierr != 0) return ierr;
2273  }
2274  else
2275  {
2276  sscanf(s,"%s %d",cc,&l);
2277  if(*(cc+1)=='P') /* NPAT */
2278  {
2279  np=l;
2280  printf("Number of patterns %d\n",np);
2281  }
2282  else if(*(cc+1)=='I') /* NINP */
2283  {
2284  nin=l;
2285  PAT.Nin = nin;
2286  printf("Number of inputs %d\n",nin);
2287  }
2288  else if(*(cc+1)=='O' && *(cc+2)=='U') /* NOUT */
2289  {
2290  nout=l;
2291  PAT.Nout = nout;
2292  printf("Number of outputs %d\n",nout);
2293  }
2294  else if(*(cc+1)=='O' && *(cc+2)=='R') /* NORM */
2295  {
2296  DIVERS.Norm=l;
2297  if(l==1) printf("Normalize inputs\n");
2298  }
2299 /* obsolete datacard NLAY */
2300  else if(*(cc+1)=='L')
2301  {
2302  printf("NLAY datacard is no longer needed\n");
2303  }
2304  else if(*(cc+1)=='E') /* NEPO */
2305  {
2306  LEARN.Nepoch=l;
2307  printf("Number of epochs %d\n",l);
2308  }
2309  else if(*(cc+1)=='R') /* NRES */
2310  {
2311  LEARN.Nreset=l;
2312  printf(
2313  "Reset to steepest descent every %d epochs\n",
2314  l);
2315  }
2316  }
2317  }
2318  else if(*s=='L')
2319  {
2320  if(*(s+1)=='P') /* LPAR */
2321  {
2322  sscanf(s,"%s %le",cc,&p);
2323  printf("Learning parameter %f\n",p);
2324  LEARN.eta = (dbl) p;
2325  }
2326  else if(*(s+1)=='M') /* LMET */
2327  {
2328  *ilearn = 1;
2329  sscanf(s,"%s %d",cc,&(LEARN.Meth));
2330  printf("Learning method = ");
2331  switch(LEARN.Meth)
2332  {
2333  case 1: printf("Stochastic Minimization\n");
2334  break;
2335  case 2: printf("Steepest descent with fixed step\n");
2336  break;
2337  case 3: printf("Steepest descent with line search\n"); break;
2338  case 4: printf("Polak-Ribiere Conjugate Gradients\n"); break;
2339  case 5: printf("Fletcher-Reeves Conjugate Gradients\n");
2340  break;
2341  case 6: printf("BFGS\n");
2342  break;
2343  case 7: printf("Hybrid BFGS-linear\n");
2344  break;
2345  default: printf("Error: unknown method\n"); break;
2346  }
2347 
2348  }
2349  else if(*(s+1)=='T') /* LTAU */
2350  {
2351  sscanf(s,"%s %lf",cc,&p);
2352  printf("Tau %f\n",p);
2353  LEARN.Tau = (dbl) p;
2354  }
2355  else if(*(s+1)=='A') /* LAMB */
2356  {
2357  sscanf(s,"%s %lf",cc,&p);
2358  printf("Lambda %f\n",p);
2359  LEARN.Lambda = (dbl) p;
2360  }
2361  }
2362  else if(*s=='F')
2363  {
2364  if(*(s+1)=='S') /* FSPO */
2365  {
2366  sscanf(s,"%s %le",cc,&p);
2367  printf("Flat spot elimination parameter %f\n",p);
2368  LEARN.delta = (dbl) p;
2369  }
2370  else if(*(s+1)=='I') /* FILE */
2371  {
2372  sscanf(s,"%s %s",cc,otherfile);
2373  ierr = ReadPatterns(otherfile,ifile, inet, ilearn, iexamples);
2374  if(ierr != 0) return ierr;
2375  }
2376  }
2377  else if(*s=='M') /* momentum */
2378  {
2379  sscanf(s,"%s %le",cc,&p);
2380  printf("Momentum term %f\n",p);
2381  LEARN.epsilon = (dbl) p;
2382  }
2383  else if(*s=='O') /* OUTx */
2384  {
2385  if(*(s+3)=='W') /* OUTW */
2386  {
2387  sscanf(s,"%s %d",cc,&OutputWeights);
2388  if(OutputWeights == 0)
2389  {
2390  printf("Never write file weights.out\n");
2391  }
2392  else if(OutputWeights == -1)
2393  {
2394  printf("Write weights to output file at the end\n");
2395  }
2396  else
2397  {
2398  printf("Write weights to file every %d epochs\n",
2399  OutputWeights);
2400  }
2401  }
2402  else if(*(s+3)=='F') /* OUTF */
2403  {
2404  sscanf(s,"%s %s",cc,cc2);
2405  if(*cc2=='F' || *cc2=='C')
2406  {
2407  DIVERS.Outf = *cc2;
2408  }
2409  else
2410  {
2411  printf(" *** Error while loading file %s at line %s :",
2412  filename,s);
2413  printf(" unknown language\n");
2414  }
2415  }
2416  else
2417  {
2418  printf(" *** Error while loading file %s at line %s\n",
2419  filename,s);
2420  }
2421  }
2422  else if(*s=='R') /* RDWT */
2423  {
2424  sscanf(s,"%s %d",cc,&(NET.Rdwt));
2425  if(NET.Rdwt == 0)
2426  {
2427  printf("Random weights \n");
2428  }
2429  else
2430  {
2431  printf("Read weights from file weights.in\n");
2432  }
2433  }
2434  else if(*s=='S') /* STAT */
2435  {
2436  sscanf(s,"%s %d",cc,&(DIVERS.Stat));
2437  }
2438 /* else if(*s=='T') TFUN
2439  {
2440  sscanf(s,"%s %d %d %d",cc,&il,&in,&tf);
2441  SetTransFunc(il,in,tf);
2442  } */
2443  else if(*s=='H') /* HESS */
2444  {
2445  sscanf(s,"%s %d",cc,&(DIVERS.Ihess));
2446  }
2447  else if(*s=='D')
2448  {
2449  if(*(s+1)=='C') /* DCAY */
2450  {
2451  sscanf(s,"%s %le",cc,&p);
2452  LEARN.Decay = p;
2453  printf("Learning parameter decay %f\n",
2454  (double) LEARN.Decay);
2455  }
2456  if(*(s+1)=='B') /* DBIN */
2457  {
2458  sscanf(s,"%s %d",cc,&(DIVERS.Dbin));
2459  printf("Fill histogram every %d epochs\n",DIVERS.Dbin);
2460  }
2461  if(*(s+1)=='E') /* DEBU */
2462  {
2463  sscanf(s,"%s %d",cc,&(NET.Debug));
2464  printf("Debug mode %d\n",NET.Debug);
2465  }
2466  }
2467  else if(*s=='P') /* POND */
2468  {
2469  npon = CountLexemes(s);
2470  if(npon==2)
2471  {
2472  sscanf(s,"%s %d",cc,&(PAT.Iponde));
2473  }
2474  else
2475  {
2476  ss = (char**) malloc((npon+1)*sizeof(char*));
2477  for(i=0;i<=npon;i++)
2478  ss[i]=(char*) malloc(40*sizeof(char));
2479  getnLexemes(npon,s,ss);
2480  sscanf(ss[1],"%d",&(PAT.Iponde));
2481  for(i=2;i<npon;i++)
2482  {
2483  sscanf(ss[i],"%le",&(PAT.Ponds[i-2]));
2484  }
2485  }
2486  if(PAT.Iponde==0)
2487  {
2488  npon = 0;
2489  }
2490  else
2491  {
2492  npon = 1;
2493  }
2494  }
2495  else if(*s=='#') /* comments */
2496  {
2497  }
2498  else /* exemple itself */
2499  {
2500  if(np==0) return 1;
2501  if(nin==0) return 2;
2502  if(nout==0) return 3;
2503 
2504 
2505 /* store number of exemples and allocate memory*/
2506  if(line==0)
2507  {
2508  PAT.Npat[ifile] = np;
2509  ierr = AllocPatterns(ifile,np,nin,nout,0);
2510  if(ierr != 0) return ierr;
2511  *iexamples = 1;
2512  }
2513 
2514 /* now get exemple */
2515 
2516  line++;
2517  ll = (line-1)%2;
2518  ipat = (line-1)/2;
2519  /* printf("Loading event \t %d\r",ipat);*/
2520 /* if(ipat>NPMAX)
2521  {
2522  printf("Too many examples in file\n");
2523  printf("Loading %d examples\n",NPMAX);
2524  PAT.Npat[ifile] = NPMAX;
2525  break;
2526  }
2527 */
2528 
2529 /* allocate the number of lines */
2530 
2531  if(line==1)
2532  {
2533 
2534  nmax = nin;
2535  if(nout>nin) nmax=nout;
2536  ss = (char**) malloc((nmax+1)*sizeof(char*));
2537  if(ss == 0) return -111;
2538  for(i=0;i<=nmax;i++)
2539  {
2540  ss[i]=(char*) malloc(40*sizeof(char));
2541  if(ss[i] == 0) return -111;
2542  }
2543  }
2544 
2545  if(ll==0) /* inputs */
2546  {
2547  getnLexemes(nin,s,ss);
2548  for(i=0;i<nin;i++)
2549  {
2550  sscanf(ss[i],"%le",&p);
2551  PAT.Rin[ifile][ipat][i] = (type_pat) p;
2552  }
2553  }
2554  else /* answers */
2555  {
2556  ntot=nout+npon;
2557  getnLexemes(ntot,s,ss);
2558  for(i=0;i<ntot;i++)
2559  {
2560  sscanf(ss[i],"%le",&p);
2561  if(i<nout)
2562  {
2563  PAT.Rans[ifile][ipat][i] = (type_pat) p;
2564  }
2565  else
2566  {
2567  if(PAT.Iponde==1)
2568  {
2569  PAT.Pond[ifile][ipat] =
2570  (type_pat) p;
2571  }
2572  else
2573  {
2574  PAT.Pond[ifile][ipat] =
2575  (type_pat) PAT.Ponds[(int) p -1];
2576  }
2577  }
2578  }
2579  }
2580  }
2581  }
2582  printf("%d examples loaded \n\n",PAT.Npat[ifile]);
2583  fclose(LVQpat);
2584  return 0;
2585 }
2586 
2587 
2588 /* CountLexemes returns the number of lexemes in s separated by blanks.*/
2589 /* extern "C"Dllexport */
2590 int CountLexemes(char *s)
2591 {
2592  char tmp[1024];
2593  int i=0;
2594 
2595  strcpy(tmp,s);
2596  if (strtok(tmp," "))
2597  {
2598  i=1;
2599  while (strtok(NULL," ")) i++;
2600  }
2601  return i;
2602 }
2603 
2604 /* splits s in substrings ss separated by blanks*/
2605 /* extern "C"Dllexport */
2606 void getnLexemes(int n, char *s, char **ss)
2607 {
2608  char tmp[1024];
2609  int i;
2610  strcpy(tmp,s);
2611  if (n>0)
2612  {
2613  strcpy(ss[0],strtok(tmp," "));
2614  for (i=1;i<n;i++)
2615  strcpy(ss[i],strtok(NULL," "));
2616  }
2617 }
2618 
2619 /* extern "C"Dllexport */
2620 void getLexemes(char *s,char **ss)
2621 {
2622  char tmp[1024];
2623  int i,n;
2624 
2625  strcpy(tmp,s);
2626  n=CountLexemes(tmp);
2627  if (n>0)
2628  {
2629  strcpy(ss[0],strtok(tmp," "));
2630  for (i=1;i<n;i++)
2631  strcpy(ss[i],strtok(NULL," "));
2632  }
2633 }
2634 
2635 
2636 /***********************************************************/
2637 /* LearnFree */
2638 /* */
2639 /* frees memory allocated for learning */
2640 /* */
2641 /* Author: J.Schwindling 28-May-99 */
2642 /***********************************************************/
2643 
2644 /* extern "C"Dllexport */
2646 {
2647  int il,in;
2648  if(LearnMemory==0) return;
2649  LearnMemory = 0;
2650  for(il=0; il<NET.Nlayer; il++)
2651  {
2652  for(in=0; in<NET.Nneur[il]; in++)
2653  {
2654  free(dir[il][in]);
2655  }
2656  free(dir[il]);
2657  }
2658  free(dir);
2659  if(BFGSMemory==0) return;
2660  BFGSMemory = 0;
2661  for(il=0; il<NET.Nweights; il++)
2662  {
2663  free(BFGSH[il]);
2664  }
2665  free(BFGSH);
2666  free(Gamma);
2667  free(delta);
2668 
2669 /* if(JacobianMemory == 0) return;
2670  JacobianMemory = 0;
2671  for(il=0; il<PAT.Npat[0]; il++) free(JacobianMatrix[il]);
2672  free(JacobianMatrix); */
2673 }
2674 
2675 
2676 /***********************************************************/
2677 /* LearnAlloc */
2678 /* */
2679 /* memory allocation for vectors and matrices used in */
2680 /* conjugate gradients or BFGS like methods */
2681 /* */
2682 /* return value (int) = error code: 0 no error */
2683 /* -111 no memory */
2684 /* */
2685 /* Author: J.Schwindling 20-Apr-99 */
2686 /* Modified: J.Schwindling 31-Jan-2000 error code */
2687 /***********************************************************/
2688 
2689 /* extern "C"Dllexport */
2691 {
2692  int il,in,i;
2693  int Nweights = 0;
2694 
2695  if(LearnMemory != 0) LearnFree();
2696  LearnMemory = 1;
2697  dir = (dbl ***) malloc(NET.Nlayer*sizeof(dbl**));
2698  if(dir == 0) return -111;
2699 
2700  for(il=0; il<NET.Nlayer; il++)
2701  {
2702  dir[il] = (dbl **) malloc(NET.Nneur[il]*sizeof(dbl*));
2703  if(dir[il] == 0) return -111;
2704  for(in=0; in<NET.Nneur[il]; in++)
2705  {
2706  if(il==0)
2707  {
2708 /* TODO: understand implications of hard-coded 101 */
2709  dir[0][in] = (dbl *)
2710  malloc(101*sizeof(dbl));
2711  if(dir[0][in] == 0) return -111;
2712  }
2713  else
2714  {
2715  dir[il][in] = (dbl *)
2716  malloc((NET.Nneur[il-1]+1)*sizeof(dbl));
2717  if(dir[il][in] == 0) return -111;
2718  Nweights += NET.Nneur[il-1]+1;
2719  }
2720  }
2721  }
2722  NET.Nweights = Nweights;
2723 
2724  if(BFGSMemory==0 && LEARN.Meth>= 6)
2725  {
2726  BFGSMemory = 1;
2727  Gamma = (dbl*) malloc(Nweights*sizeof(dbl));
2728  delta = (dbl*) malloc(Nweights*sizeof(dbl));
2729  BFGSH = (dbl**) malloc(Nweights*sizeof(dbl*));
2730  if(Gamma == 0 || delta == 0 || BFGSH == 0)
2731  return -111;
2732 
2733  for(i=0; i<Nweights; i++)
2734  {
2735  BFGSH[i] = (dbl*) malloc(Nweights*sizeof(dbl));
2736  if(BFGSH[i] == 0) return -111;
2737  }
2738  }
2739 
2740 /* if(JacobianMemory==0)
2741  {
2742  JacobianMemory = 1;
2743  printf("JacobianMemory = %d\n",JacobianMemory);
2744  JacobianMatrix = (dbl **) malloc(PAT.Npat[0]*sizeof(dbl *));
2745  for(i=0; i<PAT.Npat[0]; i++)
2746  JacobianMatrix[i] =
2747  (dbl*) malloc(Nweights*sizeof(dbl));
2748  printf("end memory alloc\n");
2749  }
2750 
2751  if(DIVERS.Ihess==1) HessianAlloc(Nweights);*/
2752 
2753  return 0;
2754 }
2755 
2756 
2757 /***********************************************************/
2758 /* MLP_PrFFun */
2759 /* */
2760 /* writes the MLP function to file as a fortran function */
2761 /* */
2762 /* inputs : char *filename = name of the output file */
2763 /* */
2764 /* return value (int) = 0: no error */
2765 /* -1: could not open file */
2766 /* */
2767 /* Author: J.Schwindling 20-Apr-99 */
2768 /* Modified: J.Schwindling 05-May-99 */
2769 /* add normalization of inputs */
2770 /* J.Schwindling 30-Nov-99 return code */
2771 /***********************************************************/
2772 
2773 /* extern "C"Dllexport */
2775 {
2776  int il,in,jn;
2777  FILE *W;
2778 
2779  W=fopen(filename,"w");
2780  if(W==0) return -1;
2781  fprintf(W," SUBROUTINE RNNFUN(rin,rout)\n");
2782  fprintf(W," DIMENSION RIN(%d)\n",NET.Nneur[0]);
2783  fprintf(W," DIMENSION ROUT(%d)\n",NET.Nneur[NET.Nlayer-1]);
2784  fprintf(W,"C\n");
2785 
2786  for(in=0; in<NET.Nneur[0]; in++)
2787  {
2788  if(DIVERS.Norm==0)
2789  {
2790  fprintf(W," OUT%d = RIN(%d)\n",in+1,in+1);
2791  }
2792  else
2793  {
2794  fprintf(W," OUT%d = (RIN(%d)-%e)/%e\n",in+1,in+1,
2795  STAT.mean[in],STAT.sigma[in]);
2796  }
2797  }
2798  for(il=1; il<NET.Nlayer-1; il++)
2799  {
2800  fprintf(W,"C\n");
2801  fprintf(W,"C layer %d\n",il+1);
2802  for(in=0; in<NET.Nneur[il]; in++)
2803  {
2804  fprintf(W," RIN%d = %e\n",in+1,
2805  (double) NET.Weights[il][in][0]);
2806  for(jn=1;jn<=NET.Nneur[il-1]; jn++)
2807  fprintf(W," > +(%e) * OUT%d\n",
2808  (double) NET.Weights[il][in][jn],jn);
2809  }
2810  fprintf(W,"C\n");
2811  for(in=0; in<NET.Nneur[il]; in++)
2812  {
2813  if(NET.T_func[il][in]==0)
2814  {
2815  fprintf(W," OUT%d = 0\n",in+1);
2816  }
2817  else if(NET.T_func[il][in]==1)
2818  {
2819  fprintf(W," OUT%d = RIN%d\n",in+1,in+1);
2820  }
2821  else if(NET.T_func[il][in]==2)
2822  {
2823  fprintf(W," OUT%d = SIGMOID(RIN%d)\n",
2824  in+1,in+1);
2825  }
2826  }
2827  }
2828  il = NET.Nlayer-1;
2829  fprintf(W,"C\n");
2830  fprintf(W,"C layer %d\n",il+1);
2831  for(in=0; in<NET.Nneur[il]; in++)
2832  {
2833  fprintf(W," RIN%d = %e\n",in+1,
2834  (double) NET.Weights[il][in][0]);
2835  for(jn=1;jn<=NET.Nneur[il-1]; jn++)
2836  fprintf(W," > +(%e) * OUT%d\n",
2837  (double) NET.Weights[il][in][jn],jn);
2838  }
2839  fprintf(W,"C\n");
2840  for(in=0; in<NET.Nneur[il]; in++)
2841  {
2842  if(NET.T_func[il][in]==0)
2843  {
2844  fprintf(W," ROUT(%d) = 0\n",in+1);
2845  }
2846  else if(NET.T_func[il][in]==1)
2847  {
2848  fprintf(W," ROUT(%d) = RIN%d\n",in+1,in+1);
2849  }
2850  else if(NET.T_func[il][in]==2)
2851  {
2852  fprintf(W," ROUT(%d) = SIGMOID(RIN%d)\n",
2853  in+1,in+1);
2854  }
2855  }
2856 
2857  fprintf(W,"C\n");
2858  fprintf(W," END\n");
2859  fprintf(W," REAL FUNCTION SIGMOID(X)\n");
2860  fprintf(W," SIGMOID = 1./(1.+EXP(-X))\n");
2861  fprintf(W," END\n");
2862 
2863  fclose(W);
2864  return 0;
2865 }
2866 
2867 
2868 /***********************************************************/
2869 /* MLP_PrCFun */
2870 /* */
2871 /* writes the MLP function to file as a C function */
2872 /* */
2873 /* inputs : char *filename = name of the output file */
2874 /* */
2875 /* return value (int) = 0: no error */
2876 /* -1: could not open file */
2877 /* */
2878 /* Author: J.Schwindling 23-Apr-99 */
2879 /* Modified: J.Schwindling 30-Nov-99 return code */
2880 /***********************************************************/
2881 
2882 /* extern "C"Dllexport */
2884 {
2885  int il,in,jn;
2886  FILE *W;
2887 
2888  W=fopen(filename,"w");
2889  if(W==0) return -1;
2890 
2891  fprintf(W,"double sigmoid(double x)\n");
2892  fprintf(W,"{\n");
2893  fprintf(W,"return 1/(1+exp(-x));\n");
2894  fprintf(W,"}\n");
2895  fprintf(W,"void rnnfun(double *rin,double *rout)\n");
2896  fprintf(W,"{\n");
2897  fprintf(W," double out1[%d];\n",NET.Nneur[0]);
2898  fprintf(W," double out2[%d];\n",NET.Nneur[1]);
2899  if(NET.Nlayer>=3) fprintf(W," double out3[%d];\n",NET.Nneur[2]);
2900  if(NET.Nlayer>=4) fprintf(W," double out4[%d];\n",NET.Nneur[3]);
2901  fprintf(W,"\n");
2902 
2903  for(in=0; in<NET.Nneur[0]; in++)
2904  {
2905  if(DIVERS.Norm==0)
2906  {
2907  fprintf(W," out1[%d] = rin[%d];\n",in,in);
2908  }
2909  else
2910  {
2911  fprintf(W," out1[%d] = (rin[%d]-%e)/%e;\n",
2912  in,in,
2913  STAT.mean[in],STAT.sigma[in]);
2914  }
2915  }
2916 
2917  for(il=1; il<=NET.Nlayer-1; il++)
2918  {
2919  fprintf(W,"\n");
2920  fprintf(W,"/* layer %d */\n",il+1);
2921  for(in=0; in<NET.Nneur[il]; in++)
2922  {
2923  fprintf(W," out%d[%d] = %e\n",il+1,in,
2924  (double) NET.Weights[il][in][0]);
2925  for(jn=1;jn<=NET.Nneur[il-1]; jn++)
2926  fprintf(W," +(%e) * out%d[%d]\n",
2927  (double) NET.Weights[il][in][jn],il,jn-1);
2928  fprintf(W," ;\n");
2929  }
2930  fprintf(W,"\n");
2931  for(in=0; in<NET.Nneur[il]; in++)
2932  {
2933  if(NET.T_func[il][in]==0)
2934  {
2935  fprintf(W," out%d[%d] = 0;\n",il+1,in);
2936  }
2937  else if(NET.T_func[il][in]==1)
2938  {
2939  }
2940  else if(NET.T_func[il][in]==2)
2941  {
2942  fprintf(W," out%d[%d] = sigmoid(out%d[%d]);\n",
2943  il+1,in,il+1,in);
2944  }
2945  }
2946  }
2947  il = NET.Nlayer-1;
2948  for(in=0; in<NET.Nneur[il]; in++)
2949  {
2950  fprintf(W," rout[%d] = out%d[%d];\n",in,il+1,in);
2951  }
2952  fprintf(W,"}\n");
2953  fclose(W);
2954  return 0;
2955 }
2956 
2957 
2958 /***********************************************************/
2959 /* SaveWeights */
2960 /* */
2961 /* writes the weights to file */
2962 /* */
2963 /* inputs : char *filename = name of the output file */
2964 /* int iepoch = epoch number */
2965 /* */
2966 /* return value (int): 0 if OK */
2967 /* -1 if file could not be opened */
2968 /* */
2969 /* Author: J.Schwindling 20-Apr-99 */
2970 /* Modified: J.Schwindling 11-Jun-99 */
2971 /* print net structure in header */
2972 /* Modified: J.Schwindling 05-Nov-99 */
2973 /* return error code */
2974 /***********************************************************/
2975 
2976 /* extern "C"Dllexport */
2977 int SaveWeights(char *filename, int iepoch)
2978 {
2979  FILE *W;
2980  int ilayer,ineur,i;
2981 
2982  W=fopen(filename,"w");
2983  if(W==0) return -1;
2984 
2985  fprintf(W,"# network structure ");
2986  for(ilayer=0; ilayer<NET.Nlayer; ilayer++)
2987  {
2988  fprintf(W,"%d ",NET.Nneur[ilayer]);
2989  }
2990 
2991  fprintf(W,"\n %d\n",iepoch);
2992  for(ilayer=1; ilayer<NET.Nlayer; ilayer++)
2993  {
2994  for(ineur=0; ineur<NET.Nneur[ilayer]; ineur++)
2995  {
2996  for(i=0; i<=NET.Nneur[ilayer-1]; i++)
2997  {
2998  fprintf(W," %1.15e\n",
2999  (double) NET.Weights[ilayer][ineur][i]);
3000  }
3001  }
3002  }
3003  fclose(W);
3004  return 0;
3005 }
3006 
3007 
3008 /***********************************************************/
3009 /* LoadWeights */
3010 /* */
3011 /* reads the weights from file */
3012 /* */
3013 /* input : char *filename = name of the input file */
3014 /* output : int *iepoch = epoch number */
3015 /* */
3016 /* return value (int): 0 if OK */
3017 /* -1 if file could not be opened */
3018 /* */
3019 /* Author: J.Schwindling 20-Apr-99 */
3020 /* Modified: J.Schwindling 11-Jun-99 */
3021 /* lines starting with # are skipped */
3022 /* Modified: J.Schwindling 05-Nov-99 */
3023 /* return error code */
3024 /***********************************************************/
3025 
3026 /* extern "C"Dllexport */
3027 int LoadWeights(char *filename, int *iepoch)
3028 {
3029  FILE *W;
3030  int ilayer,ineur,i;
3031  double p;
3032  char s[80];
3033 
3034  W=fopen(filename,"r");
3035  if(W==0) return -1;
3036  do
3037  {
3038  fgets(s,80,W);
3039  }
3040  while(*s == '#');
3041  sscanf(s," %d",iepoch);
3042  for(ilayer=1; ilayer<NET.Nlayer; ilayer++)
3043  {
3044  for(ineur=0; ineur<NET.Nneur[ilayer]; ineur++)
3045  {
3046  for(i=0; i<=NET.Nneur[ilayer-1]; i++)
3047  {
3048  fscanf(W," %le",&p);
3049  NET.Weights[ilayer][ineur][i] = (dbl) p;
3050  }
3051  }
3052  }
3053 
3054  fclose(W);
3055  return 0;
3056 }
3057 
3058 
3059 /***********************************************************/
3060 /* AllocPatterns */
3061 /* */
3062 /* allocate memory for the examples */
3063 /* */
3064 /* input : int ifile = file number (0 or 1) */
3065 /* int npat = number of examples */
3066 /* int nin = number of inputs */
3067 /* int nout = number of outputs */
3068 /* int iadd = 0: new examples */
3069 /* 1: add to existing ones */
3070 /* */
3071 /* return value (int) = error code: 0 = no error */
3072 /* 1 = wrong file number */
3073 /* -111 = no memory */
3074 /* */
3075 /* Author: J.Schwindling 21-Apr-99 */
3076 /* Modified: J.Schwindling 26-Apr-99 */
3077 /* - frees memory if already booked and iadd=0 */
3078 /* (should remove need to call mlpfree) */
3079 /* - implement iadd = 1 */
3080 /***********************************************************/
3081 
3082 /* extern "C"Dllexport */int AllocPatterns(int ifile, int npat, int nin, int nout, int iadd)
3083 {
3084  int j, ierr;
3085  type_pat *tmp, *tmp3;
3086  type_pat **tmp2;
3087  int ntot;
3088 
3089  if(ifile>1 || ifile<0) return(1);
3090 /* scanf("%d",&j); */
3091  if(ExamplesMemory==0)
3092  {
3093  ExamplesMemory=1;
3094  PAT.Pond = (type_pat **) malloc(2*sizeof(dbl*));
3095  PAT.Rin = (type_pat***) malloc(2*sizeof(type_pat**));
3096  PAT.Rans = (type_pat***) malloc(2*sizeof(type_pat**));
3097  PAT.vRin = (type_pat**) malloc(2*sizeof(type_pat*));
3098  if(PAT.Pond == 0 || PAT.Rin == 0
3099  || PAT.Rans == 0 || PAT.vRin == 0) return -111;
3100  }
3101 
3102 
3103 /* if iadd=0, check that memory not already allocated. Otherwise free it */
3104  if(iadd==0 && PatMemory[ifile]!=0)
3105  {
3106  ierr = FreePatterns(ifile);
3107  }
3108 
3109 /* allocate memory and initialize ponderations */
3110  if(iadd==0 || PatMemory[ifile]==0)
3111  {
3112  PatMemory[ifile] = 1;
3113  PAT.Pond[ifile] = (type_pat*) malloc(npat*sizeof(type_pat));
3114  if(PAT.Pond[ifile] == 0) return -111;
3115  for(j=0; j<npat; j++)
3116  PAT.Pond[ifile][j] = 1;
3117 
3118  PAT.Rin[ifile] = (type_pat**) malloc(npat*sizeof(type_pat*));
3119  if(PAT.Rin[ifile] == 0) return -111;
3120  PAT.Rans[ifile] = (type_pat**) malloc(npat*sizeof(type_pat*));
3121  if(PAT.Rans[ifile] == 0) return -111;
3122 
3123  PAT.vRin[ifile] = (type_pat *) malloc(npat*(nin+1)*
3124  sizeof(type_pat));
3125  if(PAT.vRin[ifile] == 0) return -111;
3126 
3127  for(j=0; j<npat; j++)
3128  {
3129  PAT.Rin[ifile][j] = &(PAT.vRin[ifile][j*(nin+1)+1]);
3130  PAT.vRin[ifile][j*(nin+1)] = 1;
3131  }
3132  for(j=0; j<npat; j++)
3133  {
3134  PAT.Rans[ifile][j] = (type_pat*) malloc(nout*sizeof(type_pat));
3135  if(PAT.Rans[ifile][j] == 0) return -111;
3136  }
3137  PAT.Npat[ifile] = npat;
3138 
3139  if(ifile==0)
3140  {
3141  ExamplesIndex = (int *) malloc(npat*sizeof(int));
3142  if(ExamplesIndex == 0) return -111;
3143  for(j=0; j<npat; j++) ExamplesIndex[j] = j;
3144  }
3145  }
3146  else /* add examples */
3147  {
3148  ntot = PAT.Npat[ifile]+npat;
3149 
3150 /* event weighting */
3151  tmp = (type_pat *) malloc(ntot*sizeof(type_pat));
3152  if(tmp == 0) return -111;
3153 
3154  for(j=0; j<PAT.Npat[ifile]; j++)
3155  {
3156  tmp[j] = PAT.Pond[ifile][j];
3157  }
3158  for(j=PAT.Npat[ifile];j<ntot;j++)
3159  {
3160  tmp[j] = 1;
3161  }
3162  if(PatMemory[ifile]==1) free(PAT.Pond[ifile]);
3163  PAT.Pond[ifile] = tmp;
3164 
3165 /* examples */
3166 /* tmp2 = (type_pat **) malloc(ntot*sizeof(type_pat*));
3167  for(j=0; j<PAT.Npat[ifile]; j++)
3168  {
3169  tmp2[j] = PAT.Rin[ifile][j];
3170  }
3171  for(j=PAT.Npat[ifile];j<ntot;j++)
3172  {
3173  tmp2[j] = (type_pat*) malloc(nin*sizeof(type_pat));
3174  }
3175  if(PatMemory[ifile]==1) free(PAT.Rin[ifile]);
3176  PAT.Rin[ifile] = tmp2; */
3177 
3178  tmp3 = (type_pat *) malloc(ntot*(nin+1)*sizeof(type_pat));
3179  if(tmp3 == 0) return -111;
3180 
3181  for(j=0; j<PAT.Npat[ifile]*(nin+1); j++)
3182  {
3183  tmp3[j] = PAT.vRin[ifile][j];
3184  }
3185  if(PatMemory[ifile]==1) free(PAT.vRin[ifile]);
3186  PAT.vRin[ifile] = tmp3;
3187  for(j=0; j<ntot; j++)
3188  {
3189  PAT.Rin[ifile][j] = &(PAT.vRin[ifile][j*(nin+1)+1]);
3190  PAT.vRin[ifile][j*(nin+1)] = 1;
3191  }
3192 
3193  tmp2 = (type_pat **) malloc(ntot*sizeof(type_pat*));
3194  if(tmp2 == 0) return -111;
3195  for(j=0; j<PAT.Npat[ifile]; j++)
3196  {
3197  tmp2[j] = PAT.Rans[ifile][j];
3198  }
3199  for(j=PAT.Npat[ifile];j<ntot;j++)
3200  {
3201  tmp2[j] = (type_pat*) malloc(nout*sizeof(type_pat));
3202  if(tmp2[j] == 0) return -111;
3203  }
3204  if(PatMemory[ifile]==1) free(PAT.Rans[ifile]);
3205  PAT.Rans[ifile] = tmp2;
3206  PAT.Npat[ifile] = ntot;
3207  PatMemory[ifile] = 1;
3208 
3209 /* indices */
3210  if(ifile==0)
3211  {
3212  free(ExamplesIndex);
3213  ExamplesIndex = (int *) malloc(ntot*sizeof(int));
3214  if(ExamplesIndex == 0) return -111;
3215  for(j=0; j<ntot; j++) ExamplesIndex[j] = j;
3216  }
3217  }
3218 
3219  return 0;
3220 }
3221 
3222 
3223 /***********************************************************/
3224 /* FreePatterns */
3225 /* */
3226 /* frees memory for the examples */
3227 /* */
3228 /* input : int ifile = file number (0 or 1) */
3229 /* */
3230 /* return value (int) = error code: 0 = no error */
3231 /* 1 = wrong file number */
3232 /* 2 = no mem allocated */
3233 /* */
3234 /* Author: J.Schwindling 26-Apr-99 */
3235 /***********************************************************/
3236 
3237 /* extern "C"Dllexport */int FreePatterns(int ifile)
3238 {
3239  int i;
3240 
3241  if(ifile>1 || ifile<0) return 1;
3242 /* printf("%d %d \n",ifile,PatMemory[ifile]);*/
3243  if(PatMemory[ifile]==0) return 2;
3244 
3245  free(PAT.Pond[ifile]);
3246  for(i=0; i<PAT.Npat[ifile]; i++)
3247  {
3248 /* free(PAT.Rin[ifile][i]); */
3249  free(PAT.Rans[ifile][i]);
3250  }
3251  free(PAT.Rin[ifile]);
3252  free(PAT.Rans[ifile]);
3253  free(PAT.vRin[ifile]);
3254  PatMemory[ifile] = 0;
3255  PAT.Npat[ifile] = 0;
3256 
3257  return 0;
3258 }
3259 
3260 
3261 /***********************************************************/
3262 /* MLP_StatInputs */
3263 /* */
3264 /* compute statistics about the inputs: mean, RMS, min, max*/
3265 /* */
3266 /* inputs: int Nexamples = number of examples */
3267 /* int Niputs = number of quantities */
3268 /* type_pat **inputs = input values */
3269 /* */
3270 /* outputs: dbl *mean = mean value */
3271 /* dbl *sigma = RMS */
3272 /* dbl *minimum = minimum */
3273 /* dbl *maximum = maximum */
3274 /* */
3275 /* return value (int): always = 0 */
3276 /* */
3277 /* Author: J.Schwindling 11-Oct-99 */
3278 /***********************************************************/
3279 
3280 int MLP_StatInputs(int Nexamples, int Ninputs, type_pat **inputs,
3281  dbl *mean, dbl *sigma, dbl *minimum, dbl *maximum)
3282 {
3283  dbl *fmean;
3284  int j, ipat, nmax;
3285 
3286 /* first compute a fast rough mean using the first 100 events */
3287  fmean = (dbl*) malloc(Ninputs*sizeof(dbl));
3288  nmax = 100;
3289  if(Nexamples<100) nmax=Nexamples;
3290 
3291  for(j=0;j<Ninputs;j++)
3292  {
3293  fmean[j] = 0;
3294  for(ipat=0;ipat<nmax;ipat++)
3295  {
3296  fmean[j] += (dbl) inputs[ipat][j];
3297  }
3298  fmean[j] = fmean[j]/(dbl) nmax;
3299 
3300 /* now compute real mean and sigma, min and max */
3301  mean[j] = 0;
3302  sigma[j] = 0;
3303  minimum[j] = 99999;
3304  maximum[j] = -99999;
3305  for(ipat=0;ipat<Nexamples;ipat++)
3306  {
3307  mean[j] += (dbl) inputs[ipat][j];
3308  sigma[j] += ((dbl) inputs[ipat][j]-fmean[j])*
3309  ((dbl) inputs[ipat][j]-fmean[j]);
3310  if((dbl) inputs[ipat][j] > maximum[j])
3311  maximum[j]=(dbl) inputs[ipat][j];
3312  if((dbl) inputs[ipat][j] < minimum[j])
3313  minimum[j]=(dbl) inputs[ipat][j];
3314  }
3315  mean[j] = mean[j]/(dbl) Nexamples;
3316  sigma[j] = sqrt(sigma[j]/ (dbl) Nexamples -
3317  (mean[j]-fmean[j])*
3318  (mean[j]-fmean[j]));
3319  }
3320  free(fmean);
3321  return 0;
3322 }
3323 
3324 /***********************************************************/
3325 /* MLP_PrintInputStat */
3326 /* */
3327 /* prints statistics about the inputs: mean, RMS, min, max */
3328 /* */
3329 /* return value (int) = error code: 0 = OK */
3330 /* -111 = could not */
3331 /* allocate memory */
3332 /* */
3333 /* Author: J.Schwindling 11-Oct-99 */
3334 /* Modified: J.Schwindling 31-Jan-2000: return value */
3335 /***********************************************************/
3336 
3338 {
3339  int j;
3340  dbl *mean, *sigma, *minimum, *maximum;
3341 
3342 /* allocate memory */
3343  mean = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3344  sigma = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3345  minimum = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3346  maximum = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3347 
3348  if(mean == 0 || sigma == 0 || minimum == 0
3349  || maximum == 0) return -111;
3350 
3351  MLP_StatInputs(PAT.Npat[0],NET.Nneur[0],PAT.Rin[0],
3352  mean,sigma,minimum,maximum);
3353 
3354  printf("\t mean \t\t RMS \t\t min \t\t max\n");
3355  for(j=0;j<NET.Nneur[0];j++)
3356  {
3357  printf("var%d \t %e \t %e \t %e \t %e\n",j+1,
3358  mean[j],sigma[j],minimum[j],maximum[j]);
3359  }
3360 
3361  free(mean);
3362  free(sigma);
3363  free(minimum);
3364  free(maximum);
3365  printf("\n");
3366  return 0;
3367 }
3368 
3369 
3370 /***********************************************************/
3371 /* NormalizeInputs */
3372 /* */
3373 /* normalize the inputs: I' = (I - <I>) / RMS(I) */
3374 /* */
3375 /* return value (int) = error code: 0 = OK */
3376 /* -111 = could not */
3377 /* allocate memory */
3378 /* */
3379 /* Author: J.Schwindling 04-May-1999 */
3380 /* Modified: J.Schwindling 31-Jan-2000: return value */
3381 /***********************************************************/
3382 
3383 /* extern "C"Dllexport */int NormalizeInputs()
3384 {
3385  int j, ipat;
3386  dbl *mean, *sigma, *minimum, *maximum;
3387 
3388 /* allocate memory */
3389  mean = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3390  sigma = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3391  STAT.mean = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3392  STAT.sigma = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3393  minimum = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3394  maximum = (dbl *) malloc(NET.Nneur[0]*sizeof(dbl));
3395 
3396  if(mean == 0 || sigma == 0 || minimum == 0
3397  || maximum == 0 || STAT.mean == 0 ||
3398  STAT.sigma == 0) return -111;
3399 
3400  MLP_StatInputs(PAT.Npat[0],NET.Nneur[0],PAT.Rin[0],
3401  mean,sigma,minimum,maximum);
3402 
3403  if(NET.Debug>=1) printf("\t mean \t\t RMS \t\t min \t\t max\n");
3404  for(j=0;j<NET.Nneur[0];j++)
3405  {
3406  if(NET.Debug>=1)
3407  printf("var%d \t %e \t %e \t %e \t %e\n",j+1,
3408  mean[j],sigma[j],minimum[j],maximum[j]);
3409 
3410 /* store mean and sigma for output function */
3411  STAT.mean[j] = mean[j];
3412  STAT.sigma[j] = sigma[j];
3413 
3414 /* finally apply the normalization */
3415  for(ipat=0;ipat<PAT.Npat[0];ipat++)
3416  {
3417  PAT.Rin[0][ipat][j] =
3418  (PAT.Rin[0][ipat][j]-(float) mean[j])/
3419  (float) sigma[j];
3420  }
3421  for(ipat=0;ipat<PAT.Npat[1];ipat++)
3422  {
3423  PAT.Rin[1][ipat][j] =
3424  (PAT.Rin[1][ipat][j]-(float) mean[j])/
3425  (float) sigma[j];
3426  }
3427  }
3428 
3429  free(mean);
3430  free(sigma);
3431  free(minimum);
3432  free(maximum);
3433  if(NET.Debug>=1) printf("\n");
3434  return 0;
3435 }
3436 
3437 
3438 /***********************************************************/
3439 /* AllocNetwork */
3440 /* */
3441 /* memory allocation for weights, etc */
3442 /* */
3443 /* inputs: int Nlayer: number of layers */
3444 /* int *Neurons: nulber of neurons per layer */
3445 /* */
3446 /* return value (int): error = 0: no error */
3447 /* = -111: could not allocate mem*/
3448 /* */
3449 /* Author: J.Schwindling 28-Sep-99 */
3450 /***********************************************************/
3451 
3452 int AllocNetwork(int Nlayer, int *Neurons)
3453 {
3454  int i, j, k, l;
3455 
3456  if(NetMemory != 0) FreeNetwork();
3457  NetMemory = 1;
3458 
3459  NET.Nneur = (int *) malloc(Nlayer*sizeof(int));
3460  if(NET.Nneur == 0) return -111;
3461 
3462  NET.T_func = (int **) malloc(Nlayer*sizeof(int *));
3463  NET.Deriv1 = (dbl **) malloc(Nlayer*sizeof(dbl *));
3464  NET.Inn = (dbl **) malloc(Nlayer*sizeof(dbl *));
3465  NET.Outn = (dbl **) malloc(Nlayer*sizeof(dbl *));
3466  NET.Delta = (dbl **) malloc(Nlayer*sizeof(dbl *));
3467  if(NET.T_func == 0 || NET.Deriv1 == 0
3468  || NET.Inn == 0 || NET.Outn == 0
3469  || NET.Delta == 0) return -111;
3470 
3471  for(i=0; i<Nlayer; i++)
3472  {
3473  NET.T_func[i] = (int *) malloc(Neurons[i]*sizeof(int));
3474  NET.Deriv1[i] = (dbl *) malloc(Neurons[i]*sizeof(dbl));
3475  NET.Inn[i] = (dbl *) malloc(Neurons[i]*sizeof(dbl));
3476  NET.Outn[i] = (dbl *) malloc(Neurons[i]*sizeof(dbl));
3477  NET.Delta[i] = (dbl *) malloc(Neurons[i]*sizeof(dbl));
3478  if(NET.T_func[i] == 0 || NET.Deriv1[i] == 0
3479  || NET.Inn[i] == 0 || NET.Outn[i] == 0
3480  || NET.Delta[i] ==0 ) return -111;
3481  }
3482 
3483  NET.Weights = (dbl ***) malloc(Nlayer*sizeof(dbl **));
3484  NET.vWeights = (dbl **) malloc(Nlayer*sizeof(dbl *));
3485  LEARN.Odw = (dbl ***) malloc(Nlayer*sizeof(dbl **));
3486  LEARN.ODeDw = (dbl ***) malloc(Nlayer*sizeof(dbl **));
3487  LEARN.DeDw = (dbl ***) malloc(Nlayer*sizeof(dbl **));
3488  if(NET.Weights == 0 || NET.vWeights == 0
3489  || LEARN.Odw == 0 || LEARN.ODeDw == 0
3490  || LEARN.DeDw == 0) return -111;
3491 
3492  for(i=1; i<Nlayer; i++)
3493  {
3494  k = Neurons[i-1]+1;
3495  NET.vWeights[i] = (dbl *) malloc(k * Neurons[i] *
3496  sizeof(dbl));
3497  NET.Weights[i] = (dbl **) malloc(Neurons[i]*sizeof(dbl *));
3498  LEARN.Odw[i] = (dbl **) malloc(Neurons[i]*sizeof(dbl *));
3499  LEARN.ODeDw[i] = (dbl **) malloc(Neurons[i]*sizeof(dbl *));
3500  LEARN.DeDw[i] = (dbl **) malloc(Neurons[i]*sizeof(dbl *));
3501  if(NET.Weights[i] == 0 || NET.vWeights[i] == 0
3502  || LEARN.Odw[i] == 0 || LEARN.ODeDw[i] == 0
3503  || LEARN.DeDw[i] == 0) return -111;
3504 
3505  for(j=0; j<Neurons[i]; j++)
3506  {
3507  NET.Weights[i][j] = &(NET.vWeights[i][j*k]);
3508  LEARN.Odw[i][j] = (dbl *) malloc(k*sizeof(dbl));
3509  LEARN.ODeDw[i][j] = (dbl *) malloc(k*sizeof(dbl));
3510  LEARN.DeDw[i][j] = (dbl *) malloc(k*sizeof(dbl));
3511  if(LEARN.Odw[i][j] == 0
3512  || LEARN.ODeDw[i][j] == 0
3513  || LEARN.DeDw[i][j] == 0) return -111;
3514 
3515  for(l=0; l<k; l++)
3516  {
3517  LEARN.Odw[i][j][l] = 0;
3518  LEARN.ODeDw[i][j][l] = 0;
3519  }
3520  }
3521  }
3522  return 0;
3523 }
3524 
3525 
3526 /***********************************************************/
3527 /* FreeNetwork */
3528 /* */
3529 /* frees the memory allocated for the network */
3530 /* */
3531 /* Author: J.Schwindling 06-Oct-99 */
3532 /***********************************************************/
3533 
3535 {
3536  int i, j;
3537  for(i=1; i<NET.Nlayer; i++)
3538  {
3539  for(j=0; j<NET.Nneur[i]; j++)
3540  {
3541 /* free(NET.Weights[i][j]); */
3542  free(LEARN.Odw[i][j]);
3543  free(LEARN.ODeDw[i][j]);
3544  free(LEARN.DeDw[i][j]);
3545  }
3546  free(NET.vWeights[i]);
3547  free(NET.Weights[i]);
3548  free(LEARN.Odw[i]);
3549  free(LEARN.ODeDw[i]);
3550  free(LEARN.DeDw[i]);
3551  }
3552  free(NET.Weights);
3553  free(LEARN.Odw);
3554  free(LEARN.ODeDw);
3555  free(LEARN.DeDw);
3556 
3557  free(NET.Nneur);
3558 
3559  for(i=0; i<NET.Nlayer; i++)
3560  {
3561  free(NET.T_func[i]);
3562  free(NET.Deriv1[i]);
3563  free(NET.Inn[i]);
3564  free(NET.Outn[i]);
3565  free(NET.Delta[i]);
3566  }
3567  free(NET.T_func);
3568  free(NET.Deriv1);
3569  free(NET.Inn);
3570  free(NET.Outn);
3571  free(NET.Delta);
3572 
3573  NetMemory = 0;
3574 }
3575 
3576 /***********************************************************/
3577 /* GetNetStructure */
3578 /* */
3579 /* given a strinng like "3,4,1" returns the network */
3580 /* structure */
3581 /* */
3582 /* inputs: char *s: input string */
3583 /* */
3584 /* outputs: int *Nlayer: number of layers */
3585 /* int *Nneur: number of neurons per layer */
3586 /* */
3587 /* return value (int): error = 0: no error */
3588 /* = -1: s is empty */
3589 /* = -2: s is too long */
3590 /* = -3: too many layers */
3591 /* */
3592 /* Author: J.Schwindling 04-Oct-99 */
3593 /***********************************************************/
3594 
3595 int GetNetStructure(char *s, int *Nlayer, int *Nneur)
3596 {
3597  int i=0;
3598  char tmp[1024];
3599 
3600  if(strlen(s)==0) return -1;
3601  if(strlen(s)>1024) return -2;
3602 
3603  strcpy(tmp,s);
3604  if (strtok(tmp,","))
3605  {
3606  i=1;
3607  while (strtok(NULL,",")) i++;
3608  }
3609  *Nlayer = i;
3610  if(i > NLMAX) return -3;
3611 
3612  strcpy(tmp,s);
3613  if (*Nlayer>0)
3614  {
3615  sscanf(strtok(tmp,","),"%d",&(Nneur[0]));
3616  for (i=1;i<*Nlayer;i++)
3617  sscanf(strtok(NULL,","),"%d",&(Nneur[i]));
3618  }
3619 
3620  return 0;
3621 }
3622 
3623 
3624 /***********************************************************/
3625 /* MLP_SetNet */
3626 /* */
3627 /* to set the structure of a neural network */
3628 /* inputs: int *nl = number of layers */
3629 /* int *nn = number of neurons */
3630 /* */
3631 /* return value (int) = error value: */
3632 /* 0: no error */
3633 /* 1: N layers > NLMAX */
3634 /* 2: N layers < 2 */
3635 /* -111: Error allocating memory */
3636 /* */
3637 /* Author: J.Schwindling 14-Apr-99 */
3638 /* Modified: J.Schwindling 05-Oct-99 allocate memory */
3639 /* Modified: J.Schwindling 29-Nov-99 LearnFree, LearnAlloc */
3640 /***********************************************************/
3641 
3642 int MLP_SetNet(int *nl, int *nn)
3643 {
3644  int il,ierr;
3645 
3646  if((*nl)>NLMAX) return(1);
3647  if((*nl)<2) return(2);
3648 
3649 /* LearnFree(); */
3650 /* allocate memory */
3651  ierr = AllocNetwork(*nl,nn);
3652  if(ierr != 0) return ierr;
3653 
3654 /* set number of layers */
3655  NET.Nlayer = (int) *nl;
3656 
3657 /* set number of neurons */
3658  for(il=0; il<NET.Nlayer; il++) {
3659  NET.Nneur[il] = nn[il];
3660  }
3661 
3662 /* set transfer functions */
3663  SetDefaultFuncs();
3664 /* LearnAlloc(); */
3665 
3666  return(0);
3667 }
3668 
3669 
3670 /***********************************************************/
3671 /* MLP_MatrixVectorBias */
3672 /* */
3673 /* computes a Matrix-Vector product */
3674 /* r[j] = M[j][0] + Sum_i M[j][i] v[i] */
3675 /* */
3676 /* inputs: dbl *M = matrix (n lines, m+1 columns) */
3677 /* dbl *v = vector (dimension m) */
3678 /* dbl *r = resulting vector (dimension n) */
3679 /* int n */
3680 /* int m */
3681 /* */
3682 /* Author: J.Schwindling 24-Jan-00 */
3683 /***********************************************************/
3684 
3685 void MLP_MatrixVectorBias(dbl *M, dbl *v, dbl *r, int n, int m)
3686 {
3687  int i,j;
3688  register dbl a1, a2, a3, a4, c, d;
3689  dbl *pM1 = M;
3690  dbl *pM2 = &(M[m+1]);
3691  dbl *pM3 = &(M[2*(m+1)]);
3692  dbl *pM4 = &(M[3*(m+1)]);
3693  dbl *pr = r;
3694  int mp1 = m+1;
3695 
3696  for(i=0; i<n-3;
3697  i+=4, pM1 += 3*mp1, pM2 += 3*mp1, pM3 += 3*mp1, pM4 += 3*mp1,
3698  pr+=4)
3699  {
3700  a1 = *pM1;
3701  a2 = *pM2;
3702  a3 = *pM3;
3703  a4 = *pM4;
3704  pM1++; pM2++; pM3++; pM4++;
3705  for(j=0; j<m-1; j+=2, pM1+=2, pM2+=2, pM3+=2, pM4+=2)
3706  {
3707  c = v[j];
3708  d = v[j+1];
3709  a1 = a1 + *pM1 * c + *(pM1+1) * d;
3710  a2 = a2 + *pM2 * c + *(pM2+1) * d;
3711  a3 = a3 + *pM3 * c + *(pM3+1) * d;
3712  a4 = a4 + *pM4 * c + *(pM4+1) * d;
3713  }
3714  for(j=j; j<m; j++, pM1++, pM2++, pM3++, pM4++)
3715  {
3716  c = v[j];
3717  a1 = a1 + *pM1 * c;
3718  a2 = a2 + *pM2 * c;
3719  a3 = a3 + *pM3 * c;
3720  a4 = a4 + *pM4 * c;
3721  }
3722  *pr = a1; *(pr+1) = a2; *(pr+2) = a3; *(pr+3) = a4;
3723  }
3724  for(i=i; i<n; i++)
3725  {
3726  pM1 = &(M[i*(m+1)]);
3727  a1 = *pM1;
3728  pM1++;
3729  for(j=0; j<m; j++, pM1++)
3730  {
3731  a1 = a1 + *pM1 * v[j];
3732  }
3733  r[i] = a1;
3734  }
3735 }
3736 /***********************************************************/
3737 /* MLP_MatrixVector */
3738 /* */
3739 /* computes a Matrix-Vector product */
3740 /* r[j] = Sum_i M[j][i] v[i] */
3741 /* */
3742 /* inputs: dbl *M = matrix (n lines, m+1 columns) */
3743 /* dbl *v = vector (dimension m) */
3744 /* dbl *r = resulting vector (dimension n) */
3745 /* int n */
3746 /* int m */
3747 /* */
3748 /* Author: J.Schwindling 24-Jan-00 */
3749 /***********************************************************/
3750 
3751 void MLP_MatrixVector(dbl *M, type_pat *v, dbl *r, int n, int m)
3752 {
3753  int i,j;
3754  register dbl a1, a2, a3, a4, c, d;
3755  dbl *pM1 = M;
3756  dbl *pM2 = &(M[m]);
3757  dbl *pM3 = &(M[2*m]);
3758  dbl *pM4 = &(M[3*m]);
3759  dbl *pr = r;
3760  int mp1 = m;
3761 
3762  for(i=0; i<n-3;
3763  i+=4, pM1 += 3*mp1, pM2 += 3*mp1, pM3 += 3*mp1, pM4 += 3*mp1,
3764  pr+=4)
3765  {
3766  a1 = 0;
3767  a2 = 0;
3768  a3 = 0;
3769  a4 = 0;
3770  for(j=0; j<m-1; j+=2, pM1+=2, pM2+=2, pM3+=2, pM4+=2)
3771  {
3772  c = v[j];
3773  d = v[j+1];
3774  a1 = a1 + *pM1 * c + *(pM1+1) * d;
3775  a2 = a2 + *pM2 * c + *(pM2+1) * d;
3776  a3 = a3 + *pM3 * c + *(pM3+1) * d;
3777  a4 = a4 + *pM4 * c + *(pM4+1) * d;
3778  }
3779  for(j=j; j<m; j++, pM1++, pM2++, pM3++, pM4++)
3780  {
3781  c = v[j];
3782  a1 = a1 + *pM1 * c;
3783  a2 = a2 + *pM2 * c;
3784  a3 = a3 + *pM3 * c;
3785  a4 = a4 + *pM4 * c;
3786  }
3787  *pr = a1; *(pr+1) = a2; *(pr+2) = a3; *(pr+3) = a4;
3788  }
3789  for(i=i; i<n; i++)
3790  {
3791  pM1 = &(M[i*m]);
3792  a1 = 0;
3793  for(j=0; j<m; j++, pM1++)
3794  {
3795  a1 = a1 + *pM1 * v[j];
3796  }
3797  r[i] = a1;
3798  }
3799 }
3800 
3801 
3802 /***********************************************************/
3803 /* MLP_MM2rows */
3804 /* */
3805 /* computes a Matrix-Matrix product, with the first matrix */
3806 /* having 2 rows */
3807 /* */
3808 /* inputs: dbl *c = resulting matrix (Nj * Nk) */
3809 /* dbl *a = first matrix (Ni * Nj) */
3810 /* dbl *b = second matrix (Nj * Nk) */
3811 /* int Ni */
3812 /* int Nj */
3813 /* int Nk */
3814 /* int NaOffs */
3815 /* int NbOffs */
3816 /* */
3817 /* Author: J.Schwindling 24-Jan-00 */
3818 /***********************************************************/
3819 
3821  int Ni, int Nj, int Nk, int NaOffs, int NbOffs)
3822 {
3823 //int i,j,k;
3824 int j,k;
3825 dbl s00,s01,s10,s11;
3826 type_pat *pa0,*pa1;
3827 dbl *pb0,*pb1,*pc0,*pc1;
3828 
3829  for (j=0; j<=Nj-2; j+=2)
3830  {
3831  pc0 = c+j;
3832  pc1 = c+j+Nj;
3833  s00 = 0.0; s01 = 0.0; s10 = 0.0; s11 = 0.0;
3834 
3835  for (k=0,pb0=b+k+NbOffs*j,
3836  pb1=b+k+NbOffs*(j+1),
3837  pa0=a+k,
3838  pa1=a+k+NaOffs;
3839  k<Nk;
3840  k++,pa0++,
3841  pa1++,
3842  pb0++,
3843  pb1++)
3844  {
3845  s00 += (*pa0)*(*pb0);
3846  s01 += (*pa0)*(*pb1);
3847  s10 += (*pa1)*(*pb0);
3848  s11 += (*pa1)*(*pb1);
3849  }
3850  *pc0 = s00; *(pc0+1) = s01; *pc1 = s10; *(pc1+1) = s11;
3851  }
3852  for (j=j; j<Nj; j++)
3853  {
3854  pc0 = c+j;
3855  pc1 = c+j+Nj;
3856  s00 = 0.0; s10 = 0.0;
3857  for (k=0,pb0=b+k+NbOffs*j,
3858  pa0=a+k,
3859  pa1=a+k+NaOffs;
3860  k<Nk;
3861  k++,pa0++,
3862  pa1++,
3863  pb0++)
3864  {
3865  s00 += (*pa0)*(*pb0);
3866  s10 += (*pa1)*(*pb0);
3867  }
3868  *pc0 = s00; *pc1 = s10;
3869  }
3870 }
3871 
3872 #ifdef __cplusplus
3873 } // extern "C"
3874 #endif
const double beta
dbl * delta
Definition: mlp_gen.cc:36
dbl MLP_Epoch(int iepoch, dbl *alpmin, int *Ntest)
Definition: mlp_gen.cc:710
int i
Definition: DBlmapReader.cc:9
int MessLang
Definition: mlp_gen.cc:22
float alpha
Definition: AMPTWrapper.h:95
Definition: mlp_gen.h:56
void getnLexemes(int n, char *s, char **ss)
Definition: mlp_gen.cc:2606
int JacobianMemory
Definition: mlp_gen.cc:28
Definition: mlp_gen.h:28
dbl MLP_Test_MM(int ifile, dbl *tmp)
Definition: mlp_gen.cc:266
void LearnFree()
Definition: mlp_gen.cc:2645
dbl DeDwProd()
Definition: mlp_gen.cc:1025
int ReadPatterns(char *filename, int ifile, int *inet, int *ilearn, int *iexamples)
Definition: mlp_gen.cc:2234
void CGDir(dbl beta)
Definition: mlp_gen.cc:1271
int FixedStep(dbl alpha)
Definition: mlp_gen.cc:1709
void DeDwSaveZero()
Definition: mlp_gen.cc:1102
int DecreaseSearch(dbl *alpmin, int *Ntest, dbl Err0)
Definition: mlp_gen.cc:1616
void InitWeights()
Definition: mlp_gen.cc:2152
Definition: mlp_gen.h:38
int OutputWeights
Definition: mlp_gen.cc:23
int * ExamplesIndex
Definition: mlp_gen.cc:40
long int integer
Definition: mlp_lapack.h:12
int SaveWeights(char *filename, int iepoch)
Definition: mlp_gen.cc:2977
#define NULL
Definition: scimark2.h:8
int NLineSearchFail
Definition: mlp_gen.cc:33
void MLP_Out(type_pat *rrin, dbl *rrout)
Definition: mlp_gen.cc:63
double doublereal
Definition: mlp_lapack.h:17
#define NLMAX
Definition: mlp_gen.cc:14
#define STAT
Definition: mlp_gen.h:60
TRandom random
Definition: MVATrainer.cc:138
int MLP_SetNet(int *nl, int *nn)
Definition: mlp_gen.cc:3642
int LearnMemory
Definition: mlp_gen.cc:29
T eta() const
void EtaDecay()
Definition: mlp_gen.cc:2091
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
#define PAT
Definition: mlp_gen.h:45
dbl ** BFGSH
Definition: mlp_gen.cc:37
int MLP_PrCFun(char *filename)
Definition: mlp_gen.cc:2883
tuple s2
Definition: indexGen.py:106
float MLPfitVersion
Definition: mlp_gen.cc:31
#define NET
Definition: mlp_gen.h:25
int MLP_StatInputs(int Nexamples, int Ninputs, type_pat **inputs, dbl *mean, dbl *sigma, dbl *minimum, dbl *maximum)
Definition: mlp_gen.cc:3280
void SetLambda(double Wmax)
Definition: mlp_gen.cc:1954
void SteepestDir()
Definition: mlp_gen.cc:1249
int GetBFGSH(int Nweights)
Definition: mlp_gen.cc:1423
int BFGSMemory
Definition: mlp_gen.cc:27
int GetNetStructure(char *s, int *Nlayer, int *Nneur)
Definition: mlp_gen.cc:3595
void PrintWeights()
Definition: mlp_gen.cc:2173
void MLP_Out_T(type_pat *rrin)
Definition: mlp_gen.cc:115
int StochStep()
Definition: mlp_gen.cc:970
int LearnAlloc()
Definition: mlp_gen.cc:2690
dbl MLP_Stochastic()
Definition: mlp_gen.cc:521
T sqrt(T t)
Definition: SSEVec.h:28
void MLP_vSigmoideDeriv(dbl *x, dbl *dy, int n)
double type_pat
Definition: mlp_gen.h:13
int MLP_Train(int *ipat, dbl *err)
Definition: mlp_gen.cc:900
int ShuffleExamples(int n, int *index)
Definition: mlp_gen.cc:2107
double MLP_Rand(dbl mini, dbl maxi)
Definition: mlp_gen.cc:2136
dbl DeDwNorm()
Definition: mlp_gen.cc:1003
int DeDwSum(type_pat *ans, dbl *out, int ipat)
Definition: mlp_gen.cc:1128
dbl LastAlpha
Definition: mlp_gen.cc:32
int j
Definition: DBlmapReader.cc:9
int CountLexemes(char *s)
Definition: mlp_gen.cc:2590
struct net_ net_ MLP_HIDDEN
Definition: mlp_gen.cc:16
int NetMemory
Definition: mlp_gen.cc:30
int LoadWeights(char *filename, int *iepoch)
Definition: mlp_gen.cc:3027
void MLP_LineHyb(dbl ***w0, dbl alpha)
Definition: mlp_gen.cc:1928
Definition: mlp_gen.h:16
int AllocPatterns(int ifile, int npat, int nin, int nout, int iadd)
Definition: mlp_gen.cc:3082
dbl ** JacobianMatrix
Definition: mlp_gen.cc:39
#define DIVERS
Definition: mlp_gen.h:54
tuple prod
Definition: CrabTask.py:87
int k[5][pyjets_maxn]
int SetTransFunc(int layer, int neuron, int func)
Definition: mlp_gen.cc:1205
tuple out
Definition: dbtoconf.py:99
int StochStepHyb()
Definition: mlp_gen.cc:931
void InitBFGSH(int Nweights)
Definition: mlp_gen.cc:1393
void GetGammaDelta()
Definition: mlp_gen.cc:1318
int PatMemory[2]
Definition: mlp_gen.cc:26
void MLP_ResLin()
Definition: mlp_gen.cc:1976
void MLP_MM2rows(dbl *c, type_pat *a, dbl *b, int Ni, int Nj, int Nk, int NaOffs, int NbOffs)
Definition: mlp_gen.cc:3820
void MLP_MatrixVector(dbl *M, type_pat *v, dbl *r, int n, int m)
Definition: mlp_gen.cc:3751
void DeDwScale(int Nexamples)
Definition: mlp_gen.cc:1065
void MLP_Line(dbl ***w0, dbl alpha)
Definition: mlp_gen.cc:1762
void MLP_MatrixVectorBias(dbl *M, dbl *v, dbl *r, int n, int m)
Definition: mlp_gen.cc:3685
void getLexemes(char *s, char **ss)
Definition: mlp_gen.cc:2620
#define LEARN
Definition: mlp_gen.h:36
int dgels_(char *trans, integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *work, integer *lwork, integer *info)
double b
Definition: hdecay.h:120
void MLP_vSigmoide(dbl *x, int n)
void FreeNetwork()
Definition: mlp_gen.cc:3534
int NormalizeInputs()
Definition: mlp_gen.cc:3383
int ExamplesMemory
Definition: mlp_gen.cc:24
int LineSearch(dbl *alpmin, int *Ntest, dbl Err0)
Definition: mlp_gen.cc:1478
void DeDwSave()
Definition: mlp_gen.cc:1083
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
void BFGSdir(int Nweights)
Definition: mlp_gen.cc:1345
double a
Definition: hdecay.h:121
tuple filename
Definition: lut2db_cfg.py:20
dbl * Gamma
Definition: mlp_gen.cc:38
double dbl
Definition: mlp_gen.h:12
dbl MLP_Sigmoide(dbl x)
Definition: mlp_sigmoide.cc:27
int WeightsMemory
Definition: mlp_gen.cc:25
int MLP_PrintInputStat()
Definition: mlp_gen.cc:3337
#define CLEN
Definition: mlp_gen.cc:2231
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
void MLP_Out2(type_pat *rrin)
Definition: mlp_gen.cc:186
void SetDefaultFuncs()
Definition: mlp_gen.cc:1228
tuple ifile
Definition: indexGen.py:77
int FreePatterns(int ifile)
Definition: mlp_gen.cc:3237
dbl MLP_Test(int ifile, int regul)
Definition: mlp_gen.cc:449
int MLP_PrFFun(char *filename)
Definition: mlp_gen.cc:2774
void DeDwZero()
Definition: mlp_gen.cc:1046
int AllocNetwork(int Nlayer, int *Neurons)
Definition: mlp_gen.cc:3452
mathSSE::Vec4< T > v
dbl ** Hessian
Definition: mlp_gen.cc:41
int LineSearchHyb(dbl *alpmin, int *Ntest)
Definition: mlp_gen.cc:1786
dbl DerivDir()
Definition: mlp_gen.cc:1293