disp('Running MATLAB script file EGkernelPCA.m') ;
%
%    For creating graphics exploring
%    kernel PCA 
%    using toy data sets in 3d
%
%    This uses lines from Bernhard Sch\"olkopf's program:
%        kpca_toy.m
%    which is available from the Kernel Machines web site,
%    and has comment lines:
%
% Kernel PCA toy example for k(x,y)=exp(-||x-y||^2/rbf_var), cf. Fig. 4 in 
% @article{SchSmoMue98,
%   author    = "B.~{Sch\"olkopf} and A.~Smola and K.-R.~{M\"uller}",
%   title     = "Nonlinear component analysis as a kernel Eigenvalue problem",
%   journal =	 {Neural Computation},
%   volume    = 10,
%   issue     = 5,
%   pages     = "1299 -- 1319",
%   year      = 1998}
% This file can be downloaded from http://www.kernel-machines.org.
% Last modified: 4 July 2003



ipart = 12 ;    %  1 - Good h, Colored Raw Data Plot, 
               %              Colored PCA Plot, kernel Projections
               %  2 - Too small h, kernel Projections
               %  3 - Too small h, kernel Projections
               %  11 - movie through h, kernel projections
               %  12 - movie through h, kernel projections, fixed PCA
               %  13 - movie through h, kernel projections, SVD



idat = 8 ;    %  1 - 4 well separated clusters, in axis directions
              %  2 - 7 clusters, in axis directions
              %  3 - 4 clusters, different sizes and close
              %  4 - 4 clusters, different sizes and very close
              %  5 - 2 Rings
              %  6 - Surface of a sphere
              %  7 - Standard Normal
              %  8 - Same Standard Normal, cluster colored



%  Generate data set
%  and set parameters
%
if idat == 1 ;    %  4 well separated clusters in 3d

  seed = 49385739857 ;
  randn('state',seed) ;

  mu1 = [0; 0; 30] ;
  mu2 = [0; 0; -30] ;
  mu3 = [3; 15; 0] ;
  mu4 = [3; -15; 0] ;

  n1 = 20 ;
  n2 = 20 ;
  n3 = 20 ;
  n4 = 20 ;

  mdata = randn(3,n1) + vec2matSM(mu1,n1) ;
  mdata = [mdata, (randn(3,n2) + vec2matSM(mu2,n2))] ;
  mdata = [mdata, (randn(3,n3) + vec2matSM(mu3,n3))] ;
  mdata = [mdata, (randn(3,n4) + vec2matSM(mu4,n4))] ;

  datstr = 'WS4ClustAxDir' ;
  datatitlestr1 = '4 Widely Separated Clusters' ;
  datatitlestr2 = 'Oriented in Axes Directions' ;

  goodh = 2 ;
  smallh = 0.2 ;
  bigh = 100 ;

  mcolor = [ones(n1,1) * [1 0 0] ; ...
            ones(n2,1) * [0 1 0] ; ...
            ones(n3,1) * [0 0 1] ; ...
            ones(n4,1) * [0 1 1]] ;

  markerstr = [] ;
  for i = 1:n1 ;
    markerstr = strvcat(markerstr,'+') ;
  end ;
  for i = 1:n2 ;
    markerstr = strvcat(markerstr,'x') ;
  end ;
  for i = 1:n3 ;
    markerstr = strvcat(markerstr,'s') ;
  end ;
  for i = 1:n4 ;
    markerstr = strvcat(markerstr,'d') ;
  end ;


