How to do Principle Component Analysis (PCA)

Finding literature on Principle Component Analysis (PCA) is very easy but very hard to find exactly how to do or implement PCA. Hence, I have demonstrated simple MATLB code below. Following code assumes that variable in columns and numbers of records in row. Variables name are long because I want those to be self descriptive.

covResult = cov(trainingInputs); % Calculating covariance matrix
[V,D] = eig(covResult); % Calculating eigen values and eigen vector
PCATrainingInputs = trainingInputs * V(:,HowManyMostVariantPrincipleComopnenetWantToUse:NumberOfInputVariables);

K-means clustering in C++

The following program use in-line k-means algorithm to train the centroids of the K-means clustering.

Before using this code you have to change data loading function in loadTrainingData.

The beauty of this program is you do not need to worry about what is the best value of K. You could just run for range of K values (ExperimentBestKValue) and and choose the best perform K value to re-train and predict the your target (ExperimentBestK). A kind of Evolutionary algorithm... don't you agree?

// KMeanClustring.cpp : Defines the entry point for the console application.
// Author: Rudra Poudel

//include libraries
#include "stdafx.h"
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<float.h>
#include<conio.h>
#include<time.h>

//******************************
//** Setting input parameter
//******************************
#define EXPERIMENTAL_TRIALS 1000 //Numbers of trails to estimate the average value of K's success
#define EXPERIMENTAL_TRIALS_BEST_K 1 //10000 //Numbers of trails to estimate the best centroids with K centroids
#define EXPERIMENTAL_K_FROM 2 //Starting K value for experiment
#define EXPERIMENTAL_K_TO 100 //End K value for experiment
#define TESTING_SAMPLE_SIZE 50 //Testing sample size will be taken form the end of the input data
#define NUMBER_OF_TESTING_DATA 100 //Number of testing data
//#define K 20 // value of K
#define LEARNING_RATE 0.1 // Learning rate
#define MAX_TRAIN_ITERATION 50000 //Maximum nos of iteration to train the networks
#define INPUTS_COUNT 365 // Number of input patterns i.e. training sets
#define NN 6 // Number of inputs/neurons (6 Inputs )

#define PRINT_OPTION 2
//********************************

//Random generator parameters
#define IA 16807
#define IM 2147483647
#define AM (1. /IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
long idum;

enum PrintOptions
{
poNone = 0, //Dont print to DAT file
poTraining = 2 //Printing whole training process
};

//Function decleration
void ExperimentBestKValue();
void ExperimentBestK(int exp_k);
float getRand(long *idum);
bool loadTrainingData(double input[INPUTS_COUNT][NN], double targets[INPUTS_COUNT]);
void moveTowards(double centroid[NN], double datum[NN]);
int getClosestCentroid(double datum[NN], double** centroids, int maxK);
double getDistance(double datum[NN], double centroid[NN]);
bool isOldAndNewSame(double** oldCentroids, double** centroids, int maxK);
double assignClusterTarget(double inputs[INPUTS_COUNT][NN], double** centroids, double targets[INPUTS_COUNT], int maxK, int testingSampleSize);
double predictTarget(double inputs[INPUTS_COUNT][NN], double** centroids, double targets[INPUTS_COUNT], int maxK, int testingSampleSize);


///<FunctionName>_tmain</FunctionName>
///<arg name='argc'>Size of argv[]</arg>
///<arg name='argv[]'>Array inputs</arg>
///<Summary>Main entry point of the function. It will train the all possible target patterns for NN-1 inputs for 2-class problem.</Summary>
///<Output> 0 for success, 1 for fail.</Output>
int _tmain(int argc, _TCHAR* argv[])
{
//ExperimentBestKValue();
//ExperimentBestK(EXPERIMENTAL_TRIALS_BEST_K);
ExperimentBestK(23);
return 0;
}

///<FunctionName>ExperimentBestKValue</FunctionName>
///<Summary>Function to experiment.</Summary>
///<Output>SuccessPecentage.dat- which contains average success percent for different value of K</Output>
void ExperimentBestKValue(){
int i,j, epoch, higestSuccessPecntOfAllK=-1;
double sumTotalSuccessPecent,higestSuccessPecnt=-100, higestSuccessPecntForCurrentK=-100,tempSuccessPecent;
char filename[25];
FILE *pSuccessPec;

//Input patterns, last one is always 1 as it is input for bias
double inputs[INPUTS_COUNT][NN];
double targets[INPUTS_COUNT];
double **centroids;
double **oldCentroids;
//int centroidsTargets[K]; //0 or 1 which ever greater
//int clusterPoints[K][2]; //0= Numbers of points; i=miss classified

idum = -time(0);
if(!loadTrainingData(inputs,targets))
return;

if(PRINT_OPTION==poTraining){
sprintf(filename,"SuccessPecentage.dat");
pSuccessPec = fopen(filename,"w");
}

for(int exp_k=EXPERIMENTAL_K_FROM;exp_k<=EXPERIMENTAL_K_TO;exp_k++){
centroids = (double **) malloc(exp_k * sizeof(double *));
oldCentroids = (double **) malloc(exp_k * sizeof(double *) );
for(i=0;i<exp_k;i++) {
centroids[i] = (double *) malloc(NN * sizeof(double));
oldCentroids[i] = (double *) malloc(NN * sizeof(double));
}

sumTotalSuccessPecent = 0;
higestSuccessPecntForCurrentK = -100;
for(int exp_trials=0;exp_trials<EXPERIMENTAL_TRIALS;exp_trials++){

for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
centroids[i][j] = getRand(&idum);
}
}

epoch=0;
do{//START: do ... while loop
for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
oldCentroids[i][j] = centroids[i][j];
}
}

