test_mpoly_numerical.m

几个关于多小波的程序

% TEST_MPOLY_NUMERICAL -- test the subroutines in the @mpoly directory for numerical matrices

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

disp('testing routines in mw/@mpoly for numerical arguments');
disp(' ');

% test mpoly
P2 = [112,122;212,222];
P3 = [113,123;213,223];
P23 = cat(3,P2,P3);
P = mpoly({P2,P3},2);
Q = mpoly({P2,P3});
PP = mpoly(P23,2,'polyphase',2,1);
P11 = mpoly({112,113},2);
disp('OK - mpoly');

% test display
P
PP
disp('OK - display');

% test trim

Pt = trim(P,[3,4]);
if (size(Pt.coef,3) ~= 2 | Pt.coef ~= cat(3,P3,zeros(2)))
    error('failed - trim');
end

Ppm = (eye(2) + P) - eye(2);
if (size(Ppm.coef,3) ~= 2 | Ppm.min ~= P.min)
    error('failed - trim');
end
disp('OK - trim');

% test comparisons
if ((P ~= PP) | (P == Q))
    error('failed - comparisons: eq, ne');
end

% test conversion
H = symbol(PP);
if (polyphase(H) ~= PP)
    error('failed - conversion: symbol, polyphase');
end

I = mpoly(eye(2));
Z = mpoly(zeros(2));
Iz = mpoly(eye(2),1);
SI = mpoly(sym(eye(2)));
syms a b c d
S = mpoly([a,b;c,d]);
if (sym(I) ~= SI | double(SI) ~= I)
    error('failed - conversion: sym, double');
end

if (reshape(P,1,4) ~= mpoly(reshape(P.coef,1,4,2),2))
    error('failed - reshape');
end

disp('OK - conversions: symbol, polyphase, sym, double, reshape');

% test more comparisons
if (~isconstant(I) | isconstant(Iz))
    error('failed - isconstant');
end

if (~isidentity(I) | isidentity(Iz))
    error('failed - isidentity');
end

if (~ismonomial(Iz) | ismonomial(P))
    error('failed - ismonomial');
end

if (~isnumeric(I) | isnumeric(S))
    error('failed - isnumeric');
end

if (issymbolic(I) | issymbolic(SI) | ~issymbolic(S))
    error('failed - issymbolic');
end

if (iszero(I) | ~iszero(Z))
    error('failed - iszero');
end

disp('OK - comparisons: eq, ne, isconstant, isidentity');
disp('                  ismonomial, isnumeric, issymbolic, iszero');

% test concatenation
z = mpoly(1,1);
Q = [z,0;0,z];
if (Q ~= z * eye(2))
    error('failed - concatenation');
end
disp('OK - concatenation: horzcat, vertcat');

% test right-hand subscripting
if (P.min ~= 2)
    error('failed - subscripting P.min');
end

if (P.max ~= 3)
    error('failed - subscripting P.max');
end

if (P.coef ~= P23)
    error('failed - subscripting P.coef');
end

if (P.length ~= 2)
    error('failed - subscripting P.length');
end

if (P.degree ~= 1)
    error('failed - subscripting P.degree');
end

if (~iszero(P.size - [2,2]))
    error('failed - subscripting P.size');
end

if (~strcmp(PP.type,'polyphase'))
    error('failed - subscripting P.type');
end

if (PP.m ~= 2)
    error('failed - subscripting P.m');
end

if (PP.r ~= 1)
    error('failed - subscripting P.r');
end

if (~iszero(P{end} - P3))
    error('failed - subscripting P{i}');
end

if (~iszero(P{4} - zeros(2)))
    error('failed - subscripting P{i}');
end

if (P(1,1) ~= P11)
    error('failed - subscripting P(i,j)');
end

disp('OK - right-hand subscripting: get, P.field, P(i,j), P{i}')

% test left-hand subscripting

Q = P;
Q.min = Q.min + 1;
if (Q ~= P*z)
    error('failed - left-hand subscripting P.min');
end

Q = P;
Q{2} = zeros(2);
if (Q ~= mpoly(P3,3))
    error('failed - left-hand subscripting P{i}');
end

Q = P;
Q(1,1) = z^2;
if (Q(1,1) ~= z^2)
    error('failed - left-hand subscripting P(i,j)');
end
disp('OK - left-hand subscripting: set, P.field, P(i,j), P{i}')

% test size
if (iszero(P.size - size(P)))
    disp('OK - size')
else
    error('failed - size');
end

% test addition/subtraction/scalar multiplication

if (P ~= +P)
    error('failed - unary plus');
end

if (P + P ~= 2*P)
    error('failed - addition or scalar multiplication');
end

if (~iszero(-P + P))
    error('failed - unary minus');
end

if (P + Iz - P ~= Iz)
    error('failed - addition or subtraction');
end
disp('OK - addition, subtraction, scalar multiplication');

% test elementwise multiplication/division/power

Two = 2 * mpoly(ones(2,2,2),2);

if (P .* Two ~= P * 2)
    error('failed - times');
end

if (Two .\ P ~= 2 \ P)
    error('failed - ldivide');
end

if (P ./ Two ~= P / 2)
    error('failed - rdivide');
end

if (P .^ 2 ~= mpoly(P23.^2,2))
    error('failed - power');
end
disp('OK - pointwise operations: times, ldivide, rdivide, power');

% test multiplication/division/power

% matrix with nonsingular leading term
L = eye(2);
L(1,1,2) = 1;
L = mpoly(L);

