Kmeans Clustering Using Multilayer Perceptron (MLP) in CUDA

                                
// mlpcu.h
// https://forejune.co/cuda

#ifndef __MLPCU_H__
#define __MLPCU_H__

using namespace std;

#define n1 3 
#define m1 7
#define K  4
#define numSamples 600
#define eta 0.5

// Simple MLP class
class MLP
{
private:
    double a[m1];  //linear sum of products of inputs and weights
    double h[m1];  //hidden nodes h[j] = g(a[j])
    double y[K];   //predicted output 
    double z[K];   //linear sum of products of weights and hidden nodes
    double *d_inputs;
    double *d_labels;
    double *d_a;
    double *d_h;        //hidden layer
    double *d_y;        //output
    double *d_yd;       //desired output
    double *d_z;
    double *d_w;
    double *d_wo;
public:
    MLP();
    //Sigmoid activation function
    double g(double x);
    //Derivative of sigmoid function
    double gd(double x);
    //Forward propagation
    double* predicted(double x[]);
    void train(const vector< Point> &points, int epochs);
    ~MLP();
};
#endif

                                
// mlp.cu
// https://forejune.co/cuda

#include <iostream>
#include <iomanip>
#include <cmath>
#include <stdio.h>
#include "kmeanscu.h"
#include "mlpcu.h"

using namespace std;

double w[n1][m1];
double wo[m1][K];

//Sigmoid activation function
__device__ double d_g(double x) 
{
      return 1.0 / (1.0 + exp(-x));
}
    
//Derivative of sigmoid function
__device__ double d_gd(double x) 
{
    double s = d_g(x);
     
    return s * (1 - s);
}

//y[] is the output
__global__ void d_forward(int k, double *w, double *wo, double *h, double *a, 
                               double *z,  double *y, double *d_inputs)
{
  __shared__ double x[n1];
  int ii = threadIdx.x;		//indexing input nodes
  int jj, ll;
  if (ii < n1) 	//inputs
    x[ii] = d_inputs[k*n1+ii];
  __syncthreads();

  //hidden layer activation, all hidden layer threads work in parallel
  if ( ii >= n1 && ii < n1+m1) {
    jj = ii - n1;
    a[jj] = 0;
    for (int i = 0; i < n1; i++) 
       a[jj] += w[i*m1+jj] * x[i]; 
    if ( jj > 0 )
       h[jj] = d_g(a[jj]); 	//sigmoid activation function
  }
  __syncthreads();

   //output layer activation, all output layer threads work in parallel
   if (ii >= n1 + m1){
     ll = ii - n1 - m1;
     z[ll] = 0; 
     for(int j = 0; j < m1; j++) 
       z[ll] += wo[j*K+ll] * h[j];  
      y[ll] = d_g( z[ll] );
   }
   __syncthreads();
}
  
//Backward propagation, yd is the desired output
__global__  void d_backward(int k, double *w, double *wo,  double *h, double *a, 
		            double *z, double *y, double *d_inputs, double *d_labels)
{
  __shared__ double delta[m1];
  __shared__ double deltao[K];
  __shared__ double w_share[n1][m1];	//just for easier readability, not necessarily needed
  __shared__ double wo_share[m1][K];
  __shared__ double x[n1];	//inputs
  __shared__ double yd[K];	//desired outputs
  int ii = threadIdx.x;		//indexing input nodes
  int jj = threadIdx.y;		//indexing hidden nodes
  int ll  = threadIdx.z;	//indexing output nodes
  if (jj == 0 && ll == 0 )
    x[ii] = d_inputs[k*n1+ii];
  __syncthreads();
  if ( ii == 0 && jj == 0 )
    yd[ll] = d_labels[k*K+ll];
  __syncthreads();
  if (ll == 0 )
    w_share[ii][jj] = w[ii*m1+jj]; 
  __syncthreads();
  if (ii == 0 )
    wo_share[jj][ll] = wo[jj*K+ll]; 
  __syncthreads();
  //output threads work in parallel 
  if (ii == 0 && jj == 0 ) {
  double e = yd[ll] - y[ll];  //output layer error
  deltao[ll] = e * d_gd(z[ll]);
  }
  __syncthreads();

   //Compute hidden layer error, hidden threads work in parallel
  if (ii == 0 && ll == 0 ) {
   delta[jj] = 0;
   for(int l = 0; l < K; l++)
     delta[jj] += deltao[l] * wo_share[jj][l] * d_gd(a[jj]);
  }
     //delta[jj] += deltao[l] * wo[jj*K+l] * d_gd(a[jj]);
    __syncthreads();
    
    //Update weights (hidden to output), hidden and output threads work in parallel 
  if (ii == 0 ) {
    wo_share[jj][ll] += eta * deltao[ll] * h[jj];
  }
    __syncthreads();
     
    //Update weights (input to hidden), input and hidden threads work in parallel
    if ( ll==0 ) {
      w_share[ii][jj] += eta * delta[jj] * x[ii];
    } 
    syncthreads();
   if ( ll == 0 ) 
    w[ii*m1+jj] =  w_share[ii][jj]; 
    __syncthreads();

   if ( ii == 0 ) 
    wo[jj*K+ll] = wo_share[jj][ll]; 
    __syncthreads();
} 

