using a Multi Layer Perceptron
Sébastien Marcel - Lab: Solving the Canon-ball problem using a Multi Layer Perceptron
Quote:
/* This program implements a simple Multi Layer Perceptron to solve the canonball problem (regression)
Adapted by Sebastien Marcel (2003-2004) from D. Collobert
The goal is to estimate the position X of a bullet at time t (fixed) given the initial speed (v) and angle (a)
y
/ \
|
|
|
| . X
| .
| .
| .
| .
| .
| .
| .
| .
| .
|. (v,a)
--------------------------> x
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//
// Datasets
#define N_PATTERNS_TRAIN 500
#define N_PATTERNS_TEST 500
//*****************************
// Inputs
//
struct DataIn
{
float v; // speed of bullet
float a; // angle of bullet
} X[N_PATTERNS_TRAIN + N_PATTERNS_TEST];
//*****************************
// Targets
//
struct DataOut
{
float x; // x position at time t
float y; // y position at time t
} Y[N_PATTERNS_TRAIN + N_PATTERNS_TEST];
//*****************************
// to create the data
//
const float vmax = 1.0/10;
const float gravity = 1.0/200;
//*****************************
// to train the MLP
//
const float lambda = 0.01; // learning rate
const float mu = 0.6; // inertia momentum rate
const float mse_min = 0.001; // minimum Mean Squared Error
const int max_iterations = 500; // maximum number of iterations
//*****************************
// MLP data
//
#define N_HU 3
// weights between hidden neurons and inputs
float w1[N_HU + 1][3];
float w1old[N_HU + 1][3];
// weights between outputs and hidden neurons
float w2[2][N_HU + 1];
float w2old[2][N_HU + 1];
// values of integration function
float aHidden[N_HU + 1];
float aOutput[2];
// values of transfert function
float yHidden[N_HU + 1];
// constant value of bias
float xBias = 1.0;
//*****************************
// Sigmoid transfer function
//
float f(float x)
{
return (1.0 /(1.0 + exp(-x)));
}
//*****************************
// Derivative of Sigmoid
//
float f_prime(float x)
{
float z = exp(-x);
float one_plus_z = 1.0 + z;
return (z / (one_plus_z * one_plus_z));
}
#define MY_PI 3.141592
//*****************************
// Compute the x position of bullet at time t
//
float xpos(float t, float v, float teta)
{
return v * t * cos(teta*MY_PI/180.0);
}
//*****************************
// Compute the y position of bullet at time t
//
float ypos(float t, float v, float teta)
{
return -0.5 * gravity * t * t + v * t * sin(teta*MY_PI/180.0);
}
//*****************************
// Compute the MSE
//
float MSE(struct DataOut a, struct DataOut b)
{
return 0.5*((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
//*****************************
// Random
//
float Random(float inf, float sup)
{
return( ((sup - inf)*((float)(random() & 0x7fffffffL) / (float) 0x7fffffffL) ) + inf);
}
//*****************************
// Create the MLP
//
void createMLP()
{
int i, j;
for (i = 0; i<=2; i++)
{
for (j = 0; j<= N_HU; j++)
{
w1[j][i] = Random(-1.0, 1.0);
w1old[j][i] = w1[j][i];
}
}
for (i = 0; i<2; i++)
{
for (j = 0; j<= N_HU; j++)
{
w2[i][j] = Random(-1.0, 1.0);
w2old[i][j] = w2[i][j];
}
}
}
//*****************************
// Create the datasets
//
void createDatasets()
{
int i;
i = 0;
while (i < (N_PATTERNS_TRAIN + N_PATTERNS_TEST))
{
X[i].v = (Random(0.0, vmax)) / vmax;
X[i].a = (Random(0.0, 90.0)) / (float)90;
Y[i].x = xpos(10, X[i].v * vmax, (X[i].a)*90.0);
Y[i].y = ypos(10, X[i].v * vmax, (X[i].a)*90.0);
if (Y[i].y > 0.0) i = i + 1;
}
}
//*****************************
// Forward the input in the MLP
//
void forward(struct DataIn In, struct DataOut *out)
{
int i;
aHidden[0] = xBias;
yHidden[0] = xBias;
for (i = 1; i <= N_HU; i++)
{
aHidden[i] = xBias*w1[i][0] + (In.v)*w1[i][1] + (In.a)*w1[i][2];
yHidden[i] = f(aHidden[i]);
}
aOutput[0] = 0.0;
aOutput[1] = 0.0;
for (i = 0; i <= N_HU; i++)
{
aOutput[0] = aOutput[0] + yHidden[i]*w2[0][i];
aOutput[1] = aOutput[1] + yHidden[i]*w2[1][i];
}
out->x = f(aOutput[0]);
out->y = f(aOutput[1]);
}
//*****************************
// Backward (back-propagate the gradient of error)
//
void backward(struct DataOut Yestimate, struct DataOut Outwant, struct DataIn In, float lambda, float mu)
{
float wnew, Goutput[2], Ghidden[N_HU + 1];
int i,j;
//
Goutput[0] = (Yestimate.x - Outwant.x) * f_prime(aOutput[0]);
Goutput[1] = (Yestimate.y - Outwant.y) * f_prime(aOutput[1]);
//
for (i = 0; i<=N_HU; i++)
Ghidden[i] = (Goutput[0]*w2[0][i] + Goutput[1]*w2[1][i])*f_prime(aHidden[i]);
//
for (j = 0; j<2; j++)
{
for (i = 0; i<=N_HU; i++)
{
wnew = w2[j][i] - lambda * yHidden[i] * Goutput[j] + mu * (w2[j][i] - w2old[j][i]);
w2old[j][i] = w2[j][i];
w2[j][i] = wnew;
}
}
//
for (i = 0; i<=N_HU; i++)
{
wnew = w1[i][0] - lambda * xBias * Ghidden[i] + mu * (w1[i][0] - w1old[i][0]);
w1old[i][0] = w1[i][0];
w1[i][0] = wnew;
wnew = w1[i][1] - lambda * In.v * Ghidden[i] + mu * (w1[i][1] - w1old[i][1]);
w1old[i][1] = w1[i][1];
w1[i][1] = wnew;
wnew = w1[i][2] - lambda * In.a * Ghidden[i] + mu * (w1[i][2] - w1old[i][2]);
w1old[i][2] = w1[i][2];
w1[i][2] = wnew;
}
}
//*****************************
// Stochastic gradient training
//
float train(int P)
{
struct DataOut Yestimate;
float mse_;
float mse_total;
mse_total = 0.0;
// For each train patterns
for(int p = 0 ; p < P ; p++)
{
// Forward current train pattern into the MLP
forward(X[p], &Yestimate);
// Computes the MSE
mse_ = MSE(Y[p], Yestimate);
// Accumulate the MSE
mse_total = mse_total + mse_;
// Backward MSE gradient
backward(Yestimate, Y[p], X[p], lambda, mu);
}
// Return normalized train MSE
return mse_total / (float) P;
}
//
int main()
{
float mse_train;
float mse_test;
float mse_;
struct DataOut Yestimate;
//*****************************
//
createMLP();
//*****************************
//
createDatasets();
//for (int i = 0; i< N_PATTERNS_TRAIN; i++) printf(" TRN: x=[%f %f] y=[%f %f]\n", X[i].v, X[i].a, Y[i].x, Y[i].y);
//for (int i = N_PATTERNS_TRAIN; i<(N_PATTERNS_TRAIN + N_PATTERNS_TEST); i++) printf(" TST: x=[%f %f] y=[%f %f]\n", X[i].v, X[i].a, Y[i].x, Y[i].y);
printf("Stochastic gradient training:\n");
//
FILE *pf_train = fopen("mse_train.txt", "w");
FILE *pf_test = fopen("mse_test.txt", "w");
int iter;
for(iter = 1 ; iter <= max_iterations ; iter++)
{
//*****************************
//
// train
mse_train = train(N_PATTERNS_TRAIN);
fprintf(pf_train, "%f\n", mse_train);
//*****************************
//
// test
mse_test = 0.0;
// For each test pattern
for(int i = N_PATTERNS_TRAIN ; i < (N_PATTERNS_TRAIN + N_PATTERNS_TEST) ; i++)
{
Yestimate.x = 0.0;
Yestimate.y = 0.0;
// forward current pattern into the MLP
forward(X[i], &Yestimate);
// computes MSE
mse_ = MSE(Y[i], Yestimate);
// accumulate MSE
mse_test += mse_;
}
// Normalize the MSE
mse_test /= (float)N_PATTERNS_TEST;
fprintf(pf_test, "%f\n", mse_test);
printf("."); fflush(stdout);
if (mse_train < mse_min) break;
}
printf("\n");
//
fclose(pf_test);
fclose(pf_train);
//*****************************
// Print info about training
//
printf("Number of iterations = %d\n", iter);
printf("Final MSE train = %f\n", mse_train);
printf("Final MSE test = %f\n", mse_test);
printf("End of program.\n");
return 0;
}
|
Last edited by barnix; 05-13-2008 at 08:01 AM.
|