Fixed quantity.Discrete/on to keep the right order of arguments also for 3-dim domains.

+quantity/Discrete.m

@@ -231,32 +231,40 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 	end
 	
 	methods (Access = public)
-		function obj_hat = compose(obj, g, varargin)
-			% 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) )
-			%	if f(t) = obj, g is G and f_hat is OBJ_hat.
+		function [d, I, d_size] = compositionDomain(obj, g, varargin)
 			
 			assert(isscalar(g));
-			assert(nargin(obj) == 1 );
-			
-			newDomain = g.domain();
-			composeOnDomain = g.on();
+
+			d = g.on();
 			
 			% the evaluation of obj.on( compositionDomain ) is done by:
-			domainSize = size(composeOnDomain);
+			d_size = size(d);
 			
 			% 1) vectorization of the n-d-grid: compositionDomain	
-			composeOnDomain = composeOnDomain(:);
+			d = d(:);
 
 			% 2) then it is sorted in ascending order
-			[composeOnDomain, I] = sort(composeOnDomain);			
+			[d, I] = sort(d);			
 			
 			% verify the domain to be monotonical increasing
-			deltaCOD = diff(composeOnDomain);
+			deltaCOD = diff(d);
 			assert(misc.alln(deltaCOD >= 0), 'The domain for the composition f(g(.)) must be monotonically increasing');
 
+
+		end
+		
+		function obj_hat = compose(obj, g, varargin)
+			% 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) )
+			%	if f(t) = obj, g is G and f_hat is OBJ_hat.
+			
+			assert(nargin(obj) == 1 );
+			
+			[composeOnDomain, I, domainSize] = ...
+				obj.compositionDomain(g, varargin{:});
+
 			% check if the composition domain is in the range of definition
 			% of obj.
 			if( obj.domain.lower > composeOnDomain(1) || ...
@@ -264,7 +272,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			
 				warning('quantity:Discrete:compose', ....
 					'The composition domain is not a subset of obj.domain! The missing values will be extrapolated.');
-			end
+			end			
 			
 			% 3) evaluation on the new grid:
 			newValues = obj.on( composeOnDomain );
@@ -280,7 +288,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			obj_hat = quantity.Discrete( newValues, ...
 				'name', [obj.name '°' g.name], ...
 				'size', size(obj), ...
-				'domain', newDomain );
+				'domain', g.domain());
 			
 		end
 		
@@ -346,15 +354,17 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 					gridPermuteIdx = 1:length(myDomain);
 				else
 					assert(numel(myDomain) == numel(obj(1).domain), ['Wrong grid for the evaluation of the object']);
-					[myDomain, gridPermuteIdx] = obj(1).domain.permute(myDomain);
+					% compute the permutation index, in order to bring the
+					% new domain in the same order as the original one.
+					gridPermuteIdx = obj(1).domain.getPermutationIdx(myDomain);
 				end			
 				% get the valueDiscrete data for this object. Apply the
 				% permuted myDomain. Then the obj2value will be evaluated
 				% in the order of the original domain. The permuatation to
 				% the new order will be done in the next step.