for (i=0;i<(INPUTS_COUNT-NUMBER_OF_TESTING_DATA-TESTING_SAMPLE_SIZE);i++) {
moveTowards(centroids[getClosestCentroid(inputs[i],centroids,exp_k)], inputs[i]);
}
epoch++;
//printf("epoch: %d\n",epoch);

}while(!isOldAndNewSame(oldCentroids,centroids,exp_k) && epoch<MAX_TRAIN_ITERATION);
//END: do ... while loop
tempSuccessPecent = assignClusterTarget(inputs, centroids, targets,exp_k, TESTING_SAMPLE_SIZE);
if(higestSuccessPecntForCurrentK<tempSuccessPecent){
higestSuccessPecntForCurrentK = tempSuccessPecent;
}
if(higestSuccessPecnt<tempSuccessPecent){
higestSuccessPecnt = tempSuccessPecent;
higestSuccessPecntOfAllK = exp_k;
}
sumTotalSuccessPecent += tempSuccessPecent;
}//END: for ... EXPERIMENTAL_TRIALS
if(PRINT_OPTION==poTraining){
fprintf(pSuccessPec,"%d\t%f\t%f\n",exp_k, sumTotalSuccessPecent/EXPERIMENTAL_TRIALS,higestSuccessPecntForCurrentK);
}
for(i=0;i<exp_k;i++) {
free(centroids[i]);
free(oldCentroids[i]);
}
free(centroids);
free(oldCentroids);

printf("K = %d finished.\n", exp_k);
}//END: for ... EXPERIMENTAL_K_FROM


if(PRINT_OPTION==poTraining){
fclose(pSuccessPec);

sprintf(filename,"BestKandPecentage.dat");
FILE* pFile = fopen(filename,"w");
fprintf(pFile,"%d %f\n",higestSuccessPecntOfAllK,higestSuccessPecnt);
fclose(pFile);
}

/*printf("\n");
printf("\n\nPress any to to close the program...");
getche();*/
}




///<FunctionName>ExperimentBestK</FunctionName>
///<arg name='exp_K'>Value of K</arg>
///<Summary>Function to experiment with multi starting point of K and choose the best predicted K.</Summary>
///<Output>Training_Function.dat- Best K for clustring</Output>
void ExperimentBestK(int exp_k){
int i,j, epoch;
double successPecent, bestSuccessPercent=-100;
char filename[25];
FILE *pFile;

//Input patterns, last one is always 1 as it is input for bias
double inputs[INPUTS_COUNT][NN];
double targets[INPUTS_COUNT];
double **centroids;
double **oldCentroids;
double **bestCentroids;
//int centroidsTargets[K]; //0 or 1 which ever greater
//int clusterPoints[K][2]; //0= Numbers of points; i=miss classified

idum = -time(0);
if(!loadTrainingData(inputs,targets))
return;

if(PRINT_OPTION==poTraining){
sprintf(filename,"KMeanCentroids.dat");
pFile = fopen(filename,"w");
}

centroids = (double **) malloc(exp_k * sizeof(double *));
oldCentroids = (double **) malloc(exp_k * sizeof(double *) );
bestCentroids = (double **) malloc(exp_k * sizeof(double *) );
for(i=0;i<exp_k;i++) {
centroids[i] = (double *) malloc(NN * sizeof(double));
oldCentroids[i] = (double *) malloc(NN * sizeof(double));
bestCentroids[i] = (double *) malloc(NN * sizeof(double));
}

for(int exp_trials=0;exp_trials<EXPERIMENTAL_TRIALS_BEST_K;exp_trials++){

for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
centroids[i][j] = getRand(&idum);
}
}

epoch=0;
do{//START: do ... while loop
for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
oldCentroids[i][j] = centroids[i][j];
}
}

for (i=0;i<(INPUTS_COUNT-NUMBER_OF_TESTING_DATA);i++) {
moveTowards(centroids[getClosestCentroid(inputs[i],centroids,exp_k)], inputs[i]);
}
epoch++;
//printf("epoch: %d\n",epoch);

}while(!isOldAndNewSame(oldCentroids,centroids,exp_k) && epoch<MAX_TRAIN_ITERATION);
//END: do ... while loop
successPecent = assignClusterTarget(inputs, centroids, targets,exp_k,INPUTS_COUNT-NUMBER_OF_TESTING_DATA);
if(bestSuccessPercent<successPecent){
bestSuccessPercent = successPecent;
for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
bestCentroids[i][j] = centroids[i][j];
}
}

}
}//END: for ... EXPERIMENTAL_TRIALS

if(PRINT_OPTION==poTraining){

fprintf(pFile,"%d %f\n\n",exp_k, bestSuccessPercent);
//Printing Centroids
for(i=0;i<exp_k;i++){
for(j=0;j<NN;j++){
fprintf(pFile,"%f\t",bestCentroids[i][j]);
}
printf("\n");
}
}

successPecent = predictTarget(inputs, bestCentroids, targets,exp_k, INPUTS_COUNT-NUMBER_OF_TESTING_DATA);
if(PRINT_OPTION==poTraining){

fprintf(pFile,"\n%d %f",exp_k, successPecent);
printf("\n");
fclose(pFile);
}
for(i=0;i<exp_k;i++) {
free(centroids[i]);
free(oldCentroids[i]);
free(bestCentroids[i]);
}
free(centroids);
free(oldCentroids);
free(bestCentroids);

