/*

    Copyright (C) 2002 Stefan Westerfeld
                       stefan@space.twc.de

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

    */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>

struct T
{
	float x[2];
	float y;
	float alpha;

	T(float x1 = 0, float x2 = 0, float y = 0)
	{
		x[0] = x1;
		x[1] = x2;
		this->y = y;
		alpha = 0;
	}
	
	T(const T& other)
	{
		x[0] = other.x[0];
		x[1] = other.x[1];
		y    = other.y;
	}
};

double frand()
{
  return rand() / double(RAND_MAX);
}

float sqr (float x)
{
	return x * x;
}

float dotProduct (const T& i, const T& j)
{
	/*
	 * change this to get different dot products
	 */
	enum { rbf, poly, linear } style = rbf;

	if (style == rbf)
	{
		float delta_2 = (i.x[0] - j.x[0]) * (i.x[0] - j.x[0])
					  + (i.x[1] - j.x[1]) * (i.x[1] - j.x[1]);

		return exp(-delta_2 / 1);
	}
	else if (style == linear)
	{
		return (i.x[0] * j.x[0]) + (i.x[1] * j.x[1]);
	}
	else if (style == poly)
	{
		return sqr((i.x[0] * j.x[0]) + (i.x[1] * j.x[1]) + 1);
	}
	else
	{
		assert (0);
	}
}

float score (const vector<T>& data)
{
	float result = 0.0;
	vector<T>::const_iterator i,j;

	for (i = data.begin(); i != data.end(); i++)
		result += i->alpha;

	for (i = data.begin(); i != data.end(); i++)
	{
		for (j = data.begin(); j != data.end(); j++)
		{
			result -= 0.5 * (i->alpha * j->alpha * i->y * j->y * dotProduct(*i, *j));
		}
	}
	return result;
}

float dscore (const vector<T>& data, const vector<float>& delta, float k)
{
	float result = 0.0;
	int i,j,n = data.size();

	for (i = 0; i < n; i++)
		result += data[i].alpha + delta[i] * k;

	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			result -= 0.5 * ((data[i].alpha + delta[i] * k) * (data[j].alpha + delta[j] * k)
					        * data[i].y * data[j].y * dotProduct(data[i], data[j]));
		}
	}
	return result;
}

/* derive dscore in k, set to zero to attain best k */
float bestk (const vector<T>& data, const vector<float>& delta)
{
	float k = 0.0, u = 0, v = 0;
	int i,j,n = data.size();

	for (i = 0; i < n; i++)
		k += delta[i];

	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			k -= 0.5 * ((data[i].alpha * delta[j] + data[j].alpha * delta[i])
					 * data[i].y * data[j].y * dotProduct(data[i], data[j]));
					
			v += ((delta[i] * delta[j]) * data[i].y * data[j].y * dotProduct(data[i], data[j]));
		}
	}
	k /= v;

	/* bound k: alpha_i should be positive */
	for (i = 0; i < n; i++)
	{
		if (data[i].alpha + k * delta[i] < 0)
			k = -data[i].alpha / delta[i];
	}

	return k;
}

/* classifier */
float classify (const T& t, const vector<T>& data, float b)
{
	float c = b;

	for (vector<T>::const_iterator i = data.begin(); i != data.end(); i++)
		c += dotProduct(t, *i) * i->y * i->alpha;

	return c;
}

