Discrete.m 91.8 KB
Newer Older
1001
					end
Jakob Gabriel's avatar
Jakob Gabriel committed
1002
					if isequal(values{1}, gridName2Replace{1})
1003
						% replace with same variable... everything stays the
Jakob Gabriel's avatar
Jakob Gabriel committed
1004
						% same.
1005
1006
						% Do not use "return", since, later subs might need to be
						% called recursively!
Jakob Gabriel's avatar
Jakob Gabriel committed
1007
1008
						newValue = obj.on();
					elseif any(strcmp(values{1}, obj(1).gridName))
1009
						% if for a quantity f(z, zeta) this method is
1010
1011
1012
						% called with subs(f, zeta, z), then g(z) = f(z, z)
						% results, hence the dimensions z and zeta are
						% merged.
1013
						domainIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
1014
							obj(1).domain.gridIndex(values{1})];
1015
1016
1017
1018
1019
1020
1021
1022
1023
						newDomainForOn = obj(1).domain;
						if obj(1).domain(domainIndices(1)).n > obj(1).domain(domainIndices(2)).n
							newDomainForOn(domainIndices(2)) = quantity.Domain(...
								newDomainForOn(domainIndices(2)).name, ...
								newDomainForOn(domainIndices(1)).grid);
						elseif  obj(1).domain(domainIndices(1)).n < obj(1).domain(domainIndices(2)).n
							newDomainForOn(domainIndices(1)) = quantity.Domain(...
								newDomainForOn(domainIndices(1)).name, ...
								newDomainForOn(domainIndices(2)).grid);
1024
						end
1025
1026
1027
						newValue = misc.diagNd(obj.on(newDomainForOn), domainIndices);
						newDomain = [newDomainForOn(domainIndices(2)), ...
							newDomainForOn(all(1:1:numel(newDomainForOn) ~= domainIndices(:)))];
1028
					else
1029
						% this is the default case. just grid name is changed.
1030
						%newDomain = obj(1).domain;
1031
1032
1033
						newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
							quantity.Domain(values{1}, ...
							obj(1).domain(obj(1).domain.gridIndex(gridName2Replace{1})).grid);
1034
1035
1036
1037
1038
						newValue = obj.on();
					end
					
				elseif isnumeric(values{1}) && numel(values{1}) == 1
					% if values{1} is a scalar, then obj is evaluated and
1039
					% the resulting quantity looses that spatial grid and
1040
					% gridName
1041
					newDomain = newDomain(~strcmp(gridName2Replace{1}, [newDomain.name]));
1042
1043
					% newGrid is the similar to the original grid, but the
					% grid of gridName2Replace is removed.
1044
					newGridSize = newDomain.gridLength();
1045
1046
1047
					% newGridForOn is the similar to the original grid, but
					% the grid of gridName2Replace is set to values{1} for
					% evaluation of obj.on().
1048
					newGridForOn = obj(1).grid;
1049
					newGridForOn{obj(1).domain.gridIndex(gridName2Replace{1})} = values{1};
1050
1051
1052
1053
1054
					newValue = reshape(obj.on(newGridForOn), [newGridSize, size(obj)]);
					
				elseif isnumeric(values{1}) && numel(values{1}) > 1
					% if values{1} is a double vector, then the grid is
					% replaced.
1055
1056
1057
					newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
						quantity.Domain(gridName2Replace{1}, values{1});
					newValue = obj.on(newDomain);
1058
1059
				else
					error('value must specify a gridName or a gridPoint');
1060
				end % if ischar(values{1}) || isstring(values{1})
1061
				if isempty(newDomain)
1062
1063
1064
					% TODO@Jakob: Is it possible that this case occurrs?
					% If I see it right, the only case where the newDomain is not set is the case
					% where a error is thrown.
1065
1066
					solution = newValue;
				else
1067
					solution = quantity.Discrete(newValue, newDomain, ...
1068
1069
1070
1071
1072
1073
1074
						'name', obj(1).name);
				end
				if numel(gridName2Replace) > 1
					solution = solution.subs(gridName2Replace(2:end), values(2:end));
				end
			end
		end
Jakob Gabriel's avatar
Jakob Gabriel committed
1075
		
1076
		function [idx, logicalIdx] = gridIndex(obj, varargin)
1077
			warning('DEPRICATED: use quantity.Domain.gridIndex method instead')
1078
			[~, idx, logicalIdx] = obj(1).domain.find(varargin{:});
1079
1080
		end 
		
1081
		function value = at(obj, point)
1082
			% at() evaluates the object at one point and returns it as array
1083
			% with the same size as size(obj).
1084
			value = reshape(obj.on(point), size(obj));
1085
		end % at()
1086
1087
		
		function value = atIndex(obj, varargin)
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1088
1089
1090
			% ATINDEX returns the valueDiscrete at the requested index.
			% value = atIndex(obj, varargin) returns the
			% quantity.Discrete.valueDiscrete at the index defined by
1091
			% varargin.
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1092
1093
1094
			%	value = atIndex(obj, 1) returns the first element of
			%	"valueDiscrete"
			%	value = atIndex(obj, ':') returns all elements of
1095
			%	obj.valueDiscrete in vectorized form.
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1096
1097
1098
1099
1100
			%	value = atIndex(obj, 1, end) returns the obj.valueDiscrete
			%	at the index (1, end).
			%	If a range of index is requested, the result is returned
			%	with the grids as indizes. If scalar values are requested,
			%	than the grid dimensions are neglected.
1101
			if nargin == 1
1102
				