printf("\n");
printf("\n\nPress any to to close the program...");
getche();
}


///<FunctionName>getRand</FunctionName>
///<arg name='*idum'>Seed for random number generator</arg>
///<Summary>Function follows the Random Number Generators from Numerical Recipes in C, v2.0 (C) Copr. 1986-92 Numerical Recipes Software.</Summary>
///<Output>Random number between 0-1 range.</Output>
float getRand(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;

if (*idum <= 0 || !iy) {
if (-(*idum) < 1) *idum = 1;
else *idum = -(*idum);
for (j=NTAB+7; j>=0; j--) {
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum <0) *idum += IM;
if (j < NTAB) iv[j] = *idum;
}
iy = iv[0];
}
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j = iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}


///<FunctionName>loadTrainingData</FunctionName>
///<arg name='input[INPUTS_COUNT][NN]'>Training data</arg>
///<arg name='targets[INPUTS_COUNT]'>Target class for each input with same index number</arg>
///<Summary>Load the training data form txt file.</Summary>
///<Output>ture for success, false for fail.</Output>
bool loadTrainingData(double input[INPUTS_COUNT][NN], double targets[INPUTS_COUNT]){
int tempInt;
float tempFloat;
FILE *pData;

pData = fopen("ex1-data.txt","r");
if(pData==NULL) {
printf("Error: can't access data file.\n");
return false; // terminate with error
}

for(int i=0;i<INPUTS_COUNT;i++){
for(int j=0;j<NN;j++){
fscanf (pData, "%f", &tempFloat);
input[i][j]= tempFloat;
}
fscanf (pData, "%d", &tempInt);
targets[i]= tempInt;
}
fclose(pData);

return true;
}

///<FunctionName>getClosestCentroid</FunctionName>
///<arg name='datum[NN]'>One training data</arg>
///<arg name='** centroids'>K Centroids</arg>
///<Summary>Retrun the closest centroids of the given training data among K-centroids.</Summary>
///<Output>Index of closest centroids.</Output>
int getClosestCentroid(double datum[NN], double** centroids, int maxK) {
double distance;
double min = 99999999999;
int k = -1;
for (int i=0;i<maxK;i++) {
distance = getDistance(datum, centroids[i]);
if (distance < min) { k = i; min = distance; }
}
return(k);
}

///<FunctionName>getDistance</FunctionName>
///<arg name='datum[NN]'>One training data</arg>
///<arg name='centroid[NN]'>One Centroid</arg>
///<Summary>Retrun the distance betweent given training data and centroid.</Summary>
///<Output>Distance betweent given training data and centroid.</Output>
double getDistance(double datum[NN], double centroid[NN]) {
double distance = 0.0;
for (int i=0;i<NN;i++) {
distance += pow(datum[i] - centroid[i], 2);
}
return(sqrt(distance));
}


///<FunctionName>moveTowards</FunctionName>
///<arg name='centroid[NN]'>One Centroid</arg>
///<arg name='datum[NN]'>One training data</arg>
///<Summary>Change centroid postion according to distance between points and centroid and global learning rate.</Summary>
///<Output>Void. Node centroid changed by reference</Output>
void moveTowards(double centroid[NN], double datum[NN]) {
for (int i=0;i<NN;i++) {
centroid[i] += ((datum[i] - centroid[i]) * LEARNING_RATE); }
}

///<FunctionName>isOldAndNewSame</FunctionName>
///<arg name='** oldCentroids'>Old centroids</arg>
///<arg name='** centroids'>New centroids</arg>
///<arg name='maxK'>Number of K</arg>
///<Summary>Compare whether all old centroids and new centroids are same or not.</Summary>
///<Output>True if both centroids are same else false.</Output>
bool isOldAndNewSame(double** oldCentroids, double** centroids, int maxK) {
for(int i=0;i<maxK;i++){
for(int j=0;j<NN;j++){
if(oldCentroids[i][j]!=centroids[i][j])
return false;
}
}
return true;
}


