function [ c, sigma , W_output ] = SOFNN( X, d, Kd )%SOFNN Self-Organizing Fuzzy Neural Networks%Input Parameters% X(r,n) - rth traning data from nth observation% d(n) - the desired output of the network (must be a row vector) % Kd(r) - predefined distance threshold for the rth input%Output Parameters% c(IndexInputVariable,IndexNeuron)% sigma(IndexInputVariable,IndexNeuron)% W_output is a vector%Setting up Parameters for SOFNNSigmaZero=4;delta=0.12;threshold=0.1354;k_sigma=1.12;%For more accurate results uncomment the following%format long;%Implementation of a SOFNN model[size_R,size_N]=size(X);%size_R - the number of input variablesc=[];sigma=[];W_output=[];u=0; % the number of neurons in the structureQ=[];O=[];Psi=[];for n=1:size_Nx=X(:,n);if u==0 % No neuron in the structure?c=x;sigma=SigmaZero*ones(size_R,1);u=1;Psi=GetMePsi(X,c,sigma);[Q,O] = UpdateStructure(X,Psi,d);pT_n=GetMeGreatPsi(x,Psi(n,:))';else[Q,O,pT_n] = UpdateStructureRecursively(X,Psi,Q,O,d,n);end;KeepSpinning=true;while KeepSpinning%Calculate the error and if-part criteriaae=abs(d(n)-pT_n*O); %approximation error[phi,~]=GetMePhi(x,c,sigma);[maxphi,maxindex]=max(phi); % maxindex refers to the neuron's index if ae>deltaif maxphi<threshold%enlarge width[minsigma,minindex]=min(sigma(:,maxindex));sigma(minindex,maxindex)=k_sigma*minsigma;Psi=GetMePsi(X,c,sigma);[Q,O] = UpdateStructure(X,Psi,d);pT_n=GetMeGreatPsi(x,Psi(n,:))';else%Add a new neuron and update structurectemp=[];sigmatemp=[];dist=0;for r=1:size_Rdist=abs(x(r)-c(r,1));distIndex=1;for j=2:uif abs(x(r)-c(r,j))<distdistIndex=j;dist=abs(x(r)-c(r,j));end;end;if dist<=Kd(r)ctemp=[ctemp; c(r,distIndex)];sigmatemp=[sigmatemp ; sigma(r,distIndex)];elsectemp=[ctemp; x(r)];sigmatemp=[sigmatemp ; dist];end;end;c=[c ctemp];sigma=[sigma sigmatemp];Psi=GetMePsi(X,c,sigma);[Q,O] = UpdateStructure(X,Psi,d);KeepSpinning=false;u=u+1;end;elseif maxphi<threshold%enlarge width[minsigma,minindex]=min(sigma(:,maxindex));sigma(minindex,maxindex)=k_sigma*minsigma;Psi=GetMePsi(X,c,sigma);[Q,O] = UpdateStructure(X,Psi,d);pT_n=GetMeGreatPsi(x,Psi(n,:))';else%Do nothing and exit the whileKeepSpinning=false;end;end;end;end;W_output=O;endfunction [Q_next, O_next,pT_n] = UpdateStructureRecursively(X,Psi,Q,O,d,n) %O=O(t-1) O_next=O(t)p_n=GetMeGreatPsi(X(:,n),Psi(n,:));pT_n=p_n';ee=abs(d(n)-pT_n*O); %|e(t)|temp=1+pT_n*Q*p_n;ae=abs(ee/temp);if ee>=aeL=Q*p_n*(temp)^(-1);Q_next=(eye(length(Q))-L*pT_n)*Q;O_next=O + L*ee;elseQ_next=eye(length(Q))*Q;O_next=O;end;endfunction [ Q , O ] = UpdateStructure(X,Psi,d)GreatPsiBig = GetMeGreatPsi(X,Psi);%M=u*(r+1)%n - the number of observations[M,~]=size(GreatPsiBig);%Others Ways of getting Q=[P^T(t)*P(t)]^-1%********************************************************************** ****%opts.SYM = true;%Q = linsolve(GreatPsiBig*GreatPsiBig',eye(M),opts);%%Q = inv(GreatPsiBig*GreatPsiBig');%Q = pinv(GreatPsiBig*GreatPsiBig');%********************************************************************** ****Y=GreatPsiBig\eye(M);Q=GreatPsiBig'\Y;O=Q*GreatPsiBig*d';end%This function works too with x% (X=X and Psi is a Matrix) - Gets you the whole GreatPsi% (X=x and Psi is the row related to x) - Gets you just the column related with the observationfunction [GreatPsi] = GetMeGreatPsi(X,Psi)%Psi - In a row you go through the neurons and in a column you go through number of %observations **** Psi(#obs,IndexNeuron) ****GreatPsi=[];[N,U]=size(Psi);for n=1:Nx=X(:,n);GreatPsiCol=[];for u=1:UGreatPsiCol=[ GreatPsiCol ; Psi(n,u)*[1; x] ];end;GreatPsi=[GreatPsi GreatPsiCol];end;endfunction [phi, SumPhi]=GetMePhi(x,c,sigma)[r,u]=size(c);%u - the number of neurons in the structure%r - the number of input variablesphi=[];SumPhi=0;for j=1:u % moving through the neuronsS=0;for i=1:r % moving through the input variablesS = S + ((x(i) - c(i,j))^2) / (2*sigma(i,j)^2);end;phi = [phi exp(-S)];SumPhi = SumPhi + phi(j); %phi(u)=exp(-S)end;end%This function works too with x, it will give you the row related to xfunction [Psi] = GetMePsi(X,c,sigma)[~,u]=size(c);[~,size_N]=size(X);%u - the number of neurons in the structure%size_N - the number of observationsPsi=[];for n=1:size_N[phi, SumPhi]=GetMePhi(X(:,n),c,sigma);PsiTemp=[];for j=1:u%PsiTemp is a row vector ex: [1 2 3]PsiTemp(j)=phi(j)/SumPhi;end;Psi=[Psi; PsiTemp];%Psi - In a row you go through the neurons and in a column you go through number of %observations **** Psi(#obs,IndexNeuron) ****end;end。