testDiscrete.m

function [tests] = testDiscrete()
%testQuantity Summary of this function goes here
%   Detailed explanation goes here
tests = functiontests(localfunctions);
end


function setupOnce(testCase)
syms z zeta t sy
testCase.TestData.sym.z = z;
testCase.TestData.sym.zeta = zeta;
testCase.TestData.sym.t = t;
testCase.TestData.sym.sy = sy;
end

function testHashData(testCase)

	spatialDomain = quantity.Domain("z", linspace(0, 1, 101));
	L1 = quantity.Symbolic(diag([2, -3]), ...
					'domain', spatialDomain, 'name', 'Lambda');

	L2 = quantity.Symbolic(diag([1, -1]), ...
					'domain', spatialDomain, 'name', 'Lambda');
	
	testCase.verifyNotEqual(L1.hash, L2.hash);
				
end

function testIsdiag(tc)
z = quantity.Domain("z", linspace(0, 1, 11));

myOnes = quantity.Discrete(ones(11, 4, 4), 'domain', z);
tc.verifyFalse(isdiag(myOnes));

myScalar = quantity.Discrete(ones(11, 1), 'domain', z);
myScalar2 = quantity.Discrete(zeros(11, 1), 'domain', z);
myScalar3 = quantity.Discrete(z.grid, 'domain', z);
tc.verifyTrue(isdiag(myScalar));
tc.verifyTrue(isdiag(myScalar2));
tc.verifyTrue(isdiag(myScalar3));

myEye3x3 = 2*quantity.Discrete(repmat(reshape(eye(3), [1, 3, 3]), [z.n, 1]), 'domain', z);
tc.verifyTrue(isdiag(myEye3x3));
myEye3x2 = 3*quantity.Discrete(repmat(reshape(eye(3, 2), [1, 3, 2]), [z.n, 1]), 'domain', z);
tc.verifyTrue(isdiag(myEye3x2));
end % testIsdiag(tc)

function testL2Norm(tc)

Z = tc.TestData.sym.z;
T = tc.TestData.sym.t;
ZETA = tc.TestData.sym.zeta;

t = quantity.Domain("t", linspace(0, 2, 101));
z = quantity.Domain("z", linspace(0, 1, 51));
x = quantity.Symbolic([Z * T; T], 'domain', [z, t]);

xNorm = sqrt(int(x.' * x, "z"));
tc.verifyEqual(MAX(abs(xNorm - x.l2norm('z'))), 0, 'AbsTol', 1e-12);

x = quantity.Discrete(x);
xNorm = sqrt(int(x.' * x, 'z'));
tc.verifyEqual(MAX(abs(xNorm - x.l2norm('z'))), 0, 'AbsTol', 1e-12);

tc.verifyTrue( l2norm(x) == l2norm(x, {'z', 't'}) );

end % testL2Norm()

function testCompose(testCase)

t = quantity.Domain('t', linspace(0, 2));
f = quantity.Discrete(sin(t.grid * 2 * pi) , 'domain', t, 'name', 'f' );
g = quantity.Discrete( 0.5 * t.grid, 'domain', t, 'name', 'g');
h = quantity.Discrete( 2 * t.grid, 'domain', t, 'name', 'h');

% verify simple evaluation of f( g(t) )
testCase.verifyEqual( on( f.compose(g) ), sin(t.grid * pi), 'AbsTol', 3e-3 )

% verify that a warning is thrown, if the domain of definition is exceeded
testCase.verifyWarning(@()f.compose(h), 'quantity:Discrete:compose')

% verify f(tau) = sin( tau * 2pi ); g(z, t) = t + z^2
% f( g(z, t ) ) = sin( (t + z^2 ) * 2pi )
tau = quantity.Domain('tau', linspace(-1, 3));
z = quantity.Domain('z', linspace(0, 1, 51));

f = quantity.Discrete( sin( tau.grid * 2 * pi ), 'domain', tau );
g = quantity.Discrete( ( z.grid.^2 ) + t.grid', 'domain', [z, t] );

fg = f.compose(g);
F = quantity.Function(@(z,t) sin( ( t + z.^2 ) * 2 * pi ), ...
	'domain', [z, t], 'name', 'f(g) as function handle');
testCase.verifyEqual(fg.on, F.on, 'AbsTol', 1e-2)

% verify f(z, tau) °_tau g(z, t) -> f(z, g(z,t))
%	f(z, tau) = z + tau;
%	g(z, t) = z + t
%	f(z, g(z,t) ) = 2*z + t

z = quantity.Domain('z', linspace(0, 1, 11));
t = quantity.Domain('t', linspace(0, 1, 21));
tau = quantity.Domain('tau', linspace(0, 2, 31)); 

f = quantity.Discrete( z.grid + tau.grid', 'domain', [z, tau]);
g = quantity.Discrete( z.grid + t.grid', 'domain', [z, t]);

fog = f.compose(g, 'domain', tau);
FoG = quantity.Discrete( 2 * z.grid + t.grid', 'domain', [z, t]);
testCase.verifyEqual( fog.on, FoG.on, 'AbsTol', 5e-16)

% verify the matrix valued case
% #TODO: by this, an object array can be created where different quantities
% with different grids respectively empty quantities are in the same
% matrix.
M(1:2,1:2) = f;

Mog = M.compose(g, 'domain', tau);
MoG(1:2, 1:2) = FoG;

testCase.verifyEqual( MoG.on, MoG.on, 'AbsTol', 5e-16)

% verify the case if one of the domains is a quantity.EquidistantDomain
z = quantity.EquidistantDomain("z", 0, 1, 'stepNumber', 11);
tau = quantity.EquidistantDomain("tau", 0, 2, 'stepNumber', 31);
f = quantity.Discrete( z.grid + tau.grid', 'domain', [z, tau]);

fog = f.compose(g, 'domain', tau);
testCase.verifyEqual( fog.on, FoG.on, 'AbsTol', 5e-16)

end

function testCastDiscrete2Function(testCase)

z = linspace(0, 2*pi)';
d = quantity.Discrete({sin(z); cos(z)}, 'grid', z, 'gridName', 'z');
f = quantity.Function(d);

testCase.verifyTrue(all(size(f) == size(d)));

end

function testQuadraticNorm(tc)
blub = quantity.Discrete(ones(11, 4), 'grid', linspace(0, 1, 11), ...
	'gridName', 'z', 'name', 'b');
blubQuadraticNorm = blub.quadraticNorm();
blubQuadraticWeightedNorm = blub.quadraticNorm('weight', 4*eye(4));

tc.verifyEqual(blubQuadraticNorm.on(), 2*ones(11,1));
tc.verifyEqual(blubQuadraticWeightedNorm.on(), 4*ones(11,1));
end % testQuadraticNorm()

function testSort(tc)

t = linspace(0, pi, 7)';
domA = quantity.Domain('a', t);
domB = quantity.Domain('b', t);
q = quantity.Discrete(sin(t) .* cos(t'), ...
	'domain', [domA, domB]);

q1V = q.valueDiscrete;

q.sort('descend');

q2V = q.valueDiscrete;

tc.verifyEqual(q1V , q2V')


end % testSort

function testBlkdiag(tc)
z = quantity.Domain('z', linspace(0, 1, 21));
zeta = quantity.Domain('zeta', linspace(0, 1, 41));

A = quantity.Discrete(cat(4, ...
	cat(3, 1 + z.grid .* zeta.grid', - z.ones .* zeta.grid'), ...
	cat(3, -z.grid .* zeta.ones', z.grid .^2 .* zeta.ones') ), ...
	'domain', [z, zeta], 'name', "q");

B = 2*A(:, 1);
C = 3*A(1);

% 1 input
tc.verifyEqual(A.on, A.blkdiag.on());

% 2 x same input
AA = blkdiag(A, A);
tc.verifyEqual(AA(1:2, 1:2).on(), A.on());
tc.verifyEqual(AA(3:4, 3:4).on(), A.on());
tc.verifyEqual(AA(1:2, 3:4).on(), 0*A.on());
tc.verifyEqual(AA(3:4, 1:2).on(), 0*A.on());

% 3 different sized quantites
ABCB = blkdiag(A, B, C, B);
tc.verifyEqual(ABCB(1:2, 1:2).on(), A.on());
tc.verifyEqual(ABCB(3:4, 3).on(), B.on());
tc.verifyEqual(ABCB(5, 4).on(), C.on());
tc.verifyEqual(ABCB(6:7, 5).on(), B.on());

zeroElements = ~logical(blkdiag(true(size(A)), true(size(B)), true(size(C)), true(size(B))));
tc.verifyEqual(ABCB(zeroElements(:)).on(), 0*ABCB(zeroElements(:)).on());

end % testBlkdiag()

function testSetName(tc)
blub = quantity.Discrete(ones(11, 2), 'domain', quantity.Domain("z", linspace(0, 1, 11)));
blub.setName("asdf");
tc.verifyEqual(blub(1).name, "asdf");
tc.verifyEqual(blub(2).name, "asdf");
end

function testCtranspose(tc)
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

qSymbolic = quantity.Symbolic(...
	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
	'gridName', {'z', 'zeta'}, 'name', "q");
qDiscrete = quantity.Discrete(qSymbolic);
qDiscreteCtransp = qDiscrete';
tc.verifyEqual(qDiscrete(1,1).on(), qDiscreteCtransp(1,1).on());
tc.verifyEqual(qDiscrete(2,2).on(), qDiscreteCtransp(2,2).on());
tc.verifyEqual(qDiscrete(1,2).on(), qDiscreteCtransp(2,1).on());
tc.verifyEqual(qDiscrete(2,1).on(), qDiscreteCtransp(1,2).on());

qDiscrete2 = quantity.Discrete(ones(11, 1) + 2i, ...
	'grid', linspace(0, 1, 11), 'gridName', "z", 'name', "im");
qDiscrete2Ctransp = qDiscrete2';
tc.verifyEqual(qDiscrete2.real.on(), qDiscrete2Ctransp.real.on());
tc.verifyEqual(qDiscrete2.imag.on(), -qDiscrete2Ctransp.imag.on());
end % testCtranspose

function testTranspose(tc)
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

qSymbolic = quantity.Symbolic(...
	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
	'gridName', {'z', 'zeta'}, 'name', "q");
qDiscrete = quantity.Discrete(qSymbolic);
qDiscreteTransp = qDiscrete.';

tc.verifyEqual(qDiscrete(1,1).on(), qDiscreteTransp(1,1).on());
tc.verifyEqual(qDiscrete(2,2).on(), qDiscreteTransp(2,2).on());
tc.verifyEqual(qDiscrete(1,2).on(), qDiscreteTransp(2,1).on());
tc.verifyEqual(qDiscrete(2,1).on(), qDiscreteTransp(1,2).on());
end % testTranspose

function testFlipGrid(tc)
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

myGrid = linspace(0, 1, 11);
f = quantity.Discrete(quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, ...
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
% % flip one grid
fReference = quantity.Discrete(quantity.Symbolic([1+z+(1-zeta); 2*(1-zeta)+sin(z)] + zeros(2, 1)*z*zeta, ...
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
fFlipped = f.flipGrid('zeta');
tc.verifyEqual(fReference.on(), fFlipped.on(), 'AbsTol', 10*eps);

% flip both grids
fReference2 = quantity.Discrete(quantity.Symbolic([1+(1-z)+(1-zeta); 2*(1-zeta)+sin(1-z)] + zeros(2, 1)*z*zeta, ...
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
fFlipped2 = f.flipGrid({'z', 'zeta'});
fFlipped3 = f.flipGrid({'zeta', 'z'});
tc.verifyEqual(fReference2.on(), fFlipped2.on(), 'AbsTol', 10*eps);
tc.verifyEqual(fReference2.on(), fFlipped3.on(), 'AbsTol', 10*eps);
end % testFlipGrid();

function testScalarPlusMinusQuantity(testCase)

z = quantity.Domain('z', linspace(0, 1, 7));
f = quantity.Discrete( ( [2; 1] + zeros(2, 1) .* z.grid' )', ...
	'domain', z);

testCase.verifyError(@() 1-f-1, '');
testCase.verifyError(@() 1+f+1, '');

end % testScalarPlusMinusQuantity

function testNumericVectorPlusMinusQuantity(testCase)

z = quantity.Domain('z', linspace(0, 1, 7));
f = quantity.Discrete( [1+z.grid, 2+sin(z.grid)] , ...
	'domain', z);

a = ones(size(f));
testCase.verifyEqual(on(a-f), 1-f.on());
testCase.verifyEqual(on(f-a), f.on()-1);
testCase.verifyEqual(on(a+f), f.on()+1);
testCase.verifyEqual(on(f+a), f.on()+1);
end % testNumericVectorPlusMinusQuantity

function testUnitaryPluasAndMinus(testCase)
z = quantity.Domain('z', linspace(0, 1, 7));
zeta = quantity.Domain('zeta', linspace(0, 1, 7));

q = quantity.Discrete(cat(4, ...
	cat(3, 1 + z.grid .* zeta.grid', - z.grid .* zeta.ones'), ...
	cat(3, - z.ones .* zeta.grid', z.grid.^2 .* zeta.ones') ), ...
	'domain', [z, zeta]) ;

qDoubleArray = q.on();
	
testCase.verifyEqual(on(-q), -qDoubleArray);
testCase.verifyEqual(on(+q), +qDoubleArray);

end

function testConcatenate(testCase)

t = linspace(0, pi, 7)';
A = quantity.Discrete({sin(t); cos(t)}, 'grid', {t}, 'gridName', 't');
B = quantity.Discrete({tan(t); exp(t)}, 'grid', {t}, 'gridName', 't');

AB = [A, B];
AB_ = [A', B'];
ABA = [A, B, A];

testCase.verifyTrue(all(size(AB) == [2, 2]));
testCase.verifyTrue(all(size(AB_) == [1, 4]));
testCase.verifyTrue(all(size(ABA) == [2, 3]));

t = linspace(0, pi, 13)';
C = quantity.Discrete({sin(t); cos(t)}, 'grid', {t}, 'gridName', 't');
AC = [A; C];

testCase.verifyTrue(all(size(AC) == [4, 1]));

A0s = [A, zeros(2,3)];
testCase.verifyTrue(all(all(all(A0s(:, 2:end).on() == 0))))

O = quantity.Discrete.empty([5, 0]);
O_horz = [O, O];
O_vert = [O; O];

testCase.verifyEqual(size(O_horz, 1), 5);
testCase.verifyEqual(size(O_vert, 1), 10);

z = linspace(0, 1, 5);
D = quantity.Discrete({sin(t * z); cos(t * z)}, 'grid', {t, z}, 'gridName', {'t', 'z'});
E = quantity.Discrete({tan(z' * t'); cos(z' * t')}, 'grid', {z, t}, 'gridName', {'z', 't'});

DE = [D, E];
compareDE = quantity.Discrete({sin(t * z), tan(t*z); cos(t * z), cos(t*z)}, 'grid', {t, z}, 'gridName', {'t', 'z'});

testCase.verifyEqual(DE.on(), compareDE.on())

ED = [E, D];
compareED = quantity.Discrete({tan(t*z), sin(t * z); cos(t * z), cos(t*z)}, 'grid', {t, z}, 'gridName', {'t', 'z'});
testCase.verifyEqual(ED.on(), compareED.on())

% concatenation of a constant and a function:
C = [5; 2];
const = quantity.Discrete(C, 'domain', quantity.Domain.empty);
% 
constA = [const, A];
testCase.verifyTrue(all(constA(1,1).on() == C(1)));
testCase.verifyTrue(all(constA(2,1).on() == C(2)));

Aconst = [A, const];
testCase.verifyTrue(all(Aconst(1,2).on() == C(1)));
testCase.verifyTrue(all(Aconst(2,2).on() == C(2)));


end

function testExp(testCase)
% 1 spatial variable
z = quantity.Domain('z', linspace(0, 1, 3));
s1d = quantity.Discrete(1 + z.grid .* z.grid.', 'domain', z);
testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));

% diagonal matrix
zeta = quantity.Domain('zeta', linspace(0, 1, 4));
s2dDiag = quantity.Discrete(cat(4, ...
	cat(3, 1 + z.grid .* zeta.grid', zeros(z.n, zeta.n)), ....
	cat(3, zeros(z.n, zeta.n), z.grid.^2 .* zeta.ones')), 'domain', [z, zeta]);

testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
end

function testExpm(testCase)

z = quantity.Domain('z', linspace(0, 1, 3));
zeta = quantity.Domain('zeta', linspace(0, 1, 4));

mat2d = quantity.Discrete(...
	{1 + z.grid .* zeta.grid', 3 * z.ones .* zeta.grid'; ...
	 2 + 5 * z.grid + zeta.grid', z.grid .^2 .* zeta.ones'}, ...
	 'domain', [z, zeta]);

mat2dMat = mat2d.on();
mat2dExpm = 0 * mat2d.on();
for zIdx = 1 : z.n
	for zetaIdx = 1 : zeta.n
		mat2dExpm(zIdx, zetaIdx, :, :) = expm(reshape(mat2dMat(zIdx, zetaIdx, :, :), [2, 2]));
	end
end
testCase.verifyEqual(mat2d.expm.on(), mat2dExpm, 'RelTol', 100*eps);

end

function testSqrt(testCase)
% 1 spatial variable
quanScalar1d = quantity.Discrete((linspace(0,2).^2).', 'size', [1, 1], 'grid', linspace(0, 1), ...
	'gridName', 'z', 'name', "s1d");
quanScalarRoot = quanScalar1d.sqrt();
testCase.verifyEqual(quanScalarRoot.on(), linspace(0, 2).');

% 2 spatial variables
quanScalar2d = quantity.Discrete((linspace(0, 2, 21).' + linspace(0, 2, 41)).^2, 'size', [1, 1], ...
	'grid', {linspace(0, 1, 21), linspace(0, 1, 41)}, ...
	'gridName', {'z', 'zeta'}, 'name', 's2d');
quanScalar2dRoot = quanScalar2d.sqrt();
testCase.verifyEqual(quanScalar2dRoot.on(), (linspace(0, 2, 21).' + linspace(0, 2, 41)));

% diagonal matrix
quanDiag = [quanScalar1d, 0*quanScalar1d; 0*quanScalar1d, (2)*quanScalar1d];
quanDiagRoot = quanDiag.sqrt();
testCase.verifyEqual(quanDiagRoot(1,1).on(), quanScalarRoot.on());
testCase.verifyEqual(quanDiagRoot(1,2).on(), 0*quanDiagRoot(1,2).on());
testCase.verifyEqual(quanDiagRoot(2,1).on(), 0*quanDiagRoot(2,1).on());
testCase.verifyEqual(quanDiagRoot(2,2).on(), sqrt(2)*quanDiagRoot(1,1).on(), 'AbsTol', 10*eps);

%testCase.verifyEqual(1, 0);
end

function testSqrtm(testCase)
quanScalar1d = quantity.Discrete((linspace(0,2).^2).', 'size', [1, 1], 'grid', linspace(0, 1), ...
	'gridName', 'z', 'name', 's1d');
quanScalar2d = quantity.Discrete((linspace(0, 2, 21).' + linspace(0, 2, 41)).^2, 'size', [1, 1], ...
	'grid', {linspace(0, 1, 21), linspace(0, 1, 41)}, ...
	'gridName', {'z', 'zeta'}, 'name', 's2d');
% diagonal matrix - 1 spatial variable
quanDiag = [quanScalar1d, 0*quanScalar1d; 0*quanScalar1d, (2)*quanScalar1d];
quanDiagRoot = quanDiag.sqrt(); % only works for diagonal matrices
quanDiagRootMatrix = quanDiag.sqrtm();
testCase.verifyEqual(quanDiagRootMatrix.on(), quanDiagRoot.on());

% full matrix - 1 spatial variable
quanMat1d = [quanScalar1d, 4*quanScalar1d; quanScalar1d, (2)*quanScalar1d];
quanMat1dMat = quanMat1d.on();
quanMat1dReference = 0*quanMat1dMat;
for zIdx = 1:size(quanMat1dReference, 1)
	quanMat1dReference(zIdx, :, :) = sqrt(reshape(quanMat1dMat(zIdx, :), [2, 2]));
end
testCase.verifyEqual(quanMat1d.sqrtm.on(), quanMat1dReference, 'AbsTol', 10*eps);

% full matrix - 2 spatial variables
quanMat2d = [quanScalar2d, 4*quanScalar2d; quanScalar2d, (2)*quanScalar2d];
quanMat2dMat = quanMat2d.on();
quanMat2dReference = 0*quanMat2dMat;
for zIdx = 1:size(quanMat2dMat, 1)
	for zetaIdx = 1:size(quanMat2dMat, 2)
		quanMat2dReference(zIdx, zetaIdx, :, :) = ...
			sqrt(reshape(quanMat2dMat(zIdx, zetaIdx, :), [2, 2]));
	end
end
testCase.verifyEqual(quanMat2d.sqrtm.on(), quanMat2dReference, 'AbsTol', 10*eps);

end

function testDiag2Vec(testCase)
% quantity.Symbolic
z = testCase.TestData.sym.z;

do = quantity.Domain('z', linspace(0,1,3));
myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], 'domain', do);
myVectorSymbolic = diag2vec(myMatrixSymbolic);
testCase.verifyEqual(myMatrixSymbolic(1,1).valueContinuous, myVectorSymbolic(1,1).valueContinuous);
testCase.verifyEqual(myMatrixSymbolic(2,2).valueContinuous, myVectorSymbolic(2,1).valueContinuous);
testCase.verifyEqual(numel(myVectorSymbolic), size(myMatrixSymbolic, 1));

% quantity.Discrete
myMatrixDiscrete = quantity.Discrete(myMatrixSymbolic);
myVectorDiscrete = diag2vec(myMatrixDiscrete);
testCase.verifyEqual(myMatrixDiscrete(1,1).on(), myVectorDiscrete(1,1).on());
testCase.verifyEqual(myMatrixDiscrete(2,2).on(), myVectorDiscrete(2,1).on());
testCase.verifyEqual(numel(myVectorDiscrete), size(myMatrixDiscrete, 1));
end

function testVec2Diag(testCase)
% quantity.Discrete
n = 7;
z = linspace(0,1,n)';
myMatrixDiscrete = quantity.Discrete(...
	{sin(0.5*z*pi)+1, zeros(n,1); zeros(n,1), 0.9-z/2}, ...
	'grid', z, ...
	'gridName', 'z');
myVectorDiscrete = quantity.Discrete(...
	{sin(0.5*z*pi)+1; 0.9-z/2}, ...
	'grid', z, ...
	'gridName', 'z');

testCase.verifyTrue( myVectorDiscrete.vec2diag.near(myMatrixDiscrete) );
end

function testInvert(testCase)
myGrid = linspace(0, 1, 21);
% scalar
fScalar = quantity.Discrete((myGrid.').^2, 'grid', myGrid, 'gridName', 'x');
fScalarInverse = fScalar.invert(fScalar.gridName{1});
testCase.verifyEqual(fScalarInverse.on(fScalar.on()), myGrid.'),
end

function testSolveAlgebraic(testCase)
myGrid = linspace(0, 1, 21);
% scalar
fScalar = quantity.Discrete((1+myGrid.').^2, 'grid', myGrid, 'gridName', 'x', ...
	'size', [1, 1]);
solutionScalar = fScalar.solveAlgebraic(2, fScalar.gridName{1});
testCase.verifyEqual(solutionScalar, sqrt(2)-1, 'AbsTol', 1e-3);

% array
% fArray = quantity.Discrete([2*myGrid.', myGrid.' + ones(numel(myGrid), 1)], ...
% 	'grid', myGrid', 'gridName', 'x', 'size', [2, 1]);
% solution = fArray.solveAlgebraic([1; 1], fArray(1).gridName{1});
% testCase.verifyEqual(solution, 0.5);
end

function testSolveDVariableEqualQuantityConstant(testCase)
%% simple constant case
quan = quantity.Discrete(ones(3, 1), 'grid', linspace(0, 1, 3), ...
	'gridName', 'z', 'size', [1, 1], 'name', 'a');
solution = quan.solveDVariableEqualQuantity();
[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
end

function testSolveDVariableEqualQuantityNegative(testCase)
z = testCase.TestData.sym.z;
zeta = testCase.TestData.sym.zeta;

assume(z>0 & z<1); assume(zeta>0 & zeta<1);
myParameterGrid = linspace(0, 0.1, 5);
Lambda = quantity.Symbolic(-0.1-z^2, ...%, -1.2+z^2]),...1+z*sin(z)
	'grid', myParameterGrid, 'gridName', 'z', 'name', '\Lambda');

%%
myGrid = linspace(-.2, .2, 11);
sEval = -0.09;
solveThisOdeDiscrete = solveDVariableEqualQuantity(...
	diag2vec(quantity.Discrete(Lambda)), 'variableGrid', myGrid);
solveThisOdeSymbolic = solveDVariableEqualQuantity(...
	diag2vec(Lambda), 'variableGrid', myGrid);

testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.05}), solveThisOdeDiscrete.on({sEval, 0.05}), 'RelTol', 5e-3);
end

function testSolveDVariableEqualQuantityComparedToSym(testCase)
%% compare with symbolic implementation
z = testCase.TestData.sym.z;

assume(z>0 & z<1);
quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 5)}, ...
	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 5)}, ...
	'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
solutionBSym = quanBSym.solveDVariableEqualQuantity();
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
%solutionBSym.plot(); solutionBDiscrete.plot();
testCase.verifyEqual(solutionBDiscrete.on(), solutionBSym.on(), 'RelTol', 1e-6);
end

function testSolveDVariableEqualQuantityAbsolut(testCase)
%% compare with symbolic implementation
z = testCase.TestData.sym.z;

assume(z>0 & z<1);
quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 11)}, ...
	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 11)}, ...
	'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));

solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
myGrid = solutionBDiscrete.grid{1};
solutionBDiscreteDiff = solutionBDiscrete.diff('s');
quanOfSolutionOfS = zeros(length(myGrid), length(myGrid), 1);
for icIdx = 1 : length(myGrid)
	quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
end
%%
testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-2);

end

function testMtimesDifferentGridLength(testCase)
%%
a = quantity.Discrete(ones(11, 1), 'grid', linspace(0, 1, 11), 'gridName', 'z', ...
	'size', [1, 1], 'name', 'a');
b = quantity.Discrete(ones(5, 1), 'grid', linspace(0, 1, 5), 'gridName', 'z', ...
	'size', [1, 1], 'name', 'b');
ab  = a*b;
%
z = testCase.TestData.sym.z;

c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'gridName', 'z', 'name', 'c');
ac = a*c;

%%
testCase.verifyEqual(ab.on(), ones(11, 1));
testCase.verifyEqual(ac.on(), ones(21, 1));
end

function testDiff1d(testCase)
%% 1d
constant = quantity.Discrete([2*ones(11, 1), linspace(0, 1, 11).'], 'grid', linspace(0, 1, 11), ...
	'gridName', 'z', 'name', 'constant', 'size', [2, 1]);

constantDiff = diff(constant);
testCase.verifyEqual(constantDiff.on(), [zeros(11, 1), ones(11, 1)], 'AbsTol', 10*eps);

z = linspace(0,pi)';
sinfun = quantity.Discrete(sin(z), 'grid', z, 'gridName', 'z');

% do the comparison on a smaller grid, because the numerical derivative is
% very bad a the boundarys of the domain.
Z = linspace(0.1, pi-0.1)';
testCase.verifyTrue(numeric.near(sinfun.diff().on(Z), cos(Z), 1e-3));
testCase.verifyTrue(numeric.near(sinfun.diff('z', 2).on(Z), -sin(Z), 1e-3));
testCase.verifyTrue(numeric.near(sinfun.diff('z', 3).on(Z), -cos(Z), 1e-3));

end

function testDiffConstant2d(testCase)
%% 2d
[zNdgrid, zetaNdgrid] = ndgrid(linspace(0, 1, 11), linspace(0, 1, 21));
myQuantity = quantity.Discrete(cat(3, 2*ones(11, 21), zNdgrid, zetaNdgrid), ...
	'grid', {linspace(0, 1, 11), linspace(0, 1, 21)}, ...
	'gridName', {'z', 'zeta'}, 'name', 'constant', 'size', [3, 1]);
myQuantityDz = diff(myQuantity, 'z', 1);
myQuantityDzeta = diff(myQuantity, 'zeta', 1);

testCase.verifyError( @() diff(myQuantity), 'MATLAB:functionValidation:NotScalar' )
testCase.verifyError( @() diff(myQuantity, 1, {'z', 'zeta'}), 'MATLAB:functionValidation:ClassConversionError' )
testCase.verifyError( @() diff(myQuantity,  {'z', 'zeta'}), 'MATLAB:functionValidation:NotScalar' )

% constant
testCase.verifyEqual(myQuantityDz(1).on(), zeros(11, 21));
testCase.verifyEqual(myQuantityDzeta(1).on(), zeros(11, 21));


% zNdgrid
testCase.verifyEqual(myQuantityDz(2).on(), ones(11, 21), 'AbsTol', 10*eps);
testCase.verifyEqual(myQuantityDzeta(2).on(), zeros(11, 21), 'AbsTol', 10*eps);

% zetaNdgrid
testCase.verifyEqual(myQuantityDz(3).on(), zeros(11, 21), 'AbsTol', 10*eps);
testCase.verifyEqual(myQuantityDzeta(3).on(), ones(11, 21), 'AbsTol', 10*eps);

end

function testOn(testCase)
% init data
domainVecA = quantity.Domain('z', linspace(0, 1, 27));
domainVecB = quantity.Domain('zeta', linspace(0, 1, 41));
[value] = createTestData(domainVecA.grid, domainVecB.grid);
a = quantity.Discrete(value, 'size', [2, 3], ...
	'domain', [domainVecA, domainVecB], 'name', 'A');

%
testCase.verifyEqual(value, a.on());
testCase.verifyEqual(permute(value, [2, 1, 3, 4]), ...
	a.on([domainVecB, domainVecA]));
testCase.verifyEqual(createTestData(linspace(0, 1, 100), linspace(0, 1, 21)), ...
	a.on({linspace(0, 1, 100), linspace(0, 1, 21)}), 'AbsTol', 100*eps);

d24z = quantity.Domain('z', linspace(0, 1, 24));
d21zeta = quantity.Domain('zeta', linspace(0, 1, 21));

testCase.verifyEqual(createTestData(linspace(0, 1, 24), linspace(0, 1, 21)), ...
	a.on( [d24z, d21zeta] ), 'AbsTol', 24*eps);

testCase.verifyEqual(permute(createTestData(linspace(0, 1, 21), linspace(0, 1, 24)), [2, 1, 3, 4]), ...
	a.on( {linspace(0, 1, 24), linspace(0, 1, 21)}, {'zeta', 'z'}), 'AbsTol', 2e-4);

	function [value] = createTestData(gridVecA, gridVecB)
		[zGridA , zetaGridA] = ndgrid(gridVecA, gridVecB);
		value = ones(numel(gridVecA), numel(gridVecB), 2, 3);
		value(:,:,1,1) = 1+zGridA;
		value(:,:,2,1) = 1+zetaGridA;
		value(:,:,2,2) = 1+2*zGridA - zetaGridA.^2;
		value(:,:,1,2) = 1+zGridA + zetaGridA.^2;
		value(:,:,2,3) = 1+zeros(numel(gridVecA), numel(gridVecB));
	end

%% test on on 3-dim domain
z = quantity.Domain('z', 0:3);
t = quantity.Domain('t', 0:4);
x = quantity.Domain('x', 0:5);

Z = quantity.Discrete( z.grid, 'domain', z);
T = quantity.Discrete( t.grid, 'domain', t);
X = quantity.Discrete( x.grid, 'domain', x);

ZTX = Z+T+X;
XTZ = X+T+Z;

testCase.verifyEqual( ZTX.on([t z x]), XTZ.on([t z x]), 'AbsTol', 1e-12);

end

function testOn2(testCase)
zeta = linspace(0, 1);
eta = linspace(0, 1, 71)';
z = linspace(0, 1, 51);

[ZETA, ETA, Z] = ndgrid(zeta, eta, z);

a = quantity.Discrete(cat(4, sin(ZETA.*ZETA.*Z)+ETA.*ETA, ETA+cos(ZETA.*ETA.*Z)), ...
	'size', [2 1], 'grid', {zeta, eta, z}, 'gridName', {'zeta', 'eta', 'z'});

testCase.verifyEqual(a.on({zeta, eta, z}), ...
	permute(a.on({eta, zeta, z}, {'eta', 'zeta', 'z'}), [2, 1, 3, 4]));
testCase.verifyEqual([a(1).on({zeta(2), eta(3), z(4)}); a(2).on({zeta(2), eta(3), z(4)})], ...
	[sin(zeta(2)*zeta(2)*z(4))+eta(3)*eta(3); eta(3)+cos(zeta(2)*eta(3)*z(4))]);
end

function testOnDescending(tc)
domainVecA = quantity.Domain('z', linspace(0, 1, 27));
myData = quantity.Discrete(domainVecA.grid+1, 'domain', domainVecA);
tc.verifyEqual(flip(domainVecA.grid)+1, myData.on(flip(domainVecA.grid)), 'AbsTol', 1e3*eps);
tc.verifyEqual(flip(domainVecA.grid)+1, myData.on(flip(domainVecA.grid), 'z'), 'AbsTol', 1e3*eps);
tc.verifyEqual(flip(domainVecA.grid)+1, myData.on({flip(domainVecA.grid)}, {'z'}), 'AbsTol', 1e3*eps);
end % testOnDescending()

function testMtimesZZeta2x2(testCase)



%%
z11_ = linspace(0, 1, 11);
z11 = quantity.Domain('z', z11_);
zeta11 = quantity.Domain('zeta', z11_);

a = quantity.Discrete(ones([z11.n, z11.n, 2, 2]), 'size', [2, 2], ...
	'domain', [z11, zeta11], 'name', 'A');

zeta = testCase.TestData.sym.zeta;
z10_ = linspace(0, 1, 10);

b = quantity.Symbolic((eye(2, 2)), 'gridName', 'zeta', ...
	'grid', z10_);

c = a*b;

%
testCase.verifyEqual(c.on(), a.on());
end

function testMTimesPointWise(tc)

z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

Z = linspace(0, pi, 51)';
ZETA = linspace(0, pi / 2, 71)';

P = quantity.Symbolic( sin(zeta) * z, 'grid', {Z, ZETA}, 'gridName', {'z', 'zeta'});
B = quantity.Symbolic( sin(zeta), 'grid', ZETA, 'gridName', 'zeta');

PB = P*B;

p = quantity.Discrete(P);
b = quantity.Discrete(B);

pb = p*b;

tc.verifyEqual(MAX(abs(PB-pb)), 0, 'AbsTol', 10*eps);
end

function testMldivide(testCase)
%% scalar example
s = linspace(0, 1, 21);
vL = quantity.Discrete(2*ones(21, 1), ...
	'size', [1, 1], 'grid', {s}, 'gridName', {'asdf'});
vR = quantity.Discrete((linspace(1, 2, 21).').^2, ...
	'size', [1, 1], 'grid', {s}, 'gridName', {'asdf'});
vLdVR = vL \ vR;


%% matrix example
t = linspace(0, pi);
s = linspace(0, 1, 21);
[T, ~] = ndgrid(t, s);

KL = quantity.Discrete(permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), ...
	'size', [2 2], 'grid', {t, s}, 'gridName', {'t', 's'});
KR = quantity.Discrete(permute(repmat(4*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), ...
	'size', [2 2], 'grid', {t, s}, 'gridName', {'t', 's'});
KLdKR = KL \ KR;

%%
testCase.verifyEqual(vLdVR.on(), 2 .\ (linspace(1, 2, 21).').^2);

% TODO #discussion: for the following case the operation KL \ KR seems to change the order of the grid in the old version. As this is very unexpected. I would like to change them to not change the order! 
testCase.verifyEqual(KLdKR.on(), permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]));

end

function testMrdivide(testCase)
%% scalar example
s = linspace(0, 1, 21);
vL = quantity.Discrete(2*ones(21, 1), ...
	'size', [1, 1], 'grid', {s}, 'gridName', {'asdf'});
vR = quantity.Discrete((linspace(1, 2, 21).').^2, ...
	'size', [1, 1], 'grid', {s}, 'gridName', {'asdf'});
vLdVR = vL / vR;


%% matrix example
t = linspace(0, pi);
s = linspace(0, 1, 21);
[T, ~] = ndgrid(t, s);

KL = quantity.Discrete(permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), ...
	'size', [2 2], 'grid', {t, s}, 'gridName', {'t', 's'});
KR = quantity.Discrete(permute(repmat(4*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), ...
	'size', [2 2], 'grid', {t, s}, 'gridName', {'t', 's'});
KLdKR = KL / KR;

%%
testCase.verifyEqual(vLdVR.on(), 2 ./ (linspace(1, 2, 21).').^2);

% TODO #discussion: for the following case the operation KL \ KR seems to change the order of the grid in the old version. As this is very unexpected. I would like to change them to not change the order! 
testCase.verifyEqual(KLdKR.on(),  permute(repmat(0.5*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]));
end

function testInv(testCase)
%% scalar example
s = linspace(0, 1, 21);
v = quantity.Discrete((linspace(1, 2, 21).').^2, ...
	'size', [1, 1], 'grid', {s}, 'gridName', {'asdf'});
vInv = v.inv();

%% matrix example
t = linspace(0, pi);
s = linspace(0, 1, 21);
[T, ~] = ndgrid(t, s);

K = quantity.Discrete(permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), ...
	'size', [2 2], 'grid', {t, s}, 'gridName', {'t', 's'});
kInv = inv(K);

%%
testCase.verifyEqual(vInv.on(), 1./(linspace(1, 2, 21).').^2)
testCase.verifyEqual(kInv.on(), permute(repmat(0.5*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]))

end

function testCumInt(testCase)

tGrid = linspace(pi, 1.1*pi, 51)';
s = tGrid;
[T, S] = ndgrid(tGrid, tGrid);

t = testCase.TestData.sym.t;
sy = testCase.TestData.sym.sy;

a = [ 1, sy; t, 1];
b = [ sy; 2*sy];

A = zeros([size(T), size(a)]);
B = zeros([length(s), size(b)]);

for i = 1:size(a,1)
	for j = 1:size(a,2)
		A(:,:,i,j) = subs(a(i,j), {sy, t}, {S, T});
	end
end

for i = 1:size(b,1)
	B(:,i) = subs(b(i), sy, s);
end

%% int_0_t a(t,s) * b(s) ds
% compute symbolic version of the volterra integral
v = int(a*b, sy, tGrid(1), t);
V = quantity.Symbolic(v, 'grid', tGrid, 'gridName', 't');

% compute the numeric version of the volterra integral
k = quantity.Discrete(A, 'size', size(a), 'grid', {tGrid, s}, 'gridName', {'t', 's'});
x = quantity.Discrete(B, 'size', size(b), 'grid', {s}, 'gridName', 's');

f = cumInt(k * x, 's', s(1), 't');

testCase.verifyEqual(V.on(), f.on(), 'AbsTol', 1e-5);

%% int_s_t a(t,s) * b(s) ds
v = int(a*b, sy, sy, t);
V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s}, 'gridName', {'t', 's'});

f = cumInt(k * x, 's', 's', 't');

testCase.verifyEqual( f.on(f(1).grid, f(1).gridName), V.on(f(1).grid, f(1).gridName), 'AbsTol', 1e-5 );

%% int_s_t a(t,s) * c(t,s) ds
c = [1, sy+1; t+1, 1];
C = zeros([size(T), size(c)]);
for i = 1:numel(c)
	C(:,:,i) = double(subs(c(i), {t sy}, {T S}));
end
y = quantity.Discrete(C, 'size', size(c), 'grid', {tGrid s}, 'gridName', {'t' 's'});

v = int(a*c, sy, sy, t);
V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s}, 'gridName', {'t', 's'});
f = cumInt( k * y, 's', 's', 't');

testCase.verifyEqual( f.on(f(1).grid, f(1).gridName), V.on(f(1).grid, f(1).gridName), 'AbsTol', 1e-5 );

%% testCumIntWithLowerAndUpperBoundSpecified
tGrid = linspace(pi, 1.1*pi, 51)';
s = tGrid;
[T, S] = ndgrid(tGrid, tGrid);

sy = testCase.TestData.sym.sy;
t = testCase.TestData.sym.t;

a = [ 1, sy; t, 1];

A = zeros([size(T), size(a)]);

for i = 1:size(a,1)
	for j = 1:size(a,2)
		A(:,:,i,j) = subs(a(i,j), {sy, t}, {S, T});
	end
end

% compute the numeric version of the volterra integral
k = quantity.Discrete(A, 'size', size(a), 'grid', {tGrid, s}, 'gridName', {'t', 's'});

fCumInt2Bcs = cumInt(k, 's', 'zeta', 't');
fCumInt2Cum = cumInt(k, 's', tGrid(1), 't') ...
	- cumInt(k, 's', tGrid(1), 'zeta');
testCase.verifyEqual(fCumInt2Bcs.on(), fCumInt2Cum.on(), 'AbsTol', 10*eps);

%% testCumInt with one independent variable
tGrid = linspace(0, pi, 51)';

a = [ 1, sin(t); t, 1];

% compute the numeric version of the volterra integral
f = quantity.Symbolic(a, 'grid', {tGrid}, 'gridName', {'t'});
f = quantity.Discrete(f);

intK = cumInt(f, 't', 0, 't');

K = quantity.Symbolic( int(a, 0, t), 'grid', {tGrid}, 'gridName', {'t'});

testCase.verifyEqual(intK.on(), K.on(), 'AbsTol', 1e-3);


end

function testAtIndex(testCase)

z = linspace(0,1).';
a = sin(z * pi);

A = quantity.Discrete({a}, 'grid', {z}, 'gridName', {'z'});

testCase.verifyEqual(a(end), A.atIndex(end));
testCase.verifyEqual(a(1), A.atIndex(1));
testCase.verifyEqual(a(23), A.atIndex(23));

y = linspace(0,2,51);
b1 = sin(z * pi * y);
b2 = cos(z * pi * y);
B = quantity.Discrete({b1; b2}, 'grid', {z, y}, 'gridName', {'z', 'y'});

B_1_y = B.atIndex(1,1);
b_1_y = [sin(0); cos(0)];
testCase.verifyTrue(numeric.near(B_1_y, b_1_y));

B_z_1 = B.atIndex(':',1);

testCase.verifyTrue(all(B_z_1(:,1,1) == 0));
testCase.verifyTrue(all(B_z_1(:,:,2) == 1));

end

function testMTimes(testCase)
% multiplication of a(z) * a(z)'
% a(z) \in (3, 2), z \in (100)
z = linspace(0, 2*pi)';
a = cat(3, [z*0, z, sin(z)], [z.^2, z.^3, cos(z)]); % dim = (100, 3, 2)
aa = misc.multArray(a, permute(a, [1, 3, 2]), 3, 2, 1);    % dim = (100, 3, 3)

A = quantity.Discrete(a, 'size', [3, 2], 'grid', z, 'gridName', {'z'});
AA = A*A';

testCase.verifyTrue(numeric.near(aa, AA.on()));

z = linspace(0, 2*pi, 31);
t = linspace(0, 1, 13);
v = linspace(0, 1, 5);
[z, t] = ndgrid(z, t);
[v, t2] = ndgrid(v, t(1,:));