///<FunctionName>assignClusterTarget</FunctionName>
///<arg name='inputs[INPUTS_COUNT][NN]'>Input values</arg>
///<arg name='** centroids'>Centroids</arg>
///<arg name='targets[INPUTS_COUNT]'>Targets</arg>
///<arg name='maxK'>Number of K</arg>
///<arg name='testingSampleSize'>Size of testing sample</arg>
///<Summary>Calculate the centroids targets.</Summary>
///<Output>Return the correct classified percent of the last TESTING_SAMPLE_SIZE data.</Output>
double assignClusterTarget(double inputs[INPUTS_COUNT][NN], double** centroids, double targets[INPUTS_COUNT], int maxK,int testingSampleSize){
int totalCorrectClassified=0, totalMissClassified=0;
int** centroidsTargetCounter; //0= Target 0; 1= Target 1
int* centroidsTarget; //Target as eighter 0 or 1

centroidsTargetCounter = (int**) malloc(sizeof(int *) * maxK);
centroidsTarget = (int*) malloc(sizeof(int) * maxK);
for(int i=0;i<maxK;i++) {
centroidsTargetCounter[i] = (int*) malloc(sizeof(double) * 2);
}

for(int i=0;i<maxK;i++){
for(int j=0;j<2;j++){
centroidsTargetCounter[i][j] = 0;
}
}
int testingAdjust;
if((INPUTS_COUNT-NUMBER_OF_TESTING_DATA)==testingSampleSize)
testingAdjust = 0;
else
testingAdjust = testingSampleSize;

for (int i=0;i<(INPUTS_COUNT-NUMBER_OF_TESTING_DATA-testingAdjust);i++) {
int nearestCentroid = getClosestCentroid(inputs[i],centroids,maxK);
if(targets[i]==0)
centroidsTargetCounter[nearestCentroid][0]++;
else
centroidsTargetCounter[nearestCentroid][1]++;
}

for (int i=0;i<maxK;i++) {
if(centroidsTargetCounter[i][0]>centroidsTargetCounter[i][1]){
centroidsTarget[i]=0;
}
else{
centroidsTarget[i]=1;
}
}

for (int i=(INPUTS_COUNT-NUMBER_OF_TESTING_DATA-testingSampleSize);i<(INPUTS_COUNT-NUMBER_OF_TESTING_DATA);i++) {
int nearestCentroid = getClosestCentroid(inputs[i],centroids, maxK);
if(targets[i]==centroidsTarget[nearestCentroid])
totalCorrectClassified++;
else
totalMissClassified++;
}

for(int i=0;i<maxK;i++) {
free(centroidsTargetCounter[i]);
}
free(centroidsTargetCounter);
free(centroidsTarget);

//printf("%f\n",( ((double)totalCorrectClassified/testingSampleSize) * 100));

return ( ( (((double)totalCorrectClassified)/testingSampleSize) * 100));
}

///<FunctionName>predictTarget</FunctionName>
///<arg name='inputs[INPUTS_COUNT][NN]'>Input values</arg>
///<arg name='** centroids'>Centroids</arg>
///<arg name='targets[INPUTS_COUNT]'>Targets</arg>
///<arg name='maxK'>Number of K</arg>
///<arg name='testingSampleSize'>Size of testing sample</arg>
///<Summary>Calculate the centroids targets.</Summary>
///<Output>Void.</Output>
double predictTarget(double inputs[INPUTS_COUNT][NN], double** centroids, double targets[INPUTS_COUNT], int maxK, int testingSampleSize){
int totalCorrectClassified=0, totalMissClassified=0;
int** centroidsTargetCounter; //0= Target 0; 1= Target 1
int* centroidsTarget; //Target as eighter 0 or 1

centroidsTargetCounter = (int**) malloc(sizeof(int *) * maxK);
centroidsTarget = (int*) malloc(sizeof(int) * maxK);
for(int i=0;i<maxK;i++) {
centroidsTargetCounter[i] = (int*) malloc(sizeof(double) * 2);
}

for(int i=0;i<maxK;i++){
for(int j=0;j<2;j++){
centroidsTargetCounter[i][j] = 0;
}
}

int testingAdjust;
if(INPUTS_COUNT-NUMBER_OF_TESTING_DATA==testingSampleSize)
testingAdjust = 0;
else
testingAdjust = testingSampleSize;

for (int i=0;i<INPUTS_COUNT-NUMBER_OF_TESTING_DATA-testingAdjust;i++) {
int nearestCentroid = getClosestCentroid(inputs[i],centroids,maxK);
if(targets[i]==0)
centroidsTargetCounter[nearestCentroid][0]++;
else
centroidsTargetCounter[nearestCentroid][1]++;
}

for (int i=0;i<maxK;i++) {
if(centroidsTargetCounter[i][0]>centroidsTargetCounter[i][1]){
centroidsTarget[i]=0;
totalCorrectClassified +=centroidsTargetCounter[i][0];
}
else{
centroidsTarget[i]=1;
totalCorrectClassified +=centroidsTargetCounter[i][1];
}
}


char filename[20];
FILE* pFile;
sprintf(filename,"Prediction.txt");
pFile = fopen(filename,"w");
for (int i=(INPUTS_COUNT-NUMBER_OF_TESTING_DATA);i<INPUTS_COUNT;i++) {
fprintf(pFile,"%d\n", centroidsTarget[getClosestCentroid(inputs[i],centroids, maxK)]);
}
fclose(pFile);

for(int i=0;i<maxK;i++) {
free(centroidsTargetCounter[i]);
}
free(centroidsTargetCounter);
free(centroidsTarget);

//printf("%f\n",( ((double)totalCorrectClassified/testingSampleSize) * 100));

return ( ((double)totalCorrectClassified/testingSampleSize) * 100);
}

Backpropagation in Matlab

Following code is written in Matlab using back propagation algorithm for error propagation to hidden layers. It used batch update mode and has ability to work with N numbers of hidden layers also as many as required numbers of nodes in each hidden layer.

%References: Patterson, D. [Online]Available from: http://www.csee.umbc.edu/~dpatte3/nn/backprop.html [accessed 24 February 2009]
function MLPUsingBackpropagation()
[x1,x2,x3,x4,x5,x6,t] = textread('ex1-data.txt','%f %f %f %f %f %f %f');

trainingIndex=0;
testingIndex=0;
for i=1:265
if (rand > 0.2)
trainingIndex = trainingIndex + 1;
trainingInputs(trainingIndex,1) = x1(i);
trainingInputs(trainingIndex,2) = x2(i);
trainingInputs(trainingIndex,3) = x3(i);
trainingInputs(trainingIndex,4) = x4(i);
trainingInputs(trainingIndex,5) = x5(i);
trainingInputs(trainingIndex,6) = x6(i);

trainingOutputs(trainingIndex,1) = t(i);

else
testingIndex = testingIndex + 1;
testingInputs(testingIndex,1) = x1(i);
testingInputs(testingIndex,2) = x2(i);
testingInputs(testingIndex,3) = x3(i);
testingInputs(testingIndex,4) = x4(i);
testingInputs(testingIndex,5) = x5(i);
testingInputs(testingIndex,6) = x6(i);

testingOutput(testingIndex,1) = t(i);

end

end
%***** Memory Clean up **********
clear x1
clear x2
clear x3
clear x4
clear x5
clear x6
clear t
%********************************
layers = [6 23 23 1];

for i=1:50

Output = BackpropagationTraining(trainingInputs,trainingOutputs,0.1,0.05,0.01,layers);

OutputTesting = BackpropagationTesting(testingInputs, testingOutput, layers, Output.weights);
OutputTesting = BackpropagationTesting(testingInputs, testingOutput, layers, Output.bestweights);
end

%***** Memory Clean up **********
clear trainingInputs
clear testingInputs
clear layers
%********************************

end

function Results = BackpropagationTraining(inputs,targets,learningRate,momentumWeight,minimumMSE,layers)

[trainingDataCount,inputNodesCount] = size(inputs);
[targetDataCount,targetNodesCount] = size(targets);

if trainingDataCount ~= targetDataCount
error('BackpropagationTraining:TrainingAndTargetDataLengthMismatch', 'The number of input vectors and desired ouput vectors do not match');
end

if length(layers) < 3
error('BackpropagationTraining:InvalidNetworkStructure','The network must have at least 3 layers');
end

if inputNodesCount ~= layers(1)
msg = sprintf('Dimensions of input nodes (%d) does not match with numbers of input layer (%d).',inputNodesCount,layers(1));
error('BackpropagationTraining:InvalidInputLayerSize', msg);
end

if targetNodesCount ~= layers(end)
msg = sprintf('Dimensions of output nodes (%d) does not match with numbers of output layer (%d)',targetNodesCount,layers(end));
error('BackpropagationTraining:InvalidOutLayerSize', msg);
end

leayersLength = length(layers);

weights = cell(leayersLength-1,1);
for i=1:leayersLength-2
weights{i} = [-15 + 30 .* rand(layers(i+1),layers(i)+1); zeros(1,layers(i)+1)];
end
weights{end} = -15 + 30 .* rand(layers(end),layers(end-1)+1);

MSE = Inf;
epochs = 0;

activation = cell(leayersLength,1);
activation{1} = [inputs ones(trainingDataCount,1)]; % activation{1} is the input + 1 for the bias node activation
% activation{1} remains the same throught the computation
for i=2:leayersLength-1
activation{i} = ones(trainingDataCount,layers(i)+1); % inner layers include a bias node (trainingDataCount-by-Nodes+1)
end
activation{end} = ones(trainingDataCount,layers(end)); % no bias node at output layer

net = cell(leayersLength-1,1); % one net matrix for each layer exclusive input
for i=1:leayersLength-2;
net{i} = ones(trainingDataCount,layers(i+1)+1); % affix bias node
end
net{end} = ones(trainingDataCount,layers(end));

previousDeltaW = cell(leayersLength-1,1);
sumDeltaW = cell(leayersLength-1,1);
for i=1:leayersLength-1
previousDeltaW{i} = zeros(size(weights{i})); % previousDeltaW starts at 0
sumDeltaW{i} = zeros(size(weights{i}));
end

lowestsse = 999999;
bestweights = weights;

while MSE > minimumMSE && epochs < 3000

for i=1:leayersLength-1
net{i} = activation{i} * weights{i}'; % compute inputs to current layer

if i < leayersLength-1 % inner layers
%a{i+1} = [2./(1+exp(-net{i}(:,1:end-1)))-1
%ones(trainingDataCount,1)]; % for bi-polar
activation{i+1} = [1./(1+exp(-net{i}(:,1:end-1))) ones(trainingDataCount,1)]; % for sigmoid
%activation{i+1} = [(net{i}(:,1:end-1)) ones(trainingDataCount,1)]; %without sigmoid i.e for linear activation function
else % output layers
%a{i+1} = 2 ./ (1 + exp(-net{i})) - 1;
%activation{i+1} = 1 ./ (1 + exp(-net{i}));
%a{i+1} = net{i};
activation{end} = net{i} > 0;
end
end

% calculate sum squared error of all samples
err = (targets-activation{end}); % save this for later
sse = sum(sum(err.^2)); % sum of the error for all samples, and all nodes
if(lowestsse>sse)
lowestsse = sse;
bestweights = weights;
end

%delta = err .* (1 + a{end}) .* (1 - a{end});
%delta = err .* a{end} .* (1 - a{end});
%delta = err;
for
delta = err + weights{i}' weights{i};
for i=leayersLength-1:-1:1
sumDeltaW{i} = learningRate * delta' * activation{i};
if i > 1
%delta = (1+a{i}) .* (1-a{i}) .* (delta*weights{i});
delta = activation{i} .* (1-activation{i}) .* (delta*weights{i});
%delta = (delta*weights{i}); % where there is no activation
%function
end
end

