Merge branch '14-quantity-domain' into 'master'

Resolve "quantity.Domain"

Closes #14

See merge request !5

+misc/Parser.m

@@ -6,13 +6,15 @@ classdef Parser < inputParser
 		
         
     methods
-        function obj = Parser()
+        function obj = Parser(varargin)
             obj = obj@inputParser();
             
             % set default properties
             obj.KeepUnmatched = true;
             obj.CaseSensitive = true;
             obj.PartialMatching = false;
+			
+			
 		end % Parser()
         
         function i = isDefault(obj, name)

+misc/ensureIsCell.m

+function value = ensureIsCell(value)
+%ENSUREISCELL ensures that the value is a cell.
+% c = ensureIsCell(value) checks if value is a cell. If it is not a cell,
+% it is converted as a cell.
+
+if ~iscell(value)
+	value = {value};
+end	
+	
+end
+

+quantity/BasicVariable.m

@@ -135,10 +135,9 @@ function obj = BasicVariable(valueContinuous, varargin)
 					Di(j) = obj_j.derivatives(k+1);
 				end
 			end
-			Di = reshape(Di, [ 1 size(obj)]);
 			D(i,:) = Di;
 		end
-		
+		D = reshape(D, [numel(K) size(obj)]);
 	end
 
 	%% PROPERTIES

+quantity/Domain.m

