function w = nwini(xr, m, vr)
%nwini Calculates Nugyen-Widrow initial conditions.
%  adapted from NNet toolbox             8 Axril 1999
%  xr - p-1 by 2 matrix of  [xmin xmax]
%  assumes that the bias is added
%  m - Number of neurons.
%  vr - Active region of the transfer function va = [Vmin Vmax].
%  e.g. vr = [-2 -2] for tansig , [-4 4] for logsig
%  w  is  m by p

r = size(xr,1);  p = r+1 ;

% Null case
if (r == 0) | (m == 0)
  w = zeros(s,p) ;
  return
end

% Remove constant inputs that provide no useful info
R = r;
ind = find(xr(:,1) ~= xr(:,2));
r = length(ind);
xr = xr(ind,:);

% Nguyen-Widrow Method
% Assume inputs and activation potentials range in [-1 1].

% Weights
wMag = 0.7*m^(1/r); % weight vectors magnitude
% weight vectors directions: wDir are row unity vectors
a = 2*rand(m,r)-1 ; 
if r == 1  
  b = ones./abs(a); 
else
   b=sqrt(ones./(sum((a.*a)')))';
end
wDir=b(:,ones(1,r)).*a;
w = wMag*wDir;

% Biases
if (m==1)
  wb = 0;
else
  wb = wMag*[2*(0:m-2)/(m-1)-1 1]'.*sign(w(:,1));
end

% Conversion of activation potentials of [-1 1] to [Nmin Nmax]
a1 = 0.5*(vr(2)-vr(1));
a2 = 0.5*(vr(2)+vr(1));
w = a1*w;
wb = a1*wb+a2;

% Conversion of inputs of xr to [-1 1]
a1 = 2./(xr(:,2)-xr(:,1));
a2 = 1-xr(:,2).*a1;

ap = a1';
wb = w*a2+wb;
w = w.*ap(ones(1,m),:);

% Replace constant inputs
ww = w;
w = zeros(m,R);
w(:,ind) = ww;
%  combine with biasing weights
w = [w wb] ;
%===========================================================