rbf-dda

// mt4 implementation of rbf-dda as described in:
// http://citeseer.ist.psu.edu/cache/papers/cs/6810/http:zSzzSzhttp.cs.berkeley.eduzSz~bertholdzSzFTPzSzberthold.neuroc98.pdf/berthold98constructive.pdf

#define MAX_I 1000  // Max number of samples
#define MAX_E 10    // Max number of dimensions in input
#define M_PI 3.141592653589793

int init()
  {
  TestSpiral();
  return(0);
  }

int deinit()
  {
  return(0);
  }

int start()
  {
  return(0);
  }

// Warning ! this version works with only two classes: -1 and +1
void DDALearn(int Nbi, int Nbe, double Spi[MAX_I][MAX_E], int CSpi[MAX_I], double tm, double tp,
			  int& Nbm, double& Mum[MAX_I][MAX_E], double& Sigm[MAX_I], double& Am[MAX_I],
			  int& Nbp, double& Mup[MAX_I][MAX_E], double& Sigp[MAX_I], double& Ap[MAX_I])
{
int i, j, k, nj, iter;
double ltm, a, x, mnj, minx;

ltm = MathLog(tm);
Nbm = 0;
Nbp = 0;
iter = 1;
while (iter == 1)
  {
  iter = 0;
  for (i=0;i<Nbi;i++)
    {
	 if (CSpi[i] == 1)
	   {
	   minx = 1000000000.0;

// Shrink
      for (j=0;j<Nbm;j++)
	     {
		  a = 0.0;
        for (k=0;k<Nbe;k++)
		    {
		    x = Spi[i][k]-Mum[j][k];
		    a += (x*x);
		    }
		  x = MathSqrt(-a/ltm);
		  if (x < Sigm[j])
		    Sigm[j] = x;
		  if (x < minx)
		    minx = x;
		  }

// Search for a valid prototype
	  nj = -1;
	  mnj = tp;
      for (j=0;j<Nbp;j++)
	    {
		a = 0.0;
        for (k=0;k<Nbe;k++)
		  {
		  x = Spi[i][k]-Mup[j][k];
		  a += (x*x);
		  }
		x = MathExp(-a/(Sigp[j]*Sigp[j]));
		if (x >= mnj)
		  {
		  mnj = x;
		  nj = j;
		  }
	    }

// New prototype
	  if (nj == -1)
	    {
	    Ap[Nbp] = 1.0;
       for (k=0;k<Nbe;k++)
		  Mup[Nbp][k] = Spi[i][k];
		Sigp[Nbp] = minx;  // may be 1e10
		Nbp++;
		iter = 1;
	    }
	    else
	    Ap[nj]++;
	  }

	if (CSpi[i] == -1)
	  {
	  minx = 1000000000.0;

// Shrink
      for (j=0;j<Nbp;j++)
	    {
		a = 0.0;
        for (k=0;k<Nbe;k++)
		  {
		  x = Spi[i][k]-Mup[j][k];
		  a += (x*x);
		  }
		x = MathSqrt(-a/ltm);
		if (x < Sigp[j])
		  Sigp[j] = x;
		if (x < minx)
		  minx = x;
		}

// Search for a valid prototype
	  nj = -1;
	  mnj = tp;
     for (j=0;j<Nbm;j++)
	    {
		 a = 0.0;
       for (k=0;k<Nbe;k++)
		  {
		  x = Spi[i][k]-Mum[j][k];
		  a += (x*x);
		  }
		x = MathExp(-a/(Sigm[j]*Sigm[j]));
		if (x >= mnj)
		  {
		  mnj = x;
		  nj = j;
		  }
	    }

// New prototype
	  if (nj == -1)
	    {
	    Am[Nbm] = 1.0;
        for (k=0;k<Nbe;k++)
		  Mum[Nbm][k] = Spi[i][k];
		Sigm[Nbm] = minx;  // may be 1e10
		Nbm++;
		iter = 1;
	    }
	    else
	    Am[nj]++;
	  }
    }

  if (iter == 1)
	Print("Nm = ", Nbm, " Np = ", Nbp);
  }

x = 0.0;
for (j=0;j<Nbm;j++)
  x += Am[j];
for (j=0;j<Nbm;j++)
  Am[j] /= x;

x = 0.0;
for (j=0;j<Nbp;j++)
  x += Ap[j];
for (j=0;j<Nbp;j++)
  Ap[j] /= x;
}

