initwnet.m

新的神经网络算法源程序

function [t,d,w,c,b]=initwnet(initmode,nbwavelon,x,y,min_nbw,levels)
%Initwnet: initialize WaveNet
%
%	[t,d,w,c,b] = initwnet(initmode,nbwavelon,x,y,min_nbw,levels)
%	initmode: initialization mode
%	nbwavelon: desired nb of wavelons for wavenet construction
%	x: input patterns
%	y: output patterns
%	min_nbw: minimum number of input patterns each wavelon should "cover"
%	levels: number of scale levels scanned during initialization
%
%	The arguments min_nbw and levels are optional.
%
%	t: translation parameters
%	d: dilation parameters
%	w: linear weights
%	c: linear comnibation coefficients (direct connections)
%	b: bias

%	By Qinghua Zhang. March, 1994.

if nargin < 4
  error('Wrong number of input arguments.');
end

if nargout < 5
  error('Wrong number of output arguments.');
end

[xl,xc] = size(x);
[yl,yc] = size(y);

if xc ~= yc
  error('x and y sould have the same number of columns.');
end

if yl ~= 1
  error('y should be a row vector.');
end

nbvar = xl;
nbobs = xc;

disp('   Running the initialization procedure...');

base = 2^nbvar;

% ---If arguments not defined then make empty---
if nargin < 5, min_nbw = []; end;
if nargin < 6, levels = []; end;

% determine the over nb of scale levels (nblayer)
if ~isempty(levels)
  nblayer = levels;
else
  nblayer = 1;
  nboverw = 1;
  while nboverw < nbwavelon
    nboverw = nboverw + base^nblayer;
    nblayer=nblayer+1;
  end
  nboverw = nboverw + base^nblayer;
  nblayer=nblayer+1;

  if nblayer<4, nblayer=4; end %%%%%%%
end

fprintf('\n   %g scale levels to be scanned.\n\n', nblayer);


% Default value for MINIMUM NUMBER of obs within the "support" of each wavelet
%======================================================
if isempty(min_nbw)
  %min_nbw = max(ceil(0.2 * nbobs / nbwavelon), 2+nbvar);
  min_nbw = 2+nbvar;
end 

fprintf('   Minmum number of obs points "coverd" by each wavelon = %g.\n\n',... 
min_nbw);


% Determine the wavelet candidates which "cover" at least min_nbw obs.
% The scanning is performed for each scale level and oriented by the obs
%========================================================================

d = 1;			% The first wavelon (layer=0), d=1, t=0, always taken
t = zeros(xl, 1);

for layer=1:(nblayer-1)
  trans = ceil(abs(x * 2^(layer-1))) .* sign(x); % Translation numbers

  dil = 2^layer;	% dilation parameter
  
  while 1
    [trans_l trans_c] = size(trans);

    while prod(trans(:,1)) == 0		% If translation_number(1) includes 0
      if trans_c == 1			% and there are more than one cadidates
        break;
      else
        trans = trans(:, 2:trans_c);	% then remove translation_number(1).
        trans_c = trans_c-1;
      end
    end

    if prod(trans(:,1)) == 0 & trans_c == 1 % If trans_number(1) includes 0
      break;				    % and only one candidate,
    end;  				    % then end the loop.
    
    if trans_l == 1			%Find columns different from the 1st
      ind = find(trans(1) ~= trans);	%for case nbvar=1
    else
      ind = find(~ prod(trans(:,1) * ones(1, trans_c) == trans)); %for nbvar>1
    end

    if trans_c - length(ind) >= min_nbw	   %If min_nbw or more candidates = 1st
      d = [d dil];
      t = [t (trans(:,1)*2-sign(trans(:,1)))/dil];
    end

    if isempty(ind), break; end;	%If all left candidates are id, finish.
    
    trans = trans(:,ind);	% Keep cadidates diff from the 1st.
  end
  
end

d = ones(nbvar, 1) * d;		% for dimension=nbvar
[nbvar, nbw1] = size(t);	% nbw1 = nb of wavelon candidates.

fprintf('   Number of wavelon candidates = %g.\n\n', nbw1);

if nbw1 < nbwavelon
  nbwavelon = nbw1;
  fprintf('   Number of wavelons is reduced to %g\n\n', nbw1);
elseif nbwavelon > 0
  fprintf('   %g wavelons are used to construct the wavelet net.\n\n', ...
          nbwavelon); 
end

pmat = [ones(xc,1) x' wavelon(x, t, d)'];

%normalize the vectors
pm2 = pmat .* pmat;	 
if xc ~= 1
  pm2 = sqrt(sum(pm2));
end
pmat = pmat ./ (ones(xc,1)*pm2);


%Call one of the regressor selection algorithms
if initmode==1
  [indsld, wsld] = initmod1(pmat, nbwavelon, y, nbvar);
elseif initmode==2
  [indsld, wsld] = initmod2(pmat, nbwavelon, y, nbvar);
elseif initmode==3
  [indsld, wsld] = initmod3(pmat, nbwavelon, y, nbvar);
end

%Remove the indexes correspoding to the lienar terms
indsld = indsld-(nbvar+1);
ind = find(indsld>0);
indsld = indsld(ind);

t=t(:,indsld); d=d(:,indsld);

lenrecind = length(indsld);
w = y / [wavelon(x,t,d); x; ones(1,nbobs)];
c = w(lenrecind+1:lenrecind+nbvar);
b = w(lenrecind+nbvar+1);
w = w(1:lenrecind);