% update the prev_w, weight matrices, epoch count and MSE
for i=1:leayersLength-1
previousDeltaW{i} = (sumDeltaW{i} ./ trainingDataCount) + (momentumWeight * previousDeltaW{i});
weights{i} = weights{i} + previousDeltaW{i};
end
epochs = epochs + 1;
MSE = sse/(trainingDataCount*targetNodesCount); % MSE = 1/trainingDataCount * 1/targetNodesCount * summed squared error
%MSEPerEpoch(epochs) = MSE;

end

%printing error
% hold on
% axis ( [0 3000 0 1]);
% x = 1:1:3000;
% plot(x,MSEPerEpoch);
% hold off
% s=input('Press enter to continue, enter 0 to stop \n');

% return the trained network
Results.structure = layers;
Results.weights = weights;
Results.bestweights = bestweights;
Results.epochs = epochs;
Results.MSE = MSE;
sse
lowestsse

end

function Results = BackpropagationTesting(inputs,targets,layers, weights)

[trainingDataCount,inputNodesCount] = size(inputs);
[targetDataCount,targetNodesCount] = size(targets);

if trainingDataCount ~= targetDataCount
error('BackpropagationTraining:TrainingAndTargetDataLengthMismatch', 'The number of input vectors and desired ouput vectors do not match');
end

leayersLength = length(layers);

activation = cell(leayersLength,1);
activation{1} = [inputs ones(trainingDataCount,1)]; % activation{1} is the input + 1 for the bias node activation
% activation{1} remains the same throught the computation
for i=2:leayersLength-1
activation{i} = ones(trainingDataCount,layers(i)+1); % inner layers include a bias node (trainingDataCount-by-Nodes+1)
end
activation{end} = ones(trainingDataCount,layers(end)); % no bias node at output layer

net = cell(leayersLength-1,1); % one net matrix for each layer exclusive input
for i=1:leayersLength-2;
net{i} = ones(trainingDataCount,layers(i+1)+1); % affix bias node
end
net{end} = ones(trainingDataCount,layers(end));

for i=1:leayersLength-1
net{i} = activation{i} * weights{i}'; % compute inputs to current layer

if i < leayersLength-1 % inner layers
%a{i+1} = [2./(1+exp(-net{i}(:,1:end-1)))-1 ones(trainingDataCount,1)];
activation{i+1} = [1./(1+exp(-net{i}(:,1:end-1))) ones(trainingDataCount,1)];
%a{i+1} = [(-net{i}(:,1:end-1)) ones(trainingDataCount,1)];
else % output layers
%a{i+1} = 2 ./ (1 + exp(-net{i})) - 1;
%activation{i+1} = 1 ./ (1 + exp(-net{i}));
activation{end} = net{i} > 0;
%a{i+1} = net{i};
end
end

%activation{end} = activation{end} > 0.5;
% calculate sum squared error of all samples
err = (targets-activation{end}); % save this for later
sse = sum(sum(err.^2)); % sum of the error for all samples, and all nodes

MSE = sse/(trainingDataCount*targetNodesCount); % MSE = 1/trainingDataCount * 1/targetNodesCount * summed squared error

% return the trained network
Results.MSE = MSE;
Results.sse = (sse/trainingDataCount)*100;
100-((sse/trainingDataCount)*100)

end

Genetic Algorithm

A simple genetic algorithm, It follows the Inman Harvey's microbial genetic algorithm.

Problem:

You have 10 cards numbered from 1 to 10. You have to choose a way of dividing them into 2 piles, so that the cards in Pile_0 *sum* to a number as close as possible to 36, and the remaining cards in Pile_1 *multiply* to a number as close as possible to 360.


Solution:


// CardProblembyGA.cpp : Defines the entry point for the console application.
#include "stdafx.h"
//#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <conio.h>

#define POP 1000
#define LEN 10
#define REC 0.5
#define MUT 0.1
#define END 1000

#define TARGET_SUM 36
#define TARGET_PROD 360

#define RANDOM_NUM ((float)rand()/(RAND_MAX+1))

void init_pop(char gene[POP][LEN]);
double evaluate(char gene[LEN]);
void display(char gene[LEN]);

int _tmain(int argc, _TCHAR* argv[])
{
char gene[POP][LEN];
int Winner, Loser, temp1, temp2;
bool WinnerFound = false;

srand((int)time(NULL));
init_pop(gene);
for(int t=0;t<END && !WinnerFound;t++){
temp1 = (int) (POP * RANDOM_NUM);
temp2 = (int) (POP * RANDOM_NUM);
if(evaluate(gene[temp1]) < evaluate(gene[temp2])){
Winner = temp1;
Loser = temp2;
}
else {
Winner = temp2;
Loser = temp1;
}

for(int i=0;i<LEN;i++){
if(RANDOM_NUM < REC)
gene[Loser][i] = gene[Winner][i];
if(RANDOM_NUM < MUT)
gene[Loser][i] = gene[Loser][i]=='1' ? '0':'1';

if(evaluate(gene[Loser]) == 0.0){
printf("\nIdeal condtion found on %d tournaments...", t);
//WinnerFound = true;
display(gene[Loser]);
break;
}
}
}
/*if(!WinnerFound)
printf("Unable to find ideal condtion on %d tournaments...", END);*/
getche();
return 0;
}

void init_pop(char gene[POP][LEN])
{
for(int i = 0; i<POP; i++){
for(int j = 0; j<LEN; j++){
if(RANDOM_NUM > 0.5f)
gene[i][j] = '1';
else
gene[i][j] = '0';
}
}
}