+classdef Domain < handle & matlab.mixin.CustomDisplay 
+	%DOMAIN class to describes a range of values on which a function can be defined.
+	
+	% todo:
+	%	* EquidistantDomain
+	%	* multi dimensional
+	
+	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 double {mustBeReal, mustBeFinite, mustBe.ascending};
+		% TODO add mustBeVector
+		
+		% a speaking name for this domain; Should be unique, so that the
+		% domain can be identified by the name.
+		name char;
+	end
+	
+	properties (Dependent)
+		n;			% number of discretization points
+		lower;		% lower bound of the domain
+		upper;		% upper bound of the domain
+	end
+
+	properties (Dependent, Hidden)
+		ones;
+	end
+	
+	methods
+		function obj = Domain(name, grid)
+			%DOMAIN initialize the domain
+			%
+			if nargin >= 1
+				obj.grid = grid(:);
+				obj.name = name;
+			end
+		end
+		
+		function i = get.ones(obj)
+			i = ones(size(obj.grid));
+		end
+		
+		function n = get.n(obj)
+			if isempty(obj.grid)
+				n = [];
+			else
+				n = length(obj.grid);
+			end
+		end
+		
+		function lower = get.lower(obj)
+			lower = min( obj.grid );
+		end
+		
+		function upper = get.upper(obj)
+			upper = max( obj.grid );
+		end
+		
+		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
+		
+		function s = gridLength(obj)
+			%GRIDLENGTH number of discretization points for each grid in the
+			%object-array.
+			s = [obj.n];
+		end
+		
+% 		function s = eq(obj, obj2)
+% 			
+% 			if isempty(obj) && isempty(obj2)
+% 				% if both are empty -> they are equal
+% 				s = true;
+% 			elseif isempty(obj) || isempty(obj2)
+% 				% if only one domain is empty -> they are not equal
+% 				s = false;
+% 			else
+% 				% if both are not empty: check if the grids and the
+% 				% gridNames are equal
+% 				s = all(obj.gridLength == obj2.gridLength);
+% 			
+% 			for k = 1:numel(obj)	
+% 				s = s && all(obj(k).grid == obj2(k).grid);
+% 				s = s && strcmp(obj(k).name, obj2(k).name);
+% 			end
+% 			end
+% 		end
+% 		
+		function l = ne(obj, obj2)
+			l = ~(obj == obj2);
+		end
+				
+		function [joinedDomain, index] = join(domain1, domain2)
+			%% gridJoin combines the grid and gridName of two objects (obj1,
+			% obj2), 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
+			
+			[joinedGrid, index] = unique({domain1.name, domain2.name}, 'stable');
+			
+			joinedDomain = quantity.Domain.empty;
+			
+			% #todo@domain: the gird comparison seems to be very
+			% complicated
+			
+			% 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.gridIndex(currentGridName);
+				[index2, logicalIndex2] = domain2.gridIndex(currentGridName);
+				
+				%
+				% 				if ~any(logicalIndex1)
+				% 					joinedDomain(i) = obj2(index2);
+				% 				elseif ~any(logicalIndex2)
+				% 					joinedDomain(i) = obj1(index1);
+				% 				else
+				%
+				% 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(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
+					assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
+
+					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);
+				elseif any(logicalIndex2)
+					joinedDomain(i) = domain2(index2);
+				end
+				
+			end
+		end
+		
+		function [joinedDomain, index] = gridJoin(obj1, obj2)
+			[joinedDomain, index] = obj1.join(obj2);
+		end % gridJoin()
+				
+		function [idx, logicalIdx] = gridIndex(obj, names)
+			%% GRIDINDEX returns the index of the grid
+			% [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"
+			arguments
+				obj
+				names = {obj.name};
+			end
+			
+			if isa(names, 'quantity.Domain')
+				names = {names.name};
+			end
+			
+			names = misc.ensureIsCell(names);
+			
+			idx = zeros(1, length(names));
+			nArgIdx = 1:obj.ndims();
+			
+			logicalIdx = false(1, obj.ndims());
+			
+			for k = 1:length(names)
+				log = strcmp({obj.name}, names{k});
+				logicalIdx = logicalIdx | log;
+				if any(log)
+					% #todo@domain: if the index for a grid name is searched in a
+					% list with multiple grids but same names, only the first
+					% occurrence is returned. this should be changed...
+					tmpFirst = nArgIdx(log);
+					idx(k) = tmpFirst(1);
+				else
+					idx(k) = 0;
+				end
+			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
+			% result.
+			if isempty(obj)
+				matGrid = {};
+			else
+				grids = {obj.grid};
+				[matGrid{1:obj.ndims}] = ndgrid(grids{:});
+			end
+		end % ndgrid()
+		
+		function [idx, newDomain] = getPermutationIdx(obj, order)
+			
+			if isa(order, 'quantity.Domain')
+				names = {order.name};
+			elseif ischar([order{:}])
+				names = order;
+			else
+				error('the input parameter order must be a array of quantity.Domain objects or a cell-array with string')
+			end
+						
+			idx = cellfun(@(v) obj.gridIndex(v), names);
+			
+			if isa(order, 'quantity.Domain')
+				newDomain = order;
+			else			
+				newDomain = obj(idx);
+			end
+		end
+		
+		
+		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()			
+	end
+	
+	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 % 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(gridNames{k}, grids{k})];
+			end
+		end
+		
+		function g = defaultGrid(gridSize, name)
+			
+			if nargin == 1
+				% If no names are specified, chose x_1, x_2, ... as default
+				% names.
+				name = cell(1, length(gridSize));
+				for k = 1:length(gridSize)
+					name{k} = ['x_' num2str(k)];
+				end
+			end
+			
+			% generate a default gird with given sizes
+			g = quantity.Domain.empty();
+			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));
+				
+				g(k) = quantity.Domain(name{k}, O);
+			end
+		end
+		
+		function [myDomain, unmatched] = parser(varargin)
+			%% domain parser
+			domainParser = misc.Parser();
+			domainParser.addParameter('domain', {}, @(g) isa(g, 'quantity.Domain'));
+			domainParser.addParameter('gridName', '', @(g) ischar(g) || iscell(g));
+			domainParser.addParameter('grid', [], @(g) isnumeric(g) || iscell(g));
+			domainParser.parse(varargin{:});
+
+			if domainParser.isDefault('domain') && ...
+					( domainParser.isDefault('grid') || ...
+					  domainParser.isDefault('gridName') )
+				% case 1: nothing about the grid is defined
+					error('No domain is specified! A domain is obligatory for the initialization of a quantity.')
+			elseif domainParser.isDefault('domain')
+				% case 3: the gridNames and the gridValues are defined:
+				%	-> initialize quantity.Domain objects with the
+				%	specified values
+
+				myGridName = misc.ensureIsCell(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
+			else
+				% else case: the domains are specified as domain
+				% objects.
+				myDomain = domainParser.Results.domain;
+			end	
+			
+			assert(numel(myDomain) == numel(unique({myDomain.name})), ...
+				'The names of the domain must be unique');
+			
+			unmatched = domainParser.UnmatchedNameValuePair;
+		end
+		
+	end
+end
+

+quantity/EquidistantDomain.m

+classdef EquidistantDomain < quantity.Domain
+	%EQUIDISTANTDOMAIN class to handle the discretization of the range of
+	%definition of a function. The discretization points are equally
+	%distributed over the domain.
+	
+	properties
+		Property1
+	end
+	
+	methods
+		function obj = untitled(inputArg1,inputArg2)
+			%UNTITLED Construct an instance of this class
+			%   Detailed explanation goes here
+			obj.Property1 = inputArg1 + inputArg2;
+		end
+		
+		function outputArg = method1(obj,inputArg)
+			%METHOD1 Summary of this method goes here
+			%   Detailed explanation goes here
+			outputArg = obj.Property1 + inputArg;
+		end
+	end
+end
+

+quantity/Function.m

@@ -74,57 +74,26 @@ classdef Function < quantity.Discrete
 		%-----------------------------
 		% --- overloaded functions ---
 		%-----------------------------
-		function value = on(obj, myGrid, myGridName)
-			% evaluates obj.valueContinuous function and returns an array
-			% with the dimensions (n, m, ..., z_disc). n, m, ... are the
-			% array dimensions of obj.valueContinuous and z_disc is the
-			% dimensions of the grid.
-			if nargin == 1
-				gridSize = obj(1).gridSize;
-				if gridSize > 2
-					value = cat(length(gridSize),obj.valueDiscrete);
-				else % important to include this case for empty grids!
-					value = [obj.valueDiscrete]; % = cat(2,obj.valueDiscrete)
-				end
-				value = reshape(value, [gridSize, size(obj)]);
-				return
-			elseif nargin >= 2 && ~iscell(myGrid)
-				myGrid = {myGrid};
-			end
-			gridPermuteIdx = 1:obj(1).nargin;
-			if isempty(gridPermuteIdx)
-				gridPermuteIdx = 1;
-			end
-			if nargin == 3
-				if ~iscell(myGridName)
-					myGridName = {myGridName};
-				end
-				assert(numel(myGrid) == numel(myGridName), ...
-					['If on() is called by using gridNames as third input', ...
-					', then the cell-array of grid and gridName must have ', ...
-					'equal number of elements.']);
-				assert(numel(myGridName) == obj(1).nargin, ...
-					'All (or none) gridName must be specified');
-				gridPermuteIdx = cellfun(@(v) obj(1).gridIndex(v), myGridName);
-				myGrid = myGrid(gridPermuteIdx);
-			end
-			% at first the value has to be initialized by different
-			% dimensions, to simplify the evaluation of each entry:
-
-			gridSize = cellfun('length', myGrid);
-			if isempty(gridSize), gridSize = 1; end
+		function value = obj2value(obj, myDomain, recalculate)
 
-			value = nan([numel(obj), gridSize]);
-			for k = 1:numel(obj)
-				ndGrd = obj.ndgrid( myGrid );
-				tmp = obj(k).evaluateFunction( ndGrd{:} );
-				% Replaced
-				%   double(subs(obj(k).valueSymbolic, obj.variable, grid));
-				% because this is very slow!
-				value(k,:) = tmp(:);
+			% check if the domain for the evaluation has changed. If not
+			% we can use the stored values in valueDiscrete:
+			if nargin < 3
+				recalculate = false;
+			end
+			if ~recalculate && isequal(myDomain, obj(1).domain)
+				value = obj2value@quantity.Discrete(obj);
+			else
+				% otherwise the function has to be evaluated on the new
+				% domain
+				value = cell(size(obj));
+				for k = 1:numel(obj)
+					ndGrd = myDomain.ndgrid;
+					tmp = obj(k).evaluateFunction( ndGrd{:} );
+					value{k} = tmp(:);
+				end					
+				value = reshape( cell2mat(value), [ gridLength(myDomain), size(obj)]);
 			end
-			value = reshape(permute(value, [1+(gridPermuteIdx), 1]), ...
-				[gridSize(gridPermuteIdx), size(obj)]);
 		end
 		
 		function n = nargin(obj)
@@ -134,8 +103,6 @@ classdef Function < quantity.Discrete
 		end
 	end
 	
-	
-	
 	% -------------
 	% --- Mathe ---
 	%--------------

+quantity/Operator.m

@@ -19,28 +19,63 @@ classdef Operator < handle & matlab.mixin.Copyable
 	end
 	methods
 		function obj = Operator(A, varargin)
-			
+			% OPERATOR initialization of an operator object
+			%	 obj = Operator(A, varargin) initializes an operator of the
+			%	 form 
+			%		A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + ...
+			%			A2(z) dz^2 x(z,t) + Ak dz^k x(z,t)
+			%	The input parameter A can be:
+			%		cell-array of doubles 
+			%			this is good to initialize operators with constant 
+			%			coefficients
+			%			A[x](t) = a_0 x(t) + a_1 dt x(t) + ...
+			%
+			%		cell-array of quantities
+			%			this helps to initialize operators with variable
+			%			coefficients:
+			%		A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + ...
+			%			A2(z) dz^2 x(z,t) + Ak dz^k x(z,t)
+			%
+			%
+			%	
 			if nargin > 0
 				
 				prsr = misc.Parser();
 				prsr.addOptional('s', sym('s'));
 				prsr.addOptional('name', '');
+				prsr.addOptional('domain', quantity.Domain.empty());
 				prsr.parse(varargin{:});
 				
 				s = prsr.Results.s;
 				
 				%TODO assert that A has entries with the same size
+				
+				if isa(A, 'sym') && symvar(A) == prsr.Results.s
+					% if the matrix A is a symbolic expression which
+					% contains the algebraic representation of the
+					% operator, then the object can be generated if the
+					% symbolic expression is converted into the right form
+
+					% TODO: this works only for constant coefficients:
+					% catch the case if the coefficients are variable:
+					assert(length(symvar(A)) == 1, 'Not yet implemented')
+					
+					A = num2cell( double( polyMatrix.polynomial2coefficients(A, prsr.Results.s) ), ...
+						[1, 2]);
+					
+				end
+				
 				if iscell(A)
 					for k = 1 : numel(A)
 						if isnumeric(A{k})
 							obj(k).coefficient = ...
-