function phi = abskernel(center,data)
% phi = abskernel(center,data)
%
% computes the hidden layer outputs of an RBF network with 
% an absolute value kernel and RBF centers at the points
% given by center.  Each row in phi corresponds to one
% input given by a row in data.
% The centers are specified one per row of the center matrix

dis = repmat(sum(center.*center,2)',size(data,1),1) - ...
      2*data*center' + repmat(sum(data.*data,2),1,size(center,1));
phi = abs(sqrt(dis));


%-------------------------------------------------------------------
function c = allcenters(x,n)
% c = allcenters(x,n)
%
% a function for picking the centers for an RBF network.
% Given the data points x, it ignores the value of n and
% chooses all of the data points as centers
% Each data point is one row of the x matrix.  Each center
% is one row of the c matrix (thus c=x).

c = x;


%-------------------------------------------------------------------
function evalprob6(trainx,trainy,testx,testy,t,centers,n,kernel,varargin)

limits(1,:) = [1 1];
limits(2,:) = [85 50];

[c,w] = rbftrain(kernel,centers,trainx,trainy,n,varargin{:});
err = rbftest(testx,testy,c,w,kernel,varargin{:});
t = sprintf('%s (error = %f)',t,err);
figure;
plotim('rbf',kernel,c,w,varargin{:});
figure;
plotfn(limits,limits,'rbf',1,kernel,c,w,varargin{:});
hold on;
title(t);
hold off;

%-------------------------------------------------------------------
function phi = gkernel(center,data,sigma2)
% phi = gkernel(center,data,sigma2)
%
% computes the outputs for the hidden layer of an RBF network with
% Gaussian RBFs with variance of sigma2.  The centers for each RBF
% are given as rows of the center matrix.  The output matrix phi has
% one row for each input.  The inputs are given as rows in the data
% matrix.

dis = repmat(sum(center.*center,2)',size(data,1),1) - ...
      2*data*center' + repmat(sum(data.*data,2),1,size(center,1));
phi = exp(-dis/(2*sigma2));

%-------------------------------------------------------------------
function c = randcenter(x,n)
% c = randcenter(x,n)
%
% returns n rows of the x matrix chosen at random

p = randperm(size(x,1));
c = x(p(1:n),:);


%-------------------------------------------------------------------
function y = rbf(x,kernel,c,w,varargin)
% y = rbf(x,kernel,c,w,arg1,arg2,...)
%
% returns the values (one per row) of an rbf network evalutated
% on the inputs in x (one per row).  kernel is a string of the
% kernel function to evaluate on the inputs (should take a matrix
% of one input per row and returns the hidden layer values for each
% input on one row).  w is the weight matrix connecting the hidden
% layer to the output layer.  Any optional arguments after w are
% passed on verbatum to the kernel function

phi = feval(kernel,c,x,varargin{:});
y = phi*w;


%-------------------------------------------------------------------
function sqerr = rbftest(x,y,c,w,kernel,varargin)
% sqerr = rbftest(x,y,c,w,kernel,arg1,arg2,...)
%
% returns the sum of squared errors evaluting the
% rbf defined by c (the centers), w (the weights),
% and kernel (a string with the kernel function to
% which the optional parameters are passed) on
% the input x compared to the desired output y

outy = rbf(x,kernel,c,w,varargin{:});
sqerr = sum(sum((y-outy).^2));


%-------------------------------------------------------------------
function [c,w] = rbftrain(kernel,center,x,y,n,varargin)
% [c,w] = rbftrain(kernel,center,x,y,n,arg1,arg2,...)
%
% Using the kernel function given by the string in kernel and
% the center choosing function given by the string in center, 
% computes the weights and centers for an rbf with training input
% x and training output y.  n is passed to the center function to
% indicate the number of centers and the optional remaining arguments
% are all passed to the kernel function

c = feval(center,x,n);
phi = feval(kernel,c,x,varargin{:});
% w = pinv(phi)*y;
w = phi \ y;

function [z,alpha,gammas]=katrain(x,y,eta,etarate,sigma2,maxiter)
% [z,alpha,gammas]=katrain(x,y,eta,etarate,sigma2,maxiter)
% kernel adatron, no bias, soft margin

n=size(x,1);
d=size(x,2);
alpha=zeros(n,1);
K=gkernel(x,x,sigma2);
gammas=zeros(1,maxiter);
sumalphas=zeros(1,maxiter);
go=1;
z=zeros(n,1);
C=5;
for t=1:maxiter
	r=randperm(n);
	for i=1:n
		ri=r(i);
		z(ri)=sum(alpha.*y.*K(:,ri));
		del=eta*(1-z(ri).*y(ri));
		alpha(ri)=min(max(alpha(ri)+del,0),C);
	end
	sumalphas(t)=sum(alpha);
%	subplot(1,2,1);
%	plot(sumalphas(1:t));
%	drawnow;
	gammas(t)=(min(z(y==1&alpha<C))-max(z(y==-1&alpha<C)))/2;
%	subplot(1,2,2);
%	plot(gammas(1:t));
%	drawnow;
	if gammas(t)>0.99
		go=0;
		break;
	end
	eta=eta*etarate;
end

z
t
eta
sumalpha=sum(alpha)
numsv=sum(alpha~=0)
margin=gammas(t)

%-------------------------------------------------------------------
% script for ps4

% adatron
load artificial_1v23.mat;
[z1 alpha1 gammas1]=katrain(X,Y,1,1,1,10000);
svcplot(X,Y,'rbf',alpha,0);
load artificial_2v13.mat;
[z2 alpha2 gammas2]=katrain(X,Y,1,1,1,10000);
svcplot(X,Y,'rbf',alpha,0);
load artificial_3v12.mat;
[z3 alpha3 gammas3]=katrain(X,Y,1,1,1,10000);
svcplot(X,Y,'rbf',alpha,0);

[y,i]=max([z1,z2,z3],[],2);
t=zeros(n,3);
t(i==1,1)=1;
t(i==2,2)=1;
t(i==3,3)=1;

% rbfs
load ps4a;

evalprob6(trainx,trainy,testx,testy,'All centers, Gaussian (var=10)', ...
	'allcenters',0,'gkernel',10);
evalprob6(trainx,trainy,testx,testy,'All centers, absolute value', ...
	'allcenters',0,'abskernel');

evalprob6(trainx,trainy,testx,testy,'50 random pts, Gaussian (var=10)', ...
	'randcenter',50,'gkernel',10);
evalprob6(trainx,trainy,testx,testy,'50 random pts, absolute value', ...
	'randcenter',50,'abskernel');

evalprob6(trainx,trainy,testx,testy,'200 random pts, Gaussian (var=10)', ...
	'randcenter',200,'gkernel',10);
evalprob6(trainx,trainy,testx,testy,'200 random pts, absolute value', ...
	'randcenter',200,'abskernel');

evalprob6(trainx,trainy,testx,testy,'50 pts (k-means), Gaussian (var=10)', ...
	'kmeans',50,'gkernel',10);
evalprob6(trainx,trainy,testx,testy,'50 pts (k-means), absolute value', ...
	'kmeans',50,'abskernel');

evalprob6(trainx,trainy,testx,testy,'200 pts (k-means), Gaussian (var=10)', ...
	'kmeans',200,'gkernel',10);
evalprob6(trainx,trainy,testx,testy,'200 pts (k-means), absolute value', ...
	'kmeans',200,'abskernel');