double evaluate(char gene[LEN])
{
int sum = 0, prod = 1;
double sum_deviation, prod_deviation, combined_deviation;

for(int i=0;i<LEN;i++){
if(gene[i]=='1')
prod *=(i+1);
else
sum +=(i+1);
}
sum_deviation = ((double)(sum-TARGET_SUM))/TARGET_SUM;
prod_deviation = ((double)(prod-TARGET_PROD))/TARGET_PROD;
combined_deviation = abs(sum_deviation) + abs(prod_deviation);
return combined_deviation;
}

void display(char gene[LEN])
{
int x = 0;
printf("\n Product pile elements:");
for(int i=0;i<LEN;i++){
if(gene[i]=='1')
printf(" %d",i+1);
}
printf("\n Sum pile elements:");
for(int i=0;i<LEN;i++)
{
if(gene[i]=='0')
printf(" %d",i+1);
}
}

Single Layer Perceptron

Below is the code written in C++ to train the single layer perceptron (single layer network) for 2-class classification.

// SingleLayerPerceptron.cpp : Defines the entry point for the console application.
// Single Layer Perceptron training example code
// Author: Rudra P.K. Poudel

#include "stdafx.h"
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<float.h>
#include<conio.h>

#define NN 3 // Number of neurons (2 Inputs + 1 bias )
#define LEARNING_RATE 0.1 // Learning rate
#define INPUTS_COUNT 4 // Number of input patterns i.e. training sets
#define MAX_TRAIN_ITERATION 10000 //Maximum nos of iteration to train the networks

#define PRINT_OPTION 1
//Enum flag to control the output print
enum PrintOptions
{
poNone = 0, //Dont print any output
poTraining = 1 // Printing whole training process
};

//Function decleration
int PerceptronOutput(const double weights[NN], const int inputs[NN]);
double NodeOutput(const double weights[NN], const int inputs[NN]);

//Main entry point of the function
//It will train the all possible target patterns for NN (number of nodes) inputs for 2-class problem
//This program training single layer perceptron for the XAND operation
int _tmain(int argc, _TCHAR* argv[])
{
int i,j, iterations, trainFunction;
double totalError, singleError, output;
char filename[25];
FILE *pFile;

//Input patterns, last one is always 1 as it is input for bias
int inputs[][NN]= {
{0,0,1},
{0,1,1},
{1,0,1},
{1,1,1}
};

/*int targets[INPUTS_COUNT] = {1,1,1,0};*/
//Possible combination of targets class with target +1,-1
int targets[INPUTS_COUNT] = {1,1,1,-1};

if(PRINT_OPTION==poTraining){
sprintf(filename,"Training_Function.dat");
pFile = fopen(filename,"w");
}

double weights[]= {0,0,0};
iterations=0;
//Training Perceptron with adaline learning algorithm
do{
totalError=0;
for(i=0;i<INPUTS_COUNT;i++){
output = PerceptronOutput(weights,inputs[i]);
singleError = targets[i] - output;

if(PRINT_OPTION==poTraining){
//Printing input values
for(j=0;j<NN;j++){
fprintf(pFile,"%d\t",inputs[i][j]);
}
//Printing weights
for(j=0;j<NN;j++){
fprintf(pFile,"%f\t",weights[j]);
}
//Printing target output
fprintf(pFile,"%d\t",targets[i]);
//Printing node output
fprintf(pFile,"%f\t",NodeOutput(weights,inputs[i]));
//Printing Perceptron output
fprintf(pFile,"%f\t",output);
//Printing error
fprintf(pFile,"%f\t",singleError);
//Printing adjustment
fprintf(pFile,"%f\t", ((double)LEARNING_RATE) * singleError);

}
if(singleError!=0){
for(j=0;j<NN;j++){
weights[j] = weights[j] + ((double)LEARNING_RATE) * inputs[i][j] * singleError;
}
}
if(PRINT_OPTION==poTraining){
//Printing new weights
for(j=0;j<NN;j++){
fprintf(pFile,"%f\t",weights[j]);
}
fprintf(pFile,"\n");
}
totalError += fabs(singleError);
}
iterations++;
printf("Iterations: %d\tTotal Error= %f\n",iterations,totalError);
}while(totalError!=0 && iterations<MAX_TRAIN_ITERATION);
//END- While loop to train one function

if(PRINT_OPTION==poTraining)
fclose(pFile);

printf("\n\nLearned weights:\n");
for(j=0;j<NN;j++){
printf("%f\t",weights[j]);
}
printf("\n\nPress any to to close the program...");
getche();
return 0;
}

//Function Name: PerceptronOutput
//Input: Weights of the connections and inputs to the node
//Putput: 1 if <NodeOput> is >=0 else -1
int PerceptronOutput(const double weights[NN], const int inputs[NN]){

if(NodeOutput(weights,inputs)>=0)
return 1;
else
return -1;
}

//Function Name: NodeOutput
//Input: Weights of the connections and inputs to the node
//Putput: Give sumation of weight and input sensor's product
double NodeOutput(const double weights[NN], const int inputs[NN]){
double output=0;

for(int i=0;i<NN;i++)
output +=weights[i]*inputs[i];

return output;
}

Random number and Gaussian random number algorithm / code / generatior in C++

Actually, I loose few hours to search the code for random number generator which follows the Gaussian distribution with given mean and variance. Hence I am posting my final google and some hard work's result because I think this will help you to save more hours.