int main(int argc, char **argv)
{
	vector<T> data;
	bool trace = false;

	srand (time(0));

	data.push_back(T(1.0, 1.0, -1));
	data.push_back(T(1.1, 1.6, -1));
	data.push_back(T(1.3, 1.4, -1));
	data.push_back(T(1.7, 0.8, -1));
	data.push_back(T(1.9, 1.5, -1));
	data.push_back(T(1.0,-1.1, 1));

	data.push_back(T(0.1, 0.5,  -1)); // X
	data.push_back(T(2.1, 3.5,  1));
	data.push_back(T(2.3, 4.0,  1));
	data.push_back(T(2.5, 3.2,  1));
	data.push_back(T(2.6, 3.8,  1));
	data.push_back(T(2.9, 3.3,  1));

	/*
	 * svm data
	 *
	 * outputs the data
	 */
	if ((argc == 2) && (strcmp(argv[1], "data") == 0))
	{
		for (vector<T>::iterator di = data.begin(); di != data.end(); di++)
		{
			printf("%f %f %f\n", di->x[0], di->x[1], di->y);
		}
	}

	/*
	 * svm trace
	 *
	 * outputs trace of computation process
	 */
	if ((argc == 2) && (strcmp(argv[1], "trace") == 0))
	{
		trace = true;
	}


	vector<float> delta(data.size());
	float delta_sum;

	for (int n = 0; n < 10000; n++)
	{
#if 0
		/*
		 * doesn't converge too fast as the normalization for the sum constraint
		 * usually leaves us with a direction that isn't allowed by the
		 * alpha_i >= 0 lagrange multiplicator constraints
		 */

		/* compute random gradient descend vector */
		for(int i = 0; i < delta.size(); i++)
			delta[i] = frand();

		/* normalize so that the sum_i(alpha_i*y_i) constraint remains satisfied */
		delta_sum = 0.0;
		for(int i = 0; i < delta.size(); i++)
			delta_sum += delta[i] * data[i].y;
		for(int i = 0; i < delta.size(); i++)
			delta[i] -= delta_sum / delta.size() * data[i].y;
#endif
		/*
		 * compute random gradient descend vector in two dimensions
		 *
		 * question: could this possibly break correctness?
		 */
		for(int i = 0; i < delta.size(); i++)
			delta[i] = 0;

		int dim1 = 0, dim2 = 0;
		while (dim1 == dim2)
		{
			dim1 = rand() % delta.size();
			dim2 = rand() % delta.size();
		}
		delta[dim1] = 1.0;
		delta[dim2] = -data[dim1].y * data[dim2].y;

#if 0
		/* debugging: check that adding delta doesn't loose correctness  of (15) */
		delta_sum = 0.0;
		for(int i = 0; i < delta.size(); i++)
		{
			delta_sum += delta[i] * data[i].y;
			printf("delta[%d] = %f\n", i, delta[i]);
		}
		printf("%f\n", delta_sum);
#endif

#if 0
		/* debugging: plot effects of different k values */
		for(float k = -1.0; k < 1.0; k += 0.02)
			printf("%f %f\n", k, dscore(data,delta,k));
#endif
		float k =  bestk (data, delta);
		for(int i = 0; i < delta.size(); i++)
			data[i].alpha += k * delta[i];

		if (trace)
		{
			printf ("bestk = %f\n",k);
			printf("%f\n", score(data));

			printf("alpha = < ");
			for(int i = 0; i < delta.size(); i++)
				printf("%f ", data[i].alpha);
			printf(">\n");
		}
	}

	/* calculate the optimal bias b */
	float c_neg = -1000, c_pos = +1000, b = 0.0;
	for(int i = 0; i < delta.size(); i++)
	{
		float c = classify (data[i], data, b);
		if (data[i].y < 0) c_neg = max (c_neg, c);
		if (data[i].y > 0) c_pos = min (c_pos, c);
	}
	b =  -(c_neg + c_pos) / 2.0;

	/*
	 * svm
	 *
	 * print classification results
	 */
	if (argc == 1)
	{
#if 0
		printf("w = <%f %f>\n", w0, w1);
		printf("b = %f\n", b);
#endif
		printf("# x0 x1 y c alpha\n");
		for(int i = 0; i < delta.size(); i++)
			printf("%f %f %f %f %f\n", data[i].x[0], data[i].x[1], data[i].y,
					                   classify (data[i], data, b), data[i].alpha);
	}

	/*
	 * svm test
	 *
	 * classifies a few random points
	 */
	if (argc == 2 && strcmp(argv[1], "test") == 0)
	{
		for (int i = 0; i < 10000; i++)
		{
			T test;
			test.x[0] = frand() * 10.0 - 5.0;
			test.x[1] = frand() * 10.0 - 5.0;
			printf("%f %f %f\n", test.x[0], test.x[1], classify (test, data, b));
		}
	}
	/*
	 * svm bound
	 *
	 * draws the decision boundary between the classes and the support vectors
	 */
	if (argc == 2 && strcmp(argv[1], "bound") == 0)
	{
		for (int i = 0; i < 10000; i++)
		{
			T t1, t2;
			t1.x[0] = frand() * 20.0 - 10.0;
			t1.x[1] = frand() * 20.0 - 10.0;
			t2.x[0] = frand() * 20.0 - 10.0;
			t2.x[1] = frand() * 20.0 - 10.0;
			if ((classify (t1, data, b) * classify(t2, data, b)) < 0)
			{
				T m;
				for (int i = 0; i < 32; i++)
				{
					m.x[0] = (t1.x[0]+t2.x[0]) / 2.0;
					m.x[1] = (t1.x[1]+t2.x[1]) / 2.0;
				
					if (classify (m, data, b) * classify (t1, data, b) < 0)
						t2 = m;
					else
						t1 = m;
				}
				printf("%f %f\n", m.x[0], m.x[1]);
			}
		}
		for(int i = 0; i < delta.size(); i++)
			if (data[i].alpha > 0.0001)
			printf("%f %f\n", data[i].x[0], data[i].x[1]);
	}
}