elseif idat == 2 ;    %  7 clusters, in axis directions

  seed = 49385739857 ;
  randn('state',seed) ;

  mu1 = [0; 0; 0] ;
  mu2 = [20; 0; 0] ;
  mu3 = [-20; 0; 0] ;
  mu4 = [0; 20; 0] ;
  mu5 = [0; -20; 0] ;
  mu6 = [0; 0; 20] ;
  mu7 = [0; 0; -20] ;

  n1 = 20 ;
  n2 = 20 ;
  n3 = 20 ;
  n4 = 20 ;
  n5 = 20 ;
  n6 = 20 ;
  n7 = 20 ;

  mdata = randn(3,n1) + vec2matSM(mu1,n1) ;
  mdata = [mdata, (randn(3,n2) + vec2matSM(mu2,n2))] ;
  mdata = [mdata, (randn(3,n3) + vec2matSM(mu3,n3))] ;
  mdata = [mdata, (randn(3,n4) + vec2matSM(mu4,n4))] ;
  mdata = [mdata, (randn(3,n5) + vec2matSM(mu5,n5))] ;
  mdata = [mdata, (randn(3,n6) + vec2matSM(mu6,n6))] ;
  mdata = [mdata, (randn(3,n7) + vec2matSM(mu7,n7))] ;

  datstr = 'WS7ClustAxDir' ;
  datatitlestr1 = '7 Widely Separated Clusters' ;
  datatitlestr2 = 'Oriented in Axes Directions' ;

  goodh = 2 ;
  smallh = 0.3 ;
  bigh = 100 ;

  mcolor = [ones(n1,1) * [0 0 0] ; ...
            ones(n2,1) * [1 0 0] ; ...
            ones(n3,1) * [0 1 0] ; ...
            ones(n4,1) * [0 0 1] ; ...
            ones(n5,1) * [0 1 1] ; ...
            ones(n6,1) * [1 0 1] ; ...
            ones(n7,1) * [1 1 0]] ;

  markerstr = 'o' ;
%  for i = 1:n1 ;
%    markerstr = strvcat(markerstr,'+') ;
%  end ;
%  for i = 1:n2 ;
%    markerstr = strvcat(markerstr,'x') ;
%  end ;
%  for i = 1:n3 ;
%    markerstr = strvcat(markerstr,'s') ;
%  end ;
%  for i = 1:n4 ;
%    markerstr = strvcat(markerstr,'d') ;
%  end ;


elseif idat == 3 ;    %  4 clusters, different sizes and close

  seed = 49385739857 ;
  randn('state',seed) ;

  mu1 = [0; 0; 30] ;
  mu2 = [0; 0; -30] ;
  mu3 = [3; 2.5; 0] ;
  mu4 = [3; -2.5; 0] ;

  sig1 = 0.25 ;
  sig2 = 4 ;
  sig3 = 1 ;
  sig4 = 1 ;

  n1 = 20 ;
  n2 = 20 ;
  n3 = 20 ;
  n4 = 20 ;

  mdata = sig1 * randn(3,n1) + vec2matSM(mu1,n1) ;
  mdata = [mdata, (sig2 * randn(3,n2) + vec2matSM(mu2,n2))] ;
  mdata = [mdata, (sig3 * randn(3,n3) + vec2matSM(mu3,n3))] ;
  mdata = [mdata, (sig4 * randn(3,n4) + vec2matSM(mu4,n4))] ;

  datstr = 'DS4ClustAxDir' ;
  datatitlestr1 = '4 Different Size and Close' ;
  datatitlestr2 = 'Oriented in Axes Directions' ;

  goodh = 2 ;
  smallh = 0.2 ;
  bigh = 100 ;

  mcolor = [ones(n1,1) * [1 0 0] ; ...
            ones(n2,1) * [0 1 0] ; ...
            ones(n3,1) * [0 0 1] ; ...
            ones(n4,1) * [0 1 1]] ;

  markerstr = [] ;
  for i = 1:n1 ;
    markerstr = strvcat(markerstr,'+') ;
  end ;
  for i = 1:n2 ;
    markerstr = strvcat(markerstr,'x') ;
  end ;
  for i = 1:n3 ;
    markerstr = strvcat(markerstr,'s') ;
  end ;
  for i = 1:n4 ;
    markerstr = strvcat(markerstr,'d') ;
  end ;