C++ Code (I have used Microsoft Visual Studio 2005). PLEASE PAY ATTENTION IN "stdafx.h" IF YOU ARE NOT USING Microsoft Visual Studio 2005:


// GussianDistributionGeneratior.cpp : Defines the entry point for the console application.
// Author: Rudra Poudel

#include "stdafx.h"
#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<limits.h>
#include<time.h>

//Random generator parameters
#define IA 16807
#define IM 2147483647
#define AM (1. /IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
long idum;

float getRand(long *idum);
float getGaussianRand(long *idum, double mean, double sigma);

int _tmain(int argc, _TCHAR* argv[])
{
idum= -time(0);

printf("Random Numbers:\n");
for(int i=0;i<100;i++){
printf("%f\t",getRand(&idum));
}
printf("\n\nPress any key to continue...");
getche();

printf("\n\nGaussian Random Numbers with mean=0 and variance=1:\n");
for(int i=0;i<100;i++){
printf("%f\t",getGaussianRand(&idum,4,1));
}
printf("\n\nPress any key to continue...");
getche();

return 0;
}

//Below function follows the Random Number Generators from Numerical Recipes in C, v2.0 (C) Copr. 1986-92 Numerical Recipes Software

/************************************************************************
* name: getRand & getGussianRand
* description: random generators from
* input: idum: Seed for random number generator
* output: Random number between 0-1 range
************************************************************************/
float getRand(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;

if (*idum <= 0 || !iy) {
if (-(*idum) < 1) *idum = 1;
else *idum = -(*idum);
for (j=NTAB+7; j>=0; j--) {
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum <0) *idum += IM;
if (j < NTAB) iv[j] = *idum;
}
iy = iv[0];
}
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j = iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}

/************************************************************************
* name: getGussianRand
* description: Gaussian random generators with given random number generator's seed, mean and variance
* input: idum: Seed for random number generator
* mean: mean for the distribution; for standard normal distribution mean = 0
* sigma: variance (square of standard deviation) for standard normal distribution sigma=1. As must of the case people need Gaussina Random number with mean =0 and sigma=1, hence if you are not sure just use that or plot the points and observer with different values.
* output: A random value form gussian distribution with given parameters
************************************************************************/
float getGaussianRand(long *idum, double mean, double sigma)
{
float getRand(long *idum);
static int iset=0;
static float gset;
float fac,rsq,v1,v2;

if (iset == 0) {
do {
v1=2.0*getRand(idum)-1.0;
v2=2.0*getRand(idum)-1.0;
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq==0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return ( (v2*fac*sigma) + mean);
} else {
iset =0;
return ( (gset*sigma) + mean);
}
}



As I am running out of time so I am unable to write description and theory behind it. Please don't hesitate to drop comment or question if you find something strange!!!

Analysis of Continuous-Time Recurrent Neural Network

Recent year’s research shows that Continuous-Time Recurrent Neural Network (CTRNN) is important for modelling dynamic behaviour of the system. Hence detailed analysis and understanding of CTRNN is important. As such, attempt is made in this report to analyse the CTRNN equation for single node and double nodes. Basically, I would be experimenting and analysing how individual parameter of the CTRNN affects the output of the system, and I would also try to report the relationship between different parameters as well. Finally, I would be discussing the stability of the CTRNN with respect to the different values of its parameters.

Download full pdf text

Hebbian Learning using Fixed Weight Evolved Dynamical ‘Neural’ Networks

In connectionist artificial neural network, the followers of the Hebbian learning (Hebb 1949) strongly believe that activities between two nodes can be increased or decreased by changing the connection weight between them. However, we claim and prove that learning can be produced without changing the connection weight in dynamic neural network as we believe that learning is produced by interaction with dynamic environment. To prove this I have reworked on research work of Izquierdo & Harvey (2007). However, I have achieved low fitness than their and my best evolved 4-node circuit is not as good as their in task achievement. I have used evolutionary approach to synthesize Continuous-Time Recurrent Neural Network parameters. The experimental methodology and output of the experiment have been described in details.

Download full pdf text

AquaLogic Collaboration DB connection fail with SQL Server

If you are using SQL Server as a backend for AquaLogic or Plumtree, it's diagnostic test will be fail due to schema name mismatch in configuration file and SQL Server's schema name for the Collaboration DB. To correct this problem, we need to edit the database.xml files by following steps,

- find database.xml files on AquaLogic installation directory [eg. C:\bea]

- change DB_ali_Name.DB_User_Name.'s [DB_User_Name] to [dbo]

After that diagnostic should be success as well as system will work fine. Don’t forget to restart services after xml file changes.

Action Selection using Reinforcement Learning

Success of adaptive autonomous agent is evaluated by how well it could perform the desired task in dynamic environment, i.e. selection of appropriate action even if circumstances and environment changes. Agent cannot select appropriate action in random as well as in sequence because the environment in which agent acts are dynamic in nature. Hence, action selection always remains a central question, especially in the field of adaptive autonomous agents that can function robustly and efficiently in complex and dynamic environment (Blumberg 1994, p.22). Now in this situation, reinforcement learning, getting feedback (negative or positive) of performed action has very important role in decision making of future action selection with the help of present experience. Hence, though there are few difficulties in implementation of concept of reinforcement learning, it is still the first choice in adaptive autonomous agent building. Especially in this essay I will describe how reinforcement learning helps agent for efficient action selection in dynamic environment and I will also highlight the importance of reinforcement learning in action selection as well as its limitations.

Download full pdf text