disp('Running MATLAB script file gpnprt.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION gpnpr,
%    General Purpose NonParametric Regression


itest = 14 ;     %  1,2,3,4,5,6,7,8,9,10,11,12,13,14

format compact ;
format short ;


if itest == 1 ;   %  Simplest visual check of implementations of LL
  indat = [[.2, 1];...
           [.3, 2];...
           [.5, 1];...
           [.8, 3];...
           [.4, 2]];
  invh = .1 ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpnpr(indat,invh,vxg,1) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

  titstr = 'formula looped impl.'
  subplot(2,2,2) ;
    tic ;
    gpnpr(indat,invh,vxg,2) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

  titstr = 'polyfit looped impl.'
  subplot(2,2,3) ;
    tic ;
    gpnpr(indat,invh,vxg,3,1) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;


  titstr = 'binned implementation' 
  subplot(2,2,4) ;
    tic ;
    gpnpr(indat,invh,vxg,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

elseif itest == 2 ;   %  Simplest visual check of implementations of NW
  indat = [[.2, 1];...
           [.3, 2];...
           [.5, 1];...
           [.8, 3];...
           [.4, 2]];
  invh = .1 ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpnpr(indat,invh,vxg,1,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

  titstr = 'formula looped impl.'
  subplot(2,2,2) ;
    tic ;
    gpnpr(indat,invh,vxg,2,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

  titstr = 'polyfit looped impl.'
  subplot(2,2,3) ;
    tic ;
    gpnpr(indat,invh,vxg,3,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

  titstr = 'binned implementation' 
  subplot(2,2,4) ;
    tic ;
    gpnpr(indat,invh,vxg,0,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'x') ;
      hold off ;

elseif itest == 3 ;   %  Simple check of higher degree poly's
  clf ;
  indat = [[.2, 1];...
           [.3, 2];...
           [.5, 1];...
           [.8, 3];...
           [.4, 2]];
  invh = .1 ;
  vxg = [0; 1] ;
  titstr = 'Check poly degrees 0-3'
  tic ;
  gpnpr(indat,invh,vxg,1,0) ;
  toc
    title(titstr) ;
    hold on ;
    tic ;
    gpnpr(indat,invh,vxg,3,1) ;
    toc
    tic ;
    gpnpr(indat,invh,vxg,3,2) ;
    toc
    tic ;
    gpnpr(indat,invh,vxg,3,3) ;
    toc
    plot(indat(:,1),indat(:,2),'x') ;
   hold off ;

elseif itest == 4 ;   %  Check multiple bandwidths, LL
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  invh = [.01 .03 .1] ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpnpr(indat,invh,vxg,1) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'formula looped implementation'
  subplot(2,2,2) ;
    tic ;
    gpnpr(indat,invh,vxg,2) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'polyfit looped implementation'
  subplot(2,2,3) ;
    tic ;
    gpnpr(indat,invh,vxg,3,1) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'binned implementation' 
  subplot(2,2,4) ;
    tic ;
    gpnpr(indat,invh,vxg,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;

elseif itest == 5 ;   %  Check multiple bandwidths, NW
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  invh = [.01 .03 .1] ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
  subplot(2,2,1) ;
    tic ;
    gpnpr(indat,invh,vxg,1,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'formula looped implementation'
  subplot(2,2,2) ;
    tic ;
    gpnpr(indat,invh,vxg,2,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'polyfit looped implementation'
  subplot(2,2,3) ;
    tic ;
    gpnpr(indat,invh,vxg,3,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;
  titstr = 'binned implementation' 
  subplot(2,2,4) ;
    tic ;
    gpnpr(indat,invh,vxg,0,0) ;
    toc
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'.') ;
      hold off ;

elseif itest == 6 ;   %  Check implementations against each other, LL
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  invh = .03 ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
    onpr1 = gpnpr(indat,invh,vxg,1,1) ;
  titstr = 'formula looped implementation'
    onpr2 = gpnpr(indat,invh,vxg,2,1) ;
  titstr = 'polyfit looped implementation'
    onpr3 = gpnpr(indat,invh,vxg,3,1) ;
  titstr = 'binned implementation' 
    onpr0 = gpnpr(indat,invh,vxg,0,1) ;
  titstr = 'pre-binned implementation' 
    bincts = gplbinr(indat,vxg) ;
    onprm1 = gpnpr(bincts,invh,vxg,-1,1) ;

  disp('Check these are 0 (401 grid points):') ;
  disp('    matrix - formula loop:') ;
  max(abs(onpr1 - onpr2)) 
  disp('    matrix - polyfit loop:') ;
  max(abs(onpr1 - onpr3)) 
  disp('    matrix - binned:') ;
  max(abs(onpr1 - onpr0)) 
  disp('    matrix - prebinned:') ;
  max(abs(onpr1 - onprm1)) 
  disp('    binned - prebinned:') ;
  max(abs(onpr0 - onprm1)) 

elseif itest == 7 ;   %  Check implementations against each other, NW
  indat = rand(500,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(500,1)] ;
  invh = .03 ;
  vxg = [0; 1] ;
  titstr = 'matrix implementation'
    onpr1 = gpnpr(indat,invh,vxg,1,0) ;
  titstr = 'formula looped implementation'
    onpr2 = gpnpr(indat,invh,vxg,2,0) ;
  titstr = 'polyfit looped implementation'
    onpr3 = gpnpr(indat,invh,vxg,3,0) ;
  titstr = 'binned implementation' 
    onpr0 = gpnpr(indat,invh,vxg,0,0) ;
  titstr = 'pre-binned implementation' 
    bincts = gplbinr(indat,vxg) ;
    onprm1 = gpnpr(bincts,invh,vxg,-1,0) ;

  disp('Check these are 0 (401 grid points):') ;
  disp('    matrix - formula loop:') ;
  max(abs(onpr1 - onpr2)) 
  disp('    matrix - polyfit loop:') ;
  max(abs(onpr1 - onpr3)) 
  disp('    matrix - binned:') ;
  max(abs(onpr1 - onpr0)) 
  disp('    matrix - prebinned:') ;
  max(abs(onpr1 - onprm1)) 
  disp('    binned - prebinned:') ;
  max(abs(onpr0 - onprm1)) 

elseif itest == 8 ;   %  Check endpt handling
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .05*randn(100,1)] ;
  invh = .03 ;
  subplot(2,2,1) ;
    gpnpr(indat,invh,[.2,.5],0,1) ;
      title('LL - no trunc.') ;
      axis(axis) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;
  subplot(2,2,2) ;
    gpnpr(indat,invh,[.2,.5],0,1,1) ;
      title('LL - trunc.') ;
      axis(axis) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;
  subplot(2,2,3) ;
    gpnpr(indat,invh,[.2,.5],0,0,0) ;
      title('NW - no trunc.') ;
      axis(axis) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;
  subplot(2,2,4) ;
    gpnpr(indat,invh,[.2,.5],0,0,1) ;
      title('NW - trunc.') ;
      axis(axis) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

elseif itest == 9 ;   %  Check kernel plot LL
  clf ;
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  [okde,oxgrid,omker] = gpnpr(indat,[.01; .03; .3],[0,1],0,1,0) ;
  plot(oxgrid,[okde (omker - 1.5)]) ;
    axis(axis) ;
    hold on ;
      plot(indat(:,1),indat(:,2),'.w') ;
    hold off ;

elseif itest == 10 ;   %  Check kernel plot NW
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  [okde,oxgrid,omker] = gpnpr(indat,[.01; .03; .3],[0,1],0,0,0) ;
  plot(oxgrid,[okde (omker - 1.5)]) ;
    axis(axis) ;
    hold on ;
      plot(indat(:,1),indat(:,2),'.w') ;
    hold off ;

elseif itest == 11 ;   %  Check data based bandwidth selection
  indat = rand(100,1) ;
  indat = [indat, sin(2 * pi * indat) + .2*randn(100,1)] ;
  subplot(2,2,1) ;
    gpnpr(indat,-1) ;
      title('RSW ROT') ;
    hold on ;
      plot(indat(:,1),indat(:,2),'.w') ;
    hold off ;
  subplot(2,2,2) ;
    gpnpr(indat,-2) ;
      title('RSW DPI') ;
    hold on ;
      plot(indat(:,1),indat(:,2),'.w') ;
    hold off ;
 
elseif itest == 12 ;    %  Check small h part, LL

  indat = [1,2; 2,4; 2.5,2; 4.1,1; 4.2,5; 5.5, 4] ;
%  h = .1 ;
%  h = .05 ;
%  h = .03 ;
%  h = .01 ;
  h = .001 ;
  invxg = [0; 6; 31] ;
%  invxg = [0; 6; 101] ;

%Tried this one, but the replaced it by polyfit
%  titstr = 'Small h, prebinned LL'
%  subplot(2,2,1) ;
%    inbcts = gplbinr(indat,invxg) ;
%    gpnpr(inbcts,h,invxg,-1) ;
%      title(titstr) ;
%      hold on ;
%        plot(indat(:,1),indat(:,2),'or') ;
%      hold off ;

  titstr = 'Small h, Poly Fit LL'
  subplot(2,2,1) ;
    gpnpr(indat,h,invxg,3) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Binned LL'
  subplot(2,2,2) ;
    gpnpr(indat,h,invxg) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Matrix LL'
  subplot(2,2,3) ;
    gpnpr(indat,h,invxg,1) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Looped LL'
  subplot(2,2,4) ;
    gpnpr(indat,h,invxg,2) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

elseif itest == 13 ;    %  Check small h part, NW

  indat = [1,0; 2,4; 2.5,2; 4.5,1; 5,5] ;
%  h = .1 ;
%  h = .03 ;
  h = .01 ;
  invxg = [0; 6; 31] ;

  titstr = 'Small h, prebinned NW'
  subplot(2,2,1) ;
    inbcts = gplbinr(indat,invxg) ;
    gpnpr(inbcts,h,invxg,-1,0) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Binned NW'
  subplot(2,2,2) ;
    gpnpr(indat,h,invxg,0,0) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Matrix NW'
  subplot(2,2,3) ;
    gpnpr(indat,h,invxg,1,0) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

  titstr = 'Small h, Looped NW'
  subplot(2,2,4) ;
    gpnpr(indat,h,invxg,2,0) ;
      title(titstr) ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;

elseif itest == 14 ;    %  Test widely spaced data
  titstr = 'Widely Spaced Test'
  indat = [rand(20,1); (10 + rand(20,1)); (20 + rand(20,1))] ;
  indat = [indat, abs(indat - 10.5)] ;
  vh = [.1; .3; 1; 3; 10] ;

  subplot(1,1,1) ;
    gpnpr(indat,vh) ;
      title(titstr) ;
      legend('h=0.1','h=0.3','h=1','h=3','h=10') ;
      hold on ;
        plot(indat(:,1),indat(:,2),'or') ;
      hold off ;
      text(5,10,'Try Command: axis([9.5,11.5,0,1])') ;
   
end ;


