%X,y表示學習樣本
%X是N*(p-1)矩陣，y是1*N矩陣
%x是未知樣本，q*N矩陣
%k表示支持向量所在位值
function [gg,zz,a,k,d]=matlabsvm2(C,KM,o,w,X,y,x)
global KM
global o 
global w
N=size(X,1);
p=size(X,2)+1;

%%%%%%%  產生、加載學習樣本 %%%%%%%
for i=1:N
    yy(:,i)=[X(i,:),y(i)]';
end
qq=yy';
save data.dat qq -ASCII -TABS
%%%%%%%%%  優化拉格朗日乘子 %%%%%%%%%
aa=rand(N,1);        %拉格朗日乘子初始值
%lb=zeros(N,1);
ub=C*ones(N,1);
a=fmincon('myfun',aa,[],[],[],[],[],ub,'mycon');

%%%%%%%% 構造分類函數 %%%%%%%%%
k=argk(a');
q=size(x,1);
for m=1:q
    for n=1:N
        switch KM
            case 'Kmulti'
               u1(n)=(1+sum(X(n,:).*x(m,:)))^o; %Kmulti(X(n,:),x(m,:),o,w);
               u2(n)=%Kmulti(X(n,:),X(k(1),:),o,w);
            case 'Kgauss'
               u1(n)=Kgauss(X(n,:),x(m,:),o,w);
               u2(n)=Kgauss(X(n,:),X(k(1),:),o,w);
            case 'Ktanh'
               u1(n)=Ktanh(X(n,:),x(m,:),o,w);
               u2(n)=Ktanh(X(n,:),X(k(1),:),o,w);
            otherwise 
               disp('核函數輸入錯誤')
        end   
        g(n,m)=a(n)*y(n)*u1(n)-a(n)*y(n)*u2(n)+y(k(1))/N;
    end
end
gg=sum(g);

%%%%%%%% 分類 %%%%%%%%%%
for i=1:q
    if gg(i)>0
       zz(i)=1;
    elseif gg(i)<0
        zz(i)=-1;
    end
end

%%%%%%%%% 求解margin %%%%%%%%%
for i=1:N
    d1(i,:)=a(i)*y(i)*X(i,:);
end
d2=sum(d1);
d=2/sqrt(d2*d2');

%%%%%%%% 子函數 %%%%%%%
function k=argk(z)
N=size(z,2);
m=0;
for i=1:N
    if z(i)>1e-5
        m=m+1;
        k(m)=i;
    end
end
%%%%%%%%%%%
function f=myfun(a,KM,o,w)
global KM
global o 
global w
load data.dat;
N=size(data,1);
m=size(data,2)-1;
for t=1:N
    for n=1:N
        switch KM
         case 'Kmulti'
          u(t,n)=(1+sum(data(t,1:m).*data(n,1:m)))^o
         case 'Kgauss'
          u(t,n)=exp(sum(-(data(t,1:m)-data(n,1:m)).^2)/o^2)
         case 'Ktanh'
          u(t,n)=tanh(o*sum(data(t,1:m).*data(n,1:m))-w)
         otherwise 
          disp('核函數輸入錯誤')
        end
        %u(t,n)=sum(data(t,1:m).*data(n,1:m));
        q(n,t)=a(t,1)*a(n,1)*data(t,m+1)*data(n,m+1)*u(t,n);
    end
end
qq=sum(q);
f=1/2*sum(qq)-sum(a);

%%%%%%%%%%%%%%
function [c,ceq]=mycon(a)
load data.dat;
[N,n]=size(data);
for i=1:N
    s(i)=a(i,1)*data(i,n);
end
c=-a';
ceq=sum(s);

%%%%% 核函數 %%%%%%
function kk=Kmulti(x,y,o,w)
kk=(1+sum(x.*y))^o;
%%%%%%%%%%%%%%%%
function kk=Ktanh(x,y,o,w)
kk=tanh(o*sum(x.*y)-w);
%%%%%%%%%%%%%%%%
function kk=Kgauss(x,y,o,w)
kk=exp(sum(-(x-y).^2)/o^2);