BPN.C
上传用户:libin825
上传日期:2022-07-27
资源大小:230k
文件大小:17k
开发平台:

C/C++

  1. /******************************************************************************
  2.                       ====================================================
  3.         Network:      Backpropagation Network with Bias Terms and Momentum
  4.                       ====================================================
  5.         Application:  Time-Series Forecasting
  6.                       Prediction of the Annual Number of Sunspots
  7.         Author:       Karsten Kutza
  8.         Date:         17.4.96
  9.         Reference:    D.E. Rumelhart, G.E. Hinton, R.J. Williams
  10.                       Learning Internal Representations by Error Propagation
  11.                       in:
  12.                       D.E. Rumelhart, J.L. McClelland (Eds.)
  13.                       Parallel Distributed Processing, Volume 1
  14.                       MIT Press, Cambridge, MA, pp. 318-362, 1986
  15.  ******************************************************************************/
  16. /******************************************************************************
  17.                             D E C L A R A T I O N S
  18.  ******************************************************************************/
  19. #include <stdlib.h>
  20. #include <stdio.h>
  21. #include <math.h>
  22. typedef int           BOOL;
  23. typedef int           INT;
  24. typedef double        REAL;
  25. #define FALSE         0
  26. #define TRUE          1
  27. #define NOT           !
  28. #define AND           &&
  29. #define OR            ||
  30. #define MIN_REAL      -HUGE_VAL
  31. #define MAX_REAL      +HUGE_VAL
  32. #define MIN(x,y)      ((x)<(y) ? (x) : (y))
  33. #define MAX(x,y)      ((x)>(y) ? (x) : (y))
  34. #define LO            0.1
  35. #define HI            0.9
  36. #define BIAS          1
  37. #define sqr(x)        ((x)*(x))
  38. typedef struct {                     /* A LAYER OF A NET:                     */
  39.         INT           Units;         /* - number of units in this layer       */
  40.         REAL*         Output;        /* - output of ith unit                  */
  41.         REAL*         Error;         /* - error term of ith unit              */
  42.         REAL**        Weight;        /* - connection weights to ith unit      */
  43.         REAL**        WeightSave;    /* - saved weights for stopped training  */
  44.         REAL**        dWeight;       /* - last weight deltas for momentum     */
  45. } LAYER;
  46. typedef struct {                     /* A NET:                                */
  47.         LAYER**       Layer;         /* - layers of this net                  */
  48.         LAYER*        InputLayer;    /* - input layer                         */
  49.         LAYER*        OutputLayer;   /* - output layer                        */
  50.         REAL          Alpha;         /* - momentum factor                     */
  51.         REAL          Eta;           /* - learning rate                       */
  52.         REAL          Gain;          /* - gain of sigmoid function            */
  53.         REAL          Error;         /* - total net error                     */
  54. } NET;
  55. /******************************************************************************
  56.         R A N D O M S   D R A W N   F R O M   D I S T R I B U T I O N S
  57.  ******************************************************************************/
  58. void InitializeRandoms()
  59. {
  60.   srand(4711);
  61. }
  62. INT RandomEqualINT(INT Low, INT High)
  63. {
  64.   return rand() % (High-Low+1) + Low;
  65. }      
  66. REAL RandomEqualREAL(REAL Low, REAL High)
  67. {
  68.   return ((REAL) rand() / RAND_MAX) * (High-Low) + Low;
  69. }      
  70. /******************************************************************************
  71.                A P P L I C A T I O N - S P E C I F I C   C O D E
  72.  ******************************************************************************/
  73. #define NUM_LAYERS    3
  74. #define N             30
  75. #define M             1
  76. INT                   Units[NUM_LAYERS] = {N, 10, M};
  77. #define FIRST_YEAR    1700
  78. #define NUM_YEARS     280
  79. #define TRAIN_LWB     (N)
  80. #define TRAIN_UPB     (179)
  81. #define TRAIN_YEARS   (TRAIN_UPB - TRAIN_LWB + 1)
  82. #define TEST_LWB      (180)
  83. #define TEST_UPB      (259)
  84. #define TEST_YEARS    (TEST_UPB - TEST_LWB + 1)
  85. #define EVAL_LWB      (260)
  86. #define EVAL_UPB      (NUM_YEARS - 1)
  87. #define EVAL_YEARS    (EVAL_UPB - EVAL_LWB + 1)
  88. REAL                  Sunspots_[NUM_YEARS];
  89. REAL                  Sunspots [NUM_YEARS] = {
  90.                         0.0262,  0.0575,  0.0837,  0.1203,  0.1883,  0.3033,  
  91.                         0.1517,  0.1046,  0.0523,  0.0418,  0.0157,  0.0000,  
  92.                         0.0000,  0.0105,  0.0575,  0.1412,  0.2458,  0.3295,  
  93.                         0.3138,  0.2040,  0.1464,  0.1360,  0.1151,  0.0575,  
  94.                         0.1098,  0.2092,  0.4079,  0.6381,  0.5387,  0.3818,  
  95.                         0.2458,  0.1831,  0.0575,  0.0262,  0.0837,  0.1778,  
  96.                         0.3661,  0.4236,  0.5805,  0.5282,  0.3818,  0.2092,  
  97.                         0.1046,  0.0837,  0.0262,  0.0575,  0.1151,  0.2092,  
  98.                         0.3138,  0.4231,  0.4362,  0.2495,  0.2500,  0.1606,  
  99.                         0.0638,  0.0502,  0.0534,  0.1700,  0.2489,  0.2824,  
  100.                         0.3290,  0.4493,  0.3201,  0.2359,  0.1904,  0.1093,  
  101.                         0.0596,  0.1977,  0.3651,  0.5549,  0.5272,  0.4268,  
  102.                         0.3478,  0.1820,  0.1600,  0.0366,  0.1036,  0.4838,  
  103.                         0.8075,  0.6585,  0.4435,  0.3562,  0.2014,  0.1192,  
  104.                         0.0534,  0.1260,  0.4336,  0.6904,  0.6846,  0.6177,  
  105.                         0.4702,  0.3483,  0.3138,  0.2453,  0.2144,  0.1114,  
  106.                         0.0837,  0.0335,  0.0214,  0.0356,  0.0758,  0.1778,  
  107.                         0.2354,  0.2254,  0.2484,  0.2207,  0.1470,  0.0528,  
  108.                         0.0424,  0.0131,  0.0000,  0.0073,  0.0262,  0.0638,  
  109.                         0.0727,  0.1851,  0.2395,  0.2150,  0.1574,  0.1250,  
  110.                         0.0816,  0.0345,  0.0209,  0.0094,  0.0445,  0.0868,  
  111.                         0.1898,  0.2594,  0.3358,  0.3504,  0.3708,  0.2500,  
  112.                         0.1438,  0.0445,  0.0690,  0.2976,  0.6354,  0.7233,  
  113.                         0.5397,  0.4482,  0.3379,  0.1919,  0.1266,  0.0560,  
  114.                         0.0785,  0.2097,  0.3216,  0.5152,  0.6522,  0.5036,  
  115.                         0.3483,  0.3373,  0.2829,  0.2040,  0.1077,  0.0350,  
  116.                         0.0225,  0.1187,  0.2866,  0.4906,  0.5010,  0.4038,  
  117.                         0.3091,  0.2301,  0.2458,  0.1595,  0.0853,  0.0382,  
  118.                         0.1966,  0.3870,  0.7270,  0.5816,  0.5314,  0.3462,  
  119.                         0.2338,  0.0889,  0.0591,  0.0649,  0.0178,  0.0314,  
  120.                         0.1689,  0.2840,  0.3122,  0.3332,  0.3321,  0.2730,  
  121.                         0.1328,  0.0685,  0.0356,  0.0330,  0.0371,  0.1862,  
  122.                         0.3818,  0.4451,  0.4079,  0.3347,  0.2186,  0.1370,  
  123.                         0.1396,  0.0633,  0.0497,  0.0141,  0.0262,  0.1276,  
  124.                         0.2197,  0.3321,  0.2814,  0.3243,  0.2537,  0.2296,  
  125.                         0.0973,  0.0298,  0.0188,  0.0073,  0.0502,  0.2479,  
  126.                         0.2986,  0.5434,  0.4215,  0.3326,  0.1966,  0.1365,  
  127.                         0.0743,  0.0303,  0.0873,  0.2317,  0.3342,  0.3609,  
  128.                         0.4069,  0.3394,  0.1867,  0.1109,  0.0581,  0.0298,  
  129.                         0.0455,  0.1888,  0.4168,  0.5983,  0.5732,  0.4644,  
  130.                         0.3546,  0.2484,  0.1600,  0.0853,  0.0502,  0.1736,  
  131.                         0.4843,  0.7929,  0.7128,  0.7045,  0.4388,  0.3630,  
  132.                         0.1647,  0.0727,  0.0230,  0.1987,  0.7411,  0.9947,  
  133.                         0.9665,  0.8316,  0.5873,  0.2819,  0.1961,  0.1459,  
  134.                         0.0534,  0.0790,  0.2458,  0.4906,  0.5539,  0.5518,  
  135.                         0.5465,  0.3483,  0.3603,  0.1987,  0.1804,  0.0811,  
  136.                         0.0659,  0.1428,  0.4838,  0.8127 
  137.                       };
  138. REAL                  Mean;
  139. REAL                  TrainError;
  140. REAL                  TrainErrorPredictingMean;
  141. REAL                  TestError;
  142. REAL                  TestErrorPredictingMean;
  143. FILE*                 f;
  144. void NormalizeSunspots()
  145. {
  146.   INT  Year;
  147.   REAL Min, Max;
  148.   Min = MAX_REAL;
  149.   Max = MIN_REAL;
  150.   for (Year=0; Year<NUM_YEARS; Year++) {
  151.     Min = MIN(Min, Sunspots[Year]);
  152.     Max = MAX(Max, Sunspots[Year]);
  153.   }
  154.   Mean = 0;
  155.   for (Year=0; Year<NUM_YEARS; Year++) {
  156.     Sunspots_[Year] = 
  157.     Sunspots [Year] = ((Sunspots[Year]-Min) / (Max-Min)) * (HI-LO) + LO;
  158.     Mean += Sunspots[Year] / NUM_YEARS;
  159.   }
  160. }
  161. void InitializeApplication(NET* Net)
  162. {
  163.   INT  Year, i;
  164.   REAL Out, Err;
  165.   Net->Alpha = 0.5;
  166.   Net->Eta   = 0.05;
  167.   Net->Gain  = 1;
  168.   NormalizeSunspots();
  169.   TrainErrorPredictingMean = 0;
  170.   for (Year=TRAIN_LWB; Year<=TRAIN_UPB; Year++) {
  171.     for (i=0; i<M; i++) {
  172.       Out = Sunspots[Year+i];
  173.       Err = Mean - Out;
  174.       TrainErrorPredictingMean += 0.5 * sqr(Err);
  175.     }
  176.   }
  177.   TestErrorPredictingMean = 0;
  178.   for (Year=TEST_LWB; Year<=TEST_UPB; Year++) {
  179.     for (i=0; i<M; i++) {
  180.       Out = Sunspots[Year+i];
  181.       Err = Mean - Out;
  182.       TestErrorPredictingMean += 0.5 * sqr(Err);
  183.     }
  184.   }
  185.   f = fopen("BPN.txt", "w");
  186. }
  187. void FinalizeApplication(NET* Net)
  188. {
  189.   fclose(f);
  190. }
  191. /******************************************************************************
  192.                           I N I T I A L I Z A T I O N
  193.  ******************************************************************************/
  194. void GenerateNetwork(NET* Net)
  195. {
  196.   INT l,i;
  197.   Net->Layer = (LAYER**) calloc(NUM_LAYERS, sizeof(LAYER*));
  198.    
  199.   for (l=0; l<NUM_LAYERS; l++) {
  200.     Net->Layer[l] = (LAYER*) malloc(sizeof(LAYER));
  201.       
  202.     Net->Layer[l]->Units      = Units[l];
  203.     Net->Layer[l]->Output     = (REAL*)  calloc(Units[l]+1, sizeof(REAL));
  204.     Net->Layer[l]->Error      = (REAL*)  calloc(Units[l]+1, sizeof(REAL));
  205.     Net->Layer[l]->Weight     = (REAL**) calloc(Units[l]+1, sizeof(REAL*));
  206.     Net->Layer[l]->WeightSave = (REAL**) calloc(Units[l]+1, sizeof(REAL*));
  207.     Net->Layer[l]->dWeight    = (REAL**) calloc(Units[l]+1, sizeof(REAL*));
  208.     Net->Layer[l]->Output[0]  = BIAS;
  209.       
  210.     if (l != 0) {
  211.       for (i=1; i<=Units[l]; i++) {
  212.         Net->Layer[l]->Weight[i]     = (REAL*) calloc(Units[l-1]+1, sizeof(REAL));
  213.         Net->Layer[l]->WeightSave[i] = (REAL*) calloc(Units[l-1]+1, sizeof(REAL));
  214.         Net->Layer[l]->dWeight[i]    = (REAL*) calloc(Units[l-1]+1, sizeof(REAL));
  215.       }
  216.     }
  217.   }
  218.   Net->InputLayer  = Net->Layer[0];
  219.   Net->OutputLayer = Net->Layer[NUM_LAYERS - 1];
  220.   Net->Alpha       = 0.9;
  221.   Net->Eta         = 0.25;
  222.   Net->Gain        = 1;
  223. }
  224. void RandomWeights(NET* Net)
  225. {
  226.   INT l,i,j;
  227.    
  228.   for (l=1; l<NUM_LAYERS; l++) {
  229.     for (i=1; i<=Net->Layer[l]->Units; i++) {
  230.       for (j=0; j<=Net->Layer[l-1]->Units; j++) {
  231.         Net->Layer[l]->Weight[i][j] = RandomEqualREAL(-0.5, 0.5);
  232.       }
  233.     }
  234.   }
  235. }
  236. void SetInput(NET* Net, REAL* Input)
  237. {
  238.   INT i;
  239.    
  240.   for (i=1; i<=Net->InputLayer->Units; i++) {
  241.     Net->InputLayer->Output[i] = Input[i-1];
  242.   }
  243. }
  244. void GetOutput(NET* Net, REAL* Output)
  245. {
  246.   INT i;
  247.    
  248.   for (i=1; i<=Net->OutputLayer->Units; i++) {
  249.     Output[i-1] = Net->OutputLayer->Output[i];
  250.   }
  251. }
  252. /******************************************************************************
  253.             S U P P O R T   F O R   S T O P P E D   T R A I N I N G
  254.  ******************************************************************************/
  255. void SaveWeights(NET* Net)
  256. {
  257.   INT l,i,j;
  258.   for (l=1; l<NUM_LAYERS; l++) {
  259.     for (i=1; i<=Net->Layer[l]->Units; i++) {
  260.       for (j=0; j<=Net->Layer[l-1]->Units; j++) {
  261.         Net->Layer[l]->WeightSave[i][j] = Net->Layer[l]->Weight[i][j];
  262.       }
  263.     }
  264.   }
  265. }
  266. void RestoreWeights(NET* Net)
  267. {
  268.   INT l,i,j;
  269.   for (l=1; l<NUM_LAYERS; l++) {
  270.     for (i=1; i<=Net->Layer[l]->Units; i++) {
  271.       for (j=0; j<=Net->Layer[l-1]->Units; j++) {
  272.         Net->Layer[l]->Weight[i][j] = Net->Layer[l]->WeightSave[i][j];
  273.       }
  274.     }
  275.   }
  276. }
  277. /******************************************************************************
  278.                      P R O P A G A T I N G   S I G N A L S
  279.  ******************************************************************************/
  280. void PropagateLayer(NET* Net, LAYER* Lower, LAYER* Upper)
  281. {
  282.   INT  i,j;
  283.   REAL Sum;
  284.   for (i=1; i<=Upper->Units; i++) {
  285.     Sum = 0;
  286.     for (j=0; j<=Lower->Units; j++) {
  287.       Sum += Upper->Weight[i][j] * Lower->Output[j];
  288.     }
  289.     Upper->Output[i] = 1 / (1 + exp(-Net->Gain * Sum));
  290.   }
  291. }
  292. void PropagateNet(NET* Net)
  293. {
  294.   INT l;
  295.    
  296.   for (l=0; l<NUM_LAYERS-1; l++) {
  297.     PropagateLayer(Net, Net->Layer[l], Net->Layer[l+1]);
  298.   }
  299. }
  300. /******************************************************************************
  301.                   B A C K P R O P A G A T I N G   E R R O R S
  302.  ******************************************************************************/
  303. void ComputeOutputError(NET* Net, REAL* Target)
  304. {
  305.   INT  i;
  306.   REAL Out, Err;
  307.    
  308.   Net->Error = 0;
  309.   for (i=1; i<=Net->OutputLayer->Units; i++) {
  310.     Out = Net->OutputLayer->Output[i];
  311.     Err = Target[i-1]-Out;
  312.     Net->OutputLayer->Error[i] = Net->Gain * Out * (1-Out) * Err;
  313.     Net->Error += 0.5 * sqr(Err);
  314.   }
  315. }
  316. void BackpropagateLayer(NET* Net, LAYER* Upper, LAYER* Lower)
  317. {
  318.   INT  i,j;
  319.   REAL Out, Err;
  320.    
  321.   for (i=1; i<=Lower->Units; i++) {
  322.     Out = Lower->Output[i];
  323.     Err = 0;
  324.     for (j=1; j<=Upper->Units; j++) {
  325.       Err += Upper->Weight[j][i] * Upper->Error[j];
  326.     }
  327.     Lower->Error[i] = Net->Gain * Out * (1-Out) * Err;
  328.   }
  329. }
  330. void BackpropagateNet(NET* Net)
  331. {
  332.   INT l;
  333.    
  334.   for (l=NUM_LAYERS-1; l>1; l--) {
  335.     BackpropagateLayer(Net, Net->Layer[l], Net->Layer[l-1]);
  336.   }
  337. }
  338. void AdjustWeights(NET* Net)
  339. {
  340.   INT  l,i,j;
  341.   REAL Out, Err, dWeight;
  342.    
  343.   for (l=1; l<NUM_LAYERS; l++) {
  344.     for (i=1; i<=Net->Layer[l]->Units; i++) {
  345.       for (j=0; j<=Net->Layer[l-1]->Units; j++) {
  346.         Out = Net->Layer[l-1]->Output[j];
  347.         Err = Net->Layer[l]->Error[i];
  348.         dWeight = Net->Layer[l]->dWeight[i][j];
  349.         Net->Layer[l]->Weight[i][j] += Net->Eta * Err * Out + Net->Alpha * dWeight;
  350.         Net->Layer[l]->dWeight[i][j] = Net->Eta * Err * Out;
  351.       }
  352.     }
  353.   }
  354. }
  355. /******************************************************************************
  356.                       S I M U L A T I N G   T H E   N E T
  357.  ******************************************************************************/
  358. void SimulateNet(NET* Net, REAL* Input, REAL* Output, REAL* Target, BOOL Training)
  359. {
  360.   SetInput(Net, Input);
  361.   PropagateNet(Net);
  362.   GetOutput(Net, Output);
  363.    
  364.   ComputeOutputError(Net, Target);
  365.   if (Training) {
  366.     BackpropagateNet(Net);
  367.     AdjustWeights(Net);
  368.   }
  369. }
  370. void TrainNet(NET* Net, INT Epochs)
  371. {
  372.   INT  Year, n;
  373.   REAL Output[M];
  374.   for (n=0; n<Epochs*TRAIN_YEARS; n++) {
  375.     Year = RandomEqualINT(TRAIN_LWB, TRAIN_UPB);
  376.     SimulateNet(Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), TRUE);
  377.   }
  378. }
  379. void TestNet(NET* Net)
  380. {
  381.   INT  Year;
  382.   REAL Output[M];
  383.   TrainError = 0;
  384.   for (Year=TRAIN_LWB; Year<=TRAIN_UPB; Year++) {
  385.     SimulateNet(Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), FALSE);
  386.     TrainError += Net->Error;
  387.   }
  388.   TestError = 0;
  389.   for (Year=TEST_LWB; Year<=TEST_UPB; Year++) {
  390.     SimulateNet(Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), FALSE);
  391.     TestError += Net->Error;
  392.   }
  393.   fprintf(f, "nNMSE is %0.3f on Training Set and %0.3f on Test Set",
  394.              TrainError / TrainErrorPredictingMean,
  395.              TestError / TestErrorPredictingMean);
  396. }
  397. void EvaluateNet(NET* Net)
  398. {
  399.   INT  Year;
  400.   REAL Output [M];
  401.   REAL Output_[M];
  402.   fprintf(f, "nnn");
  403.   fprintf(f, "Year    Sunspots    Open-Loop Prediction    Closed-Loop Predictionn");
  404.   fprintf(f, "n");
  405.   for (Year=EVAL_LWB; Year<=EVAL_UPB; Year++) {
  406.     SimulateNet(Net, &(Sunspots [Year-N]), Output,  &(Sunspots [Year]), FALSE);
  407.     SimulateNet(Net, &(Sunspots_[Year-N]), Output_, &(Sunspots_[Year]), FALSE);
  408.     Sunspots_[Year] = Output_[0];
  409.     fprintf(f, "%d       %0.3f                   %0.3f                     %0.3fn",
  410.                FIRST_YEAR + Year,
  411.                Sunspots[Year],
  412.                Output [0],
  413.                Output_[0]);
  414.   }
  415. }
  416. /******************************************************************************
  417.                                     M A I N
  418.  ******************************************************************************/
  419. void main()
  420. {
  421.   NET  Net;
  422.   BOOL Stop;
  423.   REAL MinTestError;
  424.   InitializeRandoms();
  425.   GenerateNetwork(&Net);
  426.   RandomWeights(&Net);
  427.   InitializeApplication(&Net);
  428.   Stop = FALSE;
  429.   MinTestError = MAX_REAL;
  430.   do {
  431.     TrainNet(&Net, 10);
  432.     TestNet(&Net);
  433.     if (TestError < MinTestError) {
  434.       fprintf(f, " - saving Weights ...");
  435.       MinTestError = TestError;
  436.       SaveWeights(&Net);
  437.     }
  438.     else if (TestError > 1.2 * MinTestError) {
  439.       fprintf(f, " - stopping Training and restoring Weights ...");
  440.       Stop = TRUE;
  441.       RestoreWeights(&Net);
  442.     }
  443.   } while (NOT Stop);
  444.   TestNet(&Net);
  445.   EvaluateNet(&Net);
  446.    
  447.   FinalizeApplication(&Net);
  448. }