testDiscrete.m 46.6 KB
Newer Older
1
2
3
4
5
6
function [tests] = testDiscrete()
%testQuantity Summary of this function goes here
%   Detailed explanation goes here
tests = functiontests(localfunctions);
end

7
8
9
10
11
12
13
14
15

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

16
17
18
19
20
21
22
23
24
25
26
27
28
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

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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)

Jakob Gabriel's avatar
Jakob Gabriel committed
48
function testL2Norm(tc)
49
50
51
52
53

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

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

58
xNorm = sqrt(int(x.' * x, "z"));
Jakob Gabriel's avatar
Jakob Gabriel committed
59
tc.verifyEqual(MAX(abs(xNorm - x.l2norm('z'))), 0, 'AbsTol', 1e-12);
60

Jakob Gabriel's avatar
Jakob Gabriel committed
61
62
63
64
x = quantity.Discrete(x);
xNorm = sqrt(int(x.' * x, 'z'));
tc.verifyEqual(MAX(abs(xNorm - x.l2norm('z'))), 0, 'AbsTol', 1e-12);

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

Jakob Gabriel's avatar
Jakob Gabriel committed
67
68
end % testL2Norm()

69
70
function testCompose(testCase)

71
t = quantity.Domain('t', linspace(0, 2));
72
73
74
75
76
77
78
79
80
81
82
83
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 )
84
85
tau = quantity.Domain('tau', linspace(-1, 3));
z = quantity.Domain('z', linspace(0, 1, 51));
86
87
88
89
90
91
92
93
94

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)

95
96
97
98
99
% 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

100
101
102
z = quantity.Domain('z', linspace(0, 1, 11));
t = quantity.Domain('t', linspace(0, 1, 21));
tau = quantity.Domain('tau', linspace(0, 2, 31)); 
103
104
105
106
107
108
109
110

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)

111
112
113
114
115
116
117
118
119
120
% 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)
121
122
123
124
125
126
127
128
129

% 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)

130
131
end

132
133
134
135
136
137
138
139
140
141
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

142
143
144
145
146
147
148
149
150
151
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()

Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
152
153
154
function testSort(tc)

t = linspace(0, pi, 7)';
155
156
domA = quantity.Domain('a', t);
domB = quantity.Domain('b', t);
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
157
158
159
q = quantity.Discrete(sin(t) .* cos(t'), ...
	'domain', [domA, domB]);

160
q1V = q.valueDiscrete;
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
161
162
163

q.sort('descend');

164
q2V = q.valueDiscrete;
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
165
166
167
168
169
170

tc.verifyEqual(q1V , q2V')


end % testSort

171
function testBlkdiag(tc)
172
173
z = quantity.Domain('z', linspace(0, 1, 21));
zeta = quantity.Domain('zeta', linspace(0, 1, 41));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
174
175
176
177

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') ), ...
178
	'domain', [z, zeta], 'name', "q");
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
179

180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
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()

205
function testSetName(tc)
206
207
208
209
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");
210
end
211

212
function testCtranspose(tc)
213
214
215
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

216
217
qSymbolic = quantity.Symbolic(...
	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
218
	'gridName', {'z', 'zeta'}, 'name', "q");
219
220
221
222
223
224
225
226
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, ...
227
	'grid', linspace(0, 1, 11), 'gridName', "z", 'name', "im");
228
229
230
231
232
qDiscrete2Ctransp = qDiscrete2';
tc.verifyEqual(qDiscrete2.real.on(), qDiscrete2Ctransp.real.on());
tc.verifyEqual(qDiscrete2.imag.on(), -qDiscrete2Ctransp.imag.on());
end % testCtranspose

233
function testTranspose(tc)
234
235
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
236

237
238
qSymbolic = quantity.Symbolic(...
	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
239
	'gridName', {'z', 'zeta'}, 'name', "q");
240
241
qDiscrete = quantity.Discrete(qSymbolic);
qDiscreteTransp = qDiscrete.';
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
242

243
244
245
246
247
248
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

249
function testFlipGrid(tc)
250
251
252
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