1103
				if numel(obj(1).domain) == 1
1104
					value = zeros(obj(1).domain.gridLength, 1);
1105
				else
1106
					value = zeros(obj(1).domain.gridLength);
1107
				end
1108
1109
1110
1111
1112
1113
1114
				if isempty(value)
					value = 0;
				end
			else
				if ~iscell(varargin)
					varargin = {varargin};
				end
1115
				value = cellfun(@(v) v(varargin{:}), {obj.valueDiscrete}, ...
1116
					'UniformOutput', false);
1117
1118
				
				valueSize = size(value{1});
1119
				
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1120
1121
1122
1123
				if all(cellfun(@numel, varargin) == 1) && all(cellfun(@isnumeric, varargin))
					outputSize = [];
				else
					outputSize = valueSize(1:obj(1).nargin);
1124
1125
				end
				
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1126
				value = reshape([value{:}], [outputSize, size(obj)]);
1127
			end
1128
		end % atIndex()
1129
1130
1131
1132
		
		function n = nargin(obj)
			% FIXME: check if all funtions in this object have the same
			% number of input values.
1133
1134
1135
1136
1137
1138
			
			% FIXME: for some combinations of constant objects, it seems to be
			% possible, that the quantity has a gridName but no grid.
			% Actually this should not be allowed. This is quick and dirty
			% work around.
			n = min(numel(obj(1).gridName), numel(obj(1).grid));
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
		end
		
		function d = gridDiff(obj)
			
			% #FIXME:
			%   1) test for multidimensional grids
			%   2) check that the grid is equally spaced
			
			d = diff(obj(1).grid{:});
			d = d(1);
		end
		
Jakob Gabriel's avatar
Jakob Gabriel committed
1151
		function s = gridSize(obj, myGridName)
1152
			% GRIDSIZE returns the size of all grid entries.
1153
			% todo: this should be called gridLength
1154
			warning('DEPRICATED: use quantity.Domain.gridLength method instead')
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1155
			if isempty(obj(1).domain)
1156
1157
				s = [];
			else
Jakob Gabriel's avatar
Jakob Gabriel committed
1158
1159
1160
1161
1162
				if nargin == 1
					s = cellfun('length', obj(1).grid);
				elseif nargin == 2
					s = cellfun('length', obj(1).gridOf(myGridName));
				end
1163
1164
1165
			end
		end
		
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1166
1167
		function H = plot(obj, varargin)
			H = [];
1168
1169
			p = misc.Parser();
			p.addParameter('fig', []);
1170
1171
			p.addParameter('dock', quantity.Settings.instance().plot.dock);
			p.addParameter('showTitle', quantity.Settings.instance().plot.showTitle);
