wolffd@0: function net = gtminit(net, options, data, samp_type, varargin)
wolffd@0: %GTMINIT Initialise the weights and latent sample in a GTM.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	NET = GTMINIT(NET, OPTIONS, DATA, SAMPTYPE) takes a GTM NET and
wolffd@0: %	generates a sample of latent data points and sets the centres (and
wolffd@0: %	widths if appropriate) of NET.RBFNET.
wolffd@0: %
wolffd@0: %	If the SAMPTYPE is 'REGULAR', then regular grids of latent data
wolffd@0: %	points and RBF centres are created.  The dimension of the latent data
wolffd@0: %	space must be 1 or 2.  For one-dimensional latent space, the
wolffd@0: %	LSAMPSIZE parameter gives the number of latent points and the
wolffd@0: %	RBFSAMPSIZE parameter gives the number of RBF centres.  For a two-
wolffd@0: %	dimensional latent space, these parameters must be vectors of length
wolffd@0: %	2 with the number of points in each of the x and y directions to
wolffd@0: %	create a rectangular grid.  The widths of the RBF basis functions are
wolffd@0: %	set by a call to RBFSETFW passing OPTIONS(7) as the scaling
wolffd@0: %	parameter.
wolffd@0: %
wolffd@0: %	If the SAMPTYPE is 'UNIFORM' or 'GAUSSIAN' then the latent data is
wolffd@0: %	found by sampling from a uniform or Gaussian distribution
wolffd@0: %	correspondingly.  The RBF basis function parameters are set by a call
wolffd@0: %	to RBFSETBF with the DATA parameter as dataset and the OPTIONS
wolffd@0: %	vector.
wolffd@0: %
wolffd@0: %	Finally, the output layer weights of the RBF are initialised by
wolffd@0: %	mapping the mean of the latent variable to the mean of the target
wolffd@0: %	variable, and the L-dimensional latent variale variance to the
wolffd@0: %	variance of the targets along the first L principal components.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GTM, GTMEM, PCA, RBFSETBF, RBFSETFW
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check for consistency
wolffd@0: errstring = consist(net, 'gtm', data);
wolffd@0: if ~isempty(errstring)
wolffd@0:   error(errstring);
wolffd@0: end
wolffd@0: 
wolffd@0: % Check type of sample
wolffd@0: stypes = {'regular', 'uniform', 'gaussian'};
wolffd@0: if (strcmp(samp_type, stypes)) == 0
wolffd@0:   error('Undefined sample type.')
wolffd@0: end
wolffd@0: 
wolffd@0: if net.dim_latent > size(data, 2)
wolffd@0:   error('Latent space dimension must not be greater than data dimension')
wolffd@0: end
wolffd@0: nlatent = net.gmmnet.ncentres;
wolffd@0: nhidden = net.rbfnet.nhidden;
wolffd@0: 
wolffd@0: % Create latent data sample and set RBF centres
wolffd@0: 
wolffd@0: switch samp_type
wolffd@0: case 'regular'
wolffd@0:    if nargin ~= 6
wolffd@0:       error('Regular type must specify latent and RBF shapes');
wolffd@0:    end
wolffd@0:    l_samp_size = varargin{1};
wolffd@0:    rbf_samp_size = varargin{2};
wolffd@0:    if round(l_samp_size) ~= l_samp_size
wolffd@0:       error('Latent sample specification must contain integers')
wolffd@0:    end
wolffd@0:    % Check existence and size of rbf specification
wolffd@0:    if any(size(rbf_samp_size) ~= [1 net.dim_latent]) | ...
wolffd@0:          prod(rbf_samp_size) ~= nhidden
wolffd@0:       error('Incorrect specification of RBF centres')
wolffd@0:    end
wolffd@0:    % Check dimension and type of latent data specification
wolffd@0:    if any(size(l_samp_size) ~= [1 net.dim_latent]) | ...
wolffd@0:          prod(l_samp_size) ~= nlatent
wolffd@0:       error('Incorrect dimension of latent sample spec.')
wolffd@0:    end
wolffd@0:    if net.dim_latent == 1
wolffd@0:       net.X = [-1:2/(l_samp_size-1):1]';
wolffd@0:       net.rbfnet.c = [-1:2/(rbf_samp_size-1):1]';
wolffd@0:       net.rbfnet = rbfsetfw(net.rbfnet, options(7));
wolffd@0:    elseif net.dim_latent == 2
wolffd@0:       net.X = gtm_rctg(l_samp_size);
wolffd@0:       net.rbfnet.c = gtm_rctg(rbf_samp_size);
wolffd@0:       net.rbfnet = rbfsetfw(net.rbfnet, options(7));
wolffd@0:    else
wolffd@0:       error('For regular sample, input dimension must be 1 or 2.')
wolffd@0:    end
wolffd@0:    
wolffd@0:    
wolffd@0: case {'uniform', 'gaussian'}
wolffd@0:    if strcmp(samp_type, 'uniform')
wolffd@0:       net.X = 2 * (rand(nlatent, net.dim_latent) - 0.5);
wolffd@0:    else
wolffd@0:       % Sample from N(0, 0.25) distribution to ensure most latent 
wolffd@0:       % data is inside square
wolffd@0:       net.X = randn(nlatent, net.dim_latent)/2;
wolffd@0:    end   
wolffd@0:    net.rbfnet = rbfsetbf(net.rbfnet, options, net.X);
wolffd@0: otherwise
wolffd@0:    % Shouldn't get here
wolffd@0:    error('Invalid sample type');
wolffd@0:    
wolffd@0: end
wolffd@0: 
wolffd@0: % Latent data sample and basis function parameters chosen.
wolffd@0: % Now set output weights
wolffd@0: [PCcoeff, PCvec] = pca(data);
wolffd@0: 
wolffd@0: % Scale PCs by eigenvalues
wolffd@0: A = PCvec(:, 1:net.dim_latent)*diag(sqrt(PCcoeff(1:net.dim_latent)));
wolffd@0: 
wolffd@0: [temp, Phi] = rbffwd(net.rbfnet, net.X);
wolffd@0: % Normalise X to ensure 1:1 mapping of variances and calculate weights
wolffd@0: % as solution of Phi*W = normX*A'
wolffd@0: normX = (net.X - ones(size(net.X))*diag(mean(net.X)))*diag(1./std(net.X));
wolffd@0: net.rbfnet.w2 = Phi \ (normX*A');
wolffd@0: % Bias is mean of target data
wolffd@0: net.rbfnet.b2 = mean(data);
wolffd@0: 
wolffd@0: % Must also set initial value of variance
wolffd@0: % Find average distance between nearest centres
wolffd@0: % Ensure that distance of centre to itself is excluded by setting diagonal
wolffd@0: % entries to realmax
wolffd@0: net.gmmnet.centres = rbffwd(net.rbfnet, net.X);
wolffd@0: d = dist2(net.gmmnet.centres, net.gmmnet.centres) + ...
wolffd@0:   diag(ones(net.gmmnet.ncentres, 1)*realmax);
wolffd@0: sigma = mean(min(d))/2;
wolffd@0: 
wolffd@0: % Now set covariance to minimum of this and next largest eigenvalue
wolffd@0: if net.dim_latent < size(data, 2)
wolffd@0:   sigma = min(sigma, PCcoeff(net.dim_latent+1));
wolffd@0: end
wolffd@0: net.gmmnet.covars = sigma*ones(1, net.gmmnet.ncentres);
wolffd@0: 
wolffd@0: % Sub-function to create the sample data in 2d
wolffd@0: function sample = gtm_rctg(samp_size)
wolffd@0: 
wolffd@0: xDim = samp_size(1);
wolffd@0: yDim = samp_size(2);
wolffd@0: % Produce a grid with the right number of rows and columns
wolffd@0: [X, Y] = meshgrid([0:1:(xDim-1)], [(yDim-1):-1:0]);
wolffd@0: 
wolffd@0: % Change grid representation 
wolffd@0: sample = [X(:), Y(:)];
wolffd@0: 
wolffd@0: % Shift grid to correct position and scale it
wolffd@0: maxXY= max(sample);
wolffd@0: sample(:,1) = 2*(sample(:,1) - maxXY(1)/2)./maxXY(1);
wolffd@0: sample(:,2) = 2*(sample(:,2) - maxXY(2)/2)./maxXY(2);
wolffd@0: return;
wolffd@0: 
wolffd@0:    
wolffd@0: