compose is working now for functions with multiple independent variables.

Thus, also the TimeDelayOperator is working for such functions.

+quantity/Discrete.m

@@ -231,11 +231,11 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 	end
 	
 	methods (Access = public)
-		function [d, I, d_size] = compositionDomain(obj, g, varargin)
+		function [d, I, d_size] = compositionDomain(obj, domainName)
 			
-			assert(isscalar(g));
-
-			d = g.on();
+			assert(isscalar(obj));
+					
+			d = obj.on();
 			
 			% the evaluation of obj.on( compositionDomain ) is done by:
 			d_size = size(d);
@@ -250,45 +250,93 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			deltaCOD = diff(d);
 			assert(misc.alln(deltaCOD >= 0), 'The domain for the composition f(g(.)) must be monotonically increasing');
 
-
+			d = quantity.Domain('grid', d, 'name', domainName);
 		end
 		
-		function obj_hat = compose(obj, g, varargin)
+		function obj_hat = compose(obj, g, optionalArgs)
 			% COMPOSE compose two functions
 			%	OBJ_hat = compose(obj, G, varargin) composes the function
 			%	defined by OBJ with the function given by G. In particular,
-			%		f_hat(z,t) = f( g(z,t) )
+			%		f_hat(z,t) = f( z, g(z,t) )
 			%	if f(t) = obj, g is G and f_hat is OBJ_hat.
-			
-			assert(nargin(obj) == 1 );
-			
+			arguments
+				obj
+				g quantity.Discrete;
+				optionalArgs.domain quantity.Domain = obj(1).domain;
+			end
+			myCompositionDomain = optionalArgs.domain;
+			originalDomain = obj(1).domain;
+			
+			assert( length( myCompositionDomain ) == 1 );
+			[idx, logOfDomain] = originalDomain.gridIndex(myCompositionDomain);
+			assert( isequal( originalDomain(idx), myCompositionDomain ), ...
+				'Composition of functions: The domains for the composition must be equal. A grid join is not implemented yet.');
+			assert( any( logOfDomain )  )
+			
+			% 2) get the composition domain:
+			%	For the argument y of a function f(y) which should be
+			%	composed by y = g(z,t) a new dommain will be created on the
+			%	basis of evaluation of g(z,t).
 			[composeOnDomain, I, domainSize] = ...
-				obj.compositionDomain(g, varargin{:});
-
+				g.compositionDomain(myCompositionDomain.name);
+			
 			% check if the composition domain is in the range of definition
 			% of obj.
-			if( obj.domain.lower > composeOnDomain(1) || ...
-				obj.domain.upper < composeOnDomain(end) )
-			
+			if ~composeOnDomain.isSubDomainOf( myCompositionDomain )
 				warning('quantity:Discrete:compose', ....
 					'The composition domain is not a subset of obj.domain! The missing values will be extrapolated.');
 			end			
 			
 			% 3) evaluation on the new grid:
-			newValues = obj.on( composeOnDomain );
+			%	the order of the domains is important. At first, the
+			%	domains which will not be replaced are taken. The last
+			%	domain must be the composed domain. For example: a function
+			%	f(x, y, z, t), where y should be composed with g(z, t) will
+			%	be resorted to f_(x, z, t, y) and then evaluated with y =
+			%	g(z,t)
+			evaluationDomain = [originalDomain( ~logOfDomain ), composeOnDomain ];
+			newValues = obj.on( evaluationDomain );
+			
+			% reshape the new values into a 2-d array so that the first
+			% dimension is any domain but the composition domain and the
+			% last dimension is the composition domain
+			sizeOldDomain = prod( [originalDomain( ~logOfDomain ).n] );
+			sizeComposeDomain = composeOnDomain.gridLength;
+			newValues = reshape(newValues, [sizeOldDomain, sizeComposeDomain]);
 			
 			% 4) reorder the computed values, dependent on the sort
-			% position
-			newValues(I) = newValues;
+			% position fo the new domain
+			newValues(:,I) = newValues(:,:);
 			
 			% 5) rearrange the computed values, to have the same dimension
 			% as the required domain
-			newValues = reshape( newValues, domainSize);
+			% *) consider the domain 
+			%		f( z, g(z,t) ) = f(z, g(zeta,t) )|_{zeta = z}
+			tmpDomain = [originalDomain( ~logOfDomain ), g(1).domain ];
+			% newValues will be reshaped into the form
+			%	f(z, t, zeta)
+			newValues = reshape( newValues, [tmpDomain.gridLength, 1] );
+			% *) now the common domains, i.e., zeta = z must be merged:
+			%	for this, find the index of the common domain in list of
+			%	temporary combined domain
+			
+			intersectDomain = intersect( {originalDomain( ~logOfDomain ).name}, ...
+				{g(1).domain.name} );
+			
+			if ~isempty(intersectDomain)
+				
+				idx = 1:length(tmpDomain);
+				idxCommon = idx(strcmp({tmpDomain.name}, intersectDomain));
+
+				% take the diagonal values of the common domain, i.e., z = zeta
+				newValues = misc.diagNd( newValues, idxCommon );
+			end
 			
+			% *) build a new valueDiscrete on the correct grid.		
 			obj_hat = quantity.Discrete( newValues, ...
 				'name', [obj.name '°' g.name], ...
 				'size', size(obj), ...
-				'domain', g.domain());
+				'domain', tmpDomain.join);
 			
 		end
 		
@@ -364,7 +412,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				% the new order will be done in the next step.
 				originalOrderedDomain(gridPermuteIdx) = myDomain;
 				value = obj.obj2value(originalOrderedDomain);