Jakob Gabriel's avatar
Jakob Gabriel committed
1172
			p.addParameter('titleWithIndex', true');
1173
			p.addParameter('hold', false);
1174
			p.addParameter('export', false);
1175
			p.addParameter('currentFigure', false);
1176
1177
1178
1179
1180
			p.addParameter('exportOptions',  ...
				{'height', [num2str(0.25*size(obj, 1)), '\textwidth'], ...
				'width', '0.8\textwidth', 'externalData', false, ...
				'showWarnings', false, 'showInfo', false, ...
				'extraAxisOptions', 'every axis title/.append style={yshift=-1.5ex}, every axis x label/.append style={yshift=2mm}'});
1181
			p.parse(varargin{:});
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1182
			additionalPlotOptions = misc.struct2namevaluepair(p.Unmatched);
1183
			if prod(obj(1).domain.gridLength) == 1
1184
1185
1186
				additionalPlotOptions = [additionalPlotOptions(:)', ...
					{'x'}];
			end
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1187
			
1188
1189
1190
			fig = p.Results.fig;
			dock = p.Results.dock;
			for figureIdx = 1:size(obj, 3)
1191
1192
1193
				if p.Results.currentFigure
					h = gcf;
				elseif isempty(p.Results.fig)
1194
					h = figure();
1195
1196
				elseif p.Results.fig == 0
					h = gcf;
1197
				else
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1198
					h = figure(fig + figureIdx - 1);
1199
				end
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1200
				H = [H, h];
1201
1202
1203
1204
1205
				
				if dock
					set(h, 'WindowStyle', 'Docked');
				end
				
1206
				assert(~isempty(obj), 'Empty quantities can not be plotted');
Jakob Gabriel's avatar
Jakob Gabriel committed
1207
				assert(obj.nargin() <= 2, 'plot only supports quantities with 2 gridNames');
1208
1209
1210
1211
1212
1213
1214
1215
1216
				
				subplotRowIdx = 1:size(obj, 1);
				subpotColumnIdx = 1:size(obj, 2);
				
				i = 1: numel(obj(:,:,figureIdx));
				i = reshape(i, size(obj, 2), size(obj, 1))';
				
				for rowIdx = subplotRowIdx
					for columnIdx = subpotColumnIdx
1217
1218
1219
						if ~p.Results.currentFigure
							subplot(size(obj, 1), size(obj, 2), i(rowIdx, columnIdx));
						end
1220
1221
						if p.Results.hold
							hold on;
1222
						else
1223
1224
							hold off;
						end
1225
						
1226
1227
1228
						if isempty(obj(rowIdx, columnIdx, figureIdx))
							warning('you are trying to plot an empty quantity');
						elseif obj.nargin() == 0
1229
1230
1231
1232
							plot(0, ...
								obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
								additionalPlotOptions{:});
						elseif obj.nargin() == 1
1233
							plot(...
Jakob Gabriel's avatar
Jakob Gabriel committed
1234
								obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1235
1236
								obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
								additionalPlotOptions{:});
1237
1238
1239
						elseif obj.nargin() == 2
							misc.isurf(obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
								obj(rowIdx, columnIdx, figureIdx).grid{2}(:), ...
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1240
1241
								obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
								additionalPlotOptions{:});
1242
1243
1244
1245
1246
							ylabel(labelHelper(2), 'Interpreter','latex');
						else
							error('number inputs not supported');
						end
						xlabel(labelHelper(1), 'Interpreter','latex');
1247
1248
						
						if p.Results.showTitle
1249
							title(titleHelper(), 'Interpreter','latex');
1250
						end
1251
1252
1253
						a = gca();
						a.TickLabelInterpreter = 'latex';
						
1254
1255
					end % for columnIdx = subpotColumnIdx
				end % for rowIdx = subplotRowIdx
1256
				
1257
			end % for figureIdx = 1:size(obj, 3)
1258
			
1259
1260
			if p.Results.export
				matlab2tikz(p.Results.exportOptions{:});
1261
			end
1262
			
1263
			function myLabel = labelHelper(gridNumber)
1264
				if ~isempty(obj(rowIdx, columnIdx, figureIdx).gridName)
1265
					myLabel = "$$" + greek2tex(obj(rowIdx, columnIdx, figureIdx).gridName{gridNumber}) + "$$";
1266
				else
1267
					myLabel = "";
1268
				end
1269
			end % labelHelper()
1270
			function myTitle = titleHelper()
Jakob Gabriel's avatar
Jakob Gabriel committed
1271
				if numel(obj) <= 1 || ~p.Results.titleWithIndex
1272
					myTitle = "$${" + greek2tex(obj(rowIdx, columnIdx, figureIdx).name) + "}$$";
Jakob Gabriel's avatar
Jakob Gabriel committed
1273
				elseif ndims(obj) <= 2
1274
1275
					myTitle = "$$[{" + greek2tex(obj(rowIdx, columnIdx, figureIdx).name) ...
						+ "]}_{" + num2str(rowIdx) + num2str(columnIdx) + "}$$";
1276
				else
1277
1278
					myTitle = "$${[" + greek2tex(obj(rowIdx, columnIdx, figureIdx).name) ...
						+ "]}_{" + num2str(rowIdx) + num2str(columnIdx) + num2str(figureIdx) + "}$$";
1279
				end
1280
			end % titleHelper()
1281
			function myText = greek2tex(myText)
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
				if ~contains(myText, "\")
					myText = strrep(myText, "Lambda", "\Lambda");
					myText = strrep(myText, "lambda", "\lambda");
					myText = strrep(myText, "Zeta", "\Zeta");
					myText = strrep(myText, "zeta", "\zeta");
					myText = strrep(myText, "Gamma", "\Gamma");
					myText = strrep(myText, "gamma", "\gamma");
					myText = strrep(myText, "psi", "\psi");
					myText = strrep(myText, "phi", "\phi");
					myText = strrep(myText, "Psi", "\Psi");
					myText = strrep(myText, "Phi", "\Phi");
					myText = strrep(myText, "Delta", "\Delta");
					myText = strrep(myText, "delta", "\delta");
					if ~contains(myText, "\zeta") && ~contains(myText, "\Zeta")
						myText = strrep(myText, "eta", "\eta");
1297
					end
1298
1299
					myText = strrep(myText, "pi", "\pi");
					myText = strrep(myText, "Pi", "\Pi");
1300
				end
1301
1302
			end % greek2tex()
		end % plot()
1303
		
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
		function [s] = optArgList(obj, exclude)
			% OPTARGLIST returns optional arguments of the quantity.Discrete
			%	S = optArgList(OBJ, EXCLUDE) returns the optional arguments for the OBJ as
			%	cell-array. With EXCLUDE as string-array, properties can be excluded
			%		
			assert(numel(obj) == 1, "argList must NOT be called for an array object");
			s = struct(obj);
			s = rmfield(s, ["domain", "grid", "gridName"]);
			if nargin == 2
				s = rmfield(s, exclude);
			end
			s = misc.struct2namevaluepair(s);
		end
		
1318
		function s = nameValuePair(obj, varargin)
1319
1320
1321
1322
			% NAMEVALUEPAIR returns a name value pair list for the object
			%	S = nameValuePair(OBJ, VARARGIN) returns the properties of the object as
			%	name-value-pair.
			warning("DEPRECATED! Use quantity.Discrete/optArgList insted")
1323
			assert(numel(obj) == 1, "nameValuePair must NOT be called for an array object");
1324
1325
1326
1327
1328
			s = struct(obj);
			if ~isempty(varargin)
				s = rmfield(s, varargin{:});
			end
			s = misc.struct2namevaluepair(s);
1329
		end % nameValuePair()
1330
		
1331
		function s = struct(obj, varargin)
1332
			if (nargin == 1) && isa(obj, "quantity.Discrete")
1333
				tempProperties = fieldnames(obj);
1334
1335
1336
				si = num2cell( size(obj) );
				s(si{:}) = struct();
				for l = 1:numel(obj)
1337
1338
					for k = 1:length(tempProperties)
						s(l).(tempProperties{k}) = obj(1).(tempProperties{k});
1339
1340
					end
				end
1341
			else
1342
				% Without this else-case
1343
1344
1345
1346
1347
				% this method is in conflict with the default matlab built-in
				% calls struct(). When calling struct('myFieldName', myQuantity) this
				% quantity-method is called (and errors) and not the
				% struct-calls-constructor.
				s = builtin('struct', obj, varargin{:});
1348
			end
1349
		end % struct()
1350
1351
1352
1353
1354
1355
		
		function s = obj2struct(obj)
			warning('depricated');
			s = struct(obj);
		end
		
1356
		function b = flipGrid(a, myGridName)
1357
1358
1359
1360
			% flipGrid() implements a flip of the grids, for example with a(z)
			% b(z) = a(1-z) = a.flipGrid('z');
			% or for the multidimensional case with c(z, zeta)
			% d(z, zeta) = c(1-z, 1-zeta) = c.flipGrid({'z', 'zeta'});
1361
1362
1363
1364
			bMat = a.on();
			if ~iscell(myGridName)
				myGridName = {myGridName};
			end
1365
			gridIdx = a(1).domain.gridIndex(myGridName);
1366
1367
1368
			for it = 1 : numel(myGridName)
				bMat = flip(bMat, gridIdx(it));
			end
1369
			b = quantity.Discrete(bMat, a(1).domain, ...
1370
				'name', "flip(" + a(1).name + ")");
1371
1372
		end % flipGrid()
		
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
		function newObj = replaceGrid(obj, myNewDomain, optArgs)
			% CHANGEGRID change the grid of the quantity.
			%	newObj = REPLACEGRID(obj, MYNEWDOMAIN, "gridName", NEWGRIDNAME)
			% replace the grid of the obj quantity. The order of grid and
			% gridName in the obj properties remains unchanged, only the
			% grid points are exchanged.
			%
			% TODO:
			%	newObj = REPLACEGRID(obj, domain) changes the domain of the
			%	object specified by the name of DOMAIN into DOMAIN.
			%
			%	newObj = REPLACEGRID(obj, domain, 'gridName', NEWNAME) changes the domain of the
			%	object specified by NEWNAME into DOMAIN. 
			
			arguments
				obj
				myNewDomain quantity.Domain
				optArgs.gridName = [myNewDomain.name];
			end
			
			if isempty(obj)
				newObj = obj.copy();
				return;
			end
			
			assert( intersect(obj(1).gridName, optArgs.gridName) == optArgs.gridName );
			
			if obj(1).isConstant
				error("Not yet implemented")
			else
				gridIndexNew = obj(1).domain.gridIndex(optArgs.gridName);
				% initialization of the newDomain array as quantity.Domain
				% array. This is required in order to handle also
				% quantity.EquidistantDomains:
				newDomain(1:obj(1).nargin) = quantity.Domain();
				newDomain(:) = obj(1).domain;

				for it = 1 : length(gridIndexNew)
					newDomain(gridIndexNew(it)) = ...
						quantity.Domain(myNewDomain(it).name, myNewDomain(it).grid);
				end
			end
			
			newObj = obj.copy();
			[newObj.domain] = deal(newDomain);
		end % replaceGrid
1419
1420
1421
1422
	end
	
	%% math
	methods (Access = public)
1423
		
1424
		function [F] = stateTransitionMatrix(obj, varargin)
1425
1426
1427
1428
1429
			% STATETRANSITIONMATRIX compute the state-transition-matrix
			%	[F0] = stateTransitionMatrix(obj, varargin) computation of the
			% state-transition matrix for a ordinary differential equation
			%	d / dt x(t) = A(t) x(t)
			% with A as this object. The system matrix can be constant or
1430
			% variable. It is the solution of
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
			%	d / dt F(t,t0) = A(t) F(t,t0)
			%		   F(t,t) = I.
			% Optional parameters for the computation can be defined as
			% name-value-pairs. The options are:
			%	# 'grid' : the grid for which the state-transition matrix
			%	should be computed.
			%	# 'gridName1' : the name of the first independent variable.
			%	For F(t,t0) this is t. The default is the first gridName of
			%	this object.
			%	# 'gridName2' : the name of the second independent
			%	variable. For F(t,t0) this is t0. The default is the first
			%	gridName + '0'.
1443
1444
1445
1446
1447
1448
1449
			
			assert(size(obj,1) == size(obj,2), 'The quantity must be a quadratic matrix for the computation of state-transition matrices');
			assert(obj(1).nargin <= 1, 'Computation of state-transition matrix is only defined for quantities with one independent variable');
			
			ip = misc.Parser();
			ip.addParameter('grid', []);
			ip.addParameter('gridName1', '');
1450
			ip.addParameter('gridName2', '');
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
			ip.parse(varargin{:});
			myGrid = ip.Results.grid;
			gridName1 = ip.Results.gridName1;
			gridName2 = ip.Results.gridName2;
			
			if ip.isDefault('grid')
				assert(numel(obj(1).grid) == 1, 'No grid is defined for the computation of the state-transition matrix. Use the name-value pair option "grid" to define the grid.');
				myGrid = obj(1).grid{1};
			end
			
			if ip.isDefault('gridName1')
				assert(numel(obj(1).gridName) == 1, 'No gridName is defined. Use name-value-pairs property "gridName1" to define a grid name');
				gridName1 = obj(1).gridName{1};
			end
1465
			
1466
1467
1468
1469
1470
1471
			if ip.isDefault('gridName2')
				gridName2 = [gridName1 '0'];
			end
			
			
			assert(numel(myGrid) > 1, 'If the state transition matrix is computed for constant values, a spatial domain has to be defined!')
1472
			
1473
			myDomain = [quantity.Domain(gridName1, myGrid), ...
1474
				quantity.Domain(gridName2, myGrid)];
1475
			
1476
1477
1478
1479
1480
1481
			if obj.isConstant
				% for a constant system matrix, the matrix exponential
				% function can be used.
				z = sym(gridName1, 'real');
				zeta = sym(gridName2, 'real');
				f0 = expm(obj.atIndex(1)*(z - zeta));
1482
				F = quantity.Symbolic(f0, myDomain);
1483
				
1484
1485
			elseif isa(obj, 'quantity.Symbolic')
				f0 = misc.fundamentalMatrix.odeSolver_par(obj.function_handle, myGrid);
1486
1487
				F = quantity.Discrete(f0, myDomain, ...
					'size', [size(obj, 1), size(obj, 2)]);
1488
1489
1490
1491
			else
				f0 = misc.fundamentalMatrix.odeSolver_par( ...
					obj.on(myGrid), ...
					myGrid );
1492
1493
				F = quantity.Discrete(f0, myDomain, ...
					'size', [size(obj, 1), size(obj, 2)]);
1494
			end
1495
		end
1496
		
1497
		function s = sum(obj, dim)
1498
1499
1500
1501
1502
1503
1504
			% sum Sum of elements
			% s = sum(X) is the sum of all elements of the array X.
			% s = sum(X, DIM) is the sum along the dimensions specified by DIM.
			arguments
				obj;
				dim = 1:ndims(obj);
			end % arguments
1505
			s = quantity.Discrete(sum(obj.on(), obj.nargin + dim), obj(1).domain, ...
1506
1507
				'name', "sum(" + obj(1).name + ")");
		end % sum()
1508
		
1509
1510
1511
1512
1513
1514
1515
1516
		function P = prod(obj, dim)
			% prod Product of elements
			% P = prod(X) is the product of all elements of the vector X.
			% P = prod(X, DIM) is the product along the dimensions specified by DIM.
			arguments
				obj;
				dim = 1:ndims(obj);
			end % arguments
1517
			P = quantity.Discrete(prod(obj.on(), obj(1).nargin + dim), obj(1).domain, ...
1518
1519
1520
				'name', "prod(" + obj(1).name + ")");
		end % prod()
		
1521
1522
		function y = sqrt(x)
			% quadratic root for scalar and diagonal quantities
1523
1524
			y = quantity.Discrete(sqrt(x.on()), x(1).domain, ...
                    'name', "sqrt(" + x(1).name + ")");
1525
		end % sqrt()
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
		
		
		function y = sqrtm(x)
			% quadratic root for matrices
			if isscalar(x)
				% use sqrt(), because its faster.
				y = sqrt(x);
			elseif (size(x, 1) == size(x, 2)) && ismatrix(x)
				% implementation of quadratic root pointwise in space in a
				% simple for-loop.
				xMat = x.on();
				permuteGridAndIdx = [[-1, 0] + ndims(xMat), 1:x.nargin];
				permuteBack = [(1:x.nargin)+2, [1, 2]];
				xPermuted = permute(xMat, permuteGridAndIdx);
				yUnmuted = 0*xPermuted;
1541
				for k = 1 : prod(x(1).domain.gridLength)
1542
1543
					yUnmuted(:,:,k) = sqrt(xPermuted(:,:,k));
				end
1544
1545
				y = quantity.Discrete(permute(yUnmuted, permuteBack), x(1).domain, ...
					'size', size(x), 'name', "sqrtm(" + x(1).name + ")");
1546
1547
1548
			else
				error('sqrtm() is only implemented for quadratic matrices');
			end
1549
		end % sqrtm()
1550
		function P = power(obj, p)
Jakob Gabriel's avatar
Jakob Gabriel committed
1551
			% a.^p implemented by multiplication
1552
1553
			% #todo unittest required
			P = quantity.Discrete.zeros(size(obj), obj(1).domain);
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
			
			if round(p) == p
				% power of integer order:
				for k = 1:numel(obj)
					P(k) = obj(k)^p;
				end
			else
				% power of not integer order
				for k = 1:numel(obj)
					P(k) = obj(k).copy();
					P(k).valueDiscrete = obj(k).valueDiscrete.^p;
					P(k).name = obj(k).name + ".^" + num2str(p);
				end
1567
1568
			end
		end
1569
		function P = mpower(a, p)
1570
			% Matrix power a^p is matrix or scalar a to the power p.
1571
			if p == 0
1572
				P = setName(eye(size(a)) + 0*a, "I");
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
			elseif p < 0
				if numel(a) > 1
					warning("mpower(a, p) implements  a^p. " ...
						+ "For matrices a with negative exponent p, inv(a^(-p)) is returned. " ...
						+ "This represents a division from left, " ...
						+ "maybe division from right is needed in your case!")
				end % this warning is not important in the scalar case.
				P = inv(mpower(a, -p));
			else % p > 0
				assert(p==floor(p), "power p must be integer");
				P = a;
				for k = 1:(p-1)
					P = P * a;
				end
1587
			end
1588
		end % mpower()
1589
1590
		function s = num2str(obj)
			s = obj.name;
1591
1592
		end % num2str()
		
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
		function P = times(a, b)
			
			assert( all( size(a) == size(b), 'all' ), 'the terms a and b must have the same size');
			
			for i = 1:numel(a)
				P(i) = a(i) * b(i);
			end
			
			P = reshape(P, size(a));
			
		end
		
1605
		function P = mtimes(a, b)
1606
1607
			% TODO rewrite the selection of the special cases! the
			% if-then-cosntruct is pretty ugly!
1608
			
1609
			% numeric computation of the matrix product
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
			if isempty(b) || isempty(a)
				% TODO: actually this only holds for multiplication of
				% matrices. If higher dimensional arrays are multiplied
				% this can lead to wrong results.
				sa = size(a);
				sb = size(b);
				
				P = quantity.Discrete.empty([sa(1:end-1) sb(2:end)]);
				return
			end
1620
1621
			if isa(b, 'double')
				if numel(b) == 1
1622
					% simple multiplication in scalar case
1623
1624
					P = quantity.Discrete(a.on() * b, a(1).domain, ...
                        'size', size(a), 'name', a(1).name + num2str(b));
1625
					return
1626
				else
1627
					b = quantity.Discrete(b, quantity.Domain.empty(), 'size', size(b));
1628
1629
1630
1631
				end
			end
			
			if isa(a, 'double')
1632
				P = (b.' * a.').';
Jakob Gabriel's avatar
Jakob Gabriel committed
1633
1634
				% this recursion is safe, because isa(b, 'double') is considered in
				% the if above.
1635
1636
1637
1638
				
				if isnumeric(P)
					return
				else
1639
					P.setName("c " + b(1).name);
1640
1641
					return
				end
1642
			end
Jakob Gabriel's avatar
Jakob Gabriel committed
1643
			if a.isConstant() && ~b.isConstant()
1644
1645
1646
1647
1648
1649
				% If the first argument a is constant value, then bad
				% things will happen. To avoid this, we calculate
				%	a * b = (b' * a')'
				% instead. Thus we have to exchange both tranposed values
				% transpose the result in the end.
				P = (b' * a')';
1650
				P.setName(a(1).name + " " + b(1).name);
1651
				return
1652
1653
			elseif a.isConstant() && b.isConstant() ...
					&& isempty( a(1).domain ) && isempty(b(1).domain)
1654
1655
				P = a.on() * b.on();
				return 
1656
			end
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1657
			if isscalar(b)
1658
				% do the scalar multiplication in a loop
1659
				for k = 1:numel(a)
1660
					P(k) = innerMTimes(a(k), b);
1661
				end
1662
1663
				P = reshape(P, size(a));
				return
1664
1665
			end
			
Ferdinand Fischer's avatar
Ferdinand Fischer committed
1666
			if isscalar(a)
1667
				% do the scalar multiplication in a loop
Jakob Gabriel's avatar
Jakob Gabriel committed
1668
				for k = 1:numel(b)
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1669
					P(k) = innerMTimes(a, b(k)); %#ok<AGROW>
1670
				end
Jakob Gabriel's avatar
Jakob Gabriel committed
1671
				P = reshape(P, size(b));
Jakob Gabriel's avatar
Jakob Gabriel committed
1672
1673
				return
			end
1674
1675
			
			P = innerMTimes(a, b);
Jakob Gabriel's avatar
Jakob Gabriel committed
1676
		end
1677
1678
		
		function P = innerMTimes(a, b)
1679
1680
			assert(size(a, 2) == size(b, 1), ['For multiplication the ', ...
				'number of columns of the left array', ...
1681
				'must be equal to number of rows of right array'])
1682
1683
1684
1685
			
			% misc.multArray is very efficient, but requires that the
			% multiple dimensions of the input are correctly arranged.
			[idx, permuteGrid] = computePermutationVectors(a, b);
1686
			
1687
			parameters.name = a(1).name + " " + b(1).name;
1688
1689
1690
1691
			parameters.figureID = a(1).figureID;
			
			domainA = a(1).domain;
			domainB = b(1).domain;
1692
			
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1693
1694
1695
1696
1697
1698
			% The joined domain of quantity A and B is a mixture from both
			% domains. If the domains have the same name, it must be
			% checked if they have the same discretization. This is done by
			% the function "gridJoin". It returns the finer grid.
			% If the domain names do not coincide, they are just
			% appended. At first, the domains of qunatity A, then the
1699
			% domains of quantity B.
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1700
1701
			joinedDomain = [ gridJoin(domainA(idx.A.common), domainB(idx.B.common)), ...
				domainA(~idx.A.common), domainB(~idx.B.common) ];
1702
1703
1704
1705
			
			%% generate the new grids
			newDomainA = quantity.Domain.empty(); % for later call of a.on() with fine grid
			newDomainB = quantity.Domain.empty(); % for later call of b.on() with fine grid
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1706
			for i = 1 : numel(joinedDomain)
1707
1708
1709
1710
				
				% if there is a component of the original domain in the
				% joined domain, it can happen the length of a domain has
				% changed. Thus, it must be updated for a later evaluation.
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1711
1712
				idxA = domainA.gridIndex(joinedDomain(i).name);
				idxB = domainB.gridIndex(joinedDomain(i).name);
1713
1714
				
				if idxA
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1715
					newDomainA = [newDomainA, joinedDomain(i)]; %#ok<AGROW>
1716
				end
1717
1718
				
				if idxB
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1719
					newDomainB = [newDomainB, joinedDomain(i)]; %#ok<AGROW>
1720
				end
1721
1722
1723
			end
			parameters = misc.struct2namevaluepair(parameters);
			
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1724
			% evaluation of the quantities on their "maybe" new domain.
1725
1726
			valueA = a.on(newDomainA);
			valueB = b.on(newDomainB);
1727
			
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1728
1729
			% do the multidimensional tensor multiplication and permute the
			% values to the right order 
1730
1731
			C = misc.multArray(valueA, valueB, idx.A.value(end), idx.B.value(1), idx.common);
			C = permute(C, permuteGrid);
1732
			P = quantity.Discrete(C, joinedDomain, parameters{:});
1733
1734
1735
1736
1737
1738
1739
1740
1741
		end
		
		function y = inv(obj)
			% inv inverts the matrix obj at every point of the domain.
			assert(ismatrix(obj) && (size(obj, 1) == size(obj, 2)), ...
				'obj to be inverted must be quadratic');
			objDiscreteOriginal = obj.on();
			if isscalar(obj)
				% use ./ for scalar case
1742
				y = quantity.Discrete(1 ./ objDiscreteOriginal, obj(1).domain, ...
1743
					'name', "(" + obj(1).name + ")^{-1}");
1744
			else
1745
1746
				% reshape and permute objDiscrete such that only on for
				% loop is needed.
1747
				objDiscreteReshaped = permute(reshape(objDiscreteOriginal, ...
1748
1749
					[prod(obj(1).domain.gridLength), size(obj)]), [2, 3, 1]);
				invDiscrete = zeros([prod(obj(1).domain.gridLength), size(obj)]);
1750
				
1751
1752
1753
				parfor it = 1 : size(invDiscrete, 1)
					invDiscrete(it, :, :) = inv(objDiscreteReshaped(:, :, it));
				end
1754
				
1755
1756
				y = quantity.Discrete(reshape(invDiscrete, size(objDiscreteOriginal)), ...
                    obj(1).domain, 'name', "{(" + obj(1).name + ")}^{-1}");
1757
			end
1758
		end % inv()
1759
		
1760
1761
		function objT = transpose(obj)
			objT = builtin('transpose', copy(obj));
1762
1763
1764
1765
1766
			if ~isempty(obj)
				objT.setName("{" + obj(1).name + "}^{T}");
			else
				objT.setName("{.}^{T}");
			end
1767
1768
		end % transpose(obj)
		
1769
1770
		function objCt = ctranspose(obj)
			objT = obj.';
1771
1772
1773
1774
			
			if ~isempty(objT)
			
				objCtMat = conj(objT.on());
1775
				objCt = quantity.Discrete(objCtMat, obj(1).domain, ...
1776
1777
1778
1779
					'name', "{" + obj(1).name + "}^{H}");
			else
				objCt = objT;
			end
1780
		end % ctranspose(obj)
1781
		
1782
1783
		function y = exp(obj)
			% exp() is the exponential function using obj as the exponent.
1784
			y = quantity.Discrete(exp(obj.on()), obj(1).domain, ...
1785
				'name', "exp(" + obj(1).name + ")");
1786
		end % exp()
1787
		
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
		function y = log(obj)
		%  log    Natural logarithm.
		%     log(X) is the natural logarithm of the elements of X.
		%     Complex results are produced if X is not positive.
			y = quantity.Discrete(log(obj.on()), ...
				'domain', obj(1).domain, ...
				'name', "log(" + obj(1).name + ")");
		end % log()
		
		function y = log10(obj)
		% log10  Common (base 10) logarithm.
		% log10(X) is the base 10 logarithm of the elements of X.   
		% Complex results are produced if X is not positive.
			y = quantity.Discrete(log10(obj.on()), ...
				'domain', obj(1).domain, ...
				'name', "log10(" + obj(1).name + ")");
		end % log()
		
1806
		function xNorm = l2norm(obj, integralGridName, optArg)
Jakob Gabriel's avatar
Jakob Gabriel committed
1807
1808
1809
1810
1811
			% calculates the l2 norm, for instance
			%	xNorm = sqrt(int_0^1 x.' * x dz).
			% Optionally, a weight can be defined, such that instead
			%	xNorm = sqrt(int_0^1 x.' * weight * x dz).
			% The integral domain is specified by integralGrid.
1812
1813
			arguments
				obj;
1814
				integralGridName {mustBe.gridName} = [obj(1).domain.name];
1815
1816
1817
				optArg.weight = eye(size(obj, 1));
			end
			
1818
			integralGridName = misc.ensureString(integralGridName);
Jakob Gabriel's avatar
Jakob Gabriel committed
1819
			
1820
			if obj.nargin == 1 && all(strcmp(obj(1).gridName, integralGridName))
1821
1822
1823
1824
1825
1826
1827
1828
				
				integratedValue = int(obj.' * optArg.weight * obj);
				
				if isnumeric(integratedValue)
					xNorm = sqrtm( integratedValue );
				else
					xNorm = sqrtm(on(integratedValue, integralGridName));
				end
1829
			else
1830
				xNorm = sqrtm( int(obj.' * optArg.weight * obj, integralGridName) );
1831
				if isa(xNorm, 'quantity.Discrete')
1832
					xNorm = xNorm.setName("||" + obj(1).name + "||_{L2}");
1833
				end
1834
			end
1835
1836
1837
1838
1839
1840
1841
		end % l2norm()
		
		function xNorm = quadraticNorm(obj, varargin)
			% calculates the quadratic norm, i.e,
			%	xNorm = sqrt(x.' * x).
			% Optionally, a weight can be defined, such that instead
			%	xNorm = sqrt(x.' * weight * x).
1842
			
1843
1844
1845
1846
			myParser = misc.Parser();
			myParser.addParameter('weight', eye(size(obj, 1)));
			myParser.parse(varargin{:});
			xNorm = sqrtm(obj.' * myParser.Results.weight * obj);
1847
			xNorm.setName("||" + obj(1).name + "||_{2}");
1848
		end % quadraticNorm()
Jakob Gabriel's avatar
Jakob Gabriel committed
1849
		
1850
		function y = expm(x)
1851
1852
			% exp() is the matrix-exponential function using obj as the
			% exponent.
1853
1854
1855
1856
			if isscalar(x)
				% use exp(), because its faster.
				y = exp(x);
			elseif ismatrix(x)
1857
1858
				% implementation of expm pointwise in space in a simple
				% for-loop.
1859
1860
1861
1862
1863
				xMat = x.on();
				permuteGridAndIdx = [[-1, 0] + ndims(xMat), 1:x.nargin];
				permuteBack = [(1:x.nargin)+2, [1, 2]];
				xPermuted = permute(xMat, permuteGridAndIdx);
				yUnmuted = 0*xPermuted;
1864
				for k = 1 : prod(x(1).domain.gridLength)
1865
1866
					yUnmuted(:,:,k) = expm(xPermuted(:,:,k));
				end
1867
1868
				y = quantity.Discrete(permute(yUnmuted, permuteBack), x(1).domain, ...
					'size', size(x), 'name', "expm(" + x(1).name + ")");
1869
1870
1871
1872
			end
		end
		
		function x = mldivide(A, B)
1873
1874
1875
			% mldivide see doc mldivide of matlab: "x = A\B is the solution
			% to the equation Ax = B. Matrices A and B must have the same
			% number of rows."
1876
			x = inv(A) * B;
1877
		end
1878
		
1879
		function x = mrdivide(B, A)
1880
1881
1882
			% mRdivide see doc mrdivide of matlab: "x = B/A solves the
			% system of linear equations x*A = B for x. The matrices A and
			% B must contain the same number of columns"
1883
			x = B * inv(A);
1884
1885
1886
1887
		end
		
		function x = rdivide(A, B)
			if isnumeric(A)	
1888
1889
				x = quantity.Discrete( A ./ B.on(), B(1).domain, ...
                    'size', size(B), 'name', "1 ./ " + B(1).name);
1890
1891
1892
			else
				error('Not yet implemented')
			end
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
		end
		
		function P = matTimes(a, b)
			
			assert(size(a,2) == size(b,1));
			
			p = a(1).gridSize();
			q = p(2);
			p = p(1);
			A = a.on();
			B = b.on();
			
			% dimensions
			n = size(a, 1);
			m = size(b, 2);
			o = size(b, 1);
			
			P = zeros(p, q, n, m);
			
			for k = 1 : n % rows of P
				for l = 1 : m % columns of P
					for r = 1 : o % rows of B or columns of A
						P(:, :, k, l) = P( :, :, k, l ) + A( :, :, k, r) .* B( :, :, r, l );
					end
				end
			end
			
1920
1921
1922
			objOptArgs = a(1).optArgList();
			P = quantity.Discrete(quantity.Discrete.value2cell(P, [p, q], [n, m]), a(1).domain,...
				objOptArgs{:});
1923
1924
1925
1926
1927
		end
		function empty = isempty(obj)
			% ISEMPTY checks if the quantity object is empty
			%   empty = isempty(obj)
			
1928
			% Check if there is any dimension which is zero
1929
1930
			empty = any(size(obj) == 0);
			
Ferdinand Fischer's avatar
tmp    
Ferdinand Fischer committed
1931
1932
1933
1934
1935
1936
1937
			% If the constructor is called without any arguments, an object
			% is still created ( this is required to allow object-arrays),
			% but the object should be recognized as an empty object. Thus,
			% we can test if the domain has been set with a domain object.
			% This happens only in the part of the constructor which is
			% entered if a valid quantity is initialized.
			empty = empty || ~isa(obj(1).domain, 'quantity.Domain');
1938
			
1939
		end % isempty()
1940
		
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
		function result = isdiag(obj)
			% ISDIAG returns true if all non-diagonal elements of the quantity array
			% are zero, i.e. it is a diagonal array. Otherwise, returns false.
			selector = ~logical(eye(size(obj)));
			if (numel(selector) == 1) ||(MAX(abs(obj(selector))) == 0)
				result = true;
			else
				result = false;
			end
		end % isdiag()
		
Jakob Gabriel's avatar
Jakob Gabriel committed
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
		function result = isnan(obj, anyNan)
			%  isnan  True for Not-a-Number.
			%     isnan(X) returns an array that contains true's where
			%     the elements of X are NaN's and false's where they are not.
			%
			%     isnan(X, true) returns the same as any(isnan(X), 'all')
			
			arguments
				obj
				anyNan = false;
			end
			
			if anyNan
				result = any(isnan(obj.on()), 'all');
			else
				result = true(size(obj));
				for it = 1 : numel(obj)
					result(it) = any(isnan(obj(it).on()), 'all');
				end
			end	
		end % isnan()
		
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
		function result = iszero(obj, allZero)
			%  iszero  True for obj == 0.
			%     isnan(X) returns an array that contains true's where
			%     the elements of X == 0 and false's where they are not.
			%
			%     iszero(X, true) returns the same as all(iszero(X), 'all')
			
			arguments
				obj
				allZero = false;
			end
			
			if allZero
				result = all(obj.on() == 0, 'all');
			else
				result = reshape(all(obj.on() == 0, 1:1:obj(1).nargin), size(obj));
			end	
		end % iszero()
		
1993
1994
1995
1996
1997
1998
1999
2000
		function P = ztzTimes(a, b)
			
			assert(size(a,2) == size(b,1))
			
			p = a(1).gridSize();
			q = b(1).gridSize();
			
			assert(p(1) == q(1));
For faster browsing, not all history is shown. View entire blame