-				value = obj.obj2value(myDomain(gridPermuteIdx));
-				
-				value = permute(reshape(value, [cellfun(@(v) numel(v), {myDomain(gridPermuteIdx).grid}), size(obj)]), ...
+				originalOrderedDomain(gridPermuteIdx) = myDomain;
+				value = obj.obj2value(originalOrderedDomain);
+				value = permute(reshape(value, [cellfun(@(v) numel(v), {originalOrderedDomain.grid}), size(obj)]), ...
 					[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
 			end
 		end
@@ -1978,36 +1988,6 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				'name', [A(1).name, '+', B(1).name]);
 		end
 		
-		function [valDiscrete] = expandValueDiscrete(obj, newDomain)
-			% EXPANDVALUEDISCRETE
-			%	[valDiscrete, gridNameSorted] = ...
-			%       expandValueDiscrete(obj, gridIndex, valDiscrete)
-			%	Expanses the value of obj, so that
-			
-			
-			gridNameJoined  = {newDomain.name};
-			gridJoinedLength = newDomain.gridLength;
-			
-			% get the index of obj.grid in the joined grid
-			[~, logicalIdx] = newDomain.gridIndex({obj(1).domain.name});
-			
-			valDiscrete = obj.on( newDomain(logicalIdx) );
-			oldDim = ndims(valDiscrete);
-			valDiscrete = permute(valDiscrete, [(1:sum(~logicalIdx)) + oldDim, 1:oldDim] );
-			valDiscrete = repmat(valDiscrete, [gridJoinedLength(~logicalIdx), ones(1, ndims(valDiscrete))]);
-			%
-			valDiscrete = reshape(valDiscrete, ...
-				[gridJoinedLength(~logicalIdx), gridJoinedLength(logicalIdx), size(obj)]);
-			
-			% permute valDiscrete such that grids are in the order specified
-			% by gridNameJoined.
-			gridIndex = 1:numel(logicalIdx);
-			gridOrder = [gridIndex(~logicalIdx), gridIndex(logicalIdx)];
-			valDiscrete = permute(valDiscrete, [gridOrder, numel(logicalIdx)+(1:ndims(obj))]);
-			
-		end
-		
-		
 		function C = minus(A, B)
 			% minus uses plus()
 			C = A + (-B);
@@ -2208,25 +2188,25 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 	%%
 	methods (Static)
 		
-		function P = zeros(valueSize, grid, varargin)
+		function P = zeros(valueSize, domain, varargin)
 			%ZEROS initializes an zero quantity.Discrete object
-			%   P = zeros(VALUESIZE, GRID) returns a quantity.Discrete
-			%   object that has only zero entries on a grid which is
-			%   defined by the cell GRID and the value size defined by the
-			%   size-vector VALUESIZE
-			%       P = quantity.Discrete([n,m], {linspace(0,1)',
-			%                                       linspace(0,10)});
-			%       creates an (n times m) zero quantity.Discrete on the
-			%       grid (0,1) x (0,10)
-			%
 			%	P = zeros(VALUESIZE, DOMAIN) creates a matrix of size
 			%	VALUESIZE on the DOMAIN with zero entries.
-			if ~iscell(grid)
-				grid = {grid};
+			
+			myParser = misc.Parser();
+			myParser.addParameter('gridName', '');
+			myParser.parse(varargin{:});
+			
+			if ~isa(domain, 'quantity.Domain')
+				% if the input parameter DOMAIN is not a quantity.Domain
+				% object. It is assumed that it is a grid.
+				grids = misc.ensureIsCell(domain);
+				gridNames = misc.ensureIsCell( myParser.Results.gridName );
+				domain = quantity.Domain.gridCells2domain(grids, gridNames);
 			end
-			gridSize = cellfun('length', grid);
-			O = zeros([gridSize(:); valueSize(:)]');
-			P = quantity.Discrete(O, 'size', valueSize, 'grid', grid, varargin{:});
+			
+			O = zeros([domain.gridLength, valueSize(:)']);
+			P = quantity.Discrete(O, 'size', valueSize, 'domain', domain, varargin{:});
 		end
 		
 		function q = value2cell(value, gridSize, valueSize)
@@ -2276,6 +2256,37 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 	end %% (Static)
 	methods(Access = protected)
 		
+		function [valDiscrete] = expandValueDiscrete(obj, newDomain)
+			% EXPANDVALUEDISCRETE expand the discrete value on the
+			% newDomain
+			%	[valDiscrete] = ...
+			%       expandValueDiscrete(obj, newDomain) expands the
+			%       discrete values on a new domain. So that a function
+			%			f(z,t) = f(z) + f(t)
+			%	can be computed.
+		
+			gridJoinedLength = newDomain.gridLength;
+			
+			% get the index of obj.grid in the joined grid
+			[idx, logicalIdx] = newDomain.gridIndex({obj(1).domain.name});
+			% evaluate the 
+			valDiscrete = obj.on( newDomain(logicalIdx) );
+			oldDim = ndims(valDiscrete);
+			valDiscrete = permute(valDiscrete, [(1:sum(~logicalIdx)) + oldDim, 1:oldDim] );
+			valDiscrete = repmat(valDiscrete, [gridJoinedLength(~logicalIdx), ones(1, ndims(valDiscrete))]);
+			%
+			valDiscrete = reshape(valDiscrete, ...
+				[gridJoinedLength(~logicalIdx), gridJoinedLength(logicalIdx), size(obj)]);
+			
+			% permute valDiscrete such that grids are in the order specified
+			% by gridNameJoined.
+			gridIndex = 1:numel(logicalIdx);
+			gridOrder = [gridIndex(~logicalIdx), gridIndex(logicalIdx)];
+			gridIndex(gridOrder) = 1:numel(logicalIdx);
+			
+			valDiscrete = permute(valDiscrete, [gridIndex, numel(logicalIdx)+(1:ndims(obj))]);
+		end
+		
 		function result = diff_inner(obj, k, diffGridName)
 			gridSelector = strcmp(obj(1).gridName, diffGridName);
 			gridSelectionIndex = find(gridSelector);

+quantity/Domain.m

@@ -213,7 +213,7 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 			end
 		end % ndgrid()
 		
-		function [newDomain, idx] = permute(obj, order)
+		function [idx, newDomain] = getPermutationIdx(obj, order)
 			
 			if isa(order, 'quantity.Domain')
 				names = {order.name};
@@ -309,6 +309,20 @@ classdef Domain < handle & matlab.mixin.CustomDisplay
 	end % Access = protected
 	
 	methods (Static)
+		
+		function d = gridCells2domain(grids, gridNames)
+			grids = misc.ensureIsCell(grids);
+			gridNames = misc.ensureIsCell(gridNames);
+			
+			assert(length( grids ) == length(gridNames))
+
+			d = quantity.Domain.empty();
+			
+			for k = 1:length(grids)
+				d = [d quantity.Domain('grid', grids{k}, 'name', gridNames{k})];
+			end
+		end
+		
 		function g = defaultGrid(gridSize, name)
 			
 			if nargin == 1

+unittests/+quantity/testDiscrete.m

@@ -587,7 +587,19 @@ testCase.verifyEqual(permute(createTestData(linspace(0, 1, 21), linspace(0, 1, 2
 		value(:,:,2,3) = 1+zeros(numel(gridVecA), numel(gridVecB));
 	end
 
+%% test on on 3-dim domain
+z = quantity.Domain('grid', 0:3, 'name', 'z');
+t = quantity.Domain('grid', 0:4, 'name', 't');
+x = quantity.Domain('grid', 0:5, 'name', 'x');
 
+Z = quantity.Discrete( z.grid, 'domain', z);
+T = quantity.Discrete( t.grid, 'domain', t);
+X = quantity.Discrete( x.grid, 'domain', x);
+
+ZTX = Z+T+X;
+XTZ = X+T+Z;
+
+testCase.verifyEqual( ZTX.on([t z x]), XTZ.on([t z x]), 'AbsTol', 1e-12);
 
 end
 
@@ -1048,7 +1060,6 @@ eMatReference = zGrid + zetaGrid;
 testCase.verifyEqual(numel(eMat), numel(eMatReference));
 testCase.verifyEqual(eMat(:), eMatReference(:));
 
-
 %% addition with constant values
 testCase.verifyEqual(permute([a b], [1 3 2]), AB.on());
 
@@ -1064,6 +1075,29 @@ cAB2 = [1 2; 3 4] + AB2;
 testCase.verifyEqual(AB2c.on(), tst)
 testCase.verifyEqual(cAB2.on(), tst)
 
+%% test plus on different domains
+z = quantity.Domain('grid', 0:3, 'name', 'z');
+t = quantity.Domain('grid', 0:4, 'name', 't');
+x = quantity.Domain('grid', 0:5, 'name', 'x');
+
+T = quantity.Discrete( t.grid, 'domain', t);
+Z = quantity.Discrete( z.grid, 'domain', z);
+X = quantity.Discrete( x.grid, 'domain', x);
+
+TZX = T+Z+X;
+
+TZ1 = T + Z + 1;
+TZX1 = subs(TZX, 'x', 1);
+
+testCase.verifyEqual( TZX1.on(), TZ1.on(), 'AbsTol', 1e-12 );
+testCase.verifyEqual( on( T + 0.5 + X  ), on( subs( TZX, 'z', 0.5) ), 'AbsTol', 1e-12 );
+
+XTZ = X+T+Z;
+
+TZX + XTZ;
+
+%TODO
+
 end
 
 function testInit(testCase)