Domain.m

classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
		handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
	%DOMAIN class to describes a range of values on which a function can be defined.
	
	properties (SetAccess = protected)
		% The discrete points of the grid for the evaluation of a
		% continuous quantity. For an example, the function f(x) should be
		% considered on the domain x \in X = [0, 1]. Then, a grid can be
		% generated by X_grid = linspace(0, 1).
		grid (:, 1) double;
		
		% a speaking name for this domain; Should be unique, so that the
		% domain can be identified by the name.
		name (1,1) string = "";
		
		% number of discretization points == gridLength
		n (1,1) double;		
		
	end
	
	properties (Dependent)
		lower;		% lower bound of the domain == grid(1)
		upper;		% upper bound of the domain == grid(end)
	end

	properties (Dependent, Hidden)
		ones;
	end
	
	methods
		function obj = Domain(name, grid)
			%DOMAIN initialize the domain
			
			arguments
				name string = "";
				grid (:,1) double = [];
			end
			
			if nargin >= 1
				obj.grid = grid;
				obj.name = name;
				obj.n = length(obj.grid);
			end
		end % Domain()
		
		function it = get.ones(obj)
			it = builtin('ones', size(obj.grid));
		end
		
		function lower = get.lower(obj)
			% minimum value of grid = grid(1), since grid must be ascending
			lower = obj.grid(1);
		end % get.lower
		
		function upper = get.upper(obj)
			% maximum value of grid = grid(end), since grid must be ascending
			upper = obj.grid(end);
		end % get.upper()
		
		function d = quantity.EquidistantDomain(obj)
			d = quantity.Domain(obj.name, obj.grid);
		end
		
		function s = string(obj)
			s = obj.name;
		end		
		
		function quan = Discrete(obj)
			% DISCRETE casts the domain into a quantity.Discrete with the values of obj.grid.
			arguments
				obj (1, 1) quantity.Domain;
			end
			quan = quantity.Discrete(obj.grid, obj, 'name', obj.name);
		end % quantity.Discrete()
		
		function quan = Symbolic(obj)
			% SYMBOLIC casts the domain into a quantity.Symbolic with the grid name as symbolic
			% value.
			arguments
				obj (1, 1) quantity.Domain;
			end
			quan = quantity.Symbolic(sym(obj.name), obj, 'name', obj.name);
		end % Symbolic()
		
		function q = Function(obj)
			arguments
				obj (1,1) quantity.Domain;
			end
			
			q = quantity.Function(eval("@(" + obj.name + ")" + obj.name), obj, 'name', obj.name);
		end
		
	end
		
	methods (Access = public)
		
		function d = split(obj, splittingPoints, optArgs)
			%SPLIT splits the domain into subdomains
			% d = split(OBJ, SPLITTINGPOINTS) splits the domain into the
			% subdomains specified by SPLITTINGPOINTS. To be particular, a
			% domain (0,1) called with the SPLITTINGPOINTS 0.5 will be
			% split into two domains (0,0.5) and (0.5,1). To consider the
			% case if the SPLITTINGPOINTS are not exactly contained in the
			% grid, the closest grid point of OBJ.grid will be chosen as
			% splitting point.
			arguments
				obj
				splittingPoints double;
				optArgs.AbsTol (1,1) double = 1e-3;
				optArgs.warning (1,1) logical = true;
			end
			
			d(1:numel(splittingPoints)+1) = quantity.Domain();
			idx = zeros(numel(splittingPoints)+1,1);
			idx(1) = 1;
			
			for k = 1:numel(splittingPoints)
				
				% Search the grid point closest to the k-th splittingPoint
				[~, idx(k+1)] = min(abs((obj.grid - splittingPoints(k))));
				
				if optArgs.warning
					% verify that the distance of the chosen grid point to the
					% k-th splitting point is within some tolerances:
					delta = abs( obj.grid(idx(k+1)) - splittingPoints(k));
					if delta > optArgs.AbsTol
						warning('DOMAIN:split', ...
							'The deviation of the grid from desired splitting point %d is %d\n', splittingPoints(k), delta);
					end
				end
				
				d(k) = quantity.Domain(obj.name, obj.grid(idx(k):idx(k+1)));
			end
			
			d(k+1) = quantity.Domain(obj.name, obj.grid(idx(k+1):end));
			
		end
		
		function [i, logIdx] = contains(obj, d)
			
			logIdx = false(obj.n,1);
			
			if isa(d, 'quantity.Domain')
				d = d.grid;
			end
			
			grd = obj.ndgrid;
			grd = grd{1};
			
			for i = 1:numel(d)
				logIdx = logIdx | grd == d(i);
			end
			
			i = any(logIdx(:));
			
		end
		
		function isEq = isequal(obj, a)
			% isequal checks if two domains are equal.
			
			isEq = numel(obj) == numel(a);
			
			for k = 1:numel(obj)
				isEq = isEq && all(obj(k).name == a(k).name);
				isEq = isEq && isequal(obj(k).grid, a(k).grid);
			end
		end % isequal()
		
		function d = rename(obj, newName, oldName)
			% RENAME renames the name of the domain OBJ. 
			%
			%	d = rename(obj, newName) replaces all names of obj with the names specified by the
			%	string array newName
			%
			%	d = rename(obj, newName, oldName) replaces the names of obj specified by the string
			%	array oldName with the corresponding elements of the string array newName
			
			arguments
				obj quantity.Domain;
				newName string;
				oldName string = [obj.name];
			end % arguments
			
			assert(numel(newName) == numel(oldName), ...
				"number of new names must be equal to number of old names.");
			d = copy(obj);
			for it = 1 : numel(newName)
				d(obj.index(oldName(it))).name = newName{it};
			end % for it = 1 : numel(newName)
			assert(misc.isunique([d.name]), "the resulting domain must contain unique names");
		end % rename

		function nd = ndims(obj)
			%NDIMS number of dimensions of the domain specified by the
			%object-array.
			nd = size(obj(:), 1);
		end
		
		function n = numGridElements(obj)
			% NUMGRIDLEMENTS returns the number of the elements of the grid
			n = prod([obj.n]);
		end % numGridElements()
		
		function s = gridLength(obj)
			%GRIDLENGTH number of discretization points for each grid in the
			%object-array.
			s = [obj.n];
		end
		
		function d = find(obj, searchName)
			% FIND a domain in object-array of quantity.Domain
			%
			%	d = find(OBJ, NAME) returns the domain with the name NAME
			%	of the object-array of quantity.Domain OBJ.
			
			%	d = find(OBJ, NAME1, NAME2, ...) returns the domain array with the
			%	names NAME1 and NAME2 of the object-array of quantity.Domain OBJ.
			%
			%	d = find(OBJ, DOMAIN) find the domain with the name
			%	DOMAIN.name in the object-array of quantity.Domain OBJ. It
			%	only looks for the domain name and not the grid.
			%
			%	d = find(OBJ, { NAME or DOMAIN}) the NAME or the DOMAIN to
			%	be found can also be defined by cell-array.
			arguments
				obj quantity.Domain;
			end
			arguments (Repeating)
				searchName string;
			end
					
			idx = obj.index( searchName );
			
			if idx == 0
				d = [];
			else 
				d = obj( idx );
			end
		end
		
		function rest = remove(obj, toBeRemoved)
			% REMOVE removes the domain specified by toBeRemoved from the quantity.Domain array obj.
			%
			%	rest = remove(obj, toBeRemoved) removes the domain specified by toBeRemoved from the 
			%	quantity.Domain array obj. toBeRemoved may be a string or a domain or a array of
			%	those.
			
			[~, logicalIdxToBeRemoved] = obj.index(toBeRemoved);
			
			rest = obj(~logicalIdxToBeRemoved);
		end % remove()
		
		function obj = replace(obj, toBeReplaced)
			% REPLACE replaced a domain with another domain.
			%
			%	obj = replace(obj, toBeReplaced) replaces the element of the obj with the same name
			%		as the toBeReplaced domain with the new domain toBeRelaced.
			
			thisIdx = obj.index(toBeReplaced);
			if thisIdx > 0
				obj(thisIdx) = toBeReplaced;
			end
		end % replace()
		
		function [idx, log] = index( obj, searchName, position)
			%% GRIDINDEX returns the index of the grid
			% [idx, log] = index(obj, names) searches in the name
			% properties of obj for the "names" and returns its index as
			% "idx" and its logical index as "log"
			arguments
				obj
				searchName = [obj.name];
				position string = "all";
			end
			
			if isa(searchName, 'quantity.Domain')
				searchName = [searchName.name];
			end
			
			objNames = [obj.name];
			
			if isempty(objNames)
				idx = 0;
				log = [];
			else
				searchName = misc.ensureString( searchName );
				
				log = false(size(objNames));
				counter = 1:numel(objNames);
				idx = [];
				for k = 1 : numel(searchName)
					log_ = objNames == searchName(k);
					log = log | log_;
					idx = [idx counter( log_ )];
				end
				if isempty(idx)
					idx = 0;
				end
			end
			
			if position == "first"
				idx = idx(1);
			end
		end
		
		function [idx, log] = gridIndex(obj, varargin)
			warning("DEPRICATED: Use quantity.Domain/index instead")
			[idx, log] = obj.index(varargin{:});
		end % index
		
		function l = ne(obj, obj2)
			% ~= Not equal.
			l = ~(obj == obj2);
		end % ne()
				
		function [joinedDomain, index] = join(domain1, domain2)
			% JOIN combine two domains
			% [joinedDomain, index] = join(domain1, domain2) combines the domain1 and domain2, such
			% that every gridName only occurs once and that the finer grid of both is used.
			arguments
				domain1 quantity.Domain = quantity.Domain.empty
				domain2 quantity.Domain = quantity.Domain.empty
			end
			
			if nargin == 1 && length(domain1) == 1
				joinedDomain = domain1;
				index = 1;
				return;
			end
			
			% remove empty domains:
			emptyIdx = arrayfun(@(c) ~isempty(c), domain1);
			domain1 = domain1(emptyIdx);
			emptyIdx = arrayfun(@(c) ~isempty(c), domain2);
			domain2 = domain2(emptyIdx);
						
			[joinedGrid, index] = unique([domain1.name, domain2.name], 'stable');
			
			joinedDomain(1 : numel(joinedGrid)) = quantity.Domain();
			
			% check for each grid if it is in the domain of obj1 or obj2 or
			% both
			for i = 1 : numel(joinedGrid)
				currentGridName = joinedGrid(i);
				[index1, logicalIndex1] = domain1.index(currentGridName, "first");
				[index2, logicalIndex2] = domain2.index(currentGridName, "first");
								
				% Check if a domain is in both domains:
				%	-> then take the finer one of both
				if any(logicalIndex1) && any(logicalIndex2)
					tempDomain1 = domain1(index1);
					tempDomain2 = domain2(index2);
					
					assert( numeric.near( tempDomain1.lower, tempDomain2.lower), 'Grids must have same domain for join')
					assert( numeric.near( tempDomain1.upper, tempDomain2.upper), 'Grids must have same domain for join')

					if tempDomain1.gridLength > tempDomain2.gridLength
						joinedDomain(i) = tempDomain1;
					else
						joinedDomain(i) = tempDomain2;
					end
				% If it is not in both, -> just take the normal grid
				elseif any(logicalIndex1)
					joinedDomain(i) = domain1(index1(1));
				elseif any(logicalIndex2)
					joinedDomain(i) = domain2(index2);
				end
				
			end
		end % join()
		
		function i = isempty(obj)
			i = any(size(obj) == 0);
			i = i || any( cellfun(@isempty, {obj.grid} ) );
		end % isempty()
		
		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 generate a mesh for this domain
			% matGrid = ndgrid(obj) calles ndgrid for the default grid, if no other grid
			% is specified. Empty grid as input returns empty cell as
			% result.
			if isempty(obj)
				matGrid = {};
			else
				grids = {obj.grid};
				[matGrid{1:obj.ndims}] = ndgrid(grids{:});
			end
		end % ndgrid()
		
		function objReal = real(obj)
			% real returns a quantity.Domain objReal whichs grid is only the real
			% part of the possibly complex grid in obj.
			objReal = quantity.Domain(obj.name, real(obj.grid));
		end % real()
		
		function objImag = imag(obj)
			% imag returns a quantity.Domain objImag whichs grid is only the
			% imaginary part of the possibly complex grid in obj.
			objImag = quantity.Domain(obj.name, imag(obj.grid));
		end % imag()
		
		function [idx, newDomain] = getPermutationIdx(obj, order)
			
			if isa(order, 'quantity.Domain')
				names = [order.name];
			elseif iscell(order) || ischar(order) || isstring(order)
				names = misc.ensureString( order );
			else
				error('the input parameter order must be a array of quantity.Domain objects or a string-array')
			end
						
			idx = cellfun(@(v) obj.index(v), names); % #todo@domainNameString
			
			if isa(order, 'quantity.Domain')
				newDomain = order;
			else			
				newDomain = obj(idx);
			end
		end % getPermutationIdx()
		
		
		function [newObj, I] = 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.
			
			if nargin == 2 && strcmp(varargin{1}, 'descend')
				descend = 1;
			elseif nargin == 1 || strcmp(varargin{1}, 'ascend')
				descend = 0;
			else
				error(['Unknown sorting order' char( varargin{1} )])
			end
			
			% only sort the grids if there is something to sort
			if obj.ndims > 1
				gridNames = [obj.name];
				
				% this is the default case for ascending alphabetical
				% order
				[sortedNames, I] = sort(gridNames);
				
				% if descending: flip the order of the entries
				if descend
					sortedNames = flip(sortedNames);
					I = flip(I);
				end
				
				% sort the grid entries
				myGrids = {obj.grid};
				myGrids = myGrids(I);
				
				newObj(1:length(myGrids)) = quantity.Domain();
				for k = 1:length(myGrids)
					newObj(k) = quantity.Domain(sortedNames{k}, myGrids{k});
				end
			else
				newObj = obj;
			end
		end% sort()
		
		function result = isequidistant(obj, AbsTol)
			% ISEQUIDISTANT(obj) returns true, if the grid is homogenious, i.e. the
			% grid points are equidistant and false else. Also works with
			% quantity.Domain arrays.
			% ISEQUIDISTANT(obj, AbsTol) allows perturbations in the spacing of
			% up to the absolute maximum defined by AbsTol. Default: 1e3*eps.
			
			arguments
				obj quantity.Domain;
				AbsTol = 1e3*eps;
			end % arguments
			
			result = true(size(obj));
			for it = 1 : numel(result)
				spacing = diff(obj(it).grid);
				if any(diff(obj(it).grid) - mean(spacing) > AbsTol)
					result = false;
				end
			end % for it = 1 : numel(result)
		end % isequidistant()
		
		function f1 = plot(obj)
			% PLOT create plot of grid points (cross makeing) and slope (dashed).
			% This can be insightful for non-equidistant grids
			arguments
				obj quantity.Domain;
			end
			
			f1 = figure();
			numDomains = numel(obj);
			for it = 1 : numDomains
				subplot(1, numDomains, it);
				equidistantSpacing = diff([obj(it).lower, obj(it).upper]) / (obj(it).n-1);
				plot(linspace(obj(it).lower, obj(it).upper, obj(it).n), obj(it).grid, 'x');
				hold on;
				plot(linspace(obj(it).lower + equidistantSpacing / 2, ...
					obj(it).upper + equidistantSpacing / 2, obj(it).n-1), ...
					diff(obj(it).grid) / equidistantSpacing, '--');
				xlim(obj(it).grid([1, end]));
				xlabel(obj(it).name);
				ylabel(obj(it).name + " (x), " ...
					+ obj(it).name + "_{" + obj(it).name + "} (- -)");
			end % for it = 1 : numDomains
		end % plot()
				
	end % methods (Public)
	
	methods (Access = protected)
		
		function s = getPropertyGroups(obj)
			% Function to display the correct values
			
			if isempty(obj)
				s = getPropertyGroups@matlab.mixin.CustomDisplay(obj);
				return;
			else
				s = getPropertyGroups@matlab.mixin.CustomDisplay(obj(1));
			end
			
			if numel(obj) ~= 1
				
				grids = {obj.grid};
				
				sizes = cellfun(@(g)size(g), grids, 'UniformOutput', false);							
								
				s.PropertyList.grid = ...
					[sprintf('%ix%i,  ', cell2mat(sizes)) sprintf('\b\b\b')];
				s.PropertyList.name = ...
					[sprintf('%s,  ', obj.name) sprintf('\b\b\b')];
				s.PropertyList.lower = ... 
					[sprintf('%d,  ', obj.lower) sprintf('\b\b\b')];
				s.PropertyList.upper = ... 
					[sprintf('%.3f,  ', obj.upper) sprintf('\b\b\b')];
				s.PropertyList.n = ...
					[sprintf('%i,  ', obj.n) sprintf('\b\b\b')];
			end
		end % getPropertyGroups
		
	end % methods (Access = protected)
	
	methods (Static)	
		function d = gridCells2domain(grids, gridNames)
			grids = misc.ensureIsCell(grids);
			if isempty(gridNames)
				gridNames = {};
			else
				gridNames = misc.ensureIsCell(gridNames);
			end
			
			assert(numel(grids) == numel(gridNames), ...
				"Number of grid vectors is not equal to number of grid names")

			d = quantity.Domain.empty();
			
			for k = 1:length(grids)
				d = [d quantity.Domain(gridNames{k}, grids{k})];
			end
		end % gridCells2domain()
		
		function d = defaultDomain(gridSize, name)
			%DEFAULTDOMAIN generate default domain
			% d = defaultDomain() will create a default domain with a grid
			%   from 0 to 1 with 100 points and the default name x_1
			% d = defaultDomain(GRIDSIZE) will create a default domain with
			% a default grid of size specified in GRIDSIZE. The domains
			% will have the default names
			%   x_1, x_2, ...
			% d = defaultDomain(GRIDSIZE, NAME) will create a default
			% domain with a default grid of size specified in GRIDSIZE. The
			% domains will have the names specified in NAME.
			arguments
				gridSize double = 100;
				name string = quantity.Domain.defaultGridNames(gridSize);
			end
			
			% generate a default gird with given sizes
			d(1:length(gridSize)) = quantity.Domain();
			for k = 1:length(gridSize)
				
				o = ones(1, length(gridSize) + 1); % + 1 is required to deal with one dimensional grids
				o(k) = gridSize(k);
				O = ones(o);
				O(:) = linspace(0, 1, gridSize(k));
				
				d(k) = quantity.Domain(name{k}, O);
			end
		end
		
		function g = defaultGrid(varargin)
			%DEFAULTGRID generate default grid for a quantity
			warning('DEPRICATED! Use quantity.Domain.defaultDomain instead.');
			g = quantity.Domain.defaultDomain(varargin{:});
		end
		
		function name = defaultGridNames(gridSize)
			%DEFAULTGRIDNAMES generates default grid names
			% name = defaultGridNames(GRIDSIZE) generates default grid
			% names for a grid of size GRIDSIZE. The default names will be
			% x_1, x_2, ...
			name = cell(1, length(gridSize));
			for k = 1:length(gridSize)
				name{k} = ['x_' num2str(k)];
			end
		end
		
		
		function [myDomain, unmatched] = parser(varargin)
			%% domain parser
			domainParser = misc.Parser();
			domainParser.addParameter('domain', {}, @(g) isa(g, 'quantity.Domain'));
			domainParser.addParameter('gridName', '');
			domainParser.addParameter('grid', [], @(g) isnumeric(g) || iscell(g));
			domainParser.parse(varargin{:});
			
			if domainParser.isDefault('domain')
				% case: the gridNames and the gridValues are defined:
				%	-> initialize quantity.Domain objects with the
				%	specified values

				myGridName = misc.ensureString(domainParser.Results.gridName);
				myGrid = misc.ensureIsCell(domainParser.Results.grid);

				assert(isequal(numel(myGrid), numel(myGridName)), ...
					'number of grids and gridNames must be equal');

				% initialize the domain objects
				myDomain = quantity.Domain.empty();
				for k = 1:numel(myGrid)
					myDomain(k) = quantity.Domain(myGridName(k), myGrid{k});
				end
			elseif ~domainParser.isDefault('domain')
				% case: the domains are specified as domain
				% objects.
				myDomain = domainParser.Results.domain;	
			else			
				error('No domain is specified! A domain is obligatory for the initialization of a quantity.')
			end	
			
			assert(misc.isunique([myDomain.name]), ...
				'The names of the domain must be unique');
			
			unmatched = domainParser.UnmatchedNameValuePair;
		end % parser()
		
	end % methods Static
end % classdef