%MPTEST
% This script contains various tests of the mp*.m files
%  Note that MPDER, MPSHIFT, MPNEXT  are tested incidentally,
%  with the major stress on MPVAL, MPMAK, MPAPI.

% cb aug97
% try trivial stuff first

nn=0;
nn=nn+1
% univariate polynomial  p(x) = 1-x;
mp = mpmak([1 -1],1);
if(any(mpval(mp,[-1:2])-[2:-1:-1]))
   fprintf('error in test1 in MVTEST.\n'), end
if any(mpbrk(mpder(mp,1),'co')-(-1))
   fprintf('error in test1der in MVTEST.\n'), end

nn=nn+1
% univariate polynomial  p(x) = (1-x)^2;
mp = mpmak([1 -2  1],1);
if(any(mpval(mp,[-1:2])-(1-[-1:2]).^2))
   fprintf('error in test2 in MVTEST.\n'), end

nn=nn+1
% univariate polynomial  p(x) = (1-x)^3;
mp = mpmak([1 -3 3 -1],1);
if(any(mpval(mp,[-1:2])-(1-[-1:2]).^3))
   fprintf('error in test3 in MVTEST.\n'), end
if any(mpbrk(mpder(mp,-2),'co')-[6 -12 6])
   fprintf('error in test3der in MVTEST.\n'), end

nn=nn+1
% a bivariate polynomial  p(x,y) = (1-x)^2;
mp = mpmak([1 0 -2 0 0 1],2);
if(any(mpval(mp,dcube(2))-[0 0 4 4]))
   fprintf('error in test4 in MVTEST.\n'), end

nn=nn+1
% a trivariate quadratic that vanishes at all the corners of the cube
mp = mpmak([-3 0 0 0 1 0 1 0 0 1],3);
if(any(mpval(mp,dcube(3))))
   fprintf('error in test5 in MVTEST.\n'), end
p = [1;-2;3]; xi = [2;1;-2]; h = 1.e-7;
if abs(mpval(mp,[p+h*xi p-h*xi])*[1;-1]/(2*h) - mpval(mpder(mp,xi),p))>h
   fprintf('error in test5der in MVTEST.\n'), end

nn=nn+1
% finally, check it out when there is less symmetry
% i.e., the polynomial (x-1)(y-1)(z-1) which vanishes at all but one corner
mp = mpmak([-1 1 1 1 0 -.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0],3);
if(abs(mpval(mp,dcube(3))-[0 0 0 0 0 0 0 -8])>1.e-15)
   fprintf('error in test6 in MVTEST.\n'), end
p = [1;-2;3]; xi = [2;1;-2]; h = 1.e-5;
if abs(mpval(mp,[p+h*xi p p-h*xi])*[1;-2;1]/(h*h) - ...
                                   mpval(mpder(mp,[xi xi]),p))>10*h
   fprintf('error in test6der in MVTEST.\n'), end
if mpbrk(mpder(mp,[xi -2*xi xi -xi 3*xi]),'co')
   fprintf('error in test6ader in MVTEST.\n'), end
% Now check the vector-valued calculation:
%  (x-1)(y-1)(z-1)
%  (x+1)(y-1)(z-1)
%  (x-1)(y+1)(z-1)
%  (x-1)(y-1)(z+1)

nn=nn+1
[-1 +1 +1 +1 0 -.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0
  1 -1 -1 +1 0 +.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0
  1 -1 +1 -1 0 -.5 0 +.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0
  1 +1 -1 -1 0 -.5 0 -.5 +.5 0 0 0 0 0 0 1/6 0 0 0 0];
mp = mpmak(ans,3);
[0 0 0 0 0 0 0 -8
 0 0 0 8 0 0 0 0
 0 0 0 0 0 8 0 0
 0 0 0 0 0 0 8 0];
if(any(abs(mpval(mp,dcube(3))-ans)>1.e-15))
   fprintf('error in test7 in MVTEST.\n'), end

nn=nn+1

nn=nn+1

nn=nn+1
% Now we are ready to tackle the least interpolant...
% first, the simplest case
mp = mpapi(0,0);
if(mpval(mp,0))
   fprintf('error in test8 in MVTEST.\n'), end

nn=nn+1
% try bivariate, linear interpolation
points = [0 1 0;0 0 1]; values = [0 1 1];
mp = mpapi(points, values);
if(mpval(mp,points) - [0 1 1])
   fprintf('error in test9 in MVTEST.\n'), end

nn=nn+1
% try bivariate, quadratic, still simple
points = [0 1 2 0 1 0; 0 0 0 1 1 2];
coefs = [1 -1 1 -1 1 -1];
mp0 = mpmak(coefs,2);
if any(abs(mpbrk(mpapi(points,mpval(mp0,points)),'co')-coefs)>1.e-14)
   fprintf('error in test10 in MVTEST.\n'), end

