Discrete.m

classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
		< handle & matlab.mixin.Copyable & matlab.mixin.CustomDisplay

	properties (SetAccess = protected)
		% Discrete evaluation of the continuous quantity
		valueDiscrete double;
	end
	
	properties (Hidden, Access = protected, Dependent)
		doNotCopy;
	end
	
	properties
		% ID of the figure handle in which the handle is plotted
		figureID double = 1;
		
		% Name of this object
		name (1,1) string;
		
		% domain
		domain;
	end
	
	properties ( Dependent )
		
		% Name of the domains that generate the grid.
		gridName {mustBe.unique};
		
		% Grid for the evaluation of the continuous quantity. For the
		% example with the function f(x,t), the grid would be
		%   {[<spatial domain>], [<temporal domain>]}
		% whereas <spatial domain> is the discret description of the
		% spatial domain and <temporal domain> the discrete description of
		% the temporal domain.
		grid; % in set.grid it is ensured that, grid is a (1,:)-cell-array
	end
	
	methods
		%--------------------
		% --- Constructor ---
		%--------------------
		function obj = Discrete(valueOriginal, varargin)
			% DISCRETE a quantity, represented by discrete values.
			%	obj = Discrete(valueOriginal, varargin) initializes a
			%	quantity. The parameters to be set are:
			% 'valueOrigin' must be
			% 1) a cell-array of double arrays with
			%	size(valueOriginal) == size(obj) and
			%	size(valueOriginal{it}) == gridSize
			%	Example: valueOrigin = { f(Z, T), g(Z, T) } is a cell array
			%	wich contains the functions f(z,t) and g(z,t) evaluated on
			%	the discrete domain (Z x T). Then, the name-value-pair
			%	parameter 'domain' must be set with quantity.Domain
			%	objects, according to the domains Z and T.
			% OR
			% 2) adouble-array with
			%	size(valueOriginal) == [gridSize, size(quantity)]
			% Furthermore, 'gridName' must be part of the name-value-pairs
			% in varargin. Additional parameters can be specified using
			% name-value-pair-syntax in varargin.
			
			% to allow the initialization as object array, the constructor
			% must be allowed to be called without arguments, i.e. nargin == 0.
			% Then no parameters are set.
			if nargin == 1
				% if nargin == 1 it can be a conversion of child-classes or an empty
				% object
				if isa(valueOriginal, 'quantity.Discrete')
					% allows the conversion of a quantity object without
					% extra check if the object is already from class
					% quantity.Discrete
					obj = valueOriginal;
				else
					% empty object. this is needed for instance, to create
					% quantity.Discrete([]), which is useful for creating default
					% values.
					obj = quantity.Discrete.empty(size(valueOriginal));
				end
			elseif nargin > 1
				
				%% input parser
				myParser = misc.Parser();
				myParser.addParameter('name', "", @mustBe.gridName);
				myParser.addParameter('figureID', 1, @isnumeric);
				myParser.parse(varargin{:});
				
				%% domain parser
				myDomain = quantity.Domain.parser(varargin{:});
												
				%% get the sizes of obj and grid
				gridLength = myDomain.gridLength;

				% convert double valued valueOriginal to cell-valued value
				% original
				if ~iscell(valueOriginal)
					valueOriginal = quantity.Discrete.value2cell(valueOriginal, gridLength);
				end
				
				% Check if the grid fits to the values. In addition, catch
				% the case if all values are empty. This is required for
				% the initialization of quantity.Function and
				% quantity.Symbolic objects
				assert( misc.alln( cellfun(@isempty, valueOriginal ) ) || ...
					numGridElements(myDomain) == numel(valueOriginal{1}), ...
					'grids do not fit to valueOriginal');				
				
				% allow initialization of empty objects
				valueOriginalSize = size(valueOriginal);
				if any(valueOriginalSize == 0)
					% If the size is specified in the arguements, it should
					% be chosen instead of the default size from the
					% valueOriginal.
					myParser = misc.Parser();
					myParser.addParameter('size', valueOriginalSize((1+ndims(myDomain)):end));
					myParser.parse(varargin{:});
					obj = quantity.Discrete.empty(myParser.Results.size);
					return;
				end
				
				% set valueDiscrete
				for k = 1:numel(valueOriginal)
					if numel(myDomain) == 1
						% for quantities on a single domain, ensure that
						% the discrete values are stored as column-vector
						% by using the (:) operator.
						obj(k).valueDiscrete = valueOriginal{k}(:); %#ok<AGROW>
					else
						obj(k).valueDiscrete = valueOriginal{k}; %#ok<AGROW>
					end
				end
				
				%% set further properties
				[obj.domain] = deal(myDomain);
				obj.setName(myParser.Results.name);
				[obj.figureID] = deal(myParser.Results.figureID);
				
				%% reshape object from vector to matrix
				obj = reshape(obj, size(valueOriginal));
			end
		end% Discrete() constructor
		
		%---------------------------
		% --- getter and setters ---
		%---------------------------
		function gridName = get.gridName(obj)
			if isempty(obj.domain)
				gridName = [];
			else
				gridName = [obj.domain.name];
			end
		end
		
		function grid = get.grid(obj)
			if isempty(obj.domain)
				grid = {};
			else
				grid = {obj.domain.grid};
			end
		end
		
		function itIs = isConstant(obj)
			% the quantity is interpreted as constant if it has no grid or
			% it has a grid that is only defined at one point.
			itIs = isempty(obj(1).domain);
		end % isConstant()
		
		function doNotCopy = get.doNotCopy(obj)
			doNotCopy = obj.doNotCopyPropertiesName();
		end
		function valueDiscrete = get.valueDiscrete(obj)
			% check if the value discrete for this object
			% has already been computed.
			empty = isempty(obj.valueDiscrete);
			if any(empty(:))
				obj.valueDiscrete = obj.obj2value(obj.domain, true);
			end
			valueDiscrete = obj.valueDiscrete;
		end
		
		%-------------------
		% --- converters ---
		%-------------------
		function exportData = exportData(obj, varargin)
			
			% make the object names:
			if obj.nargin == 1
				headers = cell(1, numel(obj) + 1);
				headers{1} = obj(1).gridName{1};
				for i= 1:numel(obj) %TODO use easier to read headers
					headers{i+1} = obj(i).name + "" + num2str(i);
				end
				exportData = export.dd(...
					'M', [obj.grid{:}, obj.valueDiscrete], ...
					'header', headers, varargin{:});
			elseif obj.nargin == 2
				error('Not yet implemented')
			else
				error('Not yet implemented')
			end
		end
		function d = double(obj)
			d = obj.on();
		end
		function o = quantity.Function(obj)
			props = nameValuePair( obj(1) );
			
			for k = 1:numel(obj)
				F = griddedInterpolant(obj(k).grid{:}', obj(k).on());
				o(k) = quantity.Function(@(varargin) F(varargin{:}), ...
					props{:});
			end
			
			o = reshape(o, size(obj));
		end
		function o = signals.PolynomialOperator(obj)
			A = cell(size(obj, 3), 1);
			for k = 1:size(obj, 3)
				A{k} = obj(:,:,k);
			end
			o = signals.PolynomialOperator(A);
		end
		
		function o = quantity.Symbolic(obj)
			if isempty(obj)
				o = quantity.Symbolic.empty(size(obj));
			else
				error('Not yet implemented')
			end
		end
		
		function obj = setName(obj, newName)
			% Function to set all names of all elements of the quantity obj to newName.
% 			if ischar(newName)
% 				warning("Depricated: use string and not char for name-property!")
% 				newName = string(newName);
% 			end
			[obj.name] = deal(newName);
		end % setName()
	end
	
	methods (Access = public)
		function d = compositionDomain(obj, domainName)
			
			assert(isscalar(obj));
					
			d = obj.on();
			
			% the evaluation of obj.on( compositionDomain ) is done by:
			d_size = size(d);
			
			% vectorization of the n-d-grid: compositionDomain	
			d = quantity.Domain(domainName, d(:));
		end
		
		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( z, g(z,t) )
			%	if f(t) = obj, g is G and f_hat is OBJ_hat.
			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 )  )
			
			% 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 = ...
				g.compositionDomain(myCompositionDomain.name);
			
			% check if the composition domain is in the range of definition
			% of obj.
			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			
			
			% evaluation on the new grid:
			%	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)
			
			% #TODO: optimize the memory consumption of this function.
			%	1) only consider the unqiue grid points in evaluationDomain
			%	2) do the conversion of the evaluationDomain directly to
			%	the target domain.			
			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]);

			%rearrange the computed values, to have the same dimension
			% as the required domain
			% 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 ), ...
				g(1).domain );
			
			if ~isempty(intersectDomain)
				
				idx = tmpDomain.gridIndex( intersectDomain );
				
				% take the diagonal values of the common domain, i.e., z = zeta		
				% use the diag_nd function because it seems to be faster
				% then the diagNd function, although the values must be
				% sorted.
				newValues = misc.diagNd(newValues, idx);
			end
			
			% *) build a new valueDiscrete on the correct grid.		
			obj_hat = quantity.Discrete( newValues, ...
				'name', obj.name + "°" + g.name, ...
				'size', size(obj), ...
				'domain', tmpDomain.join);
			
		end
		
		function value = on(obj, myDomain, gridNames)
			% ON evaluation of the quantity on a certain domain.
			%	value = on(obj) or value = obj.on(); evaluates the quantity
			%	on its standard grid. 
			%	value = obj.on( myDomain ) evalutes the quantity on the
			%	grid specified by myDomain. The order of the domains in
			%	domain, will be the same as from myDomain. 
			%	value = obj.on( grid ) evaluates the quantity specified by
			%	grid. Grid must be a cell-array with the grids as elements.
			%	value = obj.on( grid, gridName ) evaluates the quantity
			%	specified by grid. Grid must be a cell-aary with the grids
			%	as elements. By the gridName parameter the order of the
			%	grid can be specified.

			if isempty(obj)
				value = zeros(size(obj));
			else
				if nargin == 1
					% case 0: no domain was specified, hence the value is requested
					% on the default grid defined by obj(1).domain.
					value = obj.obj2value(obj(1).domain);
					
				elseif nargin == 2 && (iscell(myDomain) || isnumeric(myDomain))
					% case 1: a domain is specified by myDomain as agrid
					myDomain = misc.ensureIsCell(myDomain);
					newGrid = myDomain;

					if obj(1).isConstant()
						gridNames = repmat({''}, length(newGrid));
					else
						gridNames = {obj(1).domain.name};
					end

					% initialize the new domain
					clear('myDomain');
					myDomain(1:length(newGrid)) = quantity.Domain();					
					for k = 1:length(newGrid)
						myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
					end
					value = reshape(obj.obj2value(myDomain), ...
						           [myDomain.gridLength, size(obj)]);
				else
					% Since in the remaining cases the order of the domains is not 
					% neccessarily equal to the order in obj(1).domain, this is 
					% more involved:
					if nargin == 2
						% case 2: a domain is specified by a myDomain = domain-array
						% nothing has to be done to obtain the domain.

					elseif nargin == 3
						% case 3: a domain is specified by a grid and a grid
						% name. Then, the first input parameter is the grid,
						% i.e., myGrid = myDomain and the second is the grid
						% name.
						% Since the order of the domains is not neccessarily equal to the
						% order in obj(1).domain, this is more involved:
						myDomain = misc.ensureIsCell(myDomain);
						gridNames = misc.ensureString(gridNames);

						assert(all(cellfun(@(v)isvector(v), myDomain)), ...
							'The cell entries for a new grid have to be vectors')

						newGrid = myDomain;
						myDomain = quantity.Domain.empty();
						for k = 1:length(newGrid)
							myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
						end
					else
						error('wrong number of input arguments')
					end

					% verify the domain
					if obj(1).isConstant
						gridPermuteIdx = 1:length(myDomain);
					else
						assert(numel(myDomain) == numel(obj(1).domain), ...
							'Wrong grid for the evaluation of the object');
						% 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 permutation to
					% the new order will be done in the next step.
					originalOrderedDomain(gridPermuteIdx) = myDomain;
					value = obj.obj2value(originalOrderedDomain);
					value = permute(reshape(value, [originalOrderedDomain.gridLength, size(obj)]), ...
						[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
				end
			end % if isempty(obj)
		end % on()
		
		function interpolant = interpolant(obj)
			% get the interpolant of the obj;
			if isempty(obj)
				value = zeros(size(obj));
				indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
				interpolant = numeric.interpolant(...
					[indexGrid{:}], value);
			else
				myGrid = obj(1).grid;
				value = obj.obj2value();
				indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
				interpolant = numeric.interpolant(...
					[myGrid, indexGrid{:}], value);
			end
		end
		
		
		function assertSameGrid(a, varargin)
			% check if all elements of a have same grid and gridName. If
			% further quantites are inputs via varargin, it is verified if
			% that quantity has same grid and gridName as quantity a as
			% well.
			if isempty(a)
				if nargin > 1
					varargin{1}.assertSameGrid(varargin{2:end});
				end
				return;
			else
				referenceGridName = a(1).gridName;
				referenceGrid= a(1).grid;
			end
			for it = 1 : numel(a)
				assert(isequal(referenceGridName, a(it).gridName), ...
					'All elements of a quantity must have same gridNames');
				assert(isequal(referenceGrid, a(it).grid), ...
					'All elements of a quantity must have same grid');
			end
			if nargin > 1
				b = varargin{1};
				for it = 1 : numel(b)
					assert(isequal(referenceGridName, b(it).gridName), ...
						'All elements of a quantity must have same gridNames');
					assert(isequal(referenceGrid, b(it).grid), ...
						'All elements of a quantity must have same grid');
				end
			end
			if nargin > 2
				% if more then 1 quantity is in varargin, they are checked
				% iteratively by calling assertSameGrid() again.
				assertSameGrid(varargin{:});
			end
		end
		
		function [referenceGrid, referenceGridName] = getFinestGrid(a, varargin)
			% find the finest grid of all input quantities by comparing
			% gridSize for each iteratively.
			
			if isempty(a) || isempty(a(1).grid)
				if nargin > 1
					[referenceGrid, referenceGridName] = varargin{1}.getFinestGrid(varargin{2:end});
				else
					referenceGrid = {};
					referenceGridName = '';
				end
				return;
			else
				referenceGridName = a(1).gridName;
				referenceGrid = a(1).grid;
				referenceGridSize = [a(1).domain.n];
			end
			
			for it = 1 : numel(varargin)
				if isempty(varargin{it}) || isempty(varargin{it}(1).domain)
					continue;
				end
				assert(numel(referenceGridName) == numel(varargin{it}(1).gridName), ...
					['For getFinestGrid, the gridName of all objects must be equal', ...
					'. Maybe gridJoin() does what you want?']);
				for jt = 1 : numel(referenceGridName)
					comparisonGridSize = varargin{it}(1).domain.find(referenceGridName{jt}).n;
					comparisonGrid = varargin{it}.gridOf(referenceGridName{jt});
					assert(referenceGrid{jt}(1) == comparisonGrid(1), 'Grids must have same domain for combining them')
					assert(referenceGrid{jt}(end) == comparisonGrid(end), 'Grids must have same domain for combining them')
					if comparisonGridSize > referenceGridSize(jt)
						referenceGrid{jt} = comparisonGrid;
						referenceGridSize(jt) = comparisonGridSize;
					end
				end
			end
		end
		
		function obj = sort(obj, varargin)
			%SORT sorts the grid of the object in a desired order
			% obj = sortGrid(obj) sorts the grid in alphabetical order.
			% obj = sort(obj, 'descend') sorts the grid in descending
			% alphabetical order.
						
			% only sort the grids if there is something to sort
			if obj(1).nargin > 1
				
				[sortedDomain, I] = obj(1).domain.sort(varargin{:});
				[obj.domain] = deal(sortedDomain);
				
				for k = 1:numel(obj)
					obj(k).valueDiscrete = permute(obj(k).valueDiscrete, I);
				end
			end
		end% sort()
		
		function c = horzcat(a, varargin)
			%HORZCAT Horizontal concatenation.
			%   [A B] is the horizontal concatenation of objects A and B
			%   from the class quantity.Discrete. A and B must have the
			%   same number of rows and the same grid. [A,B] is the same
			%   thing. Any number of matrices can be concatenated within
			%   one pair of brackets. Horizontal and vertical concatenation
			%   can be combined together as in [1 2;3 4].
			%
			%   [A B; C] is allowed if the number of rows of A equals the
			%   number of rows of B and the number of columns of A plus the
			%   number of columns of B equals the number of columns of C.
			%   The matrices in a concatenation expression can themselves
			%   by formed via a concatenation as in [A B;[C D]].  These
			%   rules generalize in a hopefully obvious way to allow fairly
			%   complicated constructions.
			%
			%   N-D arrays are concatenated along the second dimension. The
			%   first and remaining dimensions must match.
			%
			%   C = HORZCAT(A,B) is called for the syntax '[A  B]' when A
			%   or B is an object.
			%
			%   Y = HORZCAT(X1,X2,X3,...) is called for the syntax '[X1 X2
			%   X3 ...]' when any of X1, X2, X3, etc. is an object.
			%
			%	See also HORZCAT, CAT.
			c = cat(2, a, varargin{:});
		end
		function c = vertcat(a, varargin)
			%VERTCAT Vertical concatenation.
			%   [A;B] is the vertical concatenation of objects A and B from
			%   the class quantity.Discrete. A and B must have the same
			%   number of columns and the same grid. Any number of matrices
			%   can be concatenated within one pair of brackets. Horizontal
			%   and vertical concatenation can be combined together as in
			%   [1 2;3 4].
			%
			%   [A B; C] is allowed if the number of rows of A equals the
			%   number of rows of B and the number of columns of A plus the
			%   number of columns of B equals the number of columns of C.
			%   The matrices in a concatenation expression can themselves
			%   by formed via a concatenation as in [A B;[C D]].  These
			%   rules generalize in a hopefully obvious way to allow fairly
			%   complicated constructions.
			%
			%   N-D arrays are concatenated along the first dimension. The
			%   remaining dimensions must match.
			%
			%   C = VERTCAT(A,B) is called for the syntax '[A; B]' when A
			%   or B is an object.
			%
			%   Y = VERTCAT(X1,X2,X3,...) is called for the syntax '[X1;
			%   X2; X3; ...]' when any of X1, X2, X3, etc. is an object.
			%
			%   See also HORZCAT, CAT.
			c = cat(1, a, varargin{:});
		end
		function c = cat(dim, a, varargin)
			%CAT Concatenate arrays.
			%   CAT(DIM,A,B) concatenates the arrays of objects A and B
			%   from the class quantity.Discrete along the dimension DIM.
			%   CAT(2,A,B) is the same as [A,B]. CAT(1,A,B) is the same as
			%   [A;B].
			%
			%   B = CAT(DIM,A1,A2,A3,A4,...) concatenates the input arrays
			%   A1, A2, etc. along the dimension DIM.
			%
			%   When used with comma separated list syntax, CAT(DIM,C{:})
			%   or CAT(DIM,C.FIELD) is a convenient way to concatenate a
			%   cell or structure array containing numeric matrices into a
			%   single matrix.
			%
			%   Examples:
			%     a = magic(3); b = pascal(3);
			%     c = cat(4,a,b)
			%   produces a 3-by-3-by-1-by-2 result and
			%     s = {a b};
			%     for i=1:length(s),
			%       siz{i} = size(s{i});
			%     end
			%     sizes = cat(1,siz{:})
			%   produces a 2-by-2 array of size vectors.
			
			if nargin == 1
				objCell = {a};
			else
				objCell = [{a}, varargin(:)'];
				
				% this function has the very special thing that it a does
				% not have to be an quantity.Discrete object. So it has to
				% be checked which of the input arguments is an
				% quantity.Discrete object. This is considered to give
				% the basic values for the initialization of new
				% quantity.Discrete values
				isAquantityDiscrete = cellfun(@(o) isa(o, 'quantity.Discrete'), objCell);
				isEmpty = cellfun(@(o) isempty(o), objCell);
				objIdx = find(isAquantityDiscrete & (~isEmpty), 1);
				
				if all(isEmpty)
					% if there are only empty entries, nothing can be
					% concatenated, so a new empty object is initialized.
					s = cellfun(@(o) size(o), objCell, 'UniformOutput', false);
					if dim == 1
						S = sum(cat(3, s{:}), 3);
					elseif dim == 2
						S = s{1};
					else
						error('Not implemented')
					end
					c = quantity.Discrete.empty(S);
					return
				else
					obj = objCell{objIdx};
				end
				
				for k = 1:numel(objCell(~isEmpty))
					
					if isa(objCell{k}, 'quantity.Discrete')
						o = objCell{k};
					else
						value = objCell{k};
						for l = 1:numel(value)
							M(:,l) = repmat(value(l), prod(obj(1).domain.gridLength), 1);
						end
						if isempty(value)
							M = zeros([prod(obj(1).domain.gridLength), size(value(l))]);
						end
						M = reshape(M, [obj(1).domain.gridLength, size(value)]);
						o = quantity.Discrete( M, ...
							'size', size(value), ...
							'domain', obj(1).domain);
					end
					
					objCell{k} = o;
				end
				
			end
			
			% sort the grid names of each quantity
			for it = 1: (numel(varargin) + 1)
				objCell{it} = objCell{it}.sort;
			end
			
			[fineGrid, fineGridName] = getFinestGrid(objCell{~isEmpty});
			for it = 1 : (numel(varargin) + 1)  % +1 because the first entry is a
				% change the grid to the finest
				objCell{it} = objCell{it}.changeGrid(fineGrid, fineGridName);
			end
			assertSameGrid(objCell{:});
			argin = [{dim}, objCell(:)'];
			c = builtin('cat', argin{:});
		end
		
		function Y = blkdiag(A, varargin)
			% blkdiag  Block diagonal concatenation of matrix input arguments.
			%									|A 0 .. 0|
			% Y = blkdiag(A,B,...)  produces	|0 B .. 0|
			%									|0 0 ..  |
			% Yet, A, B, ... must have the same gridName and grid.
			if nargin == 1
				Y = copy(A);
			else
				B = varargin{1};
				if isempty(B)
					Y = A;
				else
					assert(isequal(A(1).gridName, B(1).gridName), 'only implemented for same grid and gridName');
					assert(isequal(A(1).grid, B(1).grid), 'only implemented for same grid and gridName');
					Y = [A, zeros(size(A, 1), size(B, 2)); ...
						zeros(size(B, 1), size(A, 2)), B];
				end
				if nargin > 2
					Y = blkdiag(Y, varargin{2:end});
				end
			end
		end % blkdiag()
		
		function solution = solveAlgebraic(obj, rhs, gridName, objLimit)
			%% this method solves
			%	obj(gridName) == rhs
			% for the variable specified by gridName.
			% rhs must be of apropriate size and gridName must
			% be an gridName of obj. If the result is constant (i.e., if
			% obj only depends on variable, then a double array is
			% returned. Else the solution is of the type as obj.
			% Yet, this is only implemented for obj with one variable
			% (grid) (see quantity.invert-method).
			% The input objLimit specifies minimum and maximum of the
			% values of obj, between which the solution should be searched.
			assert(numel(obj(1).gridName) == 1);
			assert(isequal(size(obj), [1, 1]));
			
			if ~isequal(size(rhs), size(obj))
				error('rhs has not the same size as quantity');
			end
			if ~iscell(gridName)
				gridName = {gridName};
			end
			if numel(gridName) ~= 1
				error('this function can only solve for one variable');
			end
			if isempty(strcmp(obj(1).gridName, gridName{1}))
				error('quantity does not depend on variable');
			end
			
			if nargin == 4
				assert(numel(objLimit)==2, 'a lower and upper limit must be specified (or neither)');
				objValueTemp = obj.on();
				gridSelector = (objValueTemp >= objLimit(1)) & (objValueTemp <= objLimit(2));
				gridSelector([max(1, find(gridSelector, 1, 'first')-1), ...
					min(find(gridSelector, 1, 'last')+1, numel(gridSelector))]) = 1;
				limitedGrid = obj(1).grid{1}(gridSelector);
				objCopy = obj.copy();
				objCopy = objCopy.changeGrid({limitedGrid}, gridName);
				objInverseTemp = objCopy.invert(gridName);
			else
				objInverseTemp = obj.invert(gridName);
			end
			
			solution = objInverseTemp.on(rhs);
			
			% 			solution = zeros(numel(obj), 1);
			% 			for it = 1 : numel(obj)
			% 				objInverseTemp = obj(it).invert(gridName);
			% 				solution(it) = objInverseTemp.on(rhs(it));
			% 			end
			% 			solution = reshape(solution, size(obj));
		end % solveAlgebraic()
		
		function inverse = invert(obj, gridName)
			% inverse solves the function representet by the quantity for
			% its variable, for instance, if obj represents y = f(x), then
			% invert returns an object containing x = f^-1(y).
			% Yet, this is only implemented for obj with one variable
			% (grid).
			if iscell(gridName)
				% fixme: by default the first gridName is chosen as new
				% name. This works because the functions is only written
				% for quantities with one variable.
				gridName = gridName{1};
			end
			
			assert(numel(obj(1).gridName) == 1);
			assert(isequal(size(obj), [1, 1]));
			inverse = quantity.Discrete(repmat(obj(1).grid{obj(1).domain.gridIndex(gridName)}(:), [1, size(obj)]), ...
				'size', size(obj), ...
				'domain', quantity.Domain([obj(1).name], obj.on()), ...
				'name', gridName);
		end % invert()
		
		function solution = solveDVariableEqualQuantity(obj, varargin)
			% solves the first order ODE
			%	dvar / ds = obj(var(s))
			%	var(0) = ic
			% to obtain var(s, ic) depending on both the argument s and the initial 
			% condition ic. Herein, obj may only depend on one variable / gridName / ...
			% domain.
			assert(numel(obj(1).domain) == 1, ...
				'this method is only implemented for quanitities with one domain');
			
			myParser = misc.Parser();
			myParser.addParameter('initialValueGrid', obj(1).grid{1});
			myParser.addParameter('variableGrid', obj(1).grid{1});
			myParser.addParameter('newGridName', 's');
			myParser.parse(varargin{:});
			
			variableGrid = myParser.Results.variableGrid(:);
			myGridSize = [numel(variableGrid), ...
				numel(myParser.Results.initialValueGrid)];
			
			% the time (s) vector has to start at 0, to ensure the IC. If
			% variableGrid does not start with 0, it is separated in
			% negative and positive parts and later combined again.
			positiveVariableGrid = [0; variableGrid(variableGrid > 0)];
			negativeVariableGrid = [0; flip(variableGrid(variableGrid < 0))];
			
			% solve ode for every entry in obj and for every initial value
			odeSolution = zeros([myGridSize, numel(obj)]);
			for it = 1:numel(obj)
				for icIdx = 1:numel(myParser.Results.initialValueGrid)
					resultGridPositive = [];
					odeSolutionPositive = [];
					resultGridNegative = [];
					odeSolutionNegative = [];
					if numel(positiveVariableGrid) > 1
						[resultGridPositive, odeSolutionPositive] = ...
							ode45(@(y, z) obj(it).on(z), ...
							positiveVariableGrid, ...
							myParser.Results.initialValueGrid(icIdx));
					end
					if numel(negativeVariableGrid) >1
						[resultGridNegative, odeSolutionNegative] = ...
							ode45(@(y, z) obj(it).on(z), ...
							negativeVariableGrid, ...
							myParser.Results.initialValueGrid(icIdx));
					end
					if any(variableGrid == 0)
						resultGrid = [flip(resultGridNegative(2:end)); 0 ; resultGridPositive(2:end)];
						odeSolution(:, icIdx, it) = [flip(odeSolutionNegative(2:end)); ...
							myParser.Results.initialValueGrid(icIdx); odeSolutionPositive(2:end)];
					else
						resultGrid = [flip(resultGridNegative(2:end)); resultGridPositive(2:end)];
						odeSolution(:, icIdx, it) = [flip(odeSolutionNegative(2:end)); ...
							odeSolutionPositive(2:end)];
					end
					assert(isequal(resultGrid(:), variableGrid(:)));
				end
			end
			
			% return result as quantity-object
			solution = quantity.Discrete(...
				reshape(odeSolution, [myGridSize, size(obj)]), ...
				'domain', [quantity.Domain(myParser.Results.newGridName, variableGrid), ...
					quantity.Domain('ic', myParser.Results.initialValueGrid)], ...
				'size', size(obj), ...
				'name', "solve(" + obj(1).name + ")");
		end % solveDVariableEqualQuantity()
		
		function solution = subs(obj, gridName2Replace, values)
			% SUBS substitute variables of a quantity
			%	solution = SUBS(obj, NEWDOMAIN), replaces the original
			%	domain of the object with the new domain specified by
			%	NEWDOMAIN. NEWDOMAIN must have the same grid name as the
			%	original domain.
			%
			%	solution = SUBS(obj, GRIDNAMES2REPLACE, VALUES) replaces
			%	the domains which are specified by GRIDNAMES2REPLACE by
			%	VALUES. GRIDNAMES2REPLACE must be a cell-array with the
			%	names of the domains or an object-array with
			%	quantity.Domain objects which should be replaced by VALUES.
			%	VALUES must be a cell-array of the new values or new grid
			%	names.
			%
			%	Example: 
			%		q = q.subs('z', 't')
			%		will replace the domain with the name 'z' by a domain
			%		with the name 't' but with the same discretization.
			%		q = q.subs('z', linspace(0,1)')
			%		will replace the grid of domain with the name 'z' by
			%		the new grid specified by linspace.
			if nargin == 1 || isempty(gridName2Replace)
				% if gridName2Replace is empty, then nothing must be done.
				solution = obj;
			elseif isempty(obj)
				% if the object is empty, nothing must be done.
				solution = obj;
			else
				% input checks
				if nargin == 2
					assert(isa(gridName2Replace, 'quantity.Domain'), 'If only two parameters are specified, the second parameter must be a quantiy.Domain');

					values = {gridName2Replace.grid};
					gridName2Replace = {gridName2Replace.name};
					
				elseif nargin == 3
					gridName2Replace = misc.ensureIsCell(gridName2Replace);
					for k = 1:numel( gridName2Replace )
						if isa(gridName2Replace{k}, 'quantity.Domain')
							gridName2Replace{k} = gridName2Replace{k}.name;
						end
					end						
					values = misc.ensureIsCell(values);
				end
				
				assert(numel(values) == numel(gridName2Replace), ...
					'gridName2Replace and values must be of same size');
				
				% here substitution starts:
				% The first (gridName2Replace{1}, values{1})-pair is
				% replaced. If there are more cell-elements in those inputs
				% then subs() is called again for the remaining pairs
				% (gridName2Replace{2:end}, values{2:end}).
				if ischar(values{1}) || isstring(values{1})
					% if values{1} is a char-array, then the gridName is
					% replaced
					if any(strcmp(values{1}, gridName2Replace(2:end)))
						% in the case if a quantity f(z, zeta) should be
						% substituted like subs(f, {z, zeta}, {zeta, z})
						% this would cause an error, since after the first
						% substituion subs(f, z, zeta) the result would be
						% f(zeta, zeta) -> the 2nd subs(f, zeta, z) will
						% result in f(z, z) and not in f(zeta, z) as
						% intended. This is solved, by an additonal
						% substitution: 
						%	f.subs(z,zetabackUp).subs(zeta,z).subs(zetabackUp,zeta)
						values{end+1} = values{1};
						gridName2Replace{end+1} = [gridName2Replace{1}, 'backUp'];
						values{1} = gridName2Replace{end};
					end
					if isequal(values{1}, gridName2Replace{1})
						% replace with same variable... everything stays the
						% same.
						% Do not use "return", since, later subs might need to be
						% called recursively!
						newValue = obj.on();
						newDomain = obj(1).domain;
					elseif any(strcmp(values{1}, obj(1).gridName))
						% if for a quantity f(z, zeta) this method is
						% called with subs(f, zeta, z), then g(z) = f(z, z)
						% results, hence the dimensions z and zeta are
						% merged.
						domainIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
							obj(1).domain.gridIndex(values{1})];
						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);
						end
						newValue = misc.diagNd(obj.on(newDomainForOn), domainIndices);
						newDomain = [newDomainForOn(domainIndices(2)), ...
							newDomainForOn(all(1:1:numel(newDomainForOn) ~= domainIndices(:)))];
					else
						% this is the default case. just grid name is changed.
						newDomain = obj(1).domain;
						newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
							quantity.Domain(values{1}, ...
							obj(1).domain(obj(1).domain.gridIndex(gridName2Replace{1})).grid);
						newValue = obj.on();
					end
					
				elseif isnumeric(values{1}) && numel(values{1}) == 1
					% if values{1} is a scalar, then obj is evaluated and
					% the resulting quantity looses that spatial grid and
					% gridName
					newDomain = obj(1).domain;
					newDomain = newDomain(~strcmp(gridName2Replace{1}, [newDomain.name]));
					% newGrid is the similar to the original grid, but the
					% grid of gridName2Replace is removed.
					newGridSize = newDomain.gridLength();
					% newGridForOn is the similar to the original grid, but
					% the grid of gridName2Replace is set to values{1} for
					% evaluation of obj.on().
					newGridForOn = obj(1).grid;
					newGridForOn{obj(1).domain.gridIndex(gridName2Replace{1})} = values{1};
					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.
					newDomain = obj(1).domain;
					newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
						quantity.Domain(gridName2Replace{1}, values{1});
					newValue = obj.on(newDomain);
				else
					error('value must specify a gridName or a gridPoint');
				end
				if isempty(newDomain)