MLP::MLP()
{
	int wSize = n1 * m1 * sizeof(double);
	int woSize = m1 * K * sizeof(double);
	cudaMalloc((void**)&d_inputs, numSamples * n1 * sizeof(double));
	cudaMalloc((void**)&d_labels, numSamples * K * sizeof(double));
	cudaMalloc((void**)&d_a, m1 * sizeof(double));
	cudaMalloc((void**)&d_h, m1 * sizeof(double));
	cudaMalloc((void**)&d_z, K * sizeof(double));
	cudaMalloc((void**)&d_y, K * sizeof(double));
	cudaMalloc((void**)&d_yd, K * sizeof(double));
	cudaMalloc((void**)&d_w, wSize);
	cudaMalloc((void**)&d_wo, woSize);

         //random weights
        for (int i = 0; i < n1; i++)
          for (int j = 0; j < m1; j++)
            w[i][j] =  rand() % 1000 / 1000.0;

        for (int j = 0; j < m1; j++)
          for (int l = 0; l < K; l++)
            wo[j][l] = rand() % 1000 / 1000.0;

        cudaMemcpy(d_w, w, wSize, cudaMemcpyHostToDevice);
        cudaMemcpy(d_wo, wo, woSize, cudaMemcpyHostToDevice);
	h[0] = 1;	//for bias 
        cudaMemcpy(d_h, h, m1 * sizeof(double), cudaMemcpyHostToDevice);
}

//Sigmoid activation function
double MLP::g(double x) 
{
      return 1.0 / (1.0 + exp(-x));
}

//Derivative of sigmoid function
double MLP::gd(double x) 
{
       double s = g(x);
        
       return s * (1 - s);
}

//Forward propagation
double* MLP::predicted(double x[]) 
{
	//hidden layer activation
	for (int j = 0; j < m1; j++) {
          a[j] = 0;
          for (int i = 0; i < n1; i++) 
            a[j] += w[i][j] * x[i]; 
	   if ( j > 0 )
	     h[j] = g( a[j] );
         }
	 //output layer activation
	 for (int l = 0; l < K; l++){
	   z[l] = 0; 
	   for(int j = 0; j < m1; j++) 
	     z[l] += wo[j][l] * h[j];  
	   y[l] = g( z[l] );
	  }

          return y;  //return output
} 

// Train the perceptron
// eta = learning rate
// epochs = number of times of training 
// numSamples = number of different input sets (x_data)
// yd is target output
void MLP::train(const vector< Point> &points, int epochs)
{
	double inputs[numSamples][n1];
	double labels[numSamples][K];
        dim3 threads(n1, m1, K);   //input, hidden, output layers threads
	
	for (int i = 0; i < numSamples; i++) {
	  inputs[i][0] = 1;		//bias
	  inputs[i][1] = points[i].x;
	  inputs[i][2] = points[i].y;
	}
	
	for (int i = 0; i < numSamples; i++) {
          for (int k = 0; k < K; k++)
            if (points[i].cluster == k)
              labels[i][k] = 1;
            else
              labels[i][k] = 0;
	}
 
	cudaMemcpy(d_inputs, inputs, numSamples*n1*sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(d_labels, labels, numSamples*K*sizeof(double), cudaMemcpyHostToDevice);
	for (int epoch = 0; epoch < epochs; epoch++){
	  for (int i = 0; i < numSamples; i++) {
	    //for simplicity, use more threads than needed
	    d_forward<<<1,n1+m1+K>>>(i, d_w, d_wo, d_h, d_a, d_z, d_y, d_inputs);
  	    cudaDeviceSynchronize();
	    //for simplicity, use more threads than needed
	    d_backward<<<1, threads>>>(i, d_w, d_wo, d_h, d_a, d_z, d_y, d_inputs, d_labels);
  	    cudaDeviceSynchronize();
	  }
	}
	  //save the weights
	  cudaMemcpy(&w[0][0], d_w, sizeof(w), cudaMemcpyDeviceToHost);
	  cudaMemcpy(&wo[0][0], d_wo, sizeof(wo), cudaMemcpyDeviceToHost);
}

MLP::~MLP()
{
	cudaFree(d_inputs);
	cudaFree(d_labels);
	cudaFree(d_a);
	cudaFree(d_h);
	cudaFree(d_z);
	cudaFree(d_y);
}