elseif idat == 4 ;    %  4 clusters, different sizes and very close

  seed = 49385739857 ;
  randn('state',seed) ;

  mu1 = [0; 0; 30] ;
  mu2 = [0; 0; -30] ;
  mu3 = [3; 1.5; 0] ;
  mu4 = [3; -1.5; 0] ;

  sig1 = 0.25 ;
  sig2 = 4 ;
  sig3 = 1 ;
  sig4 = 1 ;

  n1 = 20 ;
  n2 = 20 ;
  n3 = 20 ;
  n4 = 20 ;

  mdata = sig1 * randn(3,n1) + vec2matSM(mu1,n1) ;
  mdata = [mdata, (sig2 * randn(3,n2) + vec2matSM(mu2,n2))] ;
  mdata = [mdata, (sig3 * randn(3,n3) + vec2matSM(mu3,n3))] ;
  mdata = [mdata, (sig4 * randn(3,n4) + vec2matSM(mu4,n4))] ;

  datstr = 'VC4ClustAxDir' ;
  datatitlestr1 = '4 Different Size and Very Close' ;
  datatitlestr2 = 'Oriented in Axes Directions' ;

  goodh = 2 ;
  smallh = 0.2 ;
  bigh = 100 ;

  mcolor = [ones(n1,1) * [1 0 0] ; ...
            ones(n2,1) * [0 1 0] ; ...
            ones(n3,1) * [0 0 1] ; ...
            ones(n4,1) * [0 1 1]] ;

  markerstr = [] ;
  for i = 1:n1 ;
    markerstr = strvcat(markerstr,'+') ;
  end ;
  for i = 1:n2 ;
    markerstr = strvcat(markerstr,'x') ;
  end ;
  for i = 1:n3 ;
    markerstr = strvcat(markerstr,'s') ;
  end ;
  for i = 1:n4 ;
    markerstr = strvcat(markerstr,'d') ;
  end ;


elseif idat == 5 ;    %  2 Rings

  seed = 49385739857 ;
  randn('state',seed) ;

  n1 = 40 ;
  n2 = 40 ;

  mdata = randn(2,n1+n2) ;
  mdata = mdata ./ vec2matSM(sqrt(ones(1,2) * mdata.^2),2) ;
      %  divide by length in 2d

  mu1 = 2 ;
  mu2 = -2 ;

  mdata = [mdata; [(ones(1,n1) * mu1) (ones(1,n2) * mu2)]] ;


  datstr = '2Rings' ;
  datatitlestr1 = '2 Rings' ;
  datatitlestr2 = 'Oriented in Axes Directions' ;

  goodh = 0.3 ;
  smallh = 0.01 ;
  bigh = 10 ;

  mcolor = [ones(n1,1) * [1 0 0] ; ...
            ones(n2,1) * [0 0 1]] ;

  markerstr = [] ;
  for i = 1:n1 ;
    markerstr = strvcat(markerstr,'+') ;
  end ;
  for i = 1:n2 ;
    markerstr = strvcat(markerstr,'x') ;
  end ;


elseif idat == 6 ;    %  Surface of a sphere

  seed = 49385739857 ;
  randn('state',seed) ;

  n = 100 ;

  mdata = randn(3,n) ;
  mdata = mdata ./ vec2matSM(sqrt(ones(1,3) * mdata.^2),3) ;
      %  divide by length in 3d


  datstr = 'Sphere' ;
  datatitlestr1 = 'Surface of a Sphere in 3d' ;
  datatitlestr2 = ' ' ;

  goodh = 0.3 ;
  smallh = 0.03 ;
  bigh = 2 ;

  mcolor = [ones(n,1) * [0 1 0]] ;

  markerstr = [] ;
  for i = 1:n ;
    markerstr = strvcat(markerstr,'o') ;
  end ;


elseif idat == 7 ;    %  Standard normal

  seed = 49385739857 ;
  randn('state',seed) ;

  n = 100 ;

  mdata = randn(3,n) ;
      %  divide by length in 3d


  datstr = 'SNorm' ;
  datatitlestr1 = 'Standard normal' ;
  datatitlestr2 = ' ' ;

%  goodh = 0.3 ;
  goodh = 0.17 ;
  smallh = 0.1 ;
  bigh = 10 ;

  mcolor = [ones(n,1) * [1 0 1]] ;

  markerstr = [] ;
  for i = 1:n ;
    markerstr = strvcat(markerstr,'+') ;
  end ;


