function [f,dd] = fgrade(x,v,p,P)
% F = FGRADE(X,V,m,P) calculates membership grade of elements in vector
% X into clusters (multi-dimensional fuzzy sets) defined by means V
% and distance norm-inducing matrix M. Both V and M are obtained from
% GKFAST. For details on meaning of V and M and for the
% algorithm, see GKFAST and the description of the Gustafson-Kessel
% method.
%
%     Input:
%       X  ... data matrix (M x N, M - number of data, N - data dimension)
%       V  ... cluster means matrix (C x N, C - number of clusters)
%       m  ... optional parameter m > 1, determines the fuzziness
%              of clusters, must be the same as in GKFAST, default 2
%       P  ... fuzzy covariance matrices
%
%     Output:
%       F  ... matrix of membership grades

% (c) Robert Babuska, 1994-96

if (nargin < 3), p = 2; elseif ~all(size(p)), p = 2; end;
if nargin < 4, P = []; end;

[mx,nx] = size(x);
[c,nc] = size(v);
nP = size(P,2);
x1 = ones(mx,1);
rho = ones(1,c);

  for j = 1 : c,                        % for all clusters
    xv = x - x1*v(j,1:nx);
if nx > 1,
% Calculate distances
    if isempty(P), 
      m = eye(nx); 
    else
%      m = (det(P((j-1)*nP+(1:nP),:))/rho(j)).^(1/nP)*inv(P((j-1)*nP+(1:nP),:));
      P1 = P((j-1)*nP+(1:nP),:);
    P1 = P1(1:nx,1:nx);
      m = (det(P1)/rho(j)).^(1/nx)*inv(P1);
%      m = M((j-1)*nx+(1:nx),:);
    end;
    d(:,j) = sum((xv*m(1:nx,1:nx).*xv)')';
else
  d(:,j) = 1./P(j)*xv.^2;
end;
end;
% Calculate "similarities" based on distance only
  dd = 1./(1+d);
%  dd = d;
% Calculate membership grades
  d0 = (d == 0);                                % find zero distances
  c0 = max(d0')';
  d(find(c0~=0),:) = zeros(sum(c0),c);                   % set d(i,:) to 0
  d(find(d0~=0)) = ones(sum(sum(d0)),1);                 % set d(d0)'s to 1
  d = d.^(-1/(p-1));
  f = (d ./ (sum(d')'*ones(1,c)));

