testSymbolic.m

function [tests ] = testSymbolic()
%TESTQUANTITY Unittests for the quantity objects
%   Add all unittests for the quantity files here
tests = functiontests(localfunctions());
end

function testSubsDomain(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
zeta = quantity.Domain("zeta", linspace(0, 1, 11));
zLr = quantity.Domain("z", linspace(0, 1, 5));
eta = quantity.Domain("eta", linspace(0, 1, 11));
f = z.Symbolic() * zeta.Symbolic();
%
fofZEta = f.subsDomain("zeta", eta);	% = f(z, eta)
tc.verifyEqual(fofZEta.on(), f.on());
tc.verifyEqual([fofZEta(1).domain.name], ["z", "eta"]);
%
fofZlrZeta = f.subsDomain("z", zLr);	% = f(z, zeta)
tc.verifyEqual(fofZlrZeta.on(), f.on([zLr, zeta]));
tc.verifyEqual([fofZlrZeta(1).domain.name], ["z", "zeta"]);

%
fofZ = f.subsDomain("zeta", zLr);		% = f(z) 
tc.verifyEqual(fofZ.on(), diag(f.on([z, quantity.Domain("zeta", z.grid)])));
tc.verifyEqual([fofZ(1).domain.name], "z");
end % testSubsDomain()

function testSubsNumeric(tc)
z = quantity.Domain("z", linspace(0, 1, 5));
zeta = quantity.Domain("zeta", linspace(0, 1, 4));
t = quantity.Domain("t", linspace(0, 1, 6));

quanScalar = z.Symbolic() + 1 + zeta.Symbolic() * t.Symbolic();
tc.verifyEqual(quanScalar.subsNumeric(["z", "t"], [0, 0]).on(), ones(zeta.n, 1));
tc.verifyEqual(quanScalar.subsNumeric(["t", "z", "zeta"], [0.1, 0.2, 0.3]), 0.2 + 1 + 0.3 * 0.1);

quanMatrix = quantity.Symbolic(...
	[sym("z") + 1 + sym("zeta") * sym("t"), 1; ...
	sym("z") + sym("t"), -sym("zeta")], [zeta, t, z]);
tc.verifyEqual(quanMatrix.subsNumeric(["z", "zeta", "t"], [0.1, 0.2, 0.3]), ...
	[0.1 + 1 + 0.2 * 0.3, 1; 0.1 + 0.3, -0.2], "AbsTol", 10*eps);
tc.verifyEqual(quanMatrix.subsNumeric(["z", "zeta", "t"], [0, 0, 0]), [1, 1; 0, 0]);
tc.verifyEqual(quanMatrix.subsNumeric(["z", "zeta", "t"], [1, 1, 1]), [3, 1; 2, -1]);
tc.verifyEqual(quanMatrix(1, 1).subsNumeric(["z", "zeta", "t"], [0.1, 0.2, 0.3]), ...
	0.1 + 1 + 0.2 * 0.3, "AbsTol", 10*eps);
end % testSubsNumeric()

function testChangeDomain(tc)
zeta = quantity.Domain("zeta", linspace(0, 1, 4));
t = quantity.Domain("t", linspace(0, 1, 4));
thisSym = sym("zeta") + sym("t");

quan = quantity.Symbolic(thisSym, [zeta, t]);
newDomain = [quantity.Domain("zeta", linspace(0, 1, 11)), quantity.Domain("t", linspace(0, 1, 21))];
newQuan = quan.changeDomain(newDomain);

tc.verifyEqual(newQuan.on(), quan.on(newDomain));
tc.verifyTrue(isa(newQuan, "quantity.Symbolic"));
end % testChangeDomain()

function testOn2(tc)
zeta = quantity.Domain("zeta", linspace(0, 1, 4));
t = quantity.Domain("t", linspace(0, 1, 4));
thisSym = sym("zeta") + sym("t");
quan = quantity.Symbolic(thisSym, [zeta, t]);

solution = 0.2 + 0.3; % this obtained by manual calculation
% verify solution
tc.verifyEqual(double(subs(subs(thisSym, "zeta", 0.2), "t", 0.3)), solution, "AbsTol", 10*eps);

% test on with numeric values only
tc.verifyEqual(quan.on({0.2, 0.3}), solution);

% test on with domain
tc.verifyEqual(quan.on([quantity.Domain("zeta", 0.2), quantity.Domain("t", 0.3)]), solution);
end % testOn2()

function testOn2b(tc)
zeta = quantity.Domain("zeta", linspace(0, 1, 4));
t = quantity.Domain("t", linspace(0, 1, 4));
thisSym = sym("zeta") + sym("t");
quan = quantity.Symbolic(thisSym, [t, zeta]);

solution = 0.2 + 0.3; % this obtained by manual calculation
% verify solution
tc.verifyEqual(double(subs(subs(thisSym, "zeta", 0.2), "t", 0.3)), solution, "AbsTol", 10*eps);

% test on with numeric values only
tc.verifyEqual(quan.on({0.3, 0.2}), solution);

% test on with domain
tc.verifyEqual(quan.on([quantity.Domain("zeta", 0.2), quantity.Domain("t", 0.3)]), solution);
end % testOn2()

function testOn3(tc)
z = quantity.Domain("z", linspace(0, 1, 4));
t = quantity.Domain("t", linspace(0, 1, 4));
thisSym = sym("z") + sym("t");
quan = quantity.Symbolic(thisSym, [z, t]);

solution = 0.2 + 0.3; % this obtained by manual calculation
% verify solution
tc.verifyEqual(double(subs(subs(thisSym, "z", 0.2), "t", 0.3)), solution, "AbsTol", 10*eps);

% test on with numeric values only
tc.verifyEqual(quan.on({0.2, 0.3}), solution);

% test on with domain
tc.verifyEqual(quan.on([quantity.Domain("z", 0.2), quantity.Domain("t", 0.3)]), solution);
end % testOn2()

function testDiffOfConstants(tc)
% unittest to fix a bug...
z = quantity.Domain("z", linspace(0, 1, 4));
K0 = quantity.Symbolic.zeros(1, z);%[z, zeta]);
tc.verifyEqual(K0.on(), zeros([4, 1]));