253
254
myGrid = linspace(0, 1, 11);
f = quantity.Discrete(quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, ...
255
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
256
% % flip one grid
257
fReference = quantity.Discrete(quantity.Symbolic([1+z+(1-zeta); 2*(1-zeta)+sin(z)] + zeros(2, 1)*z*zeta, ...
258
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
259
fFlipped = f.flipGrid('zeta');
260
261
262
263
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, ...
264
	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
265
266
267
268
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);
269
end % testFlipGrid();
270
271

function testScalarPlusMinusQuantity(testCase)
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
272

273
z = quantity.Domain('z', linspace(0, 1, 7));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
274
275
276
f = quantity.Discrete( ( [2; 1] + zeros(2, 1) .* z.grid' )', ...
	'domain', z);

277
278
testCase.verifyError(@() 1-f-1, '');
testCase.verifyError(@() 1+f+1, '');
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
279

280
end % testScalarPlusMinusQuantity
281

282
function testNumericVectorPlusMinusQuantity(testCase)
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
283

284
z = quantity.Domain('z', linspace(0, 1, 7));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
285
286
287
f = quantity.Discrete( [1+z.grid, 2+sin(z.grid)] , ...
	'domain', z);

288
289
290
291
292
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);
293
end % testNumericVectorPlusMinusQuantity
294

295
function testUnitaryPluasAndMinus(testCase)
296
297
z = quantity.Domain('z', linspace(0, 1, 7));
zeta = quantity.Domain('zeta', linspace(0, 1, 7));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
298
299
300
301
302

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]) ;
303

Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
304
305
306
307
qDoubleArray = q.on();
	
testCase.verifyEqual(on(-q), -qDoubleArray);
testCase.verifyEqual(on(+q), +qDoubleArray);
308
309
310

end

311
312
function testConcatenate(testCase)

313
t = linspace(0, pi, 7)';
314
315
316
317
318
319
320
321
322
323
324
325
326
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');
327
AC = [A; C];
328

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

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

334
335
336
337
338
339
340
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);

341
342
343
344
345
346
347
348
349
350
351
352
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())
353
354

% concatenation of a constant and a function:
355
356
C = [5; 2];
const = quantity.Discrete(C, 'domain', quantity.Domain.empty);
357
% 
358
359
360
361
362
363
364
365
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)));

366

367
368
end

369
370
function testExp(testCase)
% 1 spatial variable
371
z = quantity.Domain('z', linspace(0, 1, 3));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
372
s1d = quantity.Discrete(1 + z.grid .* z.grid.', 'domain', z);
373
374
375
testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));

% diagonal matrix
376
zeta = quantity.Domain('zeta', linspace(0, 1, 4));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
377
378
379
380
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]);

381
382
383
384
testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
end

function testExpm(testCase)
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
385

386
387
z = quantity.Domain('z', linspace(0, 1, 3));
zeta = quantity.Domain('zeta', linspace(0, 1, 4));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
388
389
390
391
392
393

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]);

394
395
mat2dMat = mat2d.on();
mat2dExpm = 0 * mat2d.on();
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
396
397
for zIdx = 1 : z.n
	for zetaIdx = 1 : zeta.n
398
399
400
401
402
403
404
405
406
407
		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), ...
408
	'gridName', 'z', 'name', "s1d");
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
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
467
468
z = testCase.TestData.sym.z;

469
do = quantity.Domain('z', linspace(0,1,3));
470
myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], 'domain', do);
471
472
473
474
475
476
477
478
479
480
481
482
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
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
483

484
485
function testVec2Diag(testCase)
% quantity.Discrete
486
487
488
489
490
491
492
493
494
495
496
497
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) );
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
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
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
525
quan = quantity.Discrete(ones(3, 1), 'grid', linspace(0, 1, 3), ...
526
527
	'gridName', 'z', 'size', [1, 1], 'name', 'a');
solution = quan.solveDVariableEqualQuantity();
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
528
[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
529
530
531
532
testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
end

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

536
assume(z>0 & z<1); assume(zeta>0 & zeta<1);
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
537
myParameterGrid = linspace(0, 0.1, 5);
538
539
540
541
Lambda = quantity.Symbolic(-0.1-z^2, ...%, -1.2+z^2]),...1+z*sin(z)
	'grid', myParameterGrid, 'gridName', 'z', 'name', '\Lambda');

%%
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
542
543
myGrid = linspace(-.2, .2, 11);
sEval = -0.09;
544
545
546
547
548
solveThisOdeDiscrete = solveDVariableEqualQuantity(...
	diag2vec(quantity.Discrete(Lambda)), 'variableGrid', myGrid);
solveThisOdeSymbolic = solveDVariableEqualQuantity(...
	diag2vec(Lambda), 'variableGrid', myGrid);

Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
549
testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.05}), solveThisOdeDiscrete.on({sEval, 0.05}), 'RelTol', 5e-3);
550
551
552
553
end

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

556
assume(z>0 & z<1);
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
557
quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 5)}, ...
558
	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
559
quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 5)}, ...
560
561
562
563
564
565
566
567
568
	'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
569
570
z = testCase.TestData.sym.z;

571
assume(z>0 & z<1);
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
572
quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 11)}, ...
573
	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
574
quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 11)}, ...
575
	'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
576