-				value = permute(reshape(value, [cellfun(@(v) numel(v), {originalOrderedDomain.grid}), size(obj)]), ...
+				value = permute(reshape(value, [originalOrderedDomain.gridLength, size(obj)]), ...
 					[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
 			end
 		end

+quantity/Domain.m

@@ -43,8 +43,6 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 				obj.grid = myParser.Results.grid(:);
 				obj.name = myParser.Results.name;
 				myParser.unmatchedWarning();
-			else
-				obj = quantity.Domain.empty();
 			end
 		end
 		
@@ -53,7 +51,11 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 		end
 		
 		function n = get.n(obj)
-			n = length(obj.grid);
+			if isempty(obj.grid)
+				n = [];
+			else
+				n = length(obj.grid);
+			end
 		end
 		
 		function lower = get.lower(obj)
@@ -174,11 +176,15 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 			% [idx, log] = gridIndex(obj, names) searches in the name
 			% properties of obj for the "names" and returns its index as
 			% "idx" and its logical index as "log"
-			
-			if nargin == 1
+			arguments
+				obj
 				names = {obj.name};
 			end
 			
+			if isa(names, 'quantity.Domain')
+				names = {names.name};
+			end
+			
 			names = misc.ensureIsCell(names);
 			
 			idx = zeros(1, length(names));
@@ -201,6 +207,20 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 			end
 		end % gridIndex
 		
+		function i = isempty(obj)
+			i = any(size(obj) == 0);
+			i = i || any( cellfun(@isempty, {obj.grid} ) );
+		end
+		
+		function i = isSubDomainOf(obj, d)
+			% ISSUBDOMAIN 
+			%	i = isSubDomainOf(OBJ, D) checks if the domain OBJ is a
+			%	sub-domain of D.
+			
+			i = obj.lower >= d.lower && ...
+				obj.upper <= d.upper;	
+		end % isSubDomainOf
+		
 		function matGrid = ndgrid(obj)
 			% ndgrid calles ndgrid for the default grid, if no other grid
 			% is specified. Empty grid as input returns empty cell as

+quantity/TimeDelayOperator.m

@@ -39,7 +39,7 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
 			obj.isZero = optionalArgs.isZero;
 		end
 		
-		function d = applyTo(obj, h)
+		function d = applyTo(obj, h, optionalArgs)
 			%APPLYON apply the operator on a function
 			%   D = applyOn(OBJ, H) computes the application of the
 			%   operator OBJ to the function H. 
@@ -47,8 +47,11 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
 			arguments
 				obj
 				h quantity.Discrete;
+ 				optionalArgs.domain quantity.Domain = h(1).domain;
 			end
 			
+			assert( length(optionalArgs.domain) == 1);
+			
 			n = size(obj, 1);
 			m = size(obj, 2); assert(m == size(h, 1));
 			l = size(h, 2);
@@ -59,7 +62,7 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
 				for j = 1 : m
 					for k = 1 : l
 						if ~obj(i, j).isZero
-							d(i, k) = d(i, k) + h(j, k).compose( obj.t + obj(i, j).coefficient );
+							d(i, k) = d(i, k) + h(j, k).compose( obj.t + obj(i, j).coefficient, 'domain', optionalArgs.domain );
 						end
 					end
 				end

+unittests/+quantity/testDiscrete.m

@@ -19,7 +19,6 @@ testCase.verifyWarning(@()f.compose(h), 'quantity:Discrete:compose')
 
 % verify f(tau) = sin( tau * 2pi ); g(z, t) = t + z^2
 % f( g(z, t ) ) = sin( (t + z^2 ) * 2pi )
-
 tau = quantity.Domain('grid', linspace(-1, 3), 'name', 'tau' );
 z = quantity.Domain('grid', linspace(0, 1, 51), 'name', 'z');
 
@@ -27,12 +26,26 @@ f = quantity.Discrete( sin( tau.grid * 2 * pi ), 'domain', tau );
 g = quantity.Discrete( ( z.grid.^2 ) + t.grid', 'domain', [z, t] );
 
 fg = f.compose(g);
-
 F = quantity.Function(@(z,t) sin( ( t + z.^2 ) * 2 * pi ), ...
 	'domain', [z, t], 'name', 'f(g) as function handle');
-
 testCase.verifyEqual(fg.on, F.on, 'AbsTol', 1e-2)
 
+% verify f(z, tau) °_tau g(z, t) -> f(z, g(z,t))
+%	f(z, tau) = z + tau;
+%	g(z, t) = z + t
+%	f(z, g(z,t) ) = 2*z + t
+
+z = quantity.Domain('grid', linspace(0, 1, 11), 'name', 'z');
+t = quantity.Domain('grid', linspace(0, 1, 21), 'name', 't');
+tau = quantity.Domain('grid', linspace(0, 2, 31), 'name', 'tau'); 
+
+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)
+
 end
 
 function testCastDiscrete2Function(testCase)
@@ -1232,10 +1245,10 @@ function testSubs(testCase)
 myTestArray = ones(21, 41, 2, 3);
 myTestArray(:,1,1,1) = linspace(0, 1, 21);
 myTestArray(:,:,2,:) = 2;
-d1 = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'z');
-d2 = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
+z = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'z');
+zeta = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
 
-quan = quantity.Discrete(myTestArray, 'domain', [d1, d2]);
+quan = quantity.Discrete(myTestArray, 'domain', [z, zeta]);
 
 quanZetaZ = quan.subs({'z', 'zeta'}, {'zeta', 'z'});
 quan10 = quan.subs({'z', 'zeta'}, {1, 0});

+unittests/+quantity/testDomain.m

@@ -71,11 +71,19 @@ testCase.verifyEqual( D.numGridElements, prod([length(Z), length(Z)]));
 
 end
 
-function testDomainEmpty(testCase)
+function testEmpty(testCase)
 
 d = quantity.Domain();
+e = quantity.Domain.empty();
+o = quantity.Domain('grid', 1, 'name', '');
 
 testCase.verifyTrue( isempty(d) )
+testCase.verifyTrue( isempty(e) )
+testCase.verifyTrue( isempty([d, e]))
+testCase.verifyTrue( isempty([d, d]))
+testCase.verifyFalse( isempty(o) )
+testCase.verifyFalse( isempty([o e]) )
+testCase.verifyTrue( isempty([o d]) )
 
 end
 

+unittests/+quantity/testTimeDelayOperator.m

@@ -76,7 +76,7 @@ function testApplyTo(testCase)
 	H2 = subs(H2, 'zeta', 0);
 	testCase.verifyEqual(H.on(), H2.on(), 'AbsTol', 5e-3);
 	
-	% test a diagonal operator
+	%% test a diagonal operator
 	h = quantity.Discrete( {sin(tau.grid * 2 * pi); ...
 		sin(tau.grid * 2 * pi)}, 'domain', tau );
 	H = quantity.Discrete( {sin(z.grid * 0 + 2 * pi + t.grid' * 2 * pi); ...
@@ -86,8 +86,23 @@ function testApplyTo(testCase)
 	
 	testCase.verifyEqual( H.on(), H_.on(), 'AbsTol', 5e-3)
 	
+	%% test application of the operator to a function h(z,t)
+	% h(z,tau) = z + sin( tau * 2 * pi)
+	%	D[h] = z + sin( ( t + g(z) ) * 2* pi )
+	%		g(z) = z
+	%	D[h] = z + sin( ( t + z ) * 2 * pi )
 	
+	h = quantity.Discrete( z.grid + sin(tau.grid' * 2 * pi), ...
+		'domain', [z, tau]);
+	
+	Dh = quantity.Discrete(	z.grid + sin( ( z.grid + t.grid') * 2 * pi ), ...
+		'domain', [z, t]);
+	
+	D = testCase.TestData.delay.variable;
+	
+	Dh2 = D.applyTo( h, 'domain', tau );
 
+	testCase.verifyEqual( Dh.on, Dh2.on, 'AbsTol', 2e-3)
 	
 end