home *** CD-ROM | disk | FTP | other *** search
/ PC World Komputer 1997 February / PCWK0297.iso / technika / nnmodel / neural.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  34.8 KB  |  1,331 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <malloc.h>
  5. #include "nndefs.h"
  6. #include "datamat.h"
  7. #include "params.h"
  8. #include "neural.h"
  9.  
  10.     float gasdev(int *idum, float dzdiv);
  11.     int sign (float x, float y);
  12.     float sign_of_arg (float x);
  13.     float Max (float a, float b);
  14.     float Random(float randz);
  15.     void Logit(const char* fmt, ...);
  16.     static char FMT_E[] = "%e ";
  17.  
  18.  
  19. NEURAL *NCreateNeural(  )
  20. {   
  21.     NEURAL *pN;
  22.     pN = (NEURAL*) malloc (sizeof(NEURAL));
  23.     
  24.     pN->m_version = REVLEVEL;
  25.     pN->m_istate = pN->m_ninputs = pN->m_noutputs = pN->m_nhidden = 0;
  26.     pN->m_params = NULL;
  27.     pN->m_sumerr2 = 0;
  28.     pN->m_cnt = 0;
  29.     pN->m_itmax = 2;
  30.     pN->m_dm=NULL;
  31.  
  32.     pN->m_ni    = NULL;
  33.     pN->m_nout    = NULL;
  34.     pN->m_hinputw    = NULL;
  35.     pN->m_houtputv    = NULL;
  36.     pN->m_htheta    = NULL;
  37.     pN->m_oinputw    = NULL;
  38.     pN->m_otheta    = NULL;
  39.  
  40.     pN->m_olastoutputv = NULL;
  41.     pN->m_hlastoutputv = NULL; 
  42.  
  43. #ifdef VERBOSE
  44.     Logit ("CNeural null constructor\n");
  45. #endif
  46.     return pN;
  47.  
  48. }
  49.  
  50. void NDeleteNeural(NEURAL *pN )
  51. {
  52. #ifdef VERBOSE
  53.     Logit ("CNeural destructor\n");
  54. #endif  
  55.  
  56.     if (pN->m_istate&NN_PARAMSLOADED) free (pN->m_params);
  57.  
  58.     if (pN->m_noutputs) {
  59. #ifdef VERBOSE
  60.         Logit ("   Destroying in=%d hid=%d out=%d\n",pN->m_ninputs,pN->m_nhidden,pN->m_noutputs);
  61. #endif
  62.         free_2d_floats (pN->m_hinputw, pN->m_nhidden);
  63.         free_2d_floats (pN->m_oinputw, pN->m_noutputs);
  64. #ifdef NN_STHRU
  65.         free_2d_floats (pN->m_iinputw, pN->m_noutputs);
  66. #endif
  67.         free (pN->m_ni);
  68.         free (pN->m_houtputv);
  69.         free (pN->m_htheta);
  70.         free (pN->m_otheta);
  71.         free (pN->m_nout);
  72.         free (pN->m_olastoutputv);
  73.     }
  74.     if (pN->m_istate&NN_DYNAMIC) {
  75.         if(pN->m_noutputs) {
  76. #ifdef VERBOSE
  77.             Logit ("   Destroy dynamic\n");
  78. #endif
  79.             free_2d_floats (pN->m_hlastdelta, pN->m_nhidden);
  80.             free_2d_floats (pN->m_olastdelta,pN->m_noutputs);
  81. #ifdef NN_STHRU
  82.             free_2d_floats (pN->m_ilastdelta,pN->m_noutputs);
  83. #endif
  84.             free (pN->m_hlastvar);
  85.             free (pN->m_hlearn);
  86.             free (pN->m_htlearn);
  87.             free (pN->m_olastvar);
  88.             free (pN->m_otraining);
  89.             free (pN->m_olearn);
  90.             free (pN->m_otlearn);
  91.             free (pN->m_hlastoutputv);
  92.         }
  93.     }
  94.     if (pN->m_istate&NN_DATAMAT) {
  95. #ifdef VERBOSE
  96.         Logit ("   Destroy DM\n");
  97. #endif
  98.         free (pN->m_dm);
  99.     }
  100. }
  101.  
  102. // add another element to a vector
  103.  
  104. float *ReallocVector(int num, float *arptr, float init) 
  105. {
  106.     float *temp;
  107.     int i;
  108.  
  109.     temp = malloc (sizeof(float)*(num+1));
  110.     for (i=0;i<num;i++) temp[i] = arptr[i];
  111.     free (arptr);
  112.     temp[num] = init;
  113.     return temp;
  114. }
  115.  
  116. // Add another column to a matrix
  117.  
  118. float **ReallocFloats1(int first, int second, float **arptr, float ran) 
  119. {
  120.     float **temp;
  121.     int i,j;
  122.  
  123.     temp = alloc_2d_floats (first+1, second);
  124.     for (i=0; i < first; i++)
  125.         for (j=0; j < second; j++)
  126.             temp[i][j] = arptr[i][j];
  127.     for (j=0; j < second; j++)
  128.         temp[first][j] = Random(ran);
  129.     free_2d_floats (arptr, first);
  130.     return temp;
  131. }
  132.  
  133. // Add another row to a matrix
  134.  
  135. float **ReallocFloats2(int first, int second, float **arptr, float ran) 
  136. {
  137.     float **temp;
  138.     int i,j;
  139.  
  140.     temp = alloc_2d_floats (first, second+1);
  141.     for (i=0; i < first; i++)
  142.         for (j=0; j < second; j++)
  143.             temp[i][j] = arptr[i][j];
  144.     for (j=0; j < first; j++)
  145.         temp[j][second] = Random(ran);
  146.     free_2d_floats (arptr, first);
  147.     return temp;
  148. }
  149.  
  150. // Add another hidden neuron to a dynamic learning neural net
  151.  
  152. void NAddHidden(NEURAL *pN)
  153. {
  154.     int oldhid;
  155.  
  156.     if (pN->m_nhidden >= pN->m_params->m_AImaxhid) return;
  157.     oldhid = pN->m_nhidden;
  158.     pN->m_nhidden++;
  159. #ifdef VERBOSE
  160.     Logit ("CNeural AddHidden Adding hidden neuron now=%d\n",pN->m_nhidden);
  161. #endif
  162.  
  163.     pN->m_hinputw    = ReallocFloats1(oldhid, pN->m_ninputs, pN->m_hinputw,pN->m_params->m_randz);
  164.     pN->m_oinputw    = ReallocFloats2(pN->m_noutputs, oldhid, pN->m_oinputw,pN->m_params->m_randz);
  165.     pN->m_hlastdelta = ReallocFloats1(oldhid, pN->m_ninputs, pN->m_hlastdelta,pN->m_params->m_randz);
  166.     pN->m_olastdelta = ReallocFloats2(pN->m_noutputs, oldhid, pN->m_olastdelta,pN->m_params->m_randz);
  167. #ifdef NN_STHRU
  168.     pN->m_ilastdelta = ReallocFloats2(pN->m_noutputs, pN->m_ninputs, pN->m_ilastdelta,pN->m_params->m_randz);
  169. #endif
  170.  
  171.     pN->m_houtputv  =  ReallocVector(oldhid, pN->m_houtputv, 0.0f);
  172.     pN->m_htheta    =  ReallocVector(oldhid, pN->m_htheta,Random(pN->m_params->m_randz) );
  173.     pN->m_hlastvar  =  ReallocVector(oldhid, pN->m_hlastvar, 0.0f);
  174.     pN->m_hlearn    =  ReallocVector(oldhid, pN->m_hlearn,pN->m_params->m_Hlearning_rate);
  175.     pN->m_htlearn   =  ReallocVector(oldhid, pN->m_htlearn,pN->m_params->m_tlearning_rate);
  176. }
  177.  
  178. float NGetROutput(NEURAL *pN, const int neuron)
  179. {
  180.     return (DRescale(pN->m_dm, pN->m_nout[neuron],'O',neuron));
  181. }
  182.  
  183. float NGetRInput(NEURAL *pN, const int neuron)
  184. {
  185.     return (DRescale(pN->m_dm, pN->m_ni[neuron],'I',neuron));
  186. }
  187.  
  188. void NSetRInput(NEURAL *pN, const int neuron, float f)
  189. {
  190.     pN->m_ni[neuron] = DScale(pN->m_dm, f,'I',neuron);
  191. }
  192.  
  193. // Attach the Design Matrix to the Dynamic Neural Net
  194.  
  195. void NSetDM( NEURAL *pN, DATAMAT *a)
  196. {
  197.     pN->m_dm = a;
  198.     pN->m_istate |= NN_DATAMAT;
  199. }
  200.  
  201.  
  202. // Back Propagate output error to weights
  203.  
  204. int NBackProp1(NEURAL *pN, int cnt) 
  205. {
  206.     int i,j,k;
  207.     static int idummy = -13;
  208.     float deltaw,ERR,VAR,sum1;
  209.  
  210. // calculate the outputs of the hidden layer 
  211.                            
  212.     for (i=0;i<pN->m_ninputs;i++) pN->m_ni[i] = pN->m_dm->m_iarray[i][cnt];
  213. // add gaussian noise
  214.     if (pN->m_params->m_inrandzdiv!=0.0) 
  215.         for (i=0;i<pN->m_ninputs;i++) pN->m_ni[i] *= (1+gasdev(&idummy,pN->m_params->m_inrandzdiv));
  216.  
  217.     for (i=0; i < pN->m_nhidden; i++) {
  218.         for (sum1=0,j=0; j < pN->m_ninputs; j++)
  219.             sum1 += (pN->m_ni[j] * pN->m_hinputw[i][j] );
  220.         sum1 += pN->m_htheta[i];
  221.         pN->m_houtputv[i] = afunc (sum1);
  222.     }
  223.  
  224. // connect hidden layer 1 to output neurons 
  225.  
  226.     for (j=0; j < pN->m_noutputs; j++) {
  227.         for (sum1=0,i=0; i < pN->m_nhidden; i++)
  228.             sum1 += (pN->m_houtputv[i] * pN->m_oinputw[j][i]);
  229. #ifdef NN_STHRU
  230.         if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) {
  231.             for (i=0; i < pN->m_ninputs; i++)
  232.                 sum1 += (pN->m_ni[i] * pN->m_iinputw[j][i]);
  233.         }
  234. #endif  
  235.         sum1 += pN->m_otheta[j];
  236.         pN->m_nout[j] = afunc (sum1);
  237.     }
  238.  
  239. // BP Output Layer to hidden
  240.     for (i=0; i < pN->m_noutputs; i++) {
  241.         if (pN->m_params->m_inrandzdiv==0.0) 
  242.             ERR = (pN->m_dm->m_oarray[i][cnt] - pN->m_nout[i]);
  243. // add gaussian noise
  244.         else
  245.             ERR =    (pN->m_dm->m_oarray[i][cnt] * 
  246.                     (1+gasdev(&idummy,pN->m_params->m_inrandzdiv))
  247.                     - pN->m_nout[i]);
  248.     
  249.         pN->m_sumerr2 += ERR*ERR;
  250.         VAR = ERR * dafunc (pN->m_nout[i]);
  251.         pN->m_olastvar[i] = VAR;
  252.         for (j=0; j < pN->m_nhidden; j++) {
  253.             deltaw = (pN->m_olearn[i] * VAR * pN->m_houtputv[j]) +
  254.                  (pN->m_params->m_alpha * pN->m_olastdelta[i][j]);
  255.             pN->m_olastdelta[i][j] = deltaw;
  256.             pN->m_oinputw[i][j] += deltaw;
  257.         }
  258. #ifdef NN_STHRU
  259.         if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) {
  260.             for (j=0; j < pN->m_ninputs; j++) {
  261.                 deltaw = (pN->m_olearn[i] * VAR * pN->m_ni[j] * pN->m_params->m_inoutlearn) +
  262.                      (pN->m_params->m_alpha * pN->m_ilastdelta[i][j]);
  263.                 pN->m_ilastdelta[i][j] = deltaw;
  264.                 pN->m_iinputw[i][j] += deltaw;
  265.             }
  266.         }
  267. #endif
  268.         pN->m_otheta[i] += (pN->m_otlearn[i] * VAR);
  269.     }
  270. // BP Hidden Layer to Input
  271.     for (j=0; j < pN->m_nhidden; j++) {
  272.         for (VAR = 0.f,i=0; i < pN->m_noutputs; i++)
  273.             VAR += (pN->m_olastvar[i] * pN->m_oinputw[i][j]);
  274.         VAR *= dafunc(pN->m_houtputv[j]);
  275.         for(k=0;k<pN->m_ninputs;k++) {
  276.             deltaw = (pN->m_hlearn[j] * VAR * pN->m_ni[k]) +
  277.                  (pN->m_params->m_alpha * pN->m_hlastdelta[j][k]);
  278.             pN->m_hlastdelta[j][k] = deltaw;
  279.             pN->m_hinputw[j][k] += deltaw;
  280.         }
  281. //        pN->m_hlastvar[j] = VAR; /* only if more hidden layers */
  282.         pN->m_htheta[j] += (pN->m_htlearn[j] * VAR);
  283.     }
  284.     return 0;
  285. }
  286.  
  287. void NFeedForward(NEURAL *pN)
  288. {
  289.     int i,j;
  290.     float sum1;
  291.  
  292. /* calculate the outputs of the hidden layer */
  293.                            
  294.     for (i=0; i < pN->m_nhidden; i++) {
  295.         for (sum1=0.f,j=0; j < pN->m_ninputs; j++)
  296.             sum1 += (pN->m_ni[j] * pN->m_hinputw[i][j] );
  297.         sum1 += pN->m_htheta[i];
  298.         pN->m_houtputv[i] = afunc (sum1);
  299.     }
  300.  
  301. /* connect hidden layer 1 to output neurons */
  302.  
  303.     for (j=0; j < pN->m_noutputs; j++) {
  304.         for (sum1=0.f,i=0; i < pN->m_nhidden; i++)
  305.             sum1 += (pN->m_houtputv[i] * pN->m_oinputw[j][i]);
  306. #ifdef NN_STHRU
  307.         if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) {
  308.             for (i=0; i < pN->m_ninputs; i++)
  309.                 sum1 += (pN->m_ni[i] * pN->m_iinputw[j][i]);
  310.         }
  311. #endif
  312.         sum1 += pN->m_otheta[j];
  313.         pN->m_nout[j] = afunc (sum1);
  314.     }
  315. }
  316.  
  317.  
  318. void NClearDelta(NEURAL *pN)
  319. {
  320.     int i,j;
  321.  
  322.     for (i=0;i<pN->m_nhidden;i++)
  323.         for (j=0;j<pN->m_ninputs;j++) pN->m_hlastdelta[i][j] = 0.f;
  324.     for (i=0;i<pN->m_noutputs;i++)
  325.         for (j=0;j<pN->m_nhidden;j++) pN->m_olastdelta[i][j] = 0.f;
  326. #ifdef NN_STHRU
  327.     for (i=0;i<pN->m_noutputs;i++)
  328.         for (j=0;j<pN->m_ninputs;j++) pN->m_ilastdelta[i][j] = 0.f;
  329. #endif
  330. }
  331.  
  332. void NInitializeNetwork(NEURAL *pN)
  333. {
  334.     int i,n;
  335.  
  336. /*  randonize weights */
  337.  
  338.     // de-allocate extra hidden units
  339.     free_2d_floats (pN->m_hinputw, pN->m_nhidden);
  340.     free_2d_floats (pN->m_oinputw, pN->m_noutputs);
  341.     free_2d_floats (pN->m_hlastdelta, pN->m_nhidden);
  342.     free_2d_floats (pN->m_olastdelta, pN->m_noutputs);
  343. #ifdef NN_STHRU
  344.     free_2d_floats (pN->m_ilastdelta, pN->m_noutputs);
  345. #endif
  346.                 
  347.     free (pN->m_houtputv);
  348.     free (pN->m_htheta);
  349.     free (pN->m_hlastvar);
  350.     free (pN->m_hlearn);
  351.     free (pN->m_htlearn);
  352.  
  353.     if (pN->m_params->m_goodness==0) pN->m_nhidden = pN->m_params->m_AImaxhid;
  354.     else pN->m_nhidden=1;
  355.             
  356.     pN->m_hinputw = alloc_2d_floats (pN->m_nhidden, pN->m_ninputs);
  357.     pN->m_oinputw = alloc_2d_floats (pN->m_noutputs, pN->m_nhidden);
  358.     pN->m_hlastdelta = alloc_2d_floats (pN->m_nhidden,pN->m_ninputs);
  359.     pN->m_olastdelta = alloc_2d_floats (pN->m_noutputs,pN->m_nhidden);
  360. #ifdef NN_STHRU
  361.     pN->m_ilastdelta = alloc_2d_floats (pN->m_noutputs,pN->m_ninputs);
  362. #endif
  363.  
  364.     pN->m_houtputv = (float *) malloc (sizeof(float)*pN->m_nhidden);
  365.     pN->m_htheta = (float *) malloc (sizeof(float)*pN->m_nhidden);
  366.     pN->m_hlastvar = (float *) malloc (sizeof(float)*pN->m_nhidden);
  367.     pN->m_hlearn = (float *) malloc (sizeof(float)*pN->m_nhidden);
  368.     pN->m_htlearn = (float *) malloc (sizeof(float)*pN->m_nhidden);
  369.               
  370.     pN->m_sumerr2=0;
  371.     pN->m_cnt = 0;
  372.  
  373.     srand(pN->m_params->m_seed);
  374.  
  375.     for (n=0; n < pN->m_ninputs; n++) pN->m_ni[n] = 0.f;
  376.     for (n=0; n < pN->m_noutputs; n++) pN->m_nout[n] = 0.f;
  377.     for (i=0; i < pN->m_nhidden; i++) {
  378.         for (n=0; n < pN->m_ninputs; n++) {
  379.             pN->m_hinputw[i][n] = Random(pN->m_params->m_randz);
  380.             pN->m_hlastdelta[i][n] = 0.f;
  381.         }
  382.         pN->m_htheta[i] = Random(pN->m_params->m_randz);
  383.         pN->m_htlearn[i] = pN->m_params->m_tlearning_rate;
  384.         pN->m_hlearn[i]  = pN->m_params->m_Hlearning_rate;
  385.         pN->m_houtputv[i]  = pN->m_hlastvar[i]  = 0.f;
  386.     }
  387.     for (i=0; i < pN->m_noutputs; i++) {
  388.         for (n=0; n < pN->m_nhidden; n++) {
  389.             pN->m_oinputw[i][n] = Random(pN->m_params->m_randz);
  390.             pN->m_olastdelta[i][n] = 0.f;
  391.         }
  392.         pN->m_otheta[i] = Random(pN->m_params->m_randz);
  393.         pN->m_otlearn[i] = pN->m_params->m_tlearning_rate;
  394.         pN->m_olearn[i] = pN->m_params->m_learning_rate;
  395.         pN->m_olastvar[i]  = 0.f;
  396. #ifdef NN_STHRU
  397.         for (n=0; n < pN->m_ninputs; n++) {
  398.             pN->m_iinputw[i][n] = Random(pN->m_params->m_randz);
  399.             pN->m_ilastdelta[i][n] = 0.f;
  400.         }
  401. #endif
  402.     }
  403. //              save_weights();
  404.     return;
  405. }
  406.                     
  407. // returns 1 when finished else 0
  408.                     
  409. int NQuickTrain(NEURAL *pN, int mode,int increment)
  410. {
  411.  
  412.     static int train_cnt;
  413.     int dograph;
  414.     int lastrow;
  415.  
  416.     if (mode==-1) { // init                                
  417.         NAI(pN,-1,0);    // init AI routine
  418.         train_cnt = 0;
  419.         NClearDelta(pN);
  420.         pN->m_sumerr2 = 0;
  421.         return 0;
  422.     }
  423.     pN->m_sumerr2 = 0;
  424.     dograph=0;
  425.     lastrow = pN->m_dm->m_numrows;
  426.  
  427. redo:
  428.  
  429. //    NewTrainingSet(train_cnt,0);
  430.  
  431.     if (train_cnt == lastrow) {
  432.         ++pN->m_cnt;
  433.         --increment;
  434.  
  435.         if (pN->m_cnt >= pN->m_params->m_cnt_max) {
  436. #ifdef VERBOSE
  437.             Logit("Done Training\n");
  438. #endif
  439.             NAI(pN, -2,0);    // deinit AI routine 
  440.             return 1;
  441.         }
  442.         if ((pN->m_cnt % pN->m_params->m_eon) == 0) {
  443.             if (NAI(pN, 0,0)) {    // adjust_learning_heuristics
  444.                 NAI(pN, -2,0);    // deinit AI routine
  445.                 return 1;
  446.             }
  447.             return 2;
  448.         }
  449.  
  450.         train_cnt = 0;
  451.         if (!increment) return dograph;
  452.         pN->m_sumerr2 = 0;
  453.         goto redo;
  454.     }
  455.     NBackProp1(pN,train_cnt); // Back Propagate error to weights
  456.     train_cnt++;
  457.     goto redo;
  458.  
  459. }
  460.  
  461. void NNewTrainingSet(NEURAL *pN, int t,int flag)
  462. {
  463.     int i,k;
  464.     static int idummy = -13;
  465.  
  466.     if (!flag) {
  467.         for (k=0;k<pN->m_ninputs;k++) pN->m_ni[k] = pN->m_dm->m_iarray[k][t];
  468.         for (i=0;i<pN->m_noutputs;i++) pN->m_otraining[i] = pN->m_dm->m_oarray[i][t];
  469.     } else {
  470.         for (k=0;k<pN->m_ninputs;k++) pN->m_ni[k] = pN->m_dm->m_itarray[k][t];
  471.         for (i=0;i<pN->m_noutputs;i++) pN->m_otraining[i] = pN->m_dm->m_otarray[i][t];
  472.     }
  473.  
  474.     return;
  475. }
  476.  
  477. /* Incremental learning adjust_learning_heuristics */
  478. int NAI(NEURAL *pN, int flag,int a) 
  479. {
  480. //      extern int keytran;
  481. //      extern int ninp,noutp;
  482.     static float RSQ,NWT,LASTRSQ;
  483.     static int nomorehidinc,nmin,nmax;
  484.     static lasthid;
  485.     static float lastmin,lastmax,min,max;
  486.     float calc_rsquare();
  487.     int i;
  488.     static int first;
  489.     static float sumdiv,pers;
  490.     static float lastsum;
  491.     static int lastnmin,lastnmax;
  492.     char cmin,cmax,cnmin,cnmax,csum;
  493.  
  494.     if (flag==-1) {  /* init  */
  495.         LASTRSQ = RSQ = lastmin = lastmax = lastsum = 0;
  496.         lastnmin = lastnmax = 0;
  497.         nomorehidinc=0;
  498.         lasthid=0;
  499.         first=1;
  500.         if (pN->m_params->m_eon<0) pN->m_params->m_eon = 10;
  501.         i = (int) ((pN->m_params->m_cnt_max-pN->m_cnt)/pN->m_params->m_eon);
  502.         if (i < 1) i = 2;
  503. //                      if (!a) InitGraph("AI TRAINING",6, i+1,0.0,1.0);
  504.  
  505. // The following is for the first graph
  506.  
  507.         RSQ = NCalcRsquare(pN);
  508.         NWT = NCheckTol(pN, &min,&max,&nmin,&nmax);
  509.         pN->m_stats   [0] = 1;
  510.         pN->m_stats   [1] = min;
  511.         pN->m_stats   [2] = max;
  512.         pN->m_stats   [3] = (float) nmin;
  513.         pN->m_stats   [4] = (float) nmax;
  514.         pN->m_stats   [5] = RSQ;
  515. #ifdef VERBOSE
  516.         Logit ("sumerr=%f min=%f max=%f nmin=%d nmax=%d rsq=%f\n",pN->m_sumerr2,min,max,nmin,nmax,RSQ);
  517. #endif
  518.         return 0;
  519.     }
  520.     if (flag==-2) {  /* deinit  */
  521. //                      if (!a) CloseGraph();
  522.         return 0;
  523.     }
  524.     lasthid++;
  525.  
  526. /* different goodness algorithms
  527.    0 = no measurement
  528.    1 = inc even thru cnt_max
  529.    2 = same a 3 manual inc
  530.    3 = inc on factors
  531. */
  532.  
  533.     switch (pN->m_params->m_goodness) {
  534.     default:
  535.         return 0;
  536.         break;
  537.     case 0: /* old training method */
  538.     case 1: /* same as 3 even thru eons */
  539.     case 2: /* same as 3 but manual hid inc */
  540.     case 3: /* number within tol */
  541.         if (first) {sumdiv = pN->m_sumerr2; first=0;}
  542.         RSQ = NCalcRsquare(pN);
  543.         NWT = NCheckTol(pN,&min,&max,&nmin,&nmax);
  544.  
  545.         pN->m_stats   [0] = pN->m_sumerr2;
  546.         pN->m_stats   [1] = min;
  547.         pN->m_stats   [2] = max;
  548.         pN->m_stats   [3] = (float) nmin;
  549.         pN->m_stats   [4] = (float) nmax;
  550.         pN->m_stats   [5] = RSQ;
  551.  
  552.         if ((pN->m_params->m_TrainFlags & TF_STOPONTOL) && (NWT == pN->m_dm->m_numrows)) return 1;
  553.         if ((pN->m_params->m_TrainFlags & TF_STOPONERR) && (pN->m_sumerr2 < pN->m_params->m_errtol)) return 1;
  554.         if ((pN->m_params->m_TrainFlags & TF_STOPONRSQ) && (RSQ > pN->m_params->m_goodrsq)) return 1;
  555.  
  556.         if ((min+pN->m_params->m_nosigninc) < lastmin) cmin  ='y'; else cmin='n';
  557.         if ((max+pN->m_params->m_nosigninc) < lastmax) cmax  ='y'; else cmax='n';
  558.         if (nmin < lastnmin)   cnmin ='y'; else cnmin='n';
  559.         if (nmax > lastnmax)   cnmax ='y'; else cnmax='n';
  560.  
  561.         if (sumdiv==0.) pers=0;
  562.         else pers= (lastsum-pN->m_sumerr2)/sumdiv;
  563.         if (pers > pN->m_params->m_nosigninc) csum  ='y'; else csum='n';
  564.  
  565.         lastmin = min;
  566.         lastmax = max;
  567.         lastnmin = nmin;
  568.         lastnmax = nmax;
  569.         lastsum = pN->m_sumerr2;
  570.  
  571.         
  572. #ifdef VERBOSE
  573.         Logit ("NWT=%6.f RSq=%6f",NWT,RSQ);
  574.         Logit (" min=%6f max=%6f",min,max);
  575.         Logit (" nmin=%4d nmax=%4d",nmin,nmax);
  576.         Logit (" %c%c%c%c%c\n",cmin,cmax,cnmin,cnmax,csum);
  577. #endif
  578.         if ((pN->m_params->m_goodness==0) || (pN->m_params->m_goodness==2)) return 0;
  579.  
  580.         if (pN->m_params->m_goodness==1) {
  581.             if (lasthid<(pN->m_params->m_cnt_max/pN->m_params->m_AImaxhid/pN->m_params->m_eon)) return 0;
  582.             lasthid=0;
  583.             goto incrnow;
  584.         }              
  585.         if (pN->m_nhidden==1) goto incrnow;
  586.         if (!((cmin|cmax|cnmin|cnmax|csum)=='n')) return 0;
  587.         if ( (pN->m_params->m_cnt_max-pN->m_cnt)/pN->m_params->m_eon < pN->m_params->m_cnt_max/pN->m_params->m_eon/5  ) nomorehidinc=1;
  588.  
  589.         break;
  590.  
  591.     } /* end of goodness switch */
  592.  
  593.     if (nomorehidinc) return 0;
  594.     if (lasthid<(pN->m_params->m_cnt_max/pN->m_params->m_eon/20)) return 0;
  595.     lasthid=0;
  596.  
  597. incrnow:
  598.     if (pN->m_params->m_goodness==0) return 0;
  599.     if (pN->m_nhidden >= pN->m_params->m_AImaxhid) return 0;
  600.  
  601.     /* decrement learning rates */
  602.  
  603.     for (i=0; i < pN->m_nhidden; i++) {
  604.         pN->m_htlearn[i] -= (pN->m_htlearn[i] * pN->m_params->m_hiddegrad);
  605.         pN->m_hlearn[i]  -= (pN->m_hlearn[i] * pN->m_params->m_hiddegrad);
  606.     }
  607.  
  608.     /* create another hidden neuron */
  609.     NAddHidden(pN);
  610.     return 0; /* continue */
  611. }
  612.  
  613. float NCalcRsquare(NEURAL *pN )
  614. {
  615.     double resid,rrq,sum;
  616.     double sumres[32],sumobs[32];
  617.     double averes[32],aveobs[32];
  618.     double rq[32];
  619.     int i,loop;
  620.     int train_cnt;
  621.  
  622.     for (i=0;i<32;i++) sumres[i] = sumobs[i] = 0.;
  623.     train_cnt = 0;
  624.     loop=1;
  625.  
  626. top:
  627.     if (train_cnt == pN->m_dm->m_numrows) {
  628.         if (loop==1) {
  629.             for(i=0;i<pN->m_noutputs;i++) {
  630.                 aveobs[i] = sumobs[i]/(float) train_cnt;
  631.                 averes[i] = sumres[i]/(float) train_cnt;
  632.                 sumres[i] = sumobs[i] = 0.;
  633.             }
  634.             train_cnt = 0;
  635.             loop=2;
  636.         }
  637.         else {
  638.             for(i=0;i<pN->m_noutputs;i++) {
  639.                 rq[i] = 1-sumres[i]/sumobs[i];
  640.                 if (rq[i] < 0.0) rq[i] = 0.0;
  641.                 if (rq[i] > 1.0) rq[i] = 1.0;
  642.             }
  643.             sum = 0.0;
  644.             for(i=0;i<pN->m_noutputs;i++) sum += rq[i];
  645.             rrq = sum/(float) pN->m_noutputs;
  646.             return (float) rrq;
  647.         }
  648.     }
  649.     NNewTrainingSet(pN, train_cnt,0);
  650.  
  651.     NFeedForward(pN);
  652.  
  653.     for(i=0;i<pN->m_noutputs;i++) {
  654.         resid = (double) DGetOutputVal(pN->m_dm, train_cnt,i)-pN->m_nout[i];
  655.         if (loop==1) {
  656.             sumobs[i] += DGetOutputVal(pN->m_dm, train_cnt,i);/*observed*/
  657.             sumres[i]  += resid; /* model error */
  658.         }
  659.         else {
  660.             sumobs[i] += pow(DGetOutputVal(pN->m_dm, train_cnt,i)-aveobs[i],2); /*observed*/
  661.             sumres[i]  += pow(resid-averes[i],2); /* model error */
  662.         }
  663.     }
  664.     train_cnt++;
  665.     goto top;
  666. }
  667.  
  668. float NCheckTol(NEURAL *pN, float *min, float *max, int *nmin, int *nmax)
  669. {
  670.     int i,t,si;
  671.     float temp;
  672.  
  673.     t = 0;  *max = 0.f; *min = 0.f;
  674.     *nmin = 0; *nmax = 0;
  675.  
  676.     for (i=0; i < pN->m_dm->m_numrows;i++) {
  677.  
  678.         NNewTrainingSet(pN, i,0);
  679.  
  680.         NFeedForward(pN);   /* Feed Forward to output */
  681.  
  682.         for (si=0;si<pN->m_noutputs;si++) {
  683.             temp = pN->m_nout[si]-pN->m_otraining[si];
  684.             if (temp > *max) *max = temp;
  685.             if (temp < -*min) *min = -temp;
  686.             temp = (float) fabs(temp);
  687.             if (temp < pN->m_params->m_tol) t++;
  688.             if (pN->m_nout[si] > (pN->m_otraining[si]+pN->m_params->m_tol)) (*nmax)++;
  689.             if (pN->m_nout[si] < (pN->m_otraining[si]-pN->m_params->m_tol)) (*nmin)++;
  690.  
  691.         }
  692.     }
  693.     return (float)t;
  694. }
  695.  
  696.  
  697.  
  698. //#define itmax 200
  699. #define EPS 1.0e-10
  700.  
  701. void Nfrprmn(NEURAL *pN, float *p, float ftol, int *iter, float *fret)
  702. {
  703.     int j,its,n;
  704.     float gg,gam,fp,dgg;
  705.     
  706. //    Logit ("frprmn\n");
  707.  
  708.     n=pN->m_NDIM;
  709.  
  710.     fp=NErrorFunction(pN,p);
  711.     NforwardDIffGradient(pN,p,pN->m_xi);
  712.     for (j=0;j<n;j++) {
  713.         pN->m_g[j] = -pN->m_xi[j];
  714.         pN->m_xi[j]=pN->m_h[j]=pN->m_g[j];
  715.     }
  716.     for (its=1;its<=pN->m_itmax;its++) {
  717.         *iter=its;
  718.         Nlinmin(pN, p,pN->m_xi,fret);
  719.         if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
  720.             //Logit("TOL\n");
  721.             //FREEALL
  722.             return;
  723.         }
  724.         fp=NErrorFunction(pN,p);
  725.         NforwardDIffGradient(pN,p,pN->m_xi);
  726.         dgg=gg=0;
  727.         for (j=0;j<n;j++) {
  728.             gg += pN->m_g[j]*pN->m_g[j];
  729. //          dgg += pN->m_xi[j]*pN->m_xi[j];   
  730.             dgg += (pN->m_xi[j]+pN->m_g[j])*pN->m_xi[j];
  731.         }
  732.         if (gg == 0.0) {
  733.             //Logit("gg=0\n");
  734.             //FREEALL
  735.             return;
  736.         }
  737.         gam=dgg/gg;
  738.         for (j=0;j<n;j++) {
  739.             pN->m_g[j] = -pN->m_xi[j];
  740.             pN->m_xi[j]=pN->m_h[j]=pN->m_g[j]+gam*pN->m_h[j];
  741.         }
  742.     }
  743.     //nrerror("Too many iterations in FRPRMN");
  744. }
  745.  
  746. #undef EPS
  747.  
  748. #define ITMAX 100
  749. #define CGOLD 0.3819660f
  750. #define ZEPS 1.0e-10
  751. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  752. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  753.  
  754. float Nbrent(NEURAL *pN,float ax, float bx, float cx, float tol,
  755.         float *xmin)
  756. {
  757.     int iter;
  758.     float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
  759.     float e=0;
  760.  
  761. //    Logit ("brent\n");
  762.  
  763.     a=((ax < cx) ? ax : cx);
  764.     b=((ax > cx) ? ax : cx);
  765.     x=w=v=bx;
  766.     fw=fv=fx=Nf1dim(pN,x);
  767.     for (iter=1;iter<=ITMAX;iter++) {
  768.         xm=0.5f*(a+b);
  769.         tol2=(float) (2.0f*(tol1=(float) (tol*fabs(x)+ZEPS)));
  770.         if (fabs(x-xm) <= (tol2-0.5f*(b-a))) {
  771.             *xmin=x;
  772.             return fx;
  773.         }
  774.         if (fabs(e) > tol1) {
  775.             r=(x-w)*(fx-fv);
  776.             q=(x-v)*(fx-fw);
  777.             p=(x-v)*q-(x-w)*r;
  778.             q=2.0f*(q-r);
  779.             if (q > 0.0) p = -p;
  780.             q=(float)fabs(q);
  781.             etemp=e;
  782.             e=d;
  783.             if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
  784.                 d=CGOLD*(e=(x >= xm ? a-x : b-x));
  785.             else {
  786.                 d=p/q;
  787.                 u=x+d;
  788.                 if (u-a < tol2 || b-u < tol2)
  789.                     d=(float)SIGN(tol1,xm-x);
  790.             }
  791.         } else {
  792.             d=CGOLD*(e=(x >= xm ? a-x : b-x));
  793.         }
  794.         u=(float)(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
  795.         fu=Nf1dim(pN,u);
  796.         if (fu <= fx) {
  797.             if (u >= x) a=x; else b=x;
  798.             SHFT(v,w,x,u)
  799.             SHFT(fv,fw,fx,fu)
  800.         } else {
  801.             if (u < x) a=u; else b=u;
  802.             if (fu <= fw || w == x) {
  803.                 v=w;
  804.                 w=u;
  805.                 fv=fw;
  806.                 fw=fu;
  807.             } else if (fu <= fv || v == x || v == w) {
  808.                 v=u;
  809.                 fv=fu;
  810.             }
  811.         }
  812.     }
  813.     Nrerror("Too many iterations in BRENT");
  814.     *xmin=x;
  815.     return fx;
  816. }
  817.  
  818. #undef ITMAX
  819. #undef CGOLD
  820. #undef ZEPS
  821. #undef SIGN
  822.  
  823. #define GOLD 1.618034f
  824. #define GLIMIT 100.0f
  825. #define TINY 1.0e-20
  826. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  827. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  828. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  829.  
  830. void Nmnbrak(NEURAL *pN, float *ax, float *bx, float *cx, float *fa, float *fb,
  831.         float *fc)
  832. {
  833.     float ulim,u,r,q,fu,dum;
  834.     long loopcount;
  835. //    Logit ("mnbrak\n");
  836.  
  837.     *fa=Nf1dim(pN,*ax);
  838.     *fb=Nf1dim(pN,*bx);
  839.     if (*fb > *fa) {
  840.         SHFT(dum,*ax,*bx,dum)
  841.         SHFT(dum,*fb,*fa,dum)
  842.     }
  843.     *cx=(*bx)+GOLD*(*bx-*ax);
  844.     *fc=Nf1dim(pN,*cx);
  845.  
  846.     loopcount = 0;
  847.     while (*fb > *fc) {
  848.         if (++loopcount > 1000) {
  849.             Logit ("mnbrak loopcount\n");
  850.             exit(1);
  851.         }
  852.         r=(*bx-*ax)*(*fb-*fc);
  853.         q=(*bx-*cx)*(*fb-*fa);
  854.         u=(float) ((*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
  855.             (2.0f*SIGN(MAX(fabs(q-r),TINY),q-r)));
  856.         ulim=(*bx)+GLIMIT*(*cx-*bx);
  857.         if ((*bx-u)*(u-*cx) > 0.0) {
  858.             fu=Nf1dim(pN,u);
  859.             if (fu < *fc) {
  860.                 *ax=(*bx);
  861.                 *bx=u;
  862.                 *fa=(*fb);
  863.                 *fb=fu;
  864.                 return;
  865.             } else if (fu > *fb) {
  866.                 *cx=u;
  867.                 *fc=fu;
  868.                 return;
  869.             }
  870.             u=(*cx)+GOLD*(*cx-*bx);
  871.             fu=Nf1dim(pN,u);
  872.         } else if ((*cx-u)*(u-ulim) > 0.0) {
  873.             fu=Nf1dim(pN,u);
  874.             if (fu < *fc) {
  875.                 SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
  876.                 SHFT(*fb,*fc,fu,Nf1dim(pN,u))
  877.             }
  878.         } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
  879.             u=ulim;
  880.             fu=Nf1dim(pN,u);
  881.         } else {
  882.             u=(*cx)+GOLD*(*cx-*bx);
  883.             fu=Nf1dim(pN,u);
  884.         }
  885.         SHFT(*ax,*bx,*cx,u)
  886.         SHFT(*fa,*fb,*fc,fu)
  887.     }
  888. }
  889.  
  890. #undef GOLD
  891. #undef GLIMIT
  892. #undef TINY
  893. #undef MAX
  894. #undef SIGN
  895. #undef SHFT
  896.  
  897. float Nf1dim(NEURAL *pN, float x)
  898. {
  899.     int j;
  900.     float f;
  901.  
  902. //Logit ("in f1dim %d\n",pN->m_NDIM);
  903.     for (j=0;j<pN->m_NDIM;j++) pN->m_xt[j]=pN->m_pcom[j]+x*pN->m_xicom[j];
  904. //Logit ("20\n");
  905.     f=NErrorFunction(pN,pN->m_xt);
  906. //Logit ("out f1dim\n");
  907.     return f;
  908. }
  909.  
  910. #define TOL 2.0e-4
  911.  
  912. void Nlinmin(NEURAL *pN, float *p, float *xi, float *fret)
  913. {
  914.     int j,n;
  915.     float xx,xmin,fx,fb,fa,bx,ax;
  916. //    Logit ("linmin\n");
  917.            
  918.     n = pN->m_NDIM;
  919.  
  920.     for (j=0;j<n;j++) {
  921.         pN->m_pcom[j]=p[j];
  922.         pN->m_xicom[j]=xi[j];
  923.     }
  924.     ax=0;
  925.     xx=1;
  926.     Nmnbrak(pN,&ax,&xx,&bx,&fa,&fx,&fb);
  927. //Logit ("returned from mnbrak\n");
  928.     *fret= Nbrent(pN,ax,xx,bx,(float) TOL,&xmin);
  929. //Logit ("returned from brent\n");
  930.     for (j=0;j<n;j++) {
  931.         xi[j] *= xmin;
  932.         p[j] += xi[j];
  933.     }
  934. }
  935.  
  936. #undef TOL
  937.  
  938.  
  939. void NforwardDIffGradient (NEURAL *pN, float *x0, float *g)
  940. {
  941.     int j;
  942.     double sqrteta, stepsizej, tempj, fj, eta, fc;
  943. //    eta = 2.22E-16;
  944.     eta = 2.22E-9;
  945.     //Logit ("forwardDIffGradient\n");
  946.  
  947.     fc = NErrorFunction(pN,x0);
  948.     sqrteta = sqrt (eta);
  949.     for (j = 0; j < pN->m_NDIM; ++j) {
  950.         stepsizej = sqrteta * Max ((float)fabs (x0[j]), sign_of_arg (x0[j]));
  951.         tempj = x0[j];
  952.         x0[j] += (float) stepsizej;
  953.  
  954.         stepsizej = x0[j] - tempj;   // reduce finite precision errors
  955.         fj = NErrorFunction (pN,x0);
  956.         g[j] = (float) ((fj - fc) / stepsizej);
  957.  
  958.         x0[j] = (float) tempj;  // reset x0[j]
  959.     }
  960.     //Logit ("end \n");
  961.  
  962. } // end forwardDIffGradient 
  963.  
  964. void Nrerror(char *error_text)
  965. {
  966.     Logit ("***ERROR %s\n",error_text);
  967.     //fprintf (pN->m_fd,"***ERROR %s\n",error_text);
  968.     //fclose(pN->m_fd);
  969.     exit(1);
  970. }
  971.  
  972. // ErrorFunction is called by the CG routine to evaluate the total system error
  973. // it takes the weight vector "x" from the CG routine and puts it into the neural
  974. // network weight structures then runs the entire training matrix through the 
  975. // Neural network to calculate the total system error
  976.  
  977. float NErrorFunction(NEURAL *pN, float* x)
  978. {     
  979.     int i,n,traincnt;
  980.     double errorsq;
  981.  
  982.  
  983.     for (i=0; i < pN->m_nhidden; i++) {
  984.         for (n=0; n < pN->m_ninputs; n++) pN->m_hinputw[i][n] = *x++;
  985.     }
  986.     for (i=0; i < pN->m_noutputs; i++) {
  987.         for (n=0; n < pN->m_nhidden; n++) pN->m_oinputw[i][n] = *x++;
  988. #ifdef NN_STHRU
  989.         if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) 
  990.             for (n=0; n < pN->m_ninputs; n++) pN->m_iinputw[i][n] = *x++;
  991. #endif
  992.     }
  993.  
  994.     
  995.     errorsq = 0.0;
  996.     for(traincnt=0;traincnt<pN->m_dm->m_numrows;traincnt++) {
  997.         for (i=0;i<pN->m_ninputs;i++) pN->m_ni[i] = pN->m_dm->m_iarray[i][traincnt];
  998.         NFeedForward(pN); // Calculate outputs
  999.         for(i=0;i<pN->m_noutputs;i++) 
  1000.             errorsq += pow ((pN->m_dm->m_oarray[i][traincnt] - pN->m_nout[i]),2);
  1001.     }
  1002. //Logit ("%.6f\n",errorsq);
  1003.     return (float) errorsq;
  1004. }
  1005.  
  1006. void NCGTrain(NEURAL *pN )
  1007. {       
  1008.     int i,j,n;
  1009.     float fret;
  1010.     int iter;
  1011.  
  1012.     if (pN->m_cnt < 201) return;
  1013.  
  1014.     pN->m_NDIM = pN->m_nhidden*pN->m_ninputs;
  1015.     pN->m_NDIM += pN->m_noutputs*pN->m_nhidden;
  1016.  
  1017. #ifdef NN_STHRU
  1018.     if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) 
  1019.         pN->m_NDIM += pN->m_noutputs*pN->m_ninputs;
  1020. #endif
  1021.  
  1022.     pN->m_startp = malloc (sizeof(float)*pN->m_NDIM);
  1023.     pN->m_pcom   = malloc (sizeof(float)*pN->m_NDIM);
  1024.     pN->m_xicom  = malloc (sizeof(float)*pN->m_NDIM);
  1025.     pN->m_xt     = malloc (sizeof(float)*pN->m_NDIM);
  1026.     pN->m_g      = malloc (sizeof(float)*pN->m_NDIM);
  1027.     pN->m_h      = malloc (sizeof(float)*pN->m_NDIM);
  1028.     pN->m_xi     = malloc (sizeof(float)*pN->m_NDIM);
  1029.  
  1030.     j=0;
  1031.     for (i=0; i < pN->m_nhidden; i++) {
  1032.         for (n=0; n < pN->m_ninputs; n++) pN->m_startp[j++] = pN->m_hinputw[i][n];
  1033.     }
  1034.     for (i=0; i < pN->m_noutputs; i++) {
  1035.         for (n=0; n < pN->m_nhidden; n++) pN->m_startp[j++] = pN->m_oinputw[i][n];
  1036. #ifdef NN_STHRU
  1037.         if (pN->m_params->m_TrainFlags & TF_NOINTOOUT) 
  1038.             for (n=0; n < pN->m_ninputs; n++) pN->m_startp[j++] = pN->m_iinputw[i][n];
  1039. #endif
  1040.     }
  1041.  
  1042. #define FTOL 1.0e-6
  1043.  
  1044. //    frprmn(pN->m_startp,pN->m_params->m_tol,&iter,&fret);
  1045.     if (pN->m_cnt < 500) pN->m_itmax = 2;
  1046.     else pN->m_itmax = (int) (pN->m_cnt * 20 / pN->m_params->m_cnt_max);
  1047.     Nfrprmn(pN, pN->m_startp,(float)FTOL,&iter,&fret);
  1048.  
  1049.     Logit ("frprmn iter=%d fret=%f\n",iter,fret);
  1050.  
  1051.     free (pN->m_startp);
  1052.     free (pN->m_pcom);
  1053.     free (pN->m_xicom);
  1054.     free (pN->m_xt);
  1055.     free (pN->m_g);
  1056.     free (pN->m_h);
  1057.     free (pN->m_xi);
  1058.  
  1059. }
  1060.  
  1061. int ntransl(char *cdummy)
  1062. {
  1063.     int val=0;
  1064.     //Logit ("ntransl <%s>\n",cdummy);
  1065.  
  1066.     if (cdummy[0] != 'N') return -1;
  1067.     sscanf(&cdummy[1],"%d",&val);
  1068. //    if (val > 12) return -1;
  1069.     return val;
  1070. }
  1071.  
  1072. int NImportNetwork (NEURAL *pN, FILE *fd)
  1073. {
  1074.     int i,j,sel,stat,h,o;
  1075.     char cdummy[80];
  1076.  
  1077.     stat=ImportParams(fd,pN->m_params);
  1078.     if (stat) return stat;
  1079.     pN->m_sumerr2 = 0;
  1080.  
  1081. top:
  1082.     stat = fscanf (fd,"%s",&cdummy);
  1083.     if (stat==EOF) goto errorexit;
  1084.     
  1085.  
  1086.     sel = ntransl(cdummy);
  1087.     switch (sel) {
  1088.     case 99:
  1089.         pN->m_dm = DCreateDataMat();
  1090.         stat=DImportDataMat(pN->m_dm,fd);
  1091. #ifdef VERBOSE
  1092.         Logit("Finished import NN\n");
  1093. #endif
  1094.         return stat;
  1095.  
  1096.         break;
  1097.     default:
  1098.         goto errorexit;
  1099.         break;
  1100.     case 1:
  1101.         fscanf (fd,"%u %d %d %d %ld ", &pN->m_istate, &pN->m_ninputs, &pN->m_nhidden, &pN->m_noutputs, &pN->m_cnt);
  1102.         pN->m_istate |= NN_PARAMSLOADED;
  1103.         pN->m_ni = (float*) malloc (sizeof(float)*pN->m_ninputs);
  1104.         pN->m_houtputv = (float*) malloc (sizeof(float)*pN->m_nhidden);
  1105.         pN->m_htheta = (float*) malloc (sizeof(float)*pN->m_nhidden);
  1106.         pN->m_otheta = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1107.         pN->m_nout = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1108.         pN->m_hinputw = alloc_2d_floats (pN->m_nhidden,pN->m_ninputs);
  1109.         pN->m_oinputw = alloc_2d_floats (pN->m_noutputs,pN->m_nhidden);
  1110. #ifdef NN_STHRU
  1111.         pN->m_iinputw = alloc_2d_floats (pN->m_noutputs,pN->m_ninputs);
  1112. #endif
  1113.  
  1114.         if (pN->m_istate&NN_DYNAMIC) {
  1115.             pN->m_hlastvar      = (float*) malloc (sizeof(float)*pN->m_nhidden);
  1116.             pN->m_hlearn        = (float*) malloc (sizeof(float)*pN->m_nhidden);
  1117.             pN->m_htlearn       = (float*) malloc (sizeof(float)*pN->m_nhidden);
  1118.             pN->m_olastvar      = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1119.             pN->m_otraining     = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1120.             pN->m_olearn        = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1121.             pN->m_otlearn       = (float*) malloc (sizeof(float)*pN->m_noutputs);
  1122.             pN->m_hlastdelta    = alloc_2d_floats (pN->m_nhidden,pN->m_ninputs);
  1123.             pN->m_olastdelta    = alloc_2d_floats (pN->m_noutputs,pN->m_nhidden);
  1124. #ifdef NN_STHRU
  1125.             pN->m_ilastdelta    = alloc_2d_floats (pN->m_noutputs,pN->m_ninputs);
  1126. #endif
  1127.  
  1128.         }
  1129.  
  1130.         for (i=0;i<pN->m_ninputs;i++) pN->m_ni[i] = 0;
  1131.         for (h=0;h<pN->m_nhidden;h++) {
  1132.             for (i=0;i<pN->m_ninputs;i++) pN->m_hinputw[h][i] = 0;
  1133.             pN->m_houtputv[h] = pN->m_htheta[h] = 0;
  1134.         }
  1135.         for (o=0;o<pN->m_noutputs;o++) {
  1136.             for (h=0;h<pN->m_nhidden;h++) pN->m_oinputw[o][h] = 0;
  1137. #ifdef NN_STHRU
  1138.             for (i=0;i<pN->m_ninputs;i++) pN->m_iinputw[o][i] = 0;
  1139. #endif
  1140.             pN->m_otheta[o] = pN->m_nout[o] = 0;
  1141.         }
  1142.         break;
  1143.     case 2:
  1144.         for (i=0;i<pN->m_ninputs;i++)
  1145.             for (j=0;j<pN->m_nhidden;j++) fscanf (fd,"%f ",&pN->m_hinputw[j][i]);
  1146.         break;
  1147.     case 3:
  1148.         for (i=0;i<pN->m_nhidden;i++) fscanf (fd,"%f ",&pN->m_htheta[i]);
  1149.         break;
  1150.     case 4:
  1151.         for (i=0;i<pN->m_noutputs;i++)
  1152.         for (j=0;j<pN->m_nhidden;j++) fscanf (fd,"%f ",&pN->m_oinputw[i][j]);
  1153.         break;
  1154.     case 5:
  1155.         for (i=0;i<pN->m_noutputs;i++) fscanf (fd,"%f ",&pN->m_otheta[i]);
  1156.         break;
  1157.     // Now load data for the dynamic part of CNeural class
  1158.     case 6:
  1159.         for (i=0;i<pN->m_nhidden;i++) fscanf (fd,"%f",&pN->m_hlastvar[i]);
  1160.         break;
  1161.     case 7:
  1162.         for (i=0;i<pN->m_nhidden;i++) fscanf (fd,"%f",&pN->m_hlearn[i]);
  1163.         break;
  1164.     case 8:
  1165.         for (i=0;i<pN->m_nhidden;i++) fscanf (fd,"%f",&pN->m_htlearn[i]);
  1166.         break;
  1167.     case 9:
  1168.         for (i=0;i<pN->m_noutputs;i++) fscanf (fd,"%f",&pN->m_olastvar[i]);
  1169.         break;
  1170.     case 10:
  1171.         for (i=0;i<pN->m_noutputs;i++) fscanf (fd,"%f",&pN->m_olearn[i]);
  1172.         break;
  1173.     case 11:
  1174.         for (i=0;i<pN->m_noutputs;i++) fscanf (fd,"%f",&pN->m_otlearn[i]);
  1175.         break;
  1176. #ifdef NN_STHRU
  1177.     case 12:
  1178.         for (i=0;i<pN->m_noutputs;i++)
  1179.             for (j=0;j<pN->m_ninputs;j++) fscanf (fd,"%f",&pN->m_iinputw[i][j]);
  1180.         break;
  1181. #endif
  1182.     }
  1183.     goto top;
  1184.  
  1185. errorexit:          
  1186.     return -1;
  1187. }
  1188.  
  1189. NEURAL *LoadNetwork (char *filename) {
  1190.     FILE *fd;
  1191.     int stat;
  1192.     PARAMS *tparams;
  1193.     NEURAL *tneural;
  1194.  
  1195.     fd = fopen(filename,"r");
  1196.     if (fd == NULL) return NULL;
  1197.  
  1198.     tparams = (PARAMS*) malloc (sizeof(PARAMS));
  1199.     tneural = (NEURAL*) malloc (sizeof(NEURAL));
  1200.     tneural->m_params = tparams;
  1201.     tneural->m_istate |= NN_PARAMSLOADED;
  1202.  
  1203.     stat = NImportNetwork(tneural,fd);
  1204.     fclose(fd);
  1205.  
  1206.  
  1207.     if (stat) {
  1208.         NDeleteNeural (tneural);
  1209. #ifdef VERBOSE
  1210.         Logit ("Error importing network\n");
  1211. #endif
  1212.         return NULL;
  1213.     }
  1214.     return tneural;
  1215. }
  1216.  
  1217. void DumpNeural(NEURAL *pN, FILE *fd)
  1218. {
  1219.     int i,j;
  1220.     static char fmt1[] = "Dump of Neural Network inp=%d hid=%d outp=%d\n";
  1221.     static char fmt2[] = "ni = ";
  1222.     static char fmt3[] = "\nhinputw = ";
  1223.     static char fmt4[] = "\nhoutputv = ";
  1224.     static char fmt5[] = "\nhtheta = ";
  1225.     static char fmt6[] = "\noinputw = ";
  1226.     static char fmt7[] = "otheta = ";
  1227.     static char fmt8[] = "\nnout = ";
  1228.  
  1229.  
  1230.  
  1231.     static char fmt21[] = "\nhlastvar = ";
  1232.     static char fmt22[] = "\nhlearn = ";
  1233.     static char fmt23[] = "\nhtlearn = ";
  1234.     static char fmt24[] = "\nhlastdelta = ";
  1235.     static char fmt25[] = "\nolastvar = ";
  1236.     static char fmt26[] = "\notraining = ";
  1237.     static char fmt27[] = "\nolearn = ";
  1238.     static char fmt28[] = "\notlearn = ";
  1239.     static char fmt29[] = "\nolastdelta = ";
  1240.     static char fmt30[] = "\nilastdelta = ";
  1241.     static char fmt31[] = "\niinputw = ";
  1242.     
  1243.  
  1244.     fprintf (fd,fmt1,pN->m_ninputs,pN->m_nhidden,pN->m_noutputs);
  1245.  
  1246.     if (!pN->m_ninputs) return;
  1247.     if (!pN->m_nhidden) return;
  1248.     if (!pN->m_noutputs) return;
  1249.     fprintf (fd,fmt2);
  1250.     for (i=0;i<pN->m_ninputs;i++)
  1251.         fprintf (fd,FMT_E,pN->m_ni[i]);
  1252.     fprintf (fd,fmt3);
  1253.     for (i=0;i<pN->m_ninputs;i++) {
  1254.         for (j=0;j<pN->m_nhidden;j++) fprintf (fd,FMT_E,pN->m_hinputw[j][i]);
  1255.     }
  1256.  
  1257.     fprintf (fd,fmt4);
  1258.     for (i=0;i<pN->m_nhidden;i++) fprintf (fd,FMT_E,pN->m_houtputv[i]);
  1259.     fprintf (fd,fmt5);
  1260.     for (i=0;i<pN->m_nhidden;i++) fprintf (fd,FMT_E,pN->m_htheta[i]);
  1261.  
  1262.     fprintf (fd,fmt6);
  1263.     for (i=0;i<pN->m_noutputs;i++) {
  1264.         for (j=0;j<pN->m_nhidden;j++) fprintf (fd,FMT_E,pN->m_oinputw[i][j]);
  1265.         fprintf (fd,"\n");
  1266.     }
  1267.     fprintf (fd,fmt7);
  1268.     for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_otheta[i]);
  1269.  
  1270. #ifdef NN_STHRU
  1271.     fprintf (fd,fmt31);
  1272.     for (i=0;i<pN->m_noutputs;i++) 
  1273.         for (j=0;j<pN->m_ninputs;j++) fprintf (fd,FMT_E,pN->m_iinputw[i][j]);
  1274. #endif
  1275.  
  1276.     fprintf (fd,fmt8);
  1277.     for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_nout[i]);
  1278.     fprintf (fd,"\n");
  1279.  
  1280.  
  1281.     if (pN->m_istate&NN_PARAMSLOADED) DumpParams(pN->m_params,fd);
  1282.  
  1283.     if (pN->m_istate&NN_DYNAMIC) {
  1284.  
  1285. // dneural part
  1286.  
  1287.         fprintf (fd,fmt21);
  1288.         for (i=0;i<pN->m_nhidden;i++) fprintf (fd,FMT_E,pN->m_hlastvar[i]);
  1289.         fprintf (fd,fmt22);
  1290.         for (i=0;i<pN->m_nhidden;i++) fprintf (fd,FMT_E,pN->m_hlearn[i]);
  1291.         fprintf (fd,fmt23);
  1292.         for (i=0;i<pN->m_nhidden;i++) fprintf (fd,FMT_E,pN->m_htlearn[i]);
  1293.         fprintf (fd,fmt24);
  1294.         for (i=0;i<pN->m_ninputs;i++) {
  1295.             for (j=0;j<pN->m_nhidden;j++) fprintf (fd,FMT_E,pN->m_hlastdelta[j][i]);
  1296.         }
  1297.     
  1298.         fprintf (fd,fmt25);
  1299.         for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_olastvar[i]);
  1300.         fprintf (fd,fmt26);
  1301.         for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_otraining[i]);
  1302.         fprintf (fd,fmt27);
  1303.         for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_olearn[i]);
  1304.         fprintf (fd,fmt28);
  1305.         for (i=0;i<pN->m_noutputs;i++) fprintf (fd,FMT_E,pN->m_otlearn[i]);
  1306.         fprintf (fd,fmt29);
  1307.         for (i=0;i<pN->m_noutputs;i++) {
  1308.             for (j=0;j<pN->m_nhidden;j++) fprintf (fd,FMT_E,pN->m_olastdelta[i][j]);
  1309.         }
  1310. #ifdef NN_STHRU
  1311.         fprintf (fd,fmt30);
  1312.         for (i=0;i<pN->m_noutputs;i++) {
  1313.             for (j=0;j<pN->m_ninputs;j++) fprintf (fd,FMT_E,pN->m_ilastdelta[i][j]);
  1314.         }
  1315. #endif
  1316.         fprintf (fd,"\n\n");
  1317.     }
  1318.     if (pN->m_dm!=NULL) DumpDataMat(pN->m_dm,fd);
  1319.     fprintf (fd,"\n\n");
  1320. }
  1321.  
  1322. void NInterrogate(NEURAL *pN,float *Ivec,float *Ovec)
  1323. {
  1324.     int i;
  1325.  
  1326.     for (i=0; i < pN->m_ninputs; i++) pN->m_ni[i] = DScale(pN->m_dm,Ivec[i],'I',i);
  1327.     NFeedForward(pN);
  1328.     for (i=0; i < pN->m_noutputs; i++) Ovec[i] = DRescale(pN->m_dm,pN->m_nout[i],'O',i);
  1329. }
  1330.  
  1331.