clear;clg;hold off;

eta=0.01;
a=0.0005;
b=0.0005;

% initial weights
%w=3/sqrt(3)*2*( rand([9 3])-0.5 ); % -3/sqrt(ki) ~ +3/sqrt(ki)
w=zeros([9 3]);
w01=w(1,1); w11=w(1,2); w21=w(1,3);
w02=w(2,1); w12=w(2,2); w22=w(2,3);
w03=w(3,1); w13=w(3,2); w23=w(3,3);
w04=w(4,1); w14=w(4,2); w24=w(4,3);
w05=w(5,1); w15=w(5,2); w25=w(5,3);
w06=w(6,1); w16=w(6,2); w26=w(6,3);
w07=w(7,1); w17=w(7,2); w27=w(7,3);
w08=w(8,1); w18=w(8,2); w28=w(8,3);
w09=w(9,1); w19=w(9,2); w29=w(9,3);

% training data
for t=1:1:11
for s=1:1:11
i=(t-1)*11+s;
X1(i)=(t-6);
X2(i)=(s-6);
D(i)=1+1*X1(i)+1*X2(i);
end
end

epoch=0;
pre_E=0;

while 1
E=0;
series=randperm(121);

for i=1:1:121
x1=X1(series(i));
x2=X2(series(i));
d=D(series(i));

% Membership
if(x1>=-5 & x1<0)
A1N=-x1/5; %x1:N
A1Z=(x1+5)/5; %x1:Z
A1P=0; %x1:P
elseif(x1>=0 & x1<=5)
A1N=0;
A1Z=-(x1-5)/5;
A1P=x1/5;
end

if(x2>=-5 & x2<0)
A2N=-x2/5; %x2:N
A2Z=(x2+5)/5; %x2:Z
A2P=0; %x2:P
elseif(x2>=0 & x2<=5)
A2N=0;
A2Z=-(x2-5)/5;
A2P=x2/5;
end

mA1=min(A1N,A2N); %Rule1: x1:N x2:N
mA2=min(A1N,A2Z); %Rule2: x1:N x2:Z
mA3=min(A1N,A2P); %Rule3: x1:N x2:P
mA4=min(A1Z,A2N); %Rule4: x1:Z x2:N
mA5=min(A1Z,A2Z); %Rule5: x1:Z x2:Z
mA6=min(A1Z,A2P); %Rule6: x1:Z x2:P
mA7=min(A1P,A2N); %Rule7: x1:P x2:N
mA8=min(A1P,A2Z); %Rule8: x1:P x2:Z
mA9=min(A1P,A2P); %Rule9: x1:P x2:P

y1=w01+w11*x1+w21*x2;
y2=w02+w12*x1+w22*x2;
y3=w03+w13*x1+w23*x2;
y4=w04+w14*x1+w24*x2;
y5=w05+w15*x1+w25*x2;
y6=w06+w16*x1+w26*x2;
y7=w07+w17*x1+w27*x2;
y8=w08+w18*x1+w28*x2;
y9=w09+w19*x1+w29*x2;

A=mA1+mA2+mA3+mA4+mA5+mA6+mA7+mA8+mA9;
y=(mA1*y1+mA2*y2+mA3*y3+mA4*y4+mA5*y5+mA6*y6+mA7*y7+mA8*y8+mA9*y9)/A;

e=d-y;

w01=w01+eta*e*mA1/A;
w02=w02+eta*e*mA2/A;
w03=w03+eta*e*mA3/A;
w04=w04+eta*e*mA4/A;
w05=w05+eta*e*mA5/A;
w06=w06+eta*e*mA6/A;
w07=w07+eta*e*mA7/A;
w08=w08+eta*e*mA8/A;
w09=w09+eta*e*mA9/A;

w11=w11+eta*e*mA1*x1/A;
w12=w12+eta*e*mA2*x1/A;
w13=w13+eta*e*mA3*x1/A;
w14=w14+eta*e*mA4*x1/A;
w15=w15+eta*e*mA5*x1/A;
w16=w16+eta*e*mA6*x1/A;
w17=w17+eta*e*mA7*x1/A;
w18=w18+eta*e*mA8*x1/A;
w19=w19+eta*e*mA9*x1/A;