void DDAclass(int Nbi, int Nbe, double Spi[MAX_I][MAX_E], double tm,
			  int Nbm, double Mum[MAX_I][MAX_E], double Sigm[MAX_I], double Am[MAX_I],
			  int Nbp, double Mup[MAX_I][MAX_E], double Sigp[MAX_I], double Ap[MAX_I],
			  double& om[MAX_I], double& op[MAX_I])
{
int i, j, k;
double a, x;

for (i=0;i<Nbi;i++)
  {
  om[i] = 0.0;
  for (j=0;j<Nbm;j++)
    {
	a = 0.0;
    for (k=0;k<Nbe;k++)
      {
	  x = Spi[i][k]-Mum[j][k];
	  a += x*x;
	  }
	x = MathExp(-a/(Sigm[j]*Sigm[j]));
	om[i] += x*Am[j];
    }

  op[i] = 0.0;
  for (j=0;j<Nbp;j++)
    {
	a = 0.0;
    for (k=0;k<Nbe;k++)
      {
	  x = Spi[i][k]-Mup[j][k];
	  a += x*x;
	  }
	x = MathExp(-a/(Sigp[j]*Sigp[j]));
	op[i] += x*Ap[j];
    }

  x = op[i] + om[i] + tm;
  op[i] /= x;
  om[i] /= x;
  }
}

void DDAStep(int Nbe, double tm, double tp,
			 int Nbi1, double Spi1[MAX_I][MAX_E], int CSpi1[MAX_I],
			 int Nbi2, double Spi2[MAX_I][MAX_E],
			 double& om[MAX_I], double& op[MAX_I])
{
static int Nbm, Nbp;
static double Mum[MAX_I][MAX_E], Sigm[MAX_I], Am[MAX_I];
static double Mup[MAX_I][MAX_E], Sigp[MAX_I], Ap[MAX_I];

DDALearn(Nbi1, Nbe, Spi1, CSpi1, tm, tp, Nbm, Mum, Sigm, Am, Nbp, Mup, Sigp, Ap);
DDAclass(Nbi2, Nbe, Spi2, tm, Nbm, Mum, Sigm, Am, Nbp, Mup, Sigp, Ap, om, op);
}

void TestSpiral()
{
int i;
double ri, spi1[98][2], spi2[98][2];
static double Spi[MAX_I][MAX_E];
static int CSpi[MAX_I];
int Nbi, Nbe;
static double Mum[MAX_I][MAX_E], Sigm[MAX_I], Am[MAX_I];
static double Mup[MAX_I][MAX_E], Sigp[MAX_I], Ap[MAX_I];
static double om[MAX_I], op[MAX_I];
double tm, tp;

Nbi = 98;
Nbe = 2;
for (i=0;i<Nbi;i++)
  {
  ri = ((140.0-i)/140.0);
  ri = ri*ri*ri;
  spi1[i][0] = ri*MathSin(i*M_PI/16.0);
  spi1[i][1] = ri*MathCos(i*M_PI/16.0);
  spi2[i][0] = -ri*MathSin(i*M_PI/16.0);
  spi2[i][1] = -ri*MathCos(i*M_PI/16.0);
  }

for (i=0;i<Nbi;i++)
  {
  ri = ((140.0-i)/140.0);
  ri = ri*ri*ri;
  Spi[2*i][0] = spi1[i][0];
  Spi[2*i][1] = spi1[i][1];
  CSpi[2*i] = 1;
  Spi[2*i+1][0] = spi2[i][0];
  Spi[2*i+1][1] = spi2[i][1];
  CSpi[2*i+1] = -1;
  }

tm = 0.2;
tp = 0.4;
// Test: learning samples and test samples are the same
DDAStep(Nbe, tm, tp, 2*Nbi, Spi, CSpi, 2*Nbi, Spi, om, op);

for (i=0;i<2*Nbi;i++)
  {
  Print(i, " ", CSpi[i], " ", op[i] - om[i]);
  }

}