nn=nn+1
% try bivariate, cubic, still simple
points = [0 1 2 3 0 1 2 0 1 0;0 0 0 0 1 1 1 2 2 3];
coefs = [1 -1 1 -1 1 -1 1 -1 1 -1];
mp0 = mpmak(coefs,2);
if any(abs(mpbrk(mpapi(points,mpval(mp0,points)),'co')-coefs)>1.e-14)
   fprintf('error in test11 in MVTEST.\n'), end

nn=nn+1
% try trivariate interpolation to a linear function
points = [0 1 0 0 1; 0 0 1 0 1; 0 0 0 1 1];
coefs = [1 -2 3 -1];
mp0 = mpmak(coefs,3);
values = mpval(mp0, points);
if any(abs(mpbrk(mpapi(points, values, 1.e-14, [0;0;0]),'co')- ...
                                              [coefs 0 0 0 0 0 0])>1.e-14)
   fprintf('error in test12 in MVTEST.\n'), end

nn=nn+1
% try trivariate, quadratic, interpolation to a quadratic function
points = [0 1 2 0 1 0 0 1 0 0;0 0 0 1 1 2 0 0 1 0;0 0 0 0 0 0 1 1 1 2];
coefs = [1 -1 1 -1 1 -1 1 -1 1 -1];
mp0 = mpmak(coefs,3);
values = mpval(mp0, points);
mp = mpapi(points, values);
mps = mpshift(mp,[0;0;0]);
if any(abs(mpbrk(mps,'co')-coefs)>1.e-14)
   fprintf('error in test13 in MVTEST.\n'), end
if any(abs(mpbrk(mp,'co') - mpbrk(mpshift(mps,mpbrk(mp,'ce')),'co'))>1.e-14)
   fprintf('error in test13shift in MVTEST.\n'), end
 
% try trivariate, cubic, interpolation to a quadratic function
points = dcube(3);
coefs = [1 -1 1 -1 1 -1 1 -1 1 -1];
mp0 = mpmak(coefs,3);
values = mpval(mp0, points);
% In this example, the coefficients are not reproduced, but the interpolant is
% in Pi_2. To check it
cc = mpbrk(mpapi(points, values),'co');
cc(1:10) = cc(1:10)-coefs;
mpd = mpmak(cc,3);
if any(abs(mpval(mpd,points))>1.e-14)
   fprintf('error in test14 in MVTEST.\n'), end

%   re-use not implemented for v4
%nn=nn+1
%% test the re-use of factorization information:
%points = rand(3,10); values = [1 -2 1]*(points.*(points - 1));
%[mp, mpfactor] = mpapi(points,values,1.e-12,[0;0;0]);
%mpf = mpapi(mpfactor,values);
%if (max(abs(mpbrk(mp,'co')-mpbrk(mpf,'co'))))  
%   fprintf('error in test15 in MVTEST.\n'), end

nn=nn+1
% test the vector capability:
points = rand(3,10); values = (points.*(points - 1));
mp = mpapi(points,values);
if any(abs(values - mpval(mp,points))>1.e-14)
   fprintf('error in test16 in MVTEST.\n'), end

nn=nn+1
% test the noise level
fprintf(['The next test checks interpolation at 10 data sites along a',...
     ' straight line in R^3 to quartic values.\n',...
  'The interpolant should formally be of degree 9 but its coefficients\n',...
  'of degree > 4 should be zero. Also, the derivative in any direction\n',...
  'perpendicular to that straight line should be zero.\n',...
'The next three numbers give the size of the supposedly zero coefficients\n',...
'relative to the size of the coefficients of the interpolant.\n'])
t = 6*rand(1,10)-3; xi = [1;2;3]; points = xi*t; values = t.^4;
mp = mpapi(points,values);
% coefficients of degree > 4 should be zero:
coefs = mpbrk(mp,'co'); dimpjd = mpbrk(mp,'p');
mc = max(abs(coefs)); mc5 = max(abs(coefs(dimpjd(6)+1:dimpjd(11))));
fprintf('noise ratio = %8g \n', mc5/mc)
% Also, the interpolant should be constant in directions orthogonal to  xi :
noise1 = max(abs(mpbrk(mpder(mp,[1;-2;1]),'co')));
noise2 = max(abs(mpbrk(mpder(mp,[3;0;-1]),'co')));
fprintf('derivative noise = %8g, %8g \n',noise1/mc, noise2/mc)

% nn=nn+1
% % Now add a speed test, to compare the old bivariate routine with the present
% % one, on a large number of points:
% nn = 8;
% for n=1:4
%    nn = nn*2
%    points = rand(2,nn); values = mpval(mpmak([1 -1 1 -1 1 -1],2),points);
%    tic
%    mp = mpapi(points,values,1.e-8);
%    toc
%    tic
%    bp = bpapi(points,values);
%    toc
% end
