function [s,R] = test_lwpr_nD(s,R)

  global lwprs;

  % test for LWPR on n-dimensional data set with 2D embedded function
  d = 5;
  n = 500;

  % a random training set using the CROSS function
  X = (rand(n,2)-.5)*2;
  Y = max([exp(-X(:,1).^2 * 10),exp(-X(:,2).^2 * 50),1.25*exp(-(X(:,1).^2+X(:,2).^2)*5)]');
  Y = Y' + randn(n,1)*0.1;

  % rotate input data into d-dimensional space
  if ~exist('R') || isempty(R)
    R = randn(d);
    R = R'*R;
    R = orth(R);
  end
  Xorg = X;
  X = [X zeros(n,d-2)]*R;

  % a systematic test set on a grid
  Xt = [];
  for i=-1:0.05:1,
    for j=-1:0.05:1,
      Xt = [Xt; i j];
    end
  end
  Yt = max([exp(-Xt(:,1).^2 * 10),exp(-Xt(:,2).^2 * 50),1.25*exp(-(Xt(:,1).^2+Xt(:,2).^2)*5)]');
  Yt = Yt';

  % rotate the test data
  Xtorg = Xt;
  Xt = [Xt zeros(length(Xt),d-2)]*R;

  % initialize LWPR, use only diagonal distance metric
  ID = 1;
  if ~exist('s') || isempty(s)
  	sprintf("noexists")%DELETE THIS%DELETE THIS%DELETE THIS%DELETE THIS
  	%lwprs(case,ID,n_in,n_out,diag_only,meta,meta_rate,penalty,init_alpha,norm,norm_out,  name)
    lwpr('Init',ID, d ,  1 ,    1,       0,    0,       1e-7,      50,  ones(d,1),[1],'lwpr_test');
  else
  	sprintf("exists")%DELETE THIS%DELETE THIS%DELETE THIS%DELETE THIS%DELETE THIS
    lwpr('Init',ID,s);
  end

  % set some parameters
  kernel = 'Gaussian';
  % kernel = 'BiSquare'; % note: the BiSquare kernel requires different values for
  %                              an initial distance metric, as in the next line
  % lwpr('Change',ID,'init_D',eye(d)*7);

  lwpr('Change',ID,'init_D',eye(d)*25); 
  lwpr('Change',ID,'init_alpha',ones(d)*100);     % this is a safe learning rate
  lwpr('Change',ID,'w_gen',0.2);                  % more overlap gives smoother surfaces
  lwpr('Change',ID,'meta',1);                     % meta learning can be faster, but numerical more dangerous
  lwpr('Change',ID,'meta_rate',100);
  lwpr('Change',ID,'init_lambda',0.995);
  lwpr('Change',ID,'final_lambda',0.9999);
  lwpr('Change',ID,'tau_lambda',0.9999);


  % train the model
  for j=1:20
    inds = randperm(n);
    mse = 0;
    for i=1:n,
    	sprintf('i,j: %d,%d', i, j)
                 %lwpr(updata,  ID,     x,         y,        e_t, alpha, composite_control)
      [yp,w,np] = lwpr('Update',ID,X(inds(i),:)',Y(inds(i),:)');
      mse = mse + (Y(inds(i),:)-yp).^2;
    end
    nMSE = mse/n/var(Y,1);	
    disp(sprintf('#Data=%d #rfs=%d nMSE=%5.3f #proj=%4.2f (TrainingSet)',lwprs(ID).n_data,length(lwprs(ID).rfs),nMSE,np));
  end

  % create predictions for the test data
  Yp = zeros(size(Yt));
  for i=1:length(Xt),
    [yp,w]=lwpr('Predict',ID,Xt(i,:)',0.001);
    Yp(i,1) = yp;
  end
  ep   = Yt-Yp;
  mse  = mean(ep.^2);
  nmse = mse/var(Y,1);
  disp(sprintf('#Data=%d #rfs=%d nMSE=%5.3f #proj=%4.2f (TestSet)',lwprs(ID).n_data,length(lwprs(ID).rfs),nmse,np));

  % get the data structure
  s = lwpr('Structure',ID);

  figure(1);
  clf;

  % plot the raw noisy data
  subplot(2,2,1);
  plot3(Xorg(:,1),Xorg(:,2),Y,'*');
  title('Noisy data samples');

  % plot the fitted surface
  axis([-1 1 -1 1 -.5 1.5]);
  subplot(2,2,2);
  [x,y,z]=makesurf([Xtorg,Yp],sqrt(length(Xtorg)));
  surfl(x,y,z);
  axis([-1 1 -1 1 -.5 1.5]);
  title(sprintf('The fitted function: nMSE=%5.3f',nmse));

  % plot the true surface
  subplot(2,2,3);
  [x,y,z]=makesurf([Xtorg,Yt],sqrt(length(Xtorg)));
  surfl(x,y,z);
  axis([-1 1 -1 1 -.5 1.5]);
  title('The true function');

  % plot the local models
  subplot(2,2,4);
  for i=1:length(lwprs(ID).rfs),
    D = R'*lwprs(ID).rfs(i).D*R;
    c = R*lwprs(ID).rfs(i).c;
    draw_ellipse(D(1:2,1:2),c(1:2),0.1,kernel);
    hold on;
  end
  plot(Xorg(:,1),Xorg(:,2),'.');
  hold off;
  axis('equal');
  title('Projected input space view of RFs');

% --------------------------------------------------------------------------------
function [X,Y,Z]=makesurf(data,nx)
% [X,Y,Z]=makesurf(data,nx) converts the 3D data file data into
% three matices as need by surf(). nx tells how long the row of the
% output matrices are

  [m,n]=size(data);

  n=0;
  for i=1:nx:m,
    n = n+1;
    X(:,n) = data(i:i+nx-1,1);
    Y(:,n) = data(i:i+nx-1,2);
    Z(:,n) = data(i:i+nx-1,3);
  end;


% --------------------------------------------------------------------------------
function []=draw_ellipse(M,C,w,kernel)
% function draw ellipse draws the ellipse corresponding to the
% eigenvalues of M at the location c.

  [V,E] = eig(M);

  E = E;
  d1 = E(1,1);
  d2 = E(2,2);

  steps = 50;
  switch kernel
   case 'Gaussian'
    start = sqrt(-2*log(w)/d1);
   case 'BiSquare'
    start = sqrt(2*(1-sqrt(w))/d1);
  end


  for i=0:steps,
    Xp(i+1,1) = -start + i*(2*start)/steps;
    switch kernel
     case 'Gaussian'
      arg = (-2*log(w)-Xp(i+1,1)^2*d1)/d2;
     case 'BiSquare'
      arg = (2*(1-sqrt(w))-Xp(i+1,1)^2*d1)/d2;
    end
    if (arg < 0), 
      arg = 0; 
    end; % should be numerical error
    Yp(i+1,1) = sqrt(arg);
  end;

  for i=1:steps+1;
    Xp(steps+1+i,1) = Xp(steps+1-i+1,1);
    Yp(steps+1+i,1) = -Yp(steps+1-i+1,1);
  end;

  % transform the rf

  M = [Xp,Yp]*V(1:2,1:2)';

  Xp = M(:,1) + C(1);
  Yp = M(:,2) + C(2);

  plot(C(1),C(2),'ro',Xp,Yp,'c');