elseif idat == 8 ;    %  Same Standard Normal, cluster colored

  seed = 49385739857 ;
  randn('state',seed) ;

  n = 100 ;

  mdata = randn(3,n) ;
      %  divide by length in 3d


  datstr = 'SNorm' ;
  datatitlestr1 = 'Standard normal' ;
  datatitlestr2 = 'KPCA Cluster Colored' ;

  goodh = 0.17 ;
  smallh = 0.1 ;
  bigh = 10 ;

  mcolor = [ones(n,1) * [0 0 0]] ;

  vindblue = [7 14 34 38 47 61 64 69 77 88 89 96] ;
  vindred = [3 4 5 9 20 25 26 27 30 32 44 55 72 81 83] ;
  vindgreen = [13 35 41 44 53 56 60 70 71 73 75 79 90] ;
      %  These were chosen by running  idat = 7  &  ipart = 1
      %  then fiddling with inequality constraints and
      %  the matrix  features

  mcolor(vindblue',:) = [ones(length(vindblue),1) * [0 0 1]] ;
  mcolor(vindred',:) = [ones(length(vindred),1) * [1 0 0]] ;
  mcolor(vindgreen',:) = [ones(length(vindgreen),1) * [0 1 0]] ;


  markerstr = [] ;
  for i = 1:n ;
    markerstr = strvcat(markerstr,'+') ;
  end ;


end ;




%  Do analysis
%
if  ipart == 1  | ...
    ipart == 2  | ...
    ipart == 3  ;


  if ipart == 1 ;    %  Good h
    h = goodh ;
  elseif ipart == 2 ;   %  Too small h
    h = smallh ;
  elseif ipart == 3 ;   %  Too big h
    h = bigh ;
  end ;
  hstr = ['h = ' num2str(h)] ; 

  %  My modifcation of Bernhard's lines
  %  
  % parameters
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %rbf_var = 0.1;
  rbf_var = h;
  %xnum = 4;
  %ynum = 2;
  %max_ev = xnum*ynum;
  max_ev = 6;
  % (extract features from the first <max_ev> Eigenvectors)
  x_test_num = 15;
  y_test_num = 15;
  %cluster_pos = [-0.5 -0.2; 0 0.6; 0.5 0];
  %cluster_size = 30;

  % generate a toy data set
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %num_clusters = size(cluster_pos,1);
  %train_num = num_clusters*cluster_size;
  %patterns = zeros(train_num, 2);
  %range = 1;
  %randn('seed', 0);
  %for i=1:num_clusters,
  %  patterns((i-1)*cluster_size+1:i*cluster_size,1) = cluster_pos(i,1)+0.1*randn(cluster_size,1);
  %  patterns((i-1)*cluster_size+1:i*cluster_size,2) = cluster_pos(i,2)+0.1*randn(cluster_size,1);
  %end
  patterns = mdata' ;
  %test_num = x_test_num*y_test_num;
  %x_range = -range:(2*range/(x_test_num - 1)):range;
  %y_offset = 0.5;
  %y_range = -range+ y_offset:(2*range/(y_test_num - 1)):range+ y_offset;
  %[xs, ys] = meshgrid(x_range, y_range);
  %test_patterns(:, 1) = xs(:);
  %test_patterns(:, 2) = ys(:);
  %cov_size = train_num;  % use all patterns to compute the covariance matrix
  cov_size = size(mdata,2) ;
  train_num = cov_size ;

  % carry out Kernel PCA
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for i=1:cov_size,
    for j=i:cov_size,
      K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
      K(j,i) = K(i,j);
    end
  end
  unit = ones(cov_size, cov_size)/cov_size;
  % centering in feature space!
  K_n = K - unit*K - K*unit + unit*K*unit;

  [evecs,evals] = eig(K_n);
  evals = real(diag(evals));
  for i=1:cov_size,
    evecs(:,i) = evecs(:,i)/(sqrt(evals(i)));
  end

  % extract features
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%  do not need the following here - only need test point features
  %   NO!, these are now used, so have uncommented 
  unit_train = ones(train_num,cov_size)/cov_size;
  for i=1:train_num,
    for j=1:cov_size,
      K_train(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
    end
  end
  K_train_n = K_train - unit_train*K - K_train*unit + unit_train*K*unit;
  features = zeros(train_num, max_ev);
  features = K_train_n * evecs(:,1:max_ev);

  %unit_test = ones(test_num,cov_size)/cov_size;
  %K_test = zeros(test_num,cov_size);
  %for i=1:test_num,
  %  for j=1:cov_size,
  %    K_test(i,j) = exp(-norm(test_patterns(i,:)-patterns(j,:))^2/rbf_var);
  %  end
  %end
  %K_test_n = K_test - unit_test*K - K_test*unit + unit_test*K*unit;
  %test_features = zeros(test_num, max_ev);
  %test_features = K_test_n * evecs(:,1:max_ev);

  % plot it
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %figure(1); clf
  %for n = 1:max_ev,
  %  subplot(ynum, xnum, n);
  %  axis([-range range -range+y_offset range+y_offset]);
  %  imag = reshape(test_features(:,n), y_test_num, x_test_num);
  %  axis('xy')
  %  colormap(gray);
  %  hold on;
  %  pcolor(x_range, y_range, imag);
  %  shading interp
  %  contour(x_range, y_range, imag, 9, 'b');
  %  plot(patterns(:,1), patterns(:,2), 'r.')
  %  text(-1,1.65,sprintf('Eigenvalue=%4.3f', evals(n)));
  %  hold off;
  %end
  


  if ipart == 1 ;

    %  Make data coordinate plot
    %
    figure(1) ;
    clf ;

    titlecellstr = {{'Data Coordinate Directions' datatitlestr1 datatitlestr2}} ;

    labelcellstr = {{'Data Coor. 1'; 'Data Coor. 2'; 'Data Coor. 3'}} ;

    paramstruct = struct('icolor',mcolor, ...
                         'npcadiradd',0, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdata,eye(3),paramstruct) ;

      orient landscape ;
      outstr = ['EGkernelPCAd' num2str(idat) datstr ...
                           'p' num2str(ipart) 'GoodhDatCoo.ps'] ;
    print('-dpsc',outstr) ;



    %  Make PCA coordinate plot
    %
    figure(2) ;
    clf ;

    titlecellstr = {{'PCA Coordinate Directions' datatitlestr1 datatitlestr2}} ;

    labelcellstr = {{'PC 1'; 'PC 2'; 'PC 3'}} ;

    paramstruct = struct('icolor',mcolor, ...
                         'npcadiradd',3, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdata,[],paramstruct) ;

      orient landscape ;
      outstr = ['EGkernelPCAd' num2str(idat) datstr ...
                           'p' num2str(ipart) 'GoodhPCACoo.ps'] ;
    print('-dpsc',outstr) ;

  end ;


  %  Make Kernel Space Coordinate plot
  %
  figure(3) ;
  clf ;

  titlecellstr = {{'Kernel Space PCA Coord. Dir''ns' ' ' datatitlestr1 ' ' datatitlestr2 hstr}} ;

  labelcellstr = {{'KS PC 1'; 'KS PC 2'; 'KS PC 3'; 'KS PC 4'; 'KS PC 5'; 'KS PC 6'}} ;

  paramstruct = struct('icolor',mcolor, ...
                       'npcadiradd',0, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'iscreenwrite',1) ;
  scatplotSM(features',eye(6),paramstruct) ;

    orient landscape ;
    if ipart == 1 ;
      outstr = ['EGkernelPCAd' num2str(idat) datstr ...
                           'p' num2str(ipart) 'GoodhKSPCACoo.ps'] ;
    elseif ipart == 2 ;
      outstr = ['EGkernelPCAd' num2str(idat) datstr ...
                           'p' num2str(ipart) 'SmallhKSPCACoo.ps'] ;
    elseif ipart == 3 ;
      outstr = ['EGkernelPCAd' num2str(idat) datstr ...
                           'p' num2str(ipart) 'GoodhKSPCACoo.ps'] ;
    end ;
  print('-dpsc',outstr) ;



elseif  ipart == 11  | ...
        ipart == 12  | ...
        ipart == 13  ;    %  Make movie through scale

  nh = 21 ;
  vh = logspace(log10(smallh),log10(bigh),nh) ;


  %  Set parameters
  %
  moviefps = 5 ;
  moviecstr = 'Cinepak' ;
  nrepeat = 2 ;
  nframe = 11 ;


  %  Initiate figure and movie
  %
  close all
  fighand = figure(1) ;
  clf ;
  vpos = get(fighand,'Position') ;
  set(fighand,'Position',[vpos(1)/4 vpos(2)/4 ...
                          (3*vpos(1)/4 + vpos(3)) ...
                          (3*vpos(2)/4 + vpos(4))]) ;
  clear moviestruct ;


  %  Do main movie loop
  %
  for ih = 1:nh ;

    h = vh(ih) ;
    hstr = ['h = ' num2str(h)] ; 

    %  My modifcation of Bernhard's lines
    %  
    % parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %rbf_var = 0.1;
    rbf_var = h;
    %xnum = 4;
    %ynum = 2;
    %max_ev = xnum*ynum;
    max_ev = 4;
    % (extract features from the first <max_ev> Eigenvectors)
    x_test_num = 15;
    y_test_num = 15;
    %cluster_pos = [-0.5 -0.2; 0 0.6; 0.5 0];
    %cluster_size = 30;

    cov_size = size(mdata,2) ;
    train_num = cov_size ;
    patterns = mdata' ;

    % carry out Kernel PCA
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:cov_size,
      for j=i:cov_size,
        K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
        K(j,i) = K(i,j);
      end
    end
    unit = ones(cov_size, cov_size)/cov_size;
    % centering in feature space!
    if  ipart == 11  |  ipart == 12 ;
      K_n = K - unit*K - K*unit + unit*K*unit;
    elseif  ipart == 13 ;
      K_n = K ;
    end ;    %  else don't center

    [evecs,evals] = eig(K_n);
    evals = real(diag(evals));
    for i=1:cov_size,
      evecs(:,i) = evecs(:,i)/(sqrt(evals(i)));
    end

    % extract features
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%  do not need the following here - only need test point features
    %   NO!, these are now used, so have uncommented 
    unit_train = ones(train_num,cov_size)/cov_size;
    for i=1:train_num,
      for j=1:cov_size,
        K_train(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
      end
    end
    K_train_n = K_train - unit_train*K - K_train*unit + unit_train*K*unit;
    features = zeros(train_num, max_ev);
    features = K_train_n * evecs(:,1:max_ev);



    if  ipart == 12  | ...
        ipart == 13  ;    %  then normalize vector directions

      for iev = 1:max_ev ;
        featev = features(:,iev) ;
        [temp,maxind] = max(abs(featev)) ;
        if featev(maxind) < 0 ;    %  then reverse direction
          features(:,iev) = -featev ;
        end ;
      end ;

    end ;



    %  Make Kernel Space Coordinate plot
    %
    titlecellstr = {{'Kernel Space PCA' datatitlestr1 datatitlestr2 hstr}} ;

    labelcellstr = {{'KS PC 1'; 'KS PC 2'; 'KS PC 3'; 'KS PC 4'}} ;

    paramstruct = struct('icolor',mcolor, ...
                         'npcadiradd',0, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'iscreenwrite',1) ;
    scatplotSM(features',eye(4),paramstruct) ;


    moviestruct(ih) = getframe(fighand) ;

  end ;    %  of iframe loop, to make movie



  %  Reorder frames, and replay
  %
  vmorder = [(1:nh),((nh-1):-1:2)] ;
  moviestruct = moviestruct(vmorder) ;
          %  reorder frames, to loop back to beginning

  movie(fighand,moviestruct,nrepeat,moviefps) ;
          %  Play movie on screen



  %  Save results
  %
  if  ipart == 11  | ...
      ipart == 12  ;
    savestr = ['EGkernelPCAd' num2str(idat) datstr ...
                          'p' num2str(ipart) 'hMovieKSPCACoo'] ;
  elseif  ipart == 13 ;
    savestr = ['EGkernelPCAd' num2str(idat) datstr ...
                          'p' num2str(ipart) 'hMovieKSSVDCoo'] ;
  end ;
  movie2avi(moviestruct,savestr,'compression',moviecstr, ...
                        'keyframe',moviefps,'fps',moviefps) ;







end ;