w21=w21+eta*e*mA1*x2/A;
w22=w22+eta*e*mA2*x2/A;
w23=w23+eta*e*mA3*x2/A;
w24=w24+eta*e*mA4*x2/A;
w25=w25+eta*e*mA5*x2/A;
w26=w26+eta*e*mA6*x2/A;
w27=w27+eta*e*mA7*x2/A;
w28=w28+eta*e*mA8*x2/A;
w29=w29+eta*e*mA9*x2/A;

E=E+1/2*e^2;
end

if (E-pre_E)<0
eta=eta+a;
elseif (E-pre_E)>0
eta=eta-b*eta;
end

pre_E=E;

epoch=epoch+1
E
if E < 0.1
break;
end
end

fprintf('************************************************\n');
fprintf('Error : %6.4f\n',E);
fprintf('Rule1 : %6.4f + %6.4f x1 + %6.4f x2\n',w01,w11,w21);
fprintf('Rule2 : %6.4f + %6.4f x1 + %6.4f x2\n',w02,w12,w22);
fprintf('Rule3 : %6.4f + %6.4f x1 + %6.4f x2\n',w03,w13,w23);
fprintf('Rule4 : %6.4f + %6.4f x1 + %6.4f x2\n',w04,w14,w24);
fprintf('Rule5 : %6.4f + %6.4f x1 + %6.4f x2\n',w05,w15,w25);
fprintf('Rule6 : %6.4f + %6.4f x1 + %6.4f x2\n',w06,w16,w26);
fprintf('Rule7 : %6.4f + %6.4f x1 + %6.4f x2\n',w07,w17,w27);
fprintf('Rule8 : %6.4f + %6.4f x1 + %6.4f x2\n',w08,w18,w28);
fprintf('Rule9 : %6.4f + %6.4f x1 + %6.4f x2\n',w09,w19,w29);
fprintf('************************************************\n');

% Check
for t=1:.5:11
for s=1:.5:11
i=2*( (t-1)*21+s );
X1(i)=(t-6);
X2(i)=(s-6);
D(i)=1+1*X1(i)+1*X2(i);
end
end

for i=1:1:(21*21+1)
x1=X1(i);
x2=X2(i);

% Membershiw
if(x1>=-5 & x1<0)
A1N=-x1/5; %x1:N
A1Z=(x1+5)/5; %x1:Z
A1P=0; %x1:P
elseif(x1>=0 & x1<=5)
A1N=0;
A1Z=-(x1-5)/5;
A1P=x1/5;
end

if(x2>=-5 & x2<0)
A2N=-x2/5; %x2:N
A2Z=(x2+5)/5; %x2:Z
A2P=0; %x2:P
elseif(x2>=0 & x2<=5)
A2N=0;
A2Z=-(x2-5)/5;
A2P=x2/5;
end

mA1=min(A1N,A2N); %Rule1: x1:N x2:N
mA2=min(A1N,A2Z); %Rule2: x1:N x2:Z
mA3=min(A1N,A2P); %Rule3: x1:N x2:P
mA4=min(A1Z,A2N); %Rule4: x1:Z x2:N
mA5=min(A1Z,A2Z); %Rule5: x1:Z x2:Z
mA6=min(A1Z,A2P); %Rule6: x1:Z x2:P
mA7=min(A1P,A2N); %Rule7: x1:P x2:N
mA8=min(A1P,A2Z); %Rule8: x1:P x2:Z
mA9=min(A1P,A2P); %Rule9: x1:P x2:P

y1=w01+w11*x1+w21*x2;
y2=w02+w12*x1+w22*x2;
y3=w03+w13*x1+w23*x2;
y4=w04+w14*x1+w24*x2;
y5=w05+w15*x1+w25*x2;
y6=w06+w16*x1+w26*x2;
y7=w07+w17*x1+w27*x2;
y8=w08+w18*x1+w28*x2;
y9=w09+w19*x1+w29*x2;

A=mA1+mA2+mA3+mA4+mA5+mA6+mA7+mA8+mA9;
y(i)=(mA1*y1+mA2*y2+mA3*y3+mA4*y4+mA5*y5+mA6*y6+mA7*y7+mA8*y8+mA9*y9)/A;
end

plot3(X1,X2,D,'o',X1,X2,y,'x');
xlabel('x1');
ylabel('x2');
zlabel('o:d x:y');
grid;