K1 = quantity.Symbolic.ones(1, z);%[z, zeta]);
tc.verifyEqual(K1.diff("z", 1).on(), zeros([4, 1]));

zeta = quantity.Domain("zeta", linspace(0, 1, 5));
K2 = quantity.Symbolic.ones([2, 2], [z, zeta]);
tc.verifyEqual(K2.diff("z", 1).on(), zeros([4, 5, 2, 2]));
end % testDiffOfConstants

function testZeroOn(tc)
% unittest to fix a bug related with isNumber() usage in on() of symbolic
z = quantity.Domain("z", linspace(0, 1, 11));
quan = 0 * quantity.Symbolic(1-0.5*sym("z"), z, "name", "Gamma");
tc.verifyEqual(quan.on({linspace(0, 1, 5)}), zeros(5, 1));
tc.verifyEqual(quan.on(linspace(0, 1, 5)), zeros(5, 1));
end % testZeroOn()

function testOnes(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(1, myDomain) - 1)), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(2, myDomain) - ones(2))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 1], myDomain) - ones([2, 1]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 2], myDomain(2)) - ones([2, 2]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([4, 4], myDomain) - ones(4))), 0);
end % testOnes()

function testEye(tc)
myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(1, myDomain) - 1)), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(2, myDomain) - eye(2))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 1], myDomain) - eye([2, 1]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 2], myDomain(2)) - eye([2, 2]))), 0);
tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([4, 4], myDomain) - eye(4))), 0);
end % testEye()

function testValueContinuousBug(tc)
z = quantity.Domain("z", linspace(0, 1, 7));
zeta = quantity.Domain("zeta", linspace(0, 1, 5));
A = quantity.Symbolic(sin(sym("z") * pi), z);

B1 = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
B2 = [zeta.Symbolic + 0*z.Symbolic, A + 0*zeta.Symbolic; ...
	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
tc.verifyEqual(B1.on([z, zeta]), B2.on([z, zeta]));

x = A * subs(A, "z", "zeta");
C1 = [0*z.Symbolic + zeta.Symbolic; x];
C2 = [zeta.Symbolic + 0*z.Symbolic; x];
tc.verifyEqual(C1.on([z, zeta]), C2.on([z, zeta]));

e = zeta.Symbolic + 0*z.Symbolic;
e = e.sortDomain();
f = e.changeDomain([z z.rename("zeta")]);


% TODO: fix this error! See issue #41
% this problem is maybe caused as input arguments of B2(1).valueContinuous is somehow ordered
% differently compared to B2(2:4).valueContinuous or all B(:).valueContinuous, 
% i.e. B2(1).valueContinuous is (zeta, z) instead of (z, zeta).
end % testValueContinuousBug()

function testPower(tc)
%% scalar case
z = quantity.Domain("z", linspace(0, 1, 11));
A = quantity.Symbolic(sin(sym("z") * pi), z);
aa = sin(z.grid * pi) .* sin(z.grid * pi);
AA = A.^2;

tc.verifyEqual( aa, AA.on());
tc.verifyEqual( on(A.^0), ones(z.n, 1));
tc.verifyEqual( on(A.^3), on(AA * A), "AbsTol", 1e-15);

%% matrix case
zeta = quantity.Domain("zeta", linspace(0, 1, 21));
B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
BBB = B.^3;
BBB_alt = B .* B .* B;
tc.verifyEqual(BBB.on(), BBB_alt.on(), "AbsTol", 1e-14);
tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, "AbsTol", 1e-15);
tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), "AbsTol", 1e-15);
end % testPower()

function testMPower(tc)
%% scalar case
z = quantity.Domain("z", linspace(0, 1, 11));
A = quantity.Symbolic(sin(sym("z") * pi), z);
aa = sin(z.grid * pi) .* sin(z.grid * pi);
AA = A^2;

tc.verifyEqual( aa, AA.on());
tc.verifyEqual( on(A^0), ones(z.n, 1));
tc.verifyEqual( on(A^3), on(AA * A));

%% matrix case
zeta = quantity.Domain("zeta", linspace(0, 1, 21));
B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
B_discrete = quantity.Discrete(B);
BBB = B^3;
BBB_alt = B * B * B;
BBB_discrete = B_discrete^3;
tc.verifyEqual(BBB.on(), BBB_alt.on(), "AbsTol", 1e-14);
tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6}))^3, "AbsTol", 1e-15);
tc.verifyEqual(BBB_discrete.at({0.5, 0.6}), (B.at({0.5, 0.6}))^3, "AbsTol", 1e-15);
tc.verifyEqual(on(B^0), permute(repmat(eye(2), [1, 1, B(1).domain(1).n, B(1).domain(2).n]), [3, 4, 1, 2]), ...
	"AbsTol", 1e-15);