577
578
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
myGrid = solutionBDiscrete.grid{1};
579
solutionBDiscreteDiff = solutionBDiscrete.diff('s');
580
581
582
583
584
quanOfSolutionOfS = zeros(length(myGrid), length(myGrid), 1);
for icIdx = 1 : length(myGrid)
	quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
end
%%
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
585
testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-2);
586
587

end
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
588

589
590
591
592
593
594
595
596
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;
%
597
598
z = testCase.TestData.sym.z;

599
c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'gridName', 'z', 'name', 'c');
600
601
602
603
604
605
606
607
608
609
610
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]);
611

612
613
constantDiff = diff(constant);
testCase.verifyEqual(constantDiff.on(), [zeros(11, 1), ones(11, 1)], 'AbsTol', 10*eps);
614
615
616
617
618
619
620
621

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));
622
623
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));
624

625
626
627
628
629
630
631
632
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]);
633
634
myQuantityDz = diff(myQuantity, 'z', 1);
myQuantityDzeta = diff(myQuantity, 'zeta', 1);
Jakob Gabriel's avatar
Jakob Gabriel committed
635

636
637
638
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' )
639
640
641
642

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

644
645
646
647
648
649
650
651

% 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);
652

653
654
655
end

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

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

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

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

676
testCase.verifyEqual(permute(createTestData(linspace(0, 1, 21), linspace(0, 1, 24)), [2, 1, 3, 4]), ...
677
	a.on( {linspace(0, 1, 24), linspace(0, 1, 21)}, {'zeta', 'z'}), 'AbsTol', 2e-4);
678
679
680
681
682
683
684
685
686
687

	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
688

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

694
695
696
697
698
699
700
701
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);
702

703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
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

721
722
723
724
725
726
727
728
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()

729
730
function testMtimesZZeta2x2(testCase)

Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
731
732
733
734


%%
z11_ = linspace(0, 1, 11);
735
736
z11 = quantity.Domain('z', z11_);
zeta11 = quantity.Domain('zeta', z11_);
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
737
738
739

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

741
zeta = testCase.TestData.sym.zeta;
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
742
743
z10_ = linspace(0, 1, 10);

744
b = quantity.Symbolic((eye(2, 2)), 'gridName', 'zeta', ...
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
745
746
	'grid', z10_);

747
748
c = a*b;

749
%
750
751
752
testCase.verifyEqual(c.on(), a.on());
end

753
function testMTimesPointWise(tc)
754

755
756
757
z = tc.TestData.sym.z;
zeta = tc.TestData.sym.zeta;

Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
758
759
760
Z = linspace(0, pi, 51)';
ZETA = linspace(0, pi / 2, 71)';

761
P = quantity.Symbolic( sin(zeta) * z, 'grid', {Z, ZETA}, 'gridName', {'z', 'zeta'});
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
762
B = quantity.Symbolic( sin(zeta), 'grid', ZETA, 'gridName', 'zeta');
763
764
765
766
767
768
769

PB = P*B;

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

pb = p*b;
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
770

771
tc.verifyEqual(MAX(abs(PB-pb)), 0, 'AbsTol', 10*eps);
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
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);
797
798
799
800

% 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]));

801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
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);
826
827
828

% 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]));
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
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)

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

859
860
t = testCase.TestData.sym.t;
sy = testCase.TestData.sym.sy;
861

862
a = [ 1, sy; t, 1];
863
864
865
866
867
868
869
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)
870
		A(:,:,i,j) = subs(a(i,j), {sy, t}, {S, T});
871
872
873
874
875
876
877
878
879
	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
880
881
v = int(a*b, sy, tGrid(1), t);
V = quantity.Symbolic(v, 'grid', tGrid, 'gridName', 't');
882
883

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

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

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

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

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

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

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

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

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

913
%% testCumIntWithLowerAndUpperBoundSpecified
914
915
916
tGrid = linspace(pi, 1.1*pi, 51)';
s = tGrid;
[T, S] = ndgrid(tGrid, tGrid);
917
918
919
920

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

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

923
924
925
926
A = zeros([size(T), size(a)]);

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

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

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

%% testCumInt with one independent variable
940
tGrid = linspace(0, pi, 51)';
941
942
943
944

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

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

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

950
K = quantity.Symbolic( int(a, 0, t), 'grid', {tGrid}, 'gridName', {'t'});
951
952
953
954

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


955
956
957
958
959
960
961
962
963
964
965
966
967
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));

968
969
970
971
972
973
974
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)];
Ferdinand Fischer's avatar
Ferdinand Fischer committed
975
testCase.verifyTrue(numeric.near(B_1_y, b_1_y));
976
977
978
979
980
981

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

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

982
983
end

984
function testMTimes(testCase)
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
% 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,:));
For faster browsing, not all history is shown. View entire blame