PL = P*L;
[Q,method] = mrdivide(PL,L);
if (Q ~= P)
    error(['failed - right multiplication or division, method ',num2str(method)]);
end
if (method ~= 1)
    warning('this should have been method 1');
end

LP = L*P;
[Q,method] = mldivide(L,LP);
if (Q ~= P)
    error(['failed - left multiplication or division, method ',num2str(method)]);
end
if (method ~= 1)
    warning('this should have been method 1');
end

% matrix with nonsingular highest term
H = [1,0;0,0];
H(:,:,2) = eye(2);
H = mpoly(H);

PH = P*H;
[Q,method] = mrdivide(PH,H);
if (Q ~= P)
    error(['failed - right multiplication or division, method ',num2str(method)]);
end
if (method ~= 1)
    warning('this should have been method 1');
end

HP = H*P;
[Q,method] = mldivide(H,HP);
if (Q ~= P)
    error(['failed - left multiplication or division, method ',num2str(method)]);
end
if (method ~= 1)
    warning('this should have been method 1');
end

% projection factor
U = [3/5;4/5];
F = mpoly({eye(size(U,1)) - U*U',U*U'},0,'polyphase',2,1);

PF = P*F;
[Q,method] = mrdivide(PF,F);
if (Q ~= P)
    error(['failed - right multiplication or division, method ',num2str(method)]);
end
if (method ~= 2)
    warning('this should have been method 2');
end

FP = F*P;
[Q,method] = mldivide(F,FP);
if (Q ~= P)
    error(['failed - left multiplication or division, method ',num2str(method)]);
end
if (method ~= 2)
    warning('this should have been method 2');
end

% upper triangular matrix
% I have to use a 3x3 example to force method 3
U = [1,1,1;0,1,1;0,0,0];
U(:,:,2) = [0,0,0;0,1,0;0,0,1];
U = mpoly(U);

R = [112,122,132;212,222,232;312,322,332];
R(:,:,2) = [113,123,133;213,223,233;313,323,333];
R = mpoly(R,2);

RU = R*U;
[Q,method] = mrdivide(RU,U);
if (Q ~= R)
    error(['failed - right multiplication or division, method ',num2str(method)]);
end
if (method ~= 3)
    warning('this should have been method 3');
end

UR = U*R;
[Q,method] = mldivide(U,UR);
if (Q ~= R)
    error(['failed - left multiplication or division, method ',num2str(method)]);
end
if (method ~= 3)
    warning('this should have been method 3');
end

% lower triangular matrix
% I have to use a 3x3 example to force method 3
U = [1,1,1;0,1,1;0,0,0]';
U(:,:,2) = [0,0,0;0,1,0;0,0,1]';
U = mpoly(U);

R = [112,122,132;212,222,232;312,322,332];
R(:,:,2) = [113,123,133;213,223,233;313,323,333];
R = mpoly(R,2);

RU = R*U;
[Q,method] = mrdivide(RU,U);
if (Q ~= R)
    error(['failed - right multiplication or division, method ',num2str(method)]);
end
if (method ~= 3)
    warning('this should have been method 3');
end

UR = U*R;
[Q,method] = mldivide(U,UR);
if (Q ~= R)
    error(['failed - left multiplication or division, method ',num2str(method)]);
end
if (method ~= 3)
    warning('this should have been method 3');
end

U = [1,1;0,0];
U(:,:,2) = [0,0;0,1];
U = mpoly(U,1);

Uinv = inv(U);
if (Uinv * U ~= eye(2))
    error('failed - inv');
end

Uinv = longinv(U);
if (Uinv * U ~= eye(2))
    error('failed - longinv');
end

if (P^2 ~= P*P)
    error('failed - mpower');
end

disp('OK - multiplication/division: mltimes, mldivide, mrtimes');
disp('                              mrdivide, mpower, inv, longinv');

% test transpose

Pt = mpoly({P3',P2'},-3);
Prt = mpoly({P2',P3'},2);

if (P' ~= Pt)
    error('failed - ctranspose');
end

if (P.' ~= Prt)
    error('failed - transpose');
end

if ((P').' ~= reverse(P))
    error('failed - reverse');
end

disp('OK - transpose, ctranspose, reverse');

% test rounding

if (ceil(P - 0.1) ~= P)
    error('failed - ceil');
end

if (floor(P + 0.1) ~= P)
    error('failed - floor');
end

if (fix(P * 1.000001) ~= P)
    error('failed - fix');
end

if (round(P * 1.000001) ~= P)
    error('failed - round');
end

disp('OK - rounding: ceil, floor, round, fix');

% test diagonal/triangle

if (P ~= (diag(P) + tril(P,-1) + triu(P,1)))
    error('failed - diag, tril or triu');
end

disp('OK - diag, tril, triu');

% test determinant

if (det(F) ~= z)
    error('failed - det')
end

disp('OK - det');

% test moment
if (moment(P) ~= P2+P3 | moment(P,1) ~= 2*P2+3*P3)
    error('failed - moment');
end
disp('OK - moment');

% test match_type
Q = P * PP;
R = PP * symbol(PP);
if (~strcmp(Q.type,'polyphase') | ~strcmp(R.type,''))
    error('failed - match_type');
end
disp('OK - match_type');

% test kron
Q = kron(P,P);
if (Q(3:4,1:2) ~= P(2,1)*P)
    error('failed - kron');
end
disp('OK - kron');