end % testMPower()

function testNorm(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
quan = Symbolic(z) * [-1; 2; 3];
tc.verifyEqual(quan.norm.on(), on(z.Discrete*sqrt(1+4+9)), "AbsTol", 1e-15);
tc.verifyEqual(quan.norm.on(), on(quan.norm(2)), "AbsTol", 1e-15);
end % testNorm()

function testLog(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
quan = Symbolic(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);

tc.verifyEqual(log(quan.on()), on(log(quan)), "AbsTol", 1e-15);
tc.verifyEqual(on(quan), on(log(exp(quan))), "AbsTol", 1e-15);
end % testLog()

function testLog10(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
quan = Symbolic(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);

tc.verifyEqual(log10(quan.on()), on(log10(quan)), "AbsTol", 1e-15);
end % testLog10()

function testDet(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
zSymbolic = quantity.Symbolic(sym("z"), z);
myMatrix2x2 = [1, 2; 3, 4] * zSymbolic;
myMatrix3x3 = [1, 2, 3; 2, -3, 4; 3, 4, 5] * zSymbolic;
myMatrix5x5 = [1, 2, 3, 4, 5; ...
				2, 3, 4, 5, 6; ...
				3, 4, 5, 6, 7; ...
				4, -6, 7, 8, 9; ...
				1, 1, 1, 1, 1] * zSymbolic;
tc.verifyEqual(zSymbolic.det.on(), zSymbolic.on(), "AbsTol", 10*eps);
tc.verifyEqual(myMatrix2x2.det.at(0.5), det(myMatrix2x2.at(0.5)), "AbsTol", 10*eps);
tc.verifyEqual(myMatrix3x3.det.at(0.5), det(myMatrix3x3.at(0.5)), "AbsTol", 10*eps);
tc.verifyEqual(myMatrix5x5.det.at(0.5), det(myMatrix5x5.at(0.5)), "AbsTol", 10*eps);
end % testDet()

function testSum(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
myOnes = ones(3, 4) * (1 + 0*quantity.Symbolic(sym("z"), z));
tc.verifyEqual(on(sum(myOnes)), 3*4*ones(11, 1), "AbsTol", 10*eps);
tc.verifyEqual(on(sum(myOnes, 2)), 4*ones(11, 3), "AbsTol", 10*eps);
tc.verifyEqual(on(sum(myOnes, 1)), 3*ones(11, 1, 4), "AbsTol", 10*eps);
tc.verifyEqual(on(sum(myOnes, [1, 2])), 3*4*ones(11, 1), "AbsTol", 10*eps);

myTwos3D = ones(2, 3, 4) * 2 * (1 + 0*quantity.Symbolic(sym("z"), z));
tc.verifyEqual(on(sum(myTwos3D)), 2*3*4*ones(11, 1)*2, "AbsTol", 10*eps);
tc.verifyEqual(on(sum(myTwos3D, 2)), 3*ones(11, 2, 1, 4)*2, "AbsTol", 10*eps);
tc.verifyEqual(on(sum(myTwos3D, [1, 3])), 2*4*ones(11, 1, 3)*2, "AbsTol", 10*eps);
end % testSum()

function testProd(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
myTwos = 2 * ones(3, 4) * (1 + 0*quantity.Symbolic(sym("z"), z));
tc.verifyEqual(on(prod(myTwos)), 2^(3*4)*ones(11, 1), "AbsTol", 10*eps);
tc.verifyEqual(on(prod(myTwos, 2)), 2^4*ones(11, 3), "AbsTol", 10*eps);
tc.verifyEqual(on(prod(myTwos, 1)), 2^3*ones(11, 1, 4), "AbsTol", 10*eps);
tc.verifyEqual(on(prod(myTwos, [1, 2])), 2^(3*4)*ones(11, 1), "AbsTol", 10*eps);

myTwos3D = quantity.Symbolic(2*(1+0*sym("z")) * ones(2, 3, 4), z);
tc.verifyEqual(on(prod(myTwos3D)), 2^(2*3*4)*ones(11, 1), "AbsTol", 10*eps);
tc.verifyEqual(on(prod(myTwos3D, 2)), 2^3*ones(11, 2, 1, 4), "AbsTol", 10*eps);
tc.verifyEqual(on(prod(myTwos3D, [1, 3])), 2^(2*4)*ones(11, 1, 3), "AbsTol", 10*eps);
end % testProd()

function testVec2Diag(testCase)
% quantity.Symbolic
syms z
Z = quantity.Domain("z", linspace(0, 1, 3));
myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], Z);
myVectorSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1; 0.9-z/2], Z);

testCase.verifyTrue( myVectorSymbolic.vec2diag.near(myMatrixSymbolic) );
end

function testCastSymbolic2Function(testCase)

z = linspace(0, 2*pi)';
Z = sym("z");
s = quantity.Symbolic([sin(Z); cos(Z)], quantity.Domain("z", z));

f = quantity.Function(s);

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

end

function testSymbolicEvaluation(tc)
syms z
myGrid = linspace(0, 1, 7);
myDomain = quantity.Domain("z", myGrid);
fS = quantity.Symbolic(z * sinh(z * 1e5) / cosh(z * 1e5), myDomain, ...
	"symbolicEvaluation", true);
fF = quantity.Symbolic(z * sinh(z * 1e5) / cosh(z * 1e5), myDomain, ...
	"symbolicEvaluation", false);

tc.verifyTrue(any(isnan(fF.on)))
tc.verifyFalse(any(isnan(fS.on)))
end

function testCtranspose(tc)
syms z zeta

dom.z = quantity.Domain("z", linspace(0, 1, 21));
dom.zeta = quantity.Domain("zeta", linspace(0, 1, 41));

qSymbolic = quantity.Symbolic( [1+zeta, -zeta; -z, z^2], [dom.z, dom.zeta], "name", "q");
qSymbolicCtransp = qSymbolic';
tc.verifyEqual(qSymbolic(1,1).on(), qSymbolicCtransp(1,1).on());
tc.verifyEqual(qSymbolic(2,2).on(), qSymbolicCtransp(2,2).on());
tc.verifyEqual(qSymbolic(1,2).on(), qSymbolicCtransp(2,1).on());
tc.verifyEqual(qSymbolic(2,1).on(), qSymbolicCtransp(1,2).on());

qSymbolic2 = quantity.Symbolic(sym("z") * 2i + sym("zeta"), [dom.z, dom.zeta], "name", "im");
qSymolic2Ctransp = qSymbolic2';
tc.verifyEqual(qSymbolic2.real.on(), qSymolic2Ctransp.real.on());
tc.verifyEqual(qSymbolic2.imag.on(), -qSymolic2Ctransp.imag.on());
end % testCtranspose

function testTranspose(tc)
syms z zeta
dom(1) = quantity.Domain("z", linspace(0, 1, 21));
dom(2) = quantity.Domain("zeta", linspace(0, 1, 41));
qSymbolic = quantity.Symbolic( [1+z*zeta, -zeta; -z, z^2], dom, "name", "q");
qSymbolicTransp = qSymbolic.';
tc.verifyEqual(qSymbolic(1,1).on(), qSymbolicTransp(1,1).on());
tc.verifyEqual(qSymbolic(2,2).on(), qSymbolicTransp(2,2).on());
tc.verifyEqual(qSymbolic(1,2).on(), qSymbolicTransp(2,1).on());
tc.verifyEqual(qSymbolic(2,1).on(), qSymbolicTransp(1,2).on());
end % testTranspose

function testFlipDomain(tc)
syms z zeta
myGrid = linspace(0, 1, 11);

myDomain(1) = quantity.Domain("z", myGrid);
myDomain(2) = quantity.Domain("zeta", myGrid);

f = quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, myDomain);
% flip zeta
fReference = quantity.Symbolic([1+z+(1-zeta); 2*(1-zeta)+sin(z)] + zeros(2, 1)*z*zeta, myDomain);
fFlipped = f.flipDomain("zeta");
tc.verifyEqual(fReference.on, fFlipped.on, "AbsTol", 10*eps);

% flip z and zeta
fReference2 = quantity.Symbolic([1+(1-z)+(1-zeta); 2*(1-zeta)+sin(1-z)] + zeros(2, 1)*z*zeta, ...
	myDomain);
fFlipped2 = f.flipDomain(["z", "zeta"]);
tc.verifyEqual(fReference2.on, fFlipped2.on, "AbsTol", 10*eps);


end % testFlipGrid();

function testInt(testCase)
t = sym("t");
T = quantity.Domain.defaultDomain(10, "t");

f = quantity.Symbolic(sin(t), T);
F = int(f);

testCase.verifyEqual(F, double( int( sin(t), t, 0, 1 ) ), "AbsTol", 6e-17);


end

function testCumInt(testCase)
tGrid = linspace(pi, 1.1*pi, 51)';
s = sym("s");
t = sym("t");

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

%% int_0_t a(t,s) * b(s) ds
% compute symbolic version of the volterra integral
myDomain(1) = quantity.Domain("t", tGrid);
myDomain(2) = quantity.Domain("s", tGrid);
integrandSymbolic = quantity.Symbolic(a*b, myDomain);
integrandDiscrete = quantity.Discrete(integrandSymbolic);
V = cumInt(integrandSymbolic, "s", tGrid(1), "t");
f = cumInt(integrandDiscrete, "s", tGrid(1), "t");

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

%% int_s_t a(t,s) * b(s) ds
V = cumInt(integrandSymbolic, "s", "s", "t");
f = cumInt(integrandDiscrete, "s", "s", "t");

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

%% testCumIntWithLowerAndUpperBoundSpecified
fCumInt2Bcs = cumInt(integrandSymbolic, "s", "zeta", "t");
fCumInt2Cum = cumInt(integrandSymbolic, "s", tGrid(1), "t") ...
	- cumInt(integrandSymbolic, "s", tGrid(1), "zeta");

testCase.verifyEqual(fCumInt2Bcs.on(fCumInt2Bcs(1).domain), ...
	fCumInt2Cum.on(fCumInt2Bcs(1).domain), "AbsTol", 100*eps);

%% test not an integrable function
f = quantity.Symbolic([ sin( sin( s ) ); cos( cos( t ) )], myDomain);

F_sym = cumInt(f, "s", 0, 1);
F_dic = cumInt( quantity.Discrete(f), "s", 0, 1);

testCase.verifyEqual( F_sym.on(), F_dic.on())

end

function testScalarPlusMinusQuantity(testCase)
syms z
myGrid = linspace(0, 1, 7);
f = quantity.Symbolic([1; 2] + zeros(2, 1)*z, quantity.Domain("z", myGrid));
testCase.verifyError(@() 1-f-1, "");
testCase.verifyError(@() 1+f+1, "");
end

function testNumericVectorPlusMinusQuantity(testCase)
syms z
myGrid = linspace(0, 1, 7);
f = quantity.Symbolic([1+z; 2+sin(z)] + zeros(2, 1)*z, quantity.Domain("z", myGrid));

a = ones(size(f));
testCase.verifyEqual(on(a-f), 1-f.on(), "RelTol", 10*eps);
testCase.verifyEqual(on(f-a), f.on()-1, "RelTol", 10*eps);
testCase.verifyEqual(on(a+f), f.on()+1, "RelTol", 10*eps);
testCase.verifyEqual(on(f+a), f.on()+1, "RelTol", 10*eps);
end

function testDiffWith2Variables(testCase)
syms z zeta
myGrid = linspace(0, 1, 7);
myDomain(1) = quantity.Domain("z", myGrid);
myDomain(2) = quantity.Domain("zeta", myGrid);

f = quantity.Symbolic([z*zeta, z; zeta, 1], myDomain);
fdz = quantity.Symbolic([zeta, 1; 0, 0], myDomain);
fdzeta = quantity.Symbolic([z, 0; 1, 0], myDomain);

testCase.verifyEqual(on(diff(f, "z")), fdz.on());
testCase.verifyEqual(on(diff(f, "zeta")), fdzeta.on());

end

function testGridName2variable(testCase)
syms z zeta eta e a
myGrid = linspace(0, 1, 7);

myDomain(1) = quantity.Domain("z", myGrid);
myDomain(2) = myDomain(1).rename("zeta");
myDomain(3) = myDomain(1).rename("eta");
myDomain(4) = myDomain(1).rename("e");
myDomain(5) = myDomain(1).rename("a");

obj = quantity.Symbolic([z*zeta, eta*z; e*a, a], myDomain);

thisGridName = {'eta', 'e', 'z'};
thisVariable = gridName2variable(obj, thisGridName);
testCase.verifyEqual(string(thisGridName).', string(thisVariable));

thisGridName = 'z';
thisVariable = gridName2variable(obj, thisGridName);
testCase.verifyEqual(thisGridName, char(thisVariable));

testCase.verifyError(@() gridName2variable(obj, 't'), "");
testCase.verifyError(@() gridName2variable(obj, 123), "");

end

function testCat(testCase)
syms z zeta
myDomain = quantity.Domain.defaultDomain(21, "z");
f1 = quantity.Symbolic(1+z*z, myDomain, "name", "f1");
f2 = quantity.Symbolic(sin(z), myDomain, "name", "f2");

% vertical concatenation
F = [f1; f2];
testCase.verifyTrue(F(1) == f1)
testCase.verifyTrue(F(2) == f2)

% horizontal concatenation
F = [f1, f2];
testCase.verifyTrue(F(1) == f1)
testCase.verifyTrue(F(2) == f2)

% combined concatenation
F = [F; f1, f2];
testCase.verifyTrue(F(1,1) == f1);
testCase.verifyTrue(F(1,2) == f2);
testCase.verifyTrue(F(2,1) == f1);
testCase.verifyTrue(F(2,2) == f2);

% concatenation on different grids
myDomain = quantity.Domain.defaultDomain(13, "z");
f3 = quantity.Symbolic(cos(z), myDomain, "name", "f1");

F = [f1, f3];
testCase.verifyEqual(F(2).domain.gridLength, f1(1).domain.gridLength)

% concatenate quantitySymbolics and other values
F = [f1, 0];
testCase.verifyEqual(F(2).valueSymbolic, sym(0))
F = [f1, sym(zeros(1, 10))];
testCase.verifyEqual([F(2:end).valueSymbolic], sym(zeros(1,10)))

% 
F = [0, f1];
testCase.verifyEqual(F(1).valueSymbolic, sym(0))
end

function testExp(testCase)
% 1 spatial variable
syms z zeta
myDomain = quantity.Domain.defaultDomain(21, "z");
s1d = quantity.Symbolic(1+z*z, myDomain, "name", "s1d");
testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));

% diagonal matrix
myDomain(2) = quantity.Domain.defaultDomain(41, "zeta");



s2dDiag = quantity.Symbolic([1+z*zeta, 0; 0, z^2], myDomain, "name", "s2dDiag");
testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
end

function testExpm(testCase)
syms z zeta

myDomain = quantity.Domain.defaultDomain(21, "z");
myDomain(2) = quantity.Domain.defaultDomain(41, "zeta");

mat2d = quantity.Symbolic(ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], myDomain, "name", "s2d");
mat2dMat = mat2d.on();
mat2dExpm = 0 * mat2d.on();
for zIdx = 1 : 21
	for zetaIdx = 1 : 41
		mat2dExpm(zIdx, zetaIdx, :, :) = expm(reshape(mat2dMat(zIdx, zetaIdx, :, :), [2, 2]));
	end
end
testCase.verifyEqual(mat2d.expm.on(), mat2dExpm, "RelTol", 100*eps);

end


function testOn(testCase)
syms z zeta
gridZ = linspace(0, 1, 40);
gridZeta = linspace(0, 1, 21);

myDomain = quantity.Domain("z", gridZ);
myDomain(2) = quantity.Domain("zeta", gridZeta);

A = quantity.Symbolic(2+[z, z^2, z+zeta; -sin(zeta), z*zeta, 1], myDomain, 'name', 'A');

%%
AOnPure = A.on();
testCase.verifyEqual(2+gridZ(:)*ones(1, numel(gridZeta)), AOnPure(:,:,1,1));
testCase.verifyEqual(2+gridZ(:)*gridZeta, AOnPure(:,:,2,2));
testCase.verifyEqual(3+0*gridZ(:)*gridZeta, AOnPure(:,:,2,3));

%%
gridZ2 = linspace(0, 1, 15);
gridZeta2 = linspace(0, 1, 31);
AOnChangedGrid = A.on({gridZ2, gridZeta2}, ["z", "zeta"]);
testCase.verifyEqual(2+gridZ2(:)*ones(1, numel(gridZeta2)), AOnChangedGrid(:,:,1,1));
testCase.verifyEqual(2+gridZ2(:)*gridZeta2, AOnChangedGrid(:,:,2,2));
testCase.verifyEqual(3+0*gridZ2(:)*gridZeta2, AOnChangedGrid(:,:,2,3));

%%
AOnChangedGrid2 = A.on({gridZeta2, gridZ2}, {'zeta', 'z'});
testCase.verifyEqual((2+gridZ2(:)*ones(1, numel(gridZeta2))).', AOnChangedGrid2(:,:,1,1));
testCase.verifyEqual((2+gridZ2(:)*gridZeta2).', AOnChangedGrid2(:,:,2,2));
testCase.verifyEqual((3+0*gridZ2(:)*gridZeta2).', AOnChangedGrid2(:,:,2,3));
% testCase.verifyEqual(ones([21, 31, 2, 3]), ...
% 	A.on({linspace(0, 1, 21), linspace(0, 1, 31)}));
% testCase.verifyEqual(ones([21, 31, 2, 3]), ...
% 	A.on({linspace(0, 1, 21), linspace(0, 1, 31)}, ["z", "zeta"]));
% testCase.verifyEqual(ones([21, 31, 2, 3]), ...
% 	A.on({linspace(0, 1, 21), linspace(0, 1, 31)}, {'zeta', 'z'}));
end

function testSqrt(testCase)
% 1 spatial variable
syms z zeta
myDomain = quantity.Domain("z", linspace(0,1,21));
s1d = quantity.Symbolic(1+z*z, myDomain, "name", "s1d");
testCase.verifyEqual(s1d.sqrt.on(), sqrt(s1d.on()));

% 2 spatial variables
myDomain(2) = quantity.Domain("zeta", linspace(0,1,41));
s2d = quantity.Symbolic(1+z*zeta, myDomain, "name", "s2d");
testCase.verifyEqual(s2d.sqrt.on(), sqrt(s2d.on()));

% diagonal matrix
s2dDiag = quantity.Symbolic([1+z*zeta, 0; 0, z^2], myDomain, "name", "s2dDiag");
testCase.verifyEqual(s2dDiag.sqrt.on(), sqrt(s2dDiag.on()));

end

function testSqrtm(testCase)
syms z zeta

myDomain = quantity.Domain("z", linspace(0,1,21));
myDomain(2) = quantity.Domain("zeta", linspace(0,1,41));

mat2d = quantity.Symbolic(ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], myDomain, "name", "s2d");
mat2dMat = mat2d.on();
mat2dSqrtm = 0 * mat2d.on();
for zIdx = 1 : 21
	for zetaIdx = 1 : 41
		mat2dSqrtm(zIdx, zetaIdx, :, :) = sqrtm(reshape(mat2dMat(zIdx, zetaIdx, :, :), [2, 2]));
	end
end
testCase.verifyEqual(mat2d.sqrtm.on(), mat2dSqrtm, "AbsTol", 10*eps);
mat2drootquad = sqrtm(mat2d)*sqrtm(mat2d);
testCase.verifyEqual(mat2drootquad.on(), mat2d.on(), "AbsTol", 100*eps);

end

function testInv(testCase)
%% scalar example
syms z
zGrid = linspace(0, 1, 21);

myDomain = quantity.Domain("z", zGrid);

v = quantity.Symbolic(1+z*z, myDomain, "name", "s1d");
vInvReference = 1 ./ (1 + zGrid.^2);
vInv = v.inv();

%% matrix example
syms zeta
zetaGrid = linspace(0, pi);
[zNdGrid, zetaNdGrid] = ndgrid(zetaGrid, zGrid);

myDomain(2) = quantity.Domain("zeta", zetaGrid);

K = quantity.Symbolic([1+z*z, 0; 0, 1+z*zeta], myDomain);
kInv = inv(K);

%%
testCase.verifyEqual(vInv.on(), vInvReference(:))
testCase.verifyEqual(kInv(1, 1).on(), repmat(vInvReference(:), [1, 100]))
testCase.verifyEqual(kInv(2, 2).on(), 1./(1+zNdGrid .* zetaNdGrid).')

end

function testZero(testCase)

x = quantity.Symbolic(0, quantity.Domain.empty());

verifyEqual(testCase, x.on(), 0);
verifyTrue(testCase, all(on(x) * magic(5) == zeros(5), "all"));

end

function testCast2Quantity(testCase)

syms z t
d = [quantity.Domain("z", linspace(0, 1, 5)), ...
	 quantity.Domain("t", linspace(0, 2, 2))];

x = quantity.Symbolic(sin(z * t * pi), d);

X = quantity.Discrete(x);

verifyTrue(testCase, numeric.near(X.on(), x.on()));
end

function testMrdivide(testCase)
syms z t
d = [quantity.Domain("z", linspace(0, 1, 5)), ...
	 quantity.Domain("t", linspace(0, 2, 2))];
x = quantity.Symbolic(sin(z * t * pi), d);
o = x / x;

verifyTrue(testCase, all(o.on() == 1, "all"));

% test constant values
invsin = 1 / ( x + 8);
ndGridNT = x(1).domain.ndgrid;
verifyTrue(testCase, numeric.near(1 ./ (sin(ndGridNT{1} .* ndGridNT{2} * pi) + 8), invsin.on()));

end

function testMTimesConstants(testCase)
syms z t
d = [quantity.Domain("z", linspace(0, 1, 5)), ...
	 quantity.Domain("t", linspace(0, 2, 2))];
x = quantity.Symbolic(sin(z * t * pi), d);

% test multiplicaiton iwth a constant scalar
x5 = x * 5;
ndGridNT = x(1).domain.ndgrid;
X = sin(ndGridNT{1} .* ndGridNT{2} * pi);

verifyTrue(testCase, numeric.near(X*5, x5.on()));

% test multiplication with a constant vector
v = [2, 3];
x11 = x * v;
X11On = x11.on();

verifyTrue(testCase, numeric.near(X*v(1), X11On(:,:,:,1)));
verifyTrue(testCase, numeric.near(X*v(2), X11On(:,:,:,2)));

% test multiplication with zero:
x0 = x*0;
testCase.verifyEqual(x0.on(), 0*ndGridNT{1});
end

function testMLRdivide2ConstantSymbolic(testCase)

syms z
assume(z>0 & z<1);
myDomain = quantity.Domain.defaultDomain(100, "z");
Lambda = quantity.Symbolic(eye(2), myDomain, "name", "Lambda");
A = quantity.Symbolic(eye(2), myDomain, "name", "A");
C1 = A / Lambda;
C2 = A \ Lambda;
testCase.verifyEqual(C1.on(), A.on());
testCase.verifyEqual(C2.on(), A.on());

end
function testGridNamesNotAllowed(testCase)

try
	myDomain(1) = quantity.Domain.defaultDomain(7, "y");
	myDomain(2) = quantity.Domain.defaultDomain(5, "t");
	myDomain(3) = quantity.Domain.defaultDomain(3, "asdf");
	quantity.Symbolic(sym("z"), myDomain);
	
	% the initialization above has to throw an exception. If the code is
	% here, the test is not valid.
	verifyTrue(testCase, false);
catch ex
	verifyTrue(testCase, isa(ex, "MException"));
end

end

function testCastSymbolic2Discrete(testCase)

syms z t
x = quantity.Symbolic(sin(z * t * pi), ...
	[quantity.Domain.defaultDomain(100, "z"), quantity.Domain("t", linspace(0, 2, 200))]);

X1 = quantity.Discrete(x);

%%
verifyEqual(testCase, X1.on, x.on);

end

function testDiff(testCase)
%%
z = sym("z", "real");
f1 = sinh(z * pi);
f2 = cosh(z * pi);

f = quantity.Symbolic([f1 ; f2 ],  quantity.Domain("z", linspace(0,1,11)), "name", "f");

d2f = f.diff("z", 2);
d1f = f.diff("z", 1);

%%
F1 = symfun([f1; f2], z);
D1F = symfun([diff(f1,1); diff(f2,1)], z);
D2F = symfun([diff(f1,2); diff(f2,2)], z);

z = f(1).domain.grid;

R1 = F1(z);
R2 = D1F(z);
R3 = D2F(z);

%%
verifyTrue(testCase, numeric.near(double([R1{:}]), f.on(), 1e-12));
verifyTrue(testCase, numeric.near(double([R2{:}]), f.diff("z", 1).on(), 1e-12));
verifyTrue(testCase, numeric.near(double([R3{:}]), f.diff("z", 2).on(), 1e-12));

end

function testInit(testCase)

z = linspace(0,1).';
t = linspace(0,1,3);

syms x y

v = [sin(x * y * pi); cos(x * y * pi)];

try
	
	myDomain(1) = quantity.Domain("z", z);
	myDomain(2) = quantity.Domain("t", t);
	
	b = quantity.Symbolic(v, myDomain);
	testCase.verifTrue(false);
end
	

myDom2(1) = quantity.Domain("x", z);
myDom2(2) = quantity.Domain("y", t);

b = quantity.Symbolic(v, myDom2);
c = quantity.Symbolic(v, myDom2);

%%
verifyTrue(testCase, all(b.on() == c.on(), "all"));
end

function testPlus(testCase)
syms z zeta

Z = quantity.Domain("z", linspace(0, 1, 3));
Zeta = quantity.Domain("zeta", linspace(0, 1, 4));

%% 1 + 2 variables
fSymVec3 = quantity.Symbolic([sinh(pi * z)], Z, "name", "sinh");
fSymVec4 = quantity.Symbolic([cosh(zeta * pi * z)], [Z Zeta], "name", "cosh");

result = fSymVec3 - fSymVec4;
resultDiscrete = quantity.Discrete(fSymVec3) - quantity.Discrete(fSymVec4);

%%
testCase.verifyEqual(result.on(), resultDiscrete.on(), "AbsTol", 10*eps);

%% quantity + constant

myQuan = quantity.Symbolic([sinh(pi * z), cosh(zeta * pi * z); 1, z^2], ...
	[Z Zeta], "name", "myQuan");
myQuanPlus1 = myQuan + ones(size(myQuan));
my1PlusQuan =  ones(size(myQuan)) + myQuan;

%%
testCase.verifyEqual(myQuan.on()+1, myQuanPlus1.on(), "AbsTol", 10*eps);
testCase.verifyEqual(myQuan.on()+1, my1PlusQuan.on(), "AbsTol", 10*eps);


%% 1 variable
fSymVec = quantity.Symbolic([sinh(z * pi) ; cosh(z * pi)], Z, "name", "sinhcosh");

F = fSymVec + fSymVec;
val = F.on();

f1fun = @(z) sinh(z * pi);
f2fun = @(z) cosh(z * pi) + 0*z;

zGrid = fSymVec(1).domain.grid;

%%
verifyTrue(testCase, numeric.near(2 * f1fun(zGrid) , val(:, 1), 1e-12));
verifyTrue(testCase, numeric.near(2 * f2fun(zGrid) , val(:, 2), 1e-12));

%% 2 variable
syms zeta
f1fun = @(z, zeta) sinh(z * pi);
f2fun = @(z, zeta) cosh(zeta * pi) + 0*z;

fSymVec = quantity.Symbolic([sinh(z * pi); cosh(zeta * pi)], [Z Zeta], "name", "sinhcosh");
fSymVec2 = [fSymVec(2); fSymVec(1)];

F = fSymVec + fSymVec;
F1Minus2 = fSymVec - fSymVec2;

zZetaNdgrid = fSymVec(1).domain.ndgrid();
zGrid = zZetaNdgrid{1}; zetaGrid = zZetaNdgrid{2};
%%
verifyTrue(testCase, numeric.near(2 * f1fun(zGrid, zetaGrid) , F(1).on(), 1e-12));
verifyTrue(testCase, numeric.near(2 * f2fun(zGrid, zetaGrid) , F(2).on(), 1e-12));
verifyTrue(testCase, numeric.near(f1fun(zGrid, zetaGrid) - f2fun(zGrid, zetaGrid), ...
	F1Minus2(1).on(), 1e-12));
verifyTrue(testCase, numeric.near(-F1Minus2(2).on(), ...
	F1Minus2(1).on(), 1e-12));
end

function testMTimes(testCase)
%%
syms z
f1 = sinh(z * pi);
f2 = cosh(z * pi);

f = quantity.Symbolic([f1 ; f2 ], quantity.Domain.defaultDomain(3, "z"), "name", "sinhcosh");

F = f * f.';
val = F.on();

f1 = matlabFunction(f1);
f2 = matlabFunction(f2);


z = f(1).domain.grid;

%%
verifyTrue(testCase, numeric.near(f1(z).^2, val(:, 1, 1), 1e-12));
verifyTrue(testCase, numeric.near(f2(z).^2, val(:, 2, 2), 1e-12));
verifyTrue(testCase, numeric.near(f1(z) .* f2(z) , val(:, 1, 2), 1e-12));


%% test multiplication with a scalar
ff = f * 2;
verifyTrue(testCase, numeric.near(ff.on(), [f1(z), f2(z)] * 2, 1e-12));


end

function testVariable(testCase)
syms z
q = quantity.Symbolic(z, quantity.Domain('z', 1));

verifyEqual(testCase, q.variable, z);

end

function testUMinusAndUPlus(testCase)
%%
syms z;
f1 = (sin(z .* pi) );
f2 =  z;
f3 = (cos(z .* pi) );
d = quantity.Domain("z", linspace(0, 1, 3).');

f = quantity.Symbolic([f1, f2; 0, f3], d);
fDisc = quantity.Discrete(f);
values = f.on();

%%
testCase.verifyEqual(on( -f ), -values);
testCase.verifyEqual(on( + ( - f ) ), -values );
testCase.verifyEqual(on( + fDisc), on(+f));
testCase.verifyEqual(on( - fDisc), on(-f));


end

function testConstantValues(testCase)
%% evaluation of a constant on the same grid
syms z
myDomain = quantity.Domain.defaultDomain(4, "z");
q = quantity.Symbolic(magic(3), myDomain);

magic7 = magic(3);
magic700 = zeros(4, 3, 3);
for k = 1:q(1).domain.gridLength()
	magic700(k, :,:) = magic7;
end

testCase.verifyEqual(magic700, q.on());

%% evaluation of a constant on a different grid
myDomain2 = quantity.Domain.defaultDomain(7, "z");
magic7 = magic(3);
magic700 = zeros(myDomain2.n, 3, 3);
for k = 1:myDomain2.n
	magic700(k, :,:) = magic7;
end
q.on(myDomain2);

testCase.verifyEqual(magic700, q.on(myDomain2));

%%
p = quantity.Symbolic([0 0 0 0 z], myDomain);
testCase.verifyFalse(p.isNumber());


%%
myDomains = [myDomain, myDomain.rename("t")];
Q = quantity.Symbolic(sym(magic(2)), myDomains);
Qon = Q.on();
magic2 = magic(2);

for k = 1:numel(magic2)