diff --git a/+misc/Parser.m b/+misc/Parser.m
index a5545a45e98306278e4c43f6e1f538e5c4bcc15d..4df1f6a5ee236e7703eec28cfcb93e6fb38fd337 100644
--- a/+misc/Parser.m
+++ b/+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)
diff --git a/+misc/ensureIsCell.m b/+misc/ensureIsCell.m
new file mode 100644
index 0000000000000000000000000000000000000000..b93e156a35a6bc0ec047105af10ddd3ab10ea612
--- /dev/null
+++ b/+misc/ensureIsCell.m
@@ -0,0 +1,11 @@
+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
+
diff --git a/+quantity/BasicVariable.m b/+quantity/BasicVariable.m
index 53f5bfaf70206700c83024c7582ae6cd8cd0dc0a..b46a9dbdbaf7dda4f7dbd1b9fd624dea08afb2ac 100644
--- a/+quantity/BasicVariable.m
+++ b/+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
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index d30e3e6489dcbbabf66d4ef52101d6392657c292..aa02851964d404afc2f720d52537343f418684bc 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -2,9 +2,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
properties (SetAccess = protected)
% Discrete evaluation of the continuous quantity
valueDiscrete double;
-
- % Name of the domains that generate the grid.
- gridName {mustBe.unique};
end
properties (Hidden, Access = protected, Dependent)
@@ -12,6 +9,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
properties
+ % ID of the figure handle in which the handle is plotted
+ figureID double = 1;
+
+ % Name of this object
+ name char;
+
+ % 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>]}
@@ -19,12 +31,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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
-
- % ID of the figure handle in which the handle is plotted
- figureID double = 1;
-
- % Name of this object
- name char;
end
methods
@@ -32,13 +38,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% --- Constructor ---
%--------------------
function obj = Discrete(valueOriginal, varargin)
- % The constructor requires valueOriginal to be
+ % 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)]
+ % 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.
@@ -56,7 +70,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
obj = valueOriginal;
else
% empty object. this is needed for instance, to create
- % quantity.Discrete([]), which is useful for creating default
+ % quantity.Discrete([]), which is useful for creating default
% values.
obj = quantity.Discrete.empty(size(valueOriginal));
end
@@ -64,106 +78,83 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
%% input parser
myParser = misc.Parser();
- myParser.addParameter('gridName', [], @(g) ischar(g) || iscell(g));
- myParser.addParameter('grid', [], @(g) isnumeric(g) || iscell(g));
myParser.addParameter('name', string(), @isstr);
myParser.addParameter('figureID', 1, @isnumeric);
myParser.parse(varargin{:});
- assert(all(~contains(myParser.UsingDefaults, 'gridName')), ...
- 'gridName is a mandatory input for quantity');
- if iscell(myParser.Results.gridName)
- myGridName = myParser.Results.gridName;
- else
- myGridName = {myParser.Results.gridName};
+ %% 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
- %% allow initialization of empty objects:
+ % 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+numel(myGridName)):end));
+ myParser.addParameter('size', valueOriginalSize((1+ndims(myDomain)):end));
myParser.parse(varargin{:});
obj = quantity.Discrete.empty(myParser.Results.size);
return;
end
- %% get the sizes of obj and grid
- if iscell(valueOriginal)
- if isempty(valueOriginal{1})
- % if valueOriginal is a cell-array with empty
- % cells, then grid must be specified as an input
- % parameter. This case is important for
- % constructing Symbolic or Function quantities
- % without discrete values.
- assert(all(~contains(myParser.UsingDefaults, 'grid')), ...
- ['grid is a mandatory input for quantity, ', ...
- 'if no discrete values are specified']);
- if ~iscell(myParser.Results.grid)
- gridSize = numel(myParser.Results.grid);
- else
- gridSize = cellfun(@(v) numel(v), myParser.Results.grid);
- end
- else
- gridSize = size(valueOriginal{1});
- end
- objSize = size(valueOriginal);
- elseif isnumeric(valueOriginal)
- gridSize = valueOriginalSize(1 : numel(myGridName));
- objSize = [valueOriginalSize(numel(myGridName)+1 : end), 1, 1];
- end
-
- %% get grid and check size
- if any(contains(myParser.UsingDefaults, 'grid'))
- myGrid = quantity.Discrete.defaultGrid(gridSize);
- else
- myGrid = myParser.Results.grid;
- end
- if ~iscell(myGrid)
- myGrid = {myGrid};
- end
- if isempty(myGridName) || isempty(myGrid)
- if ~(isempty(myGridName) && isempty(myGrid))
- error(['If one of grid and gridName is empty, ', ...
- 'then both must be empty.']);
- end
- else
- assert(isequal(size(myGrid), size(myGridName)), ...
- 'number of grids and gridNames must be equal');
- myGridSize = cellfun(@(v) numel(v), myGrid);
- assert(isequal(gridSize(gridSize>1), myGridSize(myGridSize>1)), ...
- 'grids do not fit to valueOriginal');
- end
-
- %% set valueDiscrete
- if ~iscell(valueOriginal)
- valueOriginal = quantity.Discrete.value2cell(valueOriginal, gridSize, objSize);
- end
- for k = 1:prod(objSize)
- if numel(myGrid) == 1
- obj(k).valueDiscrete = valueOriginal{k}(:);
+ % set valueDiscrete
+ for k = 1:numel(valueOriginal)
+ if numel(myDomain) == 1
+ % TODO: Which case is this? Why does it need extra
+ % treatment?
+ obj(k).valueDiscrete = valueOriginal{k}(:); %#ok<AGROW>
else
- obj(k).valueDiscrete = valueOriginal{k};
+ obj(k).valueDiscrete = valueOriginal{k}; %#ok<AGROW>
end
end
%% set further properties
- [obj.grid] = deal(myGrid);
- [obj.gridName] = deal(myGridName);
+ [obj.domain] = deal(myDomain);
[obj.name] = deal(myParser.Results.name);
[obj.figureID] = deal(myParser.Results.figureID);
%% reshape object from vector to matrix
- obj = reshape(obj, objSize);
+ 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 one point.
@@ -172,29 +163,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function doNotCopy = get.doNotCopy(obj)
doNotCopy = obj.doNotCopyPropertiesName();
end
- function set.gridName(obj, name)
- if ~iscell(name)
- name = {name};
- end
- obj.gridName = name;
- end
-
- function set.grid(obj, grid)
- if ~iscell(grid)
- grid = {grid};
- end
-
- isV = cellfun(@(v) (sum(size(v)>1) == 1) || (numel(v) == 1), grid); % also allow 1x1x1x41 grids
- assert(all(isV(:)), 'Please use vectors for the grid entries!');
-
- [obj.grid] = deal(grid);
- 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.on(obj.grid);
+ obj.valueDiscrete = obj.obj2value(obj.domain, true);
end
valueDiscrete = obj.valueDiscrete;
end
@@ -225,7 +199,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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{:}), ...
@@ -249,7 +223,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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.
[obj.name] = deal(newName);
@@ -257,46 +231,188 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
methods (Access = public)
- function value = on(obj, myGrid, myGridName)
- % TODO es sieht so aus als w�rde die Interpolation bei
- % konstanten werten ziemlichen Quatsch machen!
- % Da muss man nochmal ordentlich drauf schauen!
+ function [d, I, d_size] = compositionDomain(obj, domainName)
+
+ assert(isscalar(obj));
+
+ d = obj.on();
+
+ % the evaluation of obj.on( compositionDomain ) is done by:
+ d_size = size(d);
+
+ % 1) vectorization of the n-d-grid: compositionDomain
+ d = d(:);
+
+ % 2) then it is sorted in ascending order
+ [d, I] = sort(d);
+
+ % verify the domain to be monotonical increasing
+ deltaCOD = diff(d);
+ assert(misc.alln(deltaCOD >= 0), 'The domain for the composition f(g(.)) must be monotonically increasing');
+
+ 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 ) )
+
+ % 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] = ...
+ 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
+
+ % 3) 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)
+ 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 fo the new domain
+ newValues(:,I) = newValues(:,:);
+
+ % 5) 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 ).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', 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
- myGrid = obj(1).grid;
- myGridName = obj(1).gridName;
- elseif nargin >= 2 && ~iscell(myGrid)
- myGrid = {myGrid};
- end
- gridPermuteIdx = 1:obj(1).nargin;
- if nargin == 3
- if ~iscell(myGridName)
- myGridName = {myGridName};
+ if nargin == 2
+ % case 1: a domain is specified by myDomain or by
+ % myDomain as a cell-array with grid entries
+ if iscell(myDomain) || isnumeric(myDomain)
+ myDomain = misc.ensureIsCell(myDomain);
+% assert(all(cellfun(@(v)isvector(v), myDomain)), 'The cell entries for a new grid have to be vectors')
+ newGrid = myDomain;
+ myDomain = quantity.Domain.empty();
+
+ if obj(1).isConstant()
+ gridNames = repmat({''}, length(newGrid));
+ else
+ gridNames = {obj(1).domain.name};
+ end
+
+ for k = 1:length(newGrid)
+ myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
+ end
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
-
- value = obj.obj2value();
-
- if nargin >= 2 && (prod(obj(1).gridSize) > 1)
- indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
- tempInterpolant = numeric.interpolant(...
- [obj(1).grid, indexGrid{:}], value);
- value = tempInterpolant.evaluate(myGrid{:}, indexGrid{:});
- elseif obj.isConstant
- value = repmat(value, [cellfun(@(v) numel(v), myGrid), ones(1, length(size(obj)))]);
- gridPermuteIdx = 1:numel(myGrid);
+ elseif nargin == 3
+ % case 2: 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.
+ myDomain = misc.ensureIsCell(myDomain);
+ gridNames = misc.ensureIsCell(gridNames);
+
+ assert(all(cellfun(@(v)isvector(v), myDomain)), 'The cell entries for a new grid have to be vectors')
+ assert(iscell(gridNames), 'The gridNames parameter must be cell array')
+ assert(all(cellfun(@ischar, gridNames)), 'The gridNames must be strings')
+
+ newGrid = myDomain;
+ myDomain = quantity.Domain.empty();
+ for k = 1:length(newGrid)
+ myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
+ end
+ else
+ myDomain = obj(1).domain;
end
- value = permute(reshape(value, [cellfun(@(v) numel(v), myGrid), size(obj)]), ...
+
+ % 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 permuatation 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
@@ -307,13 +423,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
value = zeros(size(obj));
indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
interpolant = numeric.interpolant(...
- [indexGrid{:}], value);
+ [indexGrid{:}], value);
else
- myGrid = obj(1).grid;
+ myGrid = obj(1).grid;
value = obj.obj2value();
indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
interpolant = numeric.interpolant(...
- [myGrid, indexGrid{:}], value);
+ [myGrid, indexGrid{:}], value);
end
end
@@ -363,17 +479,17 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
[referenceGrid, referenceGridName] = varargin{1}.getFinestGrid(varargin{2:end});
else
referenceGrid = {};
- referenceGridName = '';
+ referenceGridName = '';
end
return;
else
referenceGridName = a(1).gridName;
referenceGrid = a(1).grid;
- referenceGridSize = a(1).gridSize(referenceGridName);
+ referenceGridSize = a(1).gridSize(referenceGridName);
end
-
+
for it = 1 : numel(varargin)
- if isempty(varargin{it})
+ if isempty(varargin{it}) || isempty(varargin{it}(1).domain)
continue;
end
assert(numel(referenceGridName) == numel(varargin{it}(1).gridName), ...
@@ -392,75 +508,24 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
- function [gridJoined, gridNameJoined] = gridJoin(obj1, obj2)
- %% 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.
-
- gridNameJoined = unique([obj1(1).gridName, obj2(1).gridName]);
- gridJoined = cell(1, numel(gridNameJoined));
- for it = 1 : numel(gridNameJoined)
- currentGridName = gridNameJoined{it};
- [index1, lolo1] = obj1.gridIndex(currentGridName);
- [index2, lolo2] = obj2.gridIndex(currentGridName);
- if ~any(lolo1)
- gridJoined{it} = obj2(1).grid{index2};
- elseif ~any(lolo2)
- gridJoined{it} = obj1(1).grid{index1};
- else
- tempGrid1 = obj1(1).grid{index1};
- tempGrid2 = obj2(1).grid{index2};
-
- if ~obj1.isConstant && ~obj2.isConstant
- assert(tempGrid1(1) == tempGrid2(1), 'Grids must have same domain for gridJoin')
- assert(tempGrid1(end) == tempGrid2(end), 'Grids must have same domain for gridJoin')
- end
- if numel(tempGrid1) > numel(tempGrid2)
- gridJoined{it} = tempGrid1;
- else
- gridJoined{it} = tempGrid2;
- end
- end
- end
- end % gridJoin()
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 = 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;
- else
- descend = 0;
- end
-
+
% only sort the grids if there is something to sort
if obj(1).nargin > 1
- gridNames = obj(1).gridName;
- % 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
- [obj.grid] = deal(obj(1).grid(I));
-
- % assign the new grid names
- [obj.gridName] = deal(sortedNames);
-
- % permute the value discrete
+ [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
end% sort()
+
function c = horzcat(a, varargin)
%HORZCAT Horizontal concatenation.
% [A B] is the horizontal concatenation of objects A and B
@@ -488,7 +553,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% X3 ...]' when any of X1, X2, X3, etc. is an object.
%
% See also HORZCAT, CAT.
- c = cat(2, a, varargin{:});
+ c = cat(2, a, varargin{:});
end
function c = vertcat(a, varargin)
%VERTCAT Vertical concatenation.
@@ -535,18 +600,18 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% single matrix.
%
% Examples:
- % a = magic(3); b = pascal(3);
+ % 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),
+ % 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};
+ objCell = {a};
else
objCell = [{a}, varargin(:)'];
@@ -565,7 +630,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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);
+ S = sum(cat(3, s{:}), 3);
elseif dim == 2
S = s{1};
else
@@ -589,7 +654,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if isempty(value)
M = zeros([prod(obj(1).gridSize), size(value(l))]);
end
- M = reshape(M, [obj(1).gridSize, size(value)]);
+ M = reshape(M, [obj(1).gridSize, size(value)]);
o = quantity.Discrete( M, ...
'size', size(value), ...
'gridName', obj(1).gridName, ...
@@ -683,13 +748,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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));
+
+ % 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
function inverse = invert(obj, gridName)
@@ -709,7 +774,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
assert(isequal(size(obj), [1, 1]));
inverse = quantity.Discrete(repmat(obj(1).grid{obj.gridIndex(gridName)}(:), [1, size(obj)]), ...
'size', size(obj), 'grid', obj.on(), 'gridName', {[obj(1).name]}, ...
- 'name', gridName);
+ 'name', gridName);
end
@@ -731,15 +796,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
myParser.addParameter('AbsTol', 1e-6);
myParser.parse(varargin{:});
- variableGrid = myParser.Results.variableGrid;
- myGridSize = [numel(variableGrid), ...
+ 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))];
+ positiveVariableGrid = [0; variableGrid(variableGrid > 0)];
+ negativeVariableGrid = [0; flip(variableGrid(variableGrid < 0))];
% solve ode for every entry in obj and for every initial value
options = odeset('RelTol', myParser.Results.RelTol, 'AbsTol', myParser.Results.AbsTol);
@@ -753,14 +818,14 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if numel(positiveVariableGrid) > 1
[resultGridPositive, odeSolutionPositive] = ...
ode45(@(y, z) obj(it).on(z), ...
- positiveVariableGrid, ...
- myParser.Results.initialValueGrid(icIdx), options);
+ positiveVariableGrid, ...
+ myParser.Results.initialValueGrid(icIdx), options);
end
if numel(negativeVariableGrid) >1
[resultGridNegative, odeSolutionNegative] = ...
ode45(@(y, z) obj(it).on(z), ...
- negativeVariableGrid, ...
- myParser.Results.initialValueGrid(icIdx), options);
+ negativeVariableGrid, ...
+ myParser.Results.initialValueGrid(icIdx), options);
end
if any(variableGrid == 0)
resultGrid = [flip(resultGridNegative(2:end)); 0 ; resultGridPositive(2:end)];
@@ -782,8 +847,29 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
{variableGrid, myParser.Results.initialValueGrid}, ...
'size', size(obj), 'name', ['solve(', obj(1).name, ')']);
end
-
+
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;
@@ -792,19 +878,23 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
solution = obj;
else
% input checks
- assert(nargin == 3, ['Wrong number of input arguments. ', ...
- 'gridName2Replace and values must be cell-arrays!']);
- if ~iscell(gridName2Replace)
- gridName2Replace = {gridName2Replace};
- end
- if ~iscell(values)
- values = {values};
+ 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);
+ 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
+ % 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}).
@@ -835,8 +925,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% called with subs(f, zeta, z), then g(z) = f(z, z)
% results, hence the dimensions z and zeta are
% merged.
- gridIndices = [obj(1).gridIndex(gridName2Replace{1}), ...
- obj(1).gridIndex(values{1})];
+ gridIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
+ obj(1).domain.gridIndex(values{1})];
newGridForOn = obj(1).grid;
if numel(obj(1).grid{gridIndices(1)}) > numel(obj(1).grid{gridIndices(2)})
newGridForOn{gridIndices(2)} = newGridForOn{gridIndices(1)};
@@ -846,17 +936,16 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
newValue = misc.diagNd(obj.on(newGridForOn), gridIndices);
newGrid = {newGridForOn{gridIndices(1)}, ...
newGridForOn{1:1:numel(newGridForOn) ~= gridIndices(1) ...
- & 1:1:numel(newGridForOn) ~= gridIndices(2)}};
+ & 1:1:numel(newGridForOn) ~= gridIndices(2)}};
newGridName = {values{1}, ...
obj(1).gridName{1:1:numel(obj(1).gridName) ~= gridIndices(1) ...
- & 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
-
+ & 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
else
% this is the default case. just grid name is
% changed.
newGrid = obj(1).grid;
newGridName = obj(1).gridName;
- newGridName{obj(1).gridIndex(gridName2Replace{1})} ...
+ newGridName{obj(1).domain.gridIndex(gridName2Replace{1})} ...
= values{1};
newValue = obj.on();
end
@@ -870,7 +959,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% newGrid is the similar to the original grid, but the
% grid of gridName2Replace is removed.
newGrid = obj(1).grid;
- newGrid = newGrid((1:1:numel(newGrid)) ~= obj.gridIndex(gridName2Replace{1}));
+ newGrid = newGrid((1:1:numel(newGrid)) ~= obj(1).domain.gridIndex(gridName2Replace{1}));
newGridSize = cellfun(@(v) numel(v), newGrid);
% newGridForOn is the similar to the original grid, but
% the grid of gridName2Replace is set to values{1} for
@@ -900,11 +989,14 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
solution = solution.subs(gridName2Replace(2:end), values(2:end));
end
end
-
end
+ function [idx, logicalIdx] = gridIndex(obj, varargin)
+ [idx, logicalIdx] = obj(1).domain.gridIndex(varargin{:});
+ end
+
function value = at(obj, point)
- % at() evaluates the object at one point and returns it as array
+ % at() evaluates the object at one point and returns it as array
% with the same size as size(obj).
value = reshape(obj.on(point), size(obj));
end
@@ -913,11 +1005,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% ATINDEX returns the valueDiscrete at the requested index.
% value = atIndex(obj, varargin) returns the
% quantity.Discrete.valueDiscrete at the index defined by
- % varargin.
+ % varargin.
% value = atIndex(obj, 1) returns the first element of
% "valueDiscrete"
% value = atIndex(obj, ':') returns all elements of
- % obj.valueDiscrete in vectorized form.
+ % obj.valueDiscrete in vectorized form.
% value = atIndex(obj, 1, end) returns the obj.valueDiscrete
% at the index (1, end).
% If a range of index is requested, the result is returned
@@ -941,15 +1033,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
'UniformOutput', false);
valueSize = size(value{1});
-
+
if all(cellfun(@numel, varargin) == 1) && all(cellfun(@isnumeric, varargin))
outputSize = [];
else
outputSize = valueSize(1:obj(1).nargin);
- end
-
+ end
+
value = reshape([value{:}], [outputSize, size(obj)]);
- end
+ end
end
function n = nargin(obj)
@@ -976,7 +1068,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function s = gridSize(obj, myGridName)
% GRIDSIZE returns the size of all grid entries.
% todo: this should be called gridLength
- if isempty(obj(1).grid)
+ if isempty(obj(1).domain)
s = [];
else
if nargin == 1
@@ -987,20 +1079,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
- function matGrid = ndgrid(obj, grid)
- % ndgrid calles ndgrid for the default grid, if no other grid
- % is specified. Empty grid as input returns empty cell as
- % result.
- if nargin == 1
- grid = obj.grid;
- end
- if isempty(grid)
- matGrid = {};
- else
- [matGrid{1:obj.nargin}] = ndgrid(grid{:});
- end
- end % ndgrid()
-
function H = plot(obj, varargin)
H = [];
p = misc.Parser();
@@ -1079,7 +1157,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
xlabel(labelHelper(1), 'Interpreter','latex');
if p.Results.showTitle
- title(titleHelper(), 'Interpreter','latex');
+ title(titleHelper(), 'Interpreter','latex');
end
a = gca();
a.TickLabelInterpreter = 'latex';
@@ -1092,10 +1170,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if p.Results.export
matlab2tikz(p.Results.exportOptions{:});
end
-
+
function myLabel = labelHelper(gridNumber)
if ~isempty(obj(rowIdx, columnIdx, figureIdx).gridName)
- myLabel = ['$$', greek2tex(obj(rowIdx, columnIdx, figureIdx).gridName{gridNumber}), '$$'];
+ myLabel = ['$$', greek2tex(obj(rowIdx, columnIdx, figureIdx).gridName{gridNumber}), '$$'];
else
myLabel = '';
end
@@ -1149,9 +1227,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
si = num2cell( size(obj) );
s(si{:}) = struct();
for l = 1:numel(obj)
-
+
doNotCopyProperties = obj(l).doNotCopy;
-
+
for k = 1:length(properties)
if ~any(strcmp(doNotCopyProperties, properties{k}))
s(l).(properties{k}) = obj(1).(properties{k});
@@ -1159,7 +1237,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
else
- % Without this else-case
+ % Without this else-case
% this method is in conflict with the default matlab built-in
% calls struct(). When calling struct('myFieldName', myQuantity) this
% quantity-method is called (and errors) and not the
@@ -1192,28 +1270,57 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end % flipGrid()
function newObj = changeGrid(obj, gridNew, gridNameNew)
+ % CHANGEGRID change the grid of the quantity.
+ % newObj = CHANGEGRID(obj, gridNew, gridNameNew)
% change the grid of the obj quantity. The order of grid and
% gridName in the obj properties remains unchanged, only the
% data points are exchanged.
+ %
+ % newObj = CHANGEGRID(obj, domain) changes the domain of the
+ % object specified by the name of DOMAIN into the
+ % corresponding domain from DOMAIN.
+ %
+ % example:
+ % q.changeGrid( linspace(0,1)', 't')
+ % will change the grid with the name 't' to the new grid linspace(0,1)'
if isempty(obj)
newObj = obj.copy();
return;
end
- gridIndexNew = obj(1).gridIndex(gridNameNew);
- myGrid = cell(1, numel(obj(1).grid));
- myGridName = cell(1, numel(obj(1).grid));
- for it = 1 : numel(myGrid)
- myGrid{gridIndexNew(it)} = gridNew{it};
- myGridName{gridIndexNew(it)} = gridNameNew{it};
+
+ if isa(gridNew, 'quantity.Domain')
+ gridNameNew = {gridNew.name};
+ gridNew = {gridNew.grid};
+ else
+ gridNew = misc.ensureIsCell(gridNew);
+ gridNameNew = misc.ensureIsCell(gridNameNew);
+ end
+
+ if obj(1).isConstant
+
+ newDomain(1:length( gridNew )) = quantity.Domain();
+ for i = 1 : length(gridNew)
+ newDomain(i) = ...
+ quantity.Domain(gridNameNew{i}, gridNew{i});
+ end
+ else
+ gridIndexNew = obj(1).domain.gridIndex(gridNameNew);
+ newDomain = obj(1).domain;
+
+ for i = 1 : length(gridIndexNew)
+ newDomain(gridIndexNew(i)) = ...
+ quantity.Domain(gridNameNew{i}, gridNew{i});
+ end
+ assert(isequal({newDomain.name}, obj(1).gridName), ...
+ 'rearranging grids failed');
end
- assert(isequal(myGridName(:), obj(1).gridName(:)), 'rearranging grids failed');
+
newObj = obj.copy();
- for it = 1 : numel(obj)
- newObj(it).valueDiscrete = obj(it).on(myGrid);
+ [newObj.domain] = deal(newDomain);
+ for i = 1 : numel(obj)
+ newObj(i).valueDiscrete = obj(i).on(newDomain);
end
- %[newObj.derivatives] = deal({});
- [newObj.grid] = deal(myGrid);
end
function [lowerBounds, upperBounds] = getBoundaryValues(obj, varargin)
@@ -1247,7 +1354,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% state-transition matrix for a ordinary differential equation
% d / dt x(t) = A(t) x(t)
% with A as this object. The system matrix can be constant or
- % variable. It is the solution of
+ % variable. It is the solution of
% d / dt F(t,t0) = A(t) F(t,t0)
% F(t,t) = I.
% Optional parameters for the computation can be defined as
@@ -1267,7 +1374,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
ip = misc.Parser();
ip.addParameter('grid', []);
ip.addParameter('gridName1', '');
- ip.addParameter('gridName2', '');
+ ip.addParameter('gridName2', '');
ip.parse(varargin{:});
myGrid = ip.Results.grid;
gridName1 = ip.Results.gridName1;
@@ -1282,38 +1389,39 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
assert(numel(obj(1).gridName) == 1, 'No gridName is defined. Use name-value-pairs property "gridName1" to define a grid name');
gridName1 = obj(1).gridName{1};
end
-
+
if ip.isDefault('gridName2')
gridName2 = [gridName1 '0'];
end
assert(numel(myGrid) > 1, 'If the state transition matrix is computed for constant values, a spatial domain has to be defined!')
-
+
+ myDomain = [quantity.Domain(gridName1, myGrid), ...
+ quantity.Domain(gridName2, myGrid)];
+
if obj.isConstant
% for a constant system matrix, the matrix exponential
% function can be used.
z = sym(gridName1, 'real');
zeta = sym(gridName2, 'real');
f0 = expm(obj.atIndex(1)*(z - zeta));
- F = quantity.Symbolic(f0, 'grid', {myGrid, myGrid});
-
+ F = quantity.Symbolic(f0, 'domain', myDomain);
+
elseif isa(obj, 'quantity.Symbolic')
f0 = misc.fundamentalMatrix.odeSolver_par(obj.function_handle, myGrid);
F = quantity.Discrete(f0, ...
'size', [size(obj, 1), size(obj, 2)], ...
- 'gridName', {gridName1, gridName2}, ...
- 'grid', {myGrid, myGrid});
+ 'domain', myDomain);
else
f0 = misc.fundamentalMatrix.odeSolver_par( ...
obj.on(myGrid), ...
myGrid );
F = quantity.Discrete(f0, ...
'size', [size(obj, 1), size(obj, 2)], ...
- 'gridName', {gridName1, gridName2}, ...
- 'grid', {myGrid, myGrid});
+ 'domain', myDomain);
end
- end
+ end
function s = sum(obj, dim)
s = quantity.Discrete(sum(obj.on(), obj.nargin + dim), ...
@@ -1367,7 +1475,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function P = mtimes(a, b)
% TODO rewrite the selection of the special cases! the
% if-then-cosntruct is pretty ugly!
-
+
% numeric computation of the matrix product
if isempty(b) || isempty(a)
% TODO: actually this only holds for multiplication of
@@ -1388,7 +1496,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
return
else
b = quantity.Discrete(b, 'size', size(b), 'grid', {}, ...
- 'gridName', {}, 'name', '');
+ 'gridName', {}, 'name', '');
end
end
@@ -1396,8 +1504,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
P = (b.' * a.').';
% this recursion is safe, because isa(b, 'double') is considered in
% the if above.
- P.setName(['c ', b(1).name]);
- return
+
+ if isnumeric(P)
+ return
+ else
+ P.setName(['c ', b(1).name]);
+ return
+ end
end
if a.isConstant() && ~b.isConstant()
% If the first argument a is constant value, then bad
@@ -1408,11 +1521,14 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
P = (b' * a')';
P.setName([a(1).name, ' ', b(1).name]);
return
+ elseif a.isConstant() && b.isConstant()
+ P = a.on() * b.on();
+ return
end
if isscalar(b)
% do the scalar multiplication in a loop
for k = 1:numel(a)
- P(k) = innerMTimes(a(k), b);
+ P(k) = innerMTimes(a(k), b);
end
P = reshape(P, size(a));
return
@@ -1421,13 +1537,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if isscalar(a)
% do the scalar multiplication in a loop
for k = 1:numel(b)
- P(k) = innerMTimes(a, b(k));
+ P(k) = innerMTimes(a, b(k)); %#ok<AGROW>
end
P = reshape(P, size(b));
return
end
-
- P = innerMTimes(a, b);
+
+ P = innerMTimes(a, b);
end
function P = innerMTimes(a, b)
@@ -1438,40 +1554,52 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% misc.multArray is very efficient, but requires that the
% multiple dimensions of the input are correctly arranged.
[idx, permuteGrid] = computePermutationVectors(a, b);
-
- parameters = struct(a(1));
+
parameters.name = [a(1).name, ' ', b(1).name];
- parameters.grid = [...
- a(1).grid(idx.A.common), ...
- b(1).grid(~idx.B.common)];
- parameters.gridName = [...
- a(1).gridName(idx.A.grid), ...
- b(1).gridName(~idx.B.common)];
-
- % select finest grid from a and b
- parameters.grid = cell(1, numel(parameters.gridName));
- [gridJoined, gridNameJoined] = gridJoin(a, b);
- % gridJoin combines the grid of a and b while using the finer
- % of both.
- aGrid = cell(1, a.nargin); % for later call of a.on() with fine grid
- bGrid = cell(1, b.nargin); % for later call of b.on() with fine grid
- for it = 1 : numel(parameters.gridName)
- parameters.grid{it} = gridJoined{...
- strcmp(parameters.gridName{it}, gridNameJoined)};
- if a.gridIndex(parameters.gridName{it})
- aGrid{a.gridIndex(parameters.gridName{it})} = ...
- parameters.grid{it};
+ parameters.figureID = a(1).figureID;
+
+ domainA = a(1).domain;
+ domainB = b(1).domain;
+
+ % The joined domain of quantity A and B is a mixture from both
+ % domains. If the domains have the same name, it must be
+ % checked if they have the same discretization. This is done by
+ % the function "gridJoin". It returns the finer grid.
+ % If the domain names do not coincide, they are just
+ % appended. At first, the domains of qunatity A, then the
+ % domains of quantity B.
+ joinedDomain = [ gridJoin(domainA(idx.A.common), domainB(idx.B.common)), ...
+ domainA(~idx.A.common), domainB(~idx.B.common) ];
+
+ %% generate the new grids
+ newDomainA = quantity.Domain.empty(); % for later call of a.on() with fine grid
+ newDomainB = quantity.Domain.empty(); % for later call of b.on() with fine grid
+ for i = 1 : numel(joinedDomain)
+
+ parameters.domain(i) = joinedDomain(i);
+
+ % if there is a component of the original domain in the
+ % joined domain, it can happen the length of a domain has
+ % changed. Thus, it must be updated for a later evaluation.
+ idxA = domainA.gridIndex(joinedDomain(i).name);
+ idxB = domainB.gridIndex(joinedDomain(i).name);
+
+ if idxA
+ newDomainA = [newDomainA, joinedDomain(i)]; %#ok<AGROW>
+ end
+
+ if idxB
+ newDomainB = [newDomainB, joinedDomain(i)]; %#ok<AGROW>
end
- if b.gridIndex(parameters.gridName{it})
- bGrid{b.gridIndex(parameters.gridName{it})} = ...
- parameters.grid{it};
- end
end
parameters = misc.struct2namevaluepair(parameters);
- valueA = permute(a.on(aGrid), idx.A.permute);
- valueB = permute(b.on(bGrid), idx.B.permute);
+ % evaluation of the quantities on their "maybe" new domain.
+ valueA = a.on(newDomainA);
+ valueB = b.on(newDomainB);
+ % do the multidimensional tensor multiplication and permute the
+ % values to the right order
C = misc.multArray(valueA, valueB, idx.A.value(end), idx.B.value(1), idx.common);
C = permute(C, permuteGrid);
P = quantity.Discrete(C, parameters{:});
@@ -1494,11 +1622,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
objDiscreteReshaped = permute(reshape(objDiscreteOriginal, ...
[prod(obj(1).gridSize), size(obj)]), [2, 3, 1]);
invDiscrete = zeros([prod(obj(1).gridSize), size(obj)]);
-
+
parfor it = 1 : size(invDiscrete, 1)
invDiscrete(it, :, :) = inv(objDiscreteReshaped(:, :, it));
end
-
+
y = quantity.Discrete(reshape(invDiscrete, size(objDiscreteOriginal)),...
'size', size(obj), ...
'name', ['{(', obj(1).name, ')}^{-1}'], ...
@@ -1515,9 +1643,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
objT = obj.';
objCtMat = conj(objT.on());
objCt = quantity.Discrete(objCtMat, 'grid', obj(1).grid, ...
- 'gridName', obj(1).gridName, 'name', ['{', obj(1).name, '}^{H}']);
+ 'gridName', obj(1).gridName, 'name', ['{', obj(1).name, '}^{H}']);
end % ctranspose(obj)
-
+
function y = exp(obj)
% exp() is the exponential function using obj as the exponent.
y = quantity.Discrete(exp(obj.on()), ...
@@ -1534,7 +1662,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% The integral domain is specified by integralGrid.
if obj.nargin > 1
- assert(ischar(integralGridName), 'integralGrid must specify a gridName as a char');
+ assert(ischar(integralGridName), 'integralGrid must specify a gridName as a char');
else
integralGridName = obj(1).gridName;
end
@@ -1555,7 +1683,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% xNorm = sqrt(x.' * x).
% Optionally, a weight can be defined, such that instead
% xNorm = sqrt(x.' * weight * x).
-
+
myParser = misc.Parser();
myParser.addParameter('weight', eye(size(obj, 1)));
myParser.parse(varargin{:});
@@ -1592,7 +1720,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% number of rows."
x = inv(A) * B;
end
-
+
function x = mrdivide(B, A)
% mRdivide see doc mrdivide of matlab: "x = B/A solves the
% system of linear equations x*A = B for x. The matrices A and
@@ -1636,11 +1764,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% Check if there is any dimension which is zero
empty = any(size(obj) == 0);
- % If the constructor is called without arguments, a
- % quantity.Discrete is created without an initialization of
- % the grid. Thus, the quantity is not initialized and empty if
- % the grid is not a cell.
- empty = empty || ~iscell(obj(1).grid);
+ % If the constructor is called without any arguments, an object
+ % is still created ( this is required to allow object-arrays),
+ % but the object should be recognized as an empty object. Thus,
+ % we can test if the domain has been set with a domain object.
+ % This happens only in the part of the constructor which is
+ % entered if a valid quantity is initialized.
+ empty = empty || ~isa(obj(1).domain, 'quantity.Domain');
end % isempty()
@@ -1682,7 +1812,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function myGrid = gridOf(obj, myGridName)
if ~iscell(myGridName)
- gridIdx = obj(1).gridIndex(myGridName);
+ gridIdx = obj(1).domain.gridIndex(myGridName);
if gridIdx > 0
myGrid = obj(1).grid{gridIdx};
else
@@ -1691,7 +1821,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
else
myGrid = cell(size(myGridName));
for it = 1 : numel(myGrid)
- gridIdx = obj(1).gridIndex(myGridName{it});
+ gridIdx = obj(1).domain.gridIndex(myGridName{it});
if gridIdx > 0
myGrid{it} = obj(1).grid{gridIdx};
else
@@ -1701,37 +1831,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
- function [idx, logicalIdx] = gridIndex(obj, names)
- %% GRIDINDEX returns the index of the grid
- % [idx, log] = gridIndex(obj, names) searches in the gridName
- % properties of obj for the "names" and returns its index as
- % "idx" and its logical index as "log"
-
- if nargin == 1
- names = obj(1).gridName;
- end
-
- if ~iscell(names)
- names = {names};
- end
-
- idx = zeros(1, length(names));
- nArgIdx = 1:obj.nargin();
-
- logicalIdx = false(1, obj.nargin());
-
- for k = 1:length(names)
- log = strcmp(obj(1).gridName, names{k});
- logicalIdx = logicalIdx | log;
- if any(log)
- idx(k) = nArgIdx(log);
- else
- idx(k) = 0;
- end
- end
-
- end
-
function result = diff(obj, k, diffGridName)
% DIFF computation of the derivative
% result = DIFF(obj, k, diffGridName) applies the
@@ -1752,11 +1851,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
return
end
- if nargin < 3 % if no diffGridName is specified, then the derivative
+ if nargin < 3 % if no diffGridName is specified, then the derivative
% w.r.t. to all gridNames is applied
diffGridName = obj(1).gridName;
end
-
+
+ if isa(diffGridName, 'quantity.Domain')
+ % a quantity.Domain is used instead of a grid name
+ % -> get the grid name from the domain object
+ %
+ % #todo@domain: make the default case to call with a
+ % quantity.Domain instead of a grid name. Then, the
+ % section about the grid selection and so can be simplified
+ %
+ diffGridName = {diffGridName.name};
+ end
% diff for each element of diffGridName (this is rather
% inefficient, but an easy implementation of the specification)
@@ -1770,37 +1879,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
else
- gridSelector = strcmp(obj(1).gridName, diffGridName);
- gridSelectionIndex = find(gridSelector);
-
- spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
- assert(numeric.near(spacing, spacing(1)), ...
- 'diff is currently only implemented for equidistant grid');
- permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
-
- objDiscrete = permute(obj.on(), ...
- [permutationVector(gridSelectionIndex), ...
- permutationVector(permutationVector ~= gridSelectionIndex)]);
-
- if size(objDiscrete, 2) == 1
- derivativeDiscrete = gradient(objDiscrete, spacing(1));
- else
- [~, derivativeDiscrete] = gradient(objDiscrete, spacing(1));
- end
-
- rePermutationVector = [2:(gridSelectionIndex), ...
- 1, (gridSelectionIndex+1):ndims(derivativeDiscrete)];
- result = quantity.Discrete(...
- permute(derivativeDiscrete, rePermutationVector), ...
- 'size', size(obj), 'grid', obj(1).grid, ...
- 'gridName', obj(1).gridName, ...
- 'name', ['(d_{', diffGridName, '}', obj(1).name, ')']);
-
- if k > 1
-% % if a higher order derivative is requested, call the function
-% % recursivly until the first-order derivative is reached
- result = result.diff(k-1, diffGridName);
- end
+ result = obj.diff_inner(k, diffGridName);
end
end
@@ -1836,15 +1915,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if iscell(intGridName) && numel(intGridName) > 1
obj = obj.int(intGridName{1}, a(1), b(1));
I = obj.int(intGridName{2:end}, a(2:end), b(2:end));
- return
+ return
end
[idxGrid, isGrid] = obj.gridIndex(intGridName);
if nargin == 4
I = cumInt( obj, varargin{:});
return;
-%
-% assert(all((varargin{2} == a(idxGrid)) & (varargin{3} == b(idxGrid))), ...
-% 'only integration from beginning to end of domain is implemented');
+ %
+ % assert(all((varargin{2} == a(idxGrid)) & (varargin{3} == b(idxGrid))), ...
+ % 'only integration from beginning to end of domain is implemented');
end
%% do the integration over one dimension
@@ -1855,12 +1934,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
[domainIdx{1:obj(1).nargin}] = deal(':');
currentGrid = obj(1).grid{ idxGrid };
domainIdx{idxGrid} = find(currentGrid >= a(idxGrid) & currentGrid <= b(idxGrid));
-
- myGrid = currentGrid(domainIdx{idxGrid});
+ myGrid = currentGrid(domainIdx{idxGrid});
+
for kObj = 1:numel(obj)
J = numeric.trapz_fast_nDim(myGrid, obj(kObj).atIndex(domainIdx{:}), idxGrid);
-
+
% If the quantity is integrated, so that the first
% dimension collapses, this dimension has to be shifted
% away, so that the valueDiscrete can be set correctly.
@@ -1868,21 +1947,20 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
J = shiftdim(J);
end
- if all(isGrid)
- grdName = '';
- newGrid = {};
- else
- grdName = obj(kObj).gridName{~isGrid};
- newGrid = obj(kObj).grid{~isGrid};
- end
-
% create result:
I(kObj).valueDiscrete = J;
- I(kObj).gridName = grdName;
- I(kObj).grid = newGrid;
I(kObj).name = ['int(' obj(kObj).name ')'];
end
-
+
+ if all(isGrid)
+ newDomain = quantity.Domain();
+ else
+ newDomain = quantity.Domain(obj(kObj).gridName{~isGrid}, ...
+ obj(kObj).grid{~isGrid});
+ end
+
+ [I.domain] = deal(newDomain);
+
end
function result = cumInt(obj, domain, lowerBound, upperBound)
@@ -1898,7 +1976,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% input parser since some default values depend on intGridName.
myParser = misc.Parser;
- myParser.addRequired('domain', @(d) obj.gridIndex(d) ~= 0);
+ myParser.addRequired('domain', @(d) obj(1).domain.gridIndex(d) ~= 0);
myParser.addRequired('lowerBound', ...
@(l) isnumeric(l) || ischar(l) );
myParser.addRequired('upperBound', ...
@@ -1907,7 +1985,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% get grid
myGrid = obj(1).grid;
- intGridIdx = obj.gridIndex(domain);
+ intGridIdx = obj(1).domain.gridIndex(domain);
% integrate
F = numeric.cumtrapz_fast_nDim(myGrid{intGridIdx}, ...
@@ -1915,7 +1993,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
result = quantity.Discrete(F, 'grid', myGrid, ...
'gridName', obj(1).gridName);
- % int_lowerBound^upperBound f(.) =
+ % int_lowerBound^upperBound f(.) =
% F(upperBound) - F(lowerBound)
result = result.subs({domain}, upperBound) ...
- result.subs({domain}, lowerBound);
@@ -1923,7 +2001,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
result.setName(deal(['int(', obj(1).name, ')']));
end
end
-
+
function C = plus(A, B)
%% PLUS is the sum of two quantities.
assert(isequal(size(A), size(B)), 'plus() not supports mismatching sizes')
@@ -1936,58 +2014,31 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% for support of numeric inputs:
if ~isa(A, 'quantity.Discrete')
if isnumeric(A)
- A = quantity.Discrete(A, 'name', 'c', 'gridName', {}, 'grid', {});
+ A = quantity.Discrete(A, 'name', 'c', 'domain', quantity.Domain.empty());
else
error('Not yet implemented')
end
elseif ~isa(B, 'quantity.Discrete')
if isnumeric(B)
- B = quantity.Discrete(B, 'name', 'c', 'gridName', {});
+ B = quantity.Discrete(B, 'name', 'c', 'domain', quantity.Domain.empty());
else
B = quantity.Discrete(B, 'name', 'c', 'gridName', A(1).gridName, 'grid', A.grid);
end
end
% combine both domains with finest grid
- [gridJoined, gridNameJoined] = gridJoin(A, B);
- gridJoinedLength = cellfun(@(v) numel(v), gridJoined);
+ joinedDomain = gridJoin(A(1).domain, B(1).domain);
- [aDiscrete] = A.expandValueDiscrete(gridJoined, gridNameJoined, gridJoinedLength);
- [bDiscrete] = B.expandValueDiscrete(gridJoined, gridNameJoined, gridJoinedLength);
-
- % create result object
- C = quantity.Discrete(aDiscrete + bDiscrete, ...
- 'grid', gridJoined, 'gridName', gridNameJoined, ...
- 'size', size(A), 'name', [A(1).name, '+', B(1).name]);
- end
-
- function [valDiscrete] = expandValueDiscrete(obj, gridJoined, gridNameJoined, gridJoinedLength)
- % EXPANDVALUEDISCRETE
- % [valDiscrete, gridNameSorted] = ...
- % expandValueDiscrete(obj, gridIndex, valDiscrete)
- % Expanses the value of obj, so that
- % todo
- % get the index of obj.grid in the joined grid
- gridLogical = logical( obj.gridIndex(gridNameJoined) );
-
- valDiscrete = obj.on(gridJoined(gridLogical > 0), gridNameJoined(gridLogical > 0));
- oldDim = ndims(valDiscrete);
- valDiscrete = permute(valDiscrete, [(1:sum(~gridLogical)) + oldDim, 1:oldDim] );
- valDiscrete = repmat(valDiscrete, [gridJoinedLength(~gridLogical), ones(1, ndims(valDiscrete))]);
-%
- valDiscrete = reshape(valDiscrete, ...
- [gridJoinedLength(~gridLogical), gridJoinedLength(gridLogical), size(obj)]);
-
-% % permute valDiscrete such that grids are in the order specified
-% % by gridNameJoined.
- gridIndex = 1:numel(gridLogical);
- gridOrder = [gridIndex(~gridLogical), gridIndex(gridLogical)];
- valDiscrete = permute(valDiscrete, [gridOrder, numel(gridLogical)+(1:ndims(obj))]);
+ [aDiscrete] = A.expandValueDiscrete(joinedDomain);
+ [bDiscrete] = B.expandValueDiscrete(joinedDomain);
+ % create result object
+ C = quantity.Discrete(aDiscrete + bDiscrete, ...
+ 'domain', joinedDomain, ...
+ 'name', [A(1).name, '+', B(1).name]);
end
-
function C = minus(A, B)
% minus uses plus()
C = A + (-B);
@@ -1997,7 +2048,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
elseif isnumeric(B)
[C.name] = deal([A(1).name, '-c']);
else
- [C.name] = deal([A(1).name, '-', B(1).name]);
+ [C.name] = deal([A(1).name, '-', B(1).name]);
end
end
@@ -2034,7 +2085,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function maxValue = max(obj)
% max returns the maximal value of all elements of A over all
% variables as a double array.
- maxValue = reshape(max(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
+ maxValue = reshape(max(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
end
function [i, m, s] = near(obj, B, varargin)
@@ -2078,7 +2129,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function minValue = min(obj)
% max returns the minimum value of all elements of A over all
% variables as a double array.
- minValue = reshape(min(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
+ minValue = reshape(min(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
end % min()
function absQuantity = abs(obj)
@@ -2086,9 +2137,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if isempty(obj)
absQuantity = quantity.Discrete.empty(size(obj));
else
- absQuantity = quantity.Discrete(abs(obj.on()), ...
- 'grid', obj(1).grid, 'gridName', obj(1).gridName, ...
- 'size', size(obj), 'name', ['|', obj(1).name, '|']);
+ absQuantity = quantity.Discrete(abs(obj.on()), ...
+ 'grid', obj(1).grid, 'gridName', obj(1).gridName, ...
+ 'size', size(obj), 'name', ['|', obj(1).name, '|']);
end
end % abs()
@@ -2111,7 +2162,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function meanValue = mean(obj)
% mean returns the mean value of all elements of A over all
% variables as a double array.
- meanValue = reshape(mean(obj.on(), 1:obj(1).nargin), size(obj));
+ meanValue = reshape(mean(obj.on(), 1:obj(1).nargin), size(obj));
end
function medianValue = median(obj)
@@ -2120,18 +2171,40 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% ATTENTION:
% median over multiple dimensions is only possible from MATLAB 2019a and above!
% TODO Update Matlab!
- medianValue = reshape(median(obj.on(), 1:obj(1).nargin), size(obj));
+ medianValue = reshape(median(obj.on(), 1:obj(1).nargin), size(obj));
end
-
- function q = obj2value(obj)
- % for a case with a 3-dimensional grid, the common cell2mat
- % with reshape did not work. So it has to be done manually:
+
+ function value = obj2value(obj, myDomain)
+ % OBJ2VALUE make the stored data in valueDiscrete available for
+ % in the output format
+ % value = obj2value(obj) returns the valueDiscrete in the
+ % form size(value) = [gridLength, objSize]
+ % obj2value(obj, myDomain) returns the valueDiscrete in the
+ % form size(value) = [myDomain.gridLength objSize]
v = zeros([numel(obj(1).valueDiscrete), numel(obj)]);
for k = 1:numel(obj)
v(:,k) = obj(k).valueDiscrete(:);
end
- q = reshape(v, [obj(1).gridSize(), size(obj)]);
+
+ value = reshape(v, [obj(1).gridSize(), size(obj)]);
+
+ if nargin == 2 && ~isequal( myDomain, obj(1).domain )
+ % if a new domain is specified for the evaluation of
+ % the quantity, ...
+ if obj.isConstant
+ % ... duplicate the constant value on the desired
+ % domain
+ value = repmat(v, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
+ else
+ %... do an interpolation based on the old data.
+ indexGrid = arrayfun(@(s) 1:s, size(obj), 'UniformOutput', false);
+ tempInterpolant = numeric.interpolant(...
+ [obj(1).grid, indexGrid{:}], value);
+ tempGrid = {myDomain.grid};
+ value = tempInterpolant.evaluate(tempGrid{:}, indexGrid{:});
+ end
+ end
end
@@ -2166,48 +2239,43 @@ 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(gridSize, valueSize, grid) returns a
- % quantity.Discrete object that has only zero entries on a
- % grid with 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)
- %
- if ~iscell(grid)
- grid = {grid};
+ % P = zeros(VALUESIZE, DOMAIN) creates a matrix of size
+ % VALUESIZE on the DOMAIN with zero entries.
+
+ 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)
% VALUE2CELL
- gs = num2cell(gridSize(:));
- vs = arrayfun(@(n) ones(n, 1), valueSize, 'UniformOutput', false);
- s = [gs(:); vs(:)]; %{gs{:}, vs{:}};%
- q = reshape(mat2cell(value, s{:}), valueSize);
- end
-
- function g = defaultGrid(gridSize)
- g = cell(1, length(gridSize));
- 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} = O;
+ fullSize = size(value);
+
+ myGridSize = num2cell( fullSize(1:length(gridSize)) );
+
+ if nargin == 2
+ valueSize = [fullSize(length(myGridSize)+1:length(fullSize)), 1, 1];
end
+
+ myValueSize = arrayfun(@(n) ones(n, 1), valueSize, 'UniformOutput', false);
+ s = [myGridSize(:); myValueSize(:)]; %{gs{:}, vs{:}};%
+ q = reshape(mat2cell(value, s{:}), valueSize);
end
-
function [p, q, parameters, gridSize] = inputNormalizer(a, b)
parameters = [];
@@ -2239,10 +2307,75 @@ 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);
+
+ spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
+ assert(numeric.near(spacing, spacing(1)), ...
+ 'diff is currently only implemented for equidistant grid');
+ permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
+
+ objDiscrete = permute(obj.on(), ...
+ [permutationVector(gridSelectionIndex), ...
+ permutationVector(permutationVector ~= gridSelectionIndex)]);
+
+ if size(objDiscrete, 2) == 1
+ derivativeDiscrete = gradient(objDiscrete, spacing(1));
+ else
+ [~, derivativeDiscrete] = gradient(objDiscrete, spacing(1));
+ end
+
+ rePermutationVector = [2:(gridSelectionIndex), ...
+ 1, (gridSelectionIndex+1):ndims(derivativeDiscrete)];
+ result = quantity.Discrete(...
+ permute(derivativeDiscrete, rePermutationVector), ...
+ 'size', size(obj), 'grid', obj(1).grid, ...
+ 'gridName', obj(1).gridName, ...
+ 'name', ['(d_{', diffGridName, '}', obj(1).name, ')']);
+
+ if k > 1
+ % % if a higher order derivative is requested, call the function
+ % % recursivly until the first-order derivative is reached
+ result = result.diff(k-1, diffGridName);
+ end
+ end
+
function [idx, permuteGrid] = computePermutationVectors(a, b)
% Computes the required permutation vectors to use
% misc.multArray for multiplication of matrices
-
+
uA = unique(a(1).gridName, 'stable');
assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
uB = unique(b(1).gridName, 'stable');
@@ -2313,7 +2446,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% quantity.Discrete class.
permuteGrid = [idxGrid(lGridA), idxGrid(lGridB), idxGrid(lGridVA), idxGrid(lGridVB)];
end
-
+
function f = setValueContinuous(obj, f)
end
@@ -2328,10 +2461,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% copied.
end
- % function s = getHeader(obj)
- % s = 'TEST'
- % end
-
function s = getPropertyGroups(obj)
% Function to display the correct values
@@ -2339,14 +2468,14 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
s = getPropertyGroups@matlab.mixin.CustomDisplay(obj);
return;
else
- s = getPropertyGroups@matlab.mixin.CustomDisplay(obj(1));
+ s = getPropertyGroups@matlab.mixin.CustomDisplay(obj(1));
end
if numel(obj) ~= 1
s.PropertyList.valueDiscrete = ...
[sprintf('%ix', gridSize(obj), size(obj)) sprintf('\b')];
end
- end
+ end % getPropertyGroups
end % methods (Access = protected)
methods ( Static, Access = protected )
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
new file mode 100644
index 0000000000000000000000000000000000000000..f9fcfbd62632982b5e632aef26a2ee109f01ee44
--- /dev/null
+++ b/+quantity/Domain.m
@@ -0,0 +1,405 @@
+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
+
diff --git a/+quantity/EquidistantDomain.m b/+quantity/EquidistantDomain.m
new file mode 100644
index 0000000000000000000000000000000000000000..ff763973dfff6ae55c52e884a40144faf8264875
--- /dev/null
+++ b/+quantity/EquidistantDomain.m
@@ -0,0 +1,24 @@
+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
+
diff --git a/+quantity/Function.m b/+quantity/Function.m
index 4b090acdc57b53388e2b9c58fb142af9540ba511..dc15eb148e7cbe49544e23c5928fbe061294f343 100644
--- a/+quantity/Function.m
+++ b/+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 ---
%--------------
diff --git a/+quantity/Operator.m b/+quantity/Operator.m
index d2ad64f5f0055cbfb8676a29f9483198cab0f4b3..197681e8f2d6b730f0553c391c2df005efc1c92d 100644
--- a/+quantity/Operator.m
+++ b/+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 = ...
- quantity.Discrete(A{k}, 'gridName', {}, 'grid', {});
+ quantity.Discrete(A{k}, 'domain', prsr.Results.domain);
else
obj(k).coefficient = A{k}; %#ok<AGROW>
end
end
elseif isnumeric(A)
- obj(1).coefficient = quantity.Discrete(A, 'gridName', {}, 'grid', {});
+ obj(1).coefficient = quantity.Discrete(A, 'domain', prsr.Results.domain);
end
[obj.s] = deal(s);
@@ -70,6 +105,29 @@ classdef Operator < handle & matlab.mixin.Copyable
I = quantity.Operator(II);
end
+ function value = applyTo(obj, y, varargin)
+ % APPLYON evaluation of the differential parameterization
+ % VALUE = applyOn(OBJ, Y, varargin) applies the operator OBJ
+ % on the function Y, i.e.,
+ % value = A0(t) y(t) + A1(t) dt y(t) + A2(t) dz^2 y(t)...
+ % is evaluated.
+
+ myParser = misc.Parser();
+ myParser.addParameter('domain', y.domain);
+ myParser.addParameter('n', length(obj));
+ myParser.parse(varargin{:});
+
+ myDomain = myParser.Results.domain;
+ n = myParser.Results.n;
+
+ value = zeros(size(obj(1).coefficient));
+
+ for k = 1:n
+ value = value + obj(k).coefficient * y.diff(k - 1, myDomain );
+ end
+
+ end
+
function C = mtimes(A, B)
% if isnumeric(A) && numel(A) == 1
@@ -105,6 +163,7 @@ classdef Operator < handle & matlab.mixin.Copyable
end
C = quantity.Operator(C);
end
+
function C = plus(A, B)
if ~isa(A, 'quantity.Operator')
@@ -130,9 +189,11 @@ classdef Operator < handle & matlab.mixin.Copyable
end
C = quantity.Operator(C);
end
+
function C = minus(A, B)
C = A + (-1)*B;
end
+
function C = adj(A)
assert(A(1).coefficient.isConstant())
@@ -145,6 +206,7 @@ classdef Operator < handle & matlab.mixin.Copyable
C = quantity.Operator(c);
end
+
function C = det(A)
assert(A(1).coefficient.isConstant())
@@ -158,13 +220,16 @@ classdef Operator < handle & matlab.mixin.Copyable
C = quantity.Operator(c);
end
+
function C = uminus(A)
C = (-1)*A;
end
+
function s = size(obj, varargin)
s = [length(obj), size(obj(1).coefficient)];
s = s(varargin{:});
end
+
function l = length(obj)
l = numel(obj);
end
@@ -415,8 +480,15 @@ classdef Operator < handle & matlab.mixin.Copyable
function newOperator = changeGrid(obj, gridNew, gridNameNew)
newOperator = quantity.Operator(obj.M.changeGrid(gridNew, gridNameNew));
end
+
function newOperator = subs(obj, varargin)
newOperator = obj.M.subs(varargin{:});
+
+ if isnumeric(newOperator)
+ newOperator = squeeze( num2cell(newOperator, [1 2]) );
+ newOperator = quantity.Operator(newOperator);
+ end
+
end
end
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index f65d413b5fa1393f0a1e10f99719deb5f60a8e09..e49270ab4594fa01bf9d5a6e9d3250296765374e 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -15,6 +15,12 @@ classdef Symbolic < quantity.Function
function obj = Symbolic(valueContinuous, varargin)
+ %TODO Contruktor überarbeiten:
+ % # bei matlabFunction( symbolischer ausdruck ) muss auf die
+ % Reihenfolge der variablen geachtet werden!
+ % # Manchmal werden mehrere grid variablen an den
+ % untergeordneten Konstruktor durchgereicht.
+
parentVarargin = {};
if nargin > 0
@@ -26,21 +32,22 @@ classdef Symbolic < quantity.Function
variableParser = misc.Parser();
variableParser.addParameter('variable', quantity.Symbolic.getVariable(valueContinuous));
variableParser.addParameter('symbolicEvaluation', false);
+ variableParser.addParameter('gridName', '');
+ variableParser.addParameter('domain', []);
variableParser.parse(varargin{:});
- variable = variableParser.Results.variable;
- numVar = numel(variable);
- defaultGrid = cell(1, numVar);
- for it = 1:numVar
- if it==1
- defaultGrid{it} = linspace(0, 1, 101).';
- else
- defaultGrid{it} = reshape(linspace(0, 1, 101), [ones(1, it-1), numel(linspace(0, 1, 101))]);
- end
+
+ if ~variableParser.isDefault('variable')
+ warning('Do not set the variable property. It will be ignored.')
end
- gridParser = misc.Parser();
- gridParser.addParameter('grid', defaultGrid);
- gridParser.parse(varargin{:});
+ if variableParser.isDefault('gridName') && variableParser.isDefault('domain')
+ warning('Grid names are generated from the symbolic variable. Please set the grid names in a quantity.Domain object explicitly')
+ variableNames = cellstr(string(variableParser.Results.variable));
+ varargin = [varargin(:)', {'gridName'}, {variableNames}];
+ end
+
+ [myDomain, varargin] = quantity.Domain.parser(varargin{:});
+
s = size(valueContinuous);
S = num2cell(s);
@@ -63,10 +70,10 @@ classdef Symbolic < quantity.Function
% variable
if isa(fun{k}, 'sym')
symb{k} = fun{k};
- fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, variable);
+ fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, {myDomain.name});
elseif isa(fun{k}, 'double')
symb{k} = sym(fun{k});
- fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, variable);
+ fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, {myDomain.name});
elseif isa(fun{k}, 'function_handle')
symb{k} = sym(fun{k});
else
@@ -74,16 +81,14 @@ classdef Symbolic < quantity.Function
end
end
- gridName = cellstr(string(variable));
- % check if the reuqired gridName fits to the current
- % variables:
- if isfield(gridParser.Unmatched, 'gridName') ...
- && ~all(strcmp(gridParser.Unmatched.gridName, gridName))
- error('If the gridName(s) are set explicitly, they have to be the same as the variable names!')
- end
+%
+% % check if the reuqired gridName fits to the current
+% % variables:
+% if ~all(strcmp(sort({myDomain.name}), sort(variableName)))
+% error('If the gridName(s) are set explicitly, they have to be the same as the variable names!')
+% end
- parentVarargin = [{fun}, varargin(:)', {'grid'}, {gridParser.Results.grid}, ...
- {'gridName'}, {gridName}];
+ parentVarargin = [{fun}, varargin(:)', {'domain'}, {myDomain}];
end
% call parent class
obj@quantity.Function(parentVarargin{:});
@@ -96,7 +101,8 @@ classdef Symbolic < quantity.Function
% special properties
for k = 1:numel(valueContinuous)
- obj(k).variable = variable; %#ok<AGROW>
+ obj(k).variable = ...
+ quantity.Symbolic.getVariable(valueContinuous); %#ok<AGROW>
obj(k).valueSymbolic = symb{k}; %#ok<AGROW>
obj(k).symbolicEvaluation = ...
variableParser.Results.symbolicEvaluation; %#ok<AGROW>
@@ -243,11 +249,9 @@ classdef Symbolic < quantity.Function
obj = objCell{objIdx};
if ~isempty(obj)
- myVariable = obj(1).variable;
- myGrid = obj(1).grid;
+ myDomain = obj(1).domain;
else
- myVariable = '';
- myGrid = {};
+ myDomain = quantity.Domain.empty();
end
for k = 1:numel(objCell)
@@ -256,8 +260,7 @@ classdef Symbolic < quantity.Function
o = objCell{k};
elseif isa(objCell{k}, 'double') || isa(objCell{k}, 'sym')
o = quantity.Symbolic( objCell{k}, ...
- 'variable', myVariable, ...
- 'grid', myGrid);
+ 'domain', myDomain);
end
objCell{k} = o;
end
@@ -267,23 +270,53 @@ classdef Symbolic < quantity.Function
function solution = subs(obj, gridName, values)
%% SUBS substitues gridName and values
- % solution = subs(obj, gridName, values) TODO
- % This function substitutes the variables specified with
- % gradName with values. It can be used to rename grids or to
- % evaluate (maybe just some) grids.
- % GridName is cell-array of char-arrays
- % chosen from obj.gridName-property. Values is a cell-array of
- % the size of gridName. Each cell can contain arrays themself,
- % but those arrays must be of same size. Values can be
- % char-arrays standing for new gridName or or numerics.
- % In contrast to on() or at(), only some, but not necessarily
- % all variables are evaluated.
- if ~iscell(gridName)
- gridName = {gridName};
+ % solution = subs(obj, DOMAIN) can be used to change the grid
+ % of this quantity. The domains with the same domain name as
+ % DOMAIN will be replaced with the corresponding grid in
+ % DOMAIN.
+ %
+ % solution = subs(obj, DOMAINNAME, GRID) can be used to
+ % change the grid of this quantity. The grid of the domains with the
+ % DOMAINNAME will be replaced by the values specified in
+ % GRID. DOMAINNAME must be a cell-array containing the names
+ % and GRID must be a cell-array containing the new grids.
+ %
+ % solution = sub(obj, DOMAINNAME, NEWDOMAIN) can be used to
+ % change the grid and/or the name of a domain. DOMAINNAME
+ % must be a cell-array of domain names, or a oboject-array of
+ % quantity.Domain objects. The NEWDOMAIN must be an
+ % object-array of quantity.Domain.
+
+ if isa(gridName, 'quantity.Domain')
+ if nargin == 2
+ gridName = {gridName.name};
+ values = {gridName.grid};
+ else
+ gridName = {gridName.name};
+ end
+ else
+ gridName = misc.ensureIsCell(gridName);
end
- if ~iscell(values)
- values = {values};
+
+ if nargin == 3 && isa(values, 'quantity.Domain')
+ % replacement of the grid AND the gridName
+ % 1) replace the grid
+ solution = obj.changeGrid( {values.grid}, gridName );
+ % 2) replace the name
+ solution = solution.subs( gridName, {values.name} );
+ return
+ else
+ values = misc.ensureIsCell(values);
end
+
+% % ensure that domains which should be replaced are known by the
+% % object
+% for k = 1:length(values)
+% if ischar(values{k})
+% assert( any( strcmp( obj(1).gridName, values{k} ) ), 'The domain name to be replaced must be a domain name of the object' );
+% end
+% end
+
isNumericValue = cellfun(@isnumeric, values);
if any((cellfun(@(v) numel(v(:)), values)>1) & isNumericValue)
error('only implemented for one value per grid');
@@ -342,7 +375,7 @@ classdef Symbolic < quantity.Function
end
end
solution = quantity.Symbolic(double(symbolicSolution), ...
- 'grid', newGrid, 'variable', newGridName, 'name', obj(1).name);
+ 'grid', newGrid, 'gridName', newGridName, 'name', obj(1).name);
else
% before creating a new quantity, it is checked that
% newGridName is unique. If there are non-unique
@@ -351,7 +384,7 @@ classdef Symbolic < quantity.Function
uniqueGridName = unique(newGridName, 'stable');
if numel(newGridName) == numel(uniqueGridName)
solution = quantity.Symbolic(symbolicSolution, ...
- 'grid', newGrid, 'variable', newGridName, 'name', obj(1).name);
+ 'grid', newGrid, 'gridName', newGridName, 'name', obj(1).name);
else
uniqueGrid = cell(1, numel(uniqueGridName));
for it = 1 : numel(uniqueGrid)
@@ -370,10 +403,9 @@ classdef Symbolic < quantity.Function
end
end
solution = quantity.Symbolic(symbolicSolution, ...
- 'grid', uniqueGrid, 'variable', uniqueGridName, 'name', obj(1).name);
+ 'grid', uniqueGrid, 'gridName', uniqueGridName, 'name', obj(1).name);
end
end
-
end
end % subs()
@@ -469,7 +501,7 @@ classdef Symbolic < quantity.Function
obj(it).variable, varNew(s)), varNew(0)==ic);
end
solution = quantity.Symbolic(symbolicSolution, ...
- 'gridName', {myParser.Results.newGridName, 'ic'}, 'variable', {s, ic}, ...
+ 'gridName', {myParser.Results.newGridName, 'ic'}, ...
'grid', {myParser.Results.variableGrid, myParser.Results.initialValueGrid}, ...
'name', ['solve(', obj(1).name, ')']);
end % solveDVariableEqualQuantity()
@@ -488,45 +520,17 @@ classdef Symbolic < quantity.Function
function y = sqrt(x)
% quadratic root for scalar and diagonal symbolic quantities
y = quantity.Symbolic(sqrt(x.sym()), ...
- 'grid', x(1).grid, 'variable', x(1).variable, ...
+ 'domain', x(1).domain, ...
'name', ['sqrt(', x(1).name, ')']);
end % sqrt()
function y = sqrtm(x)
% quadratic root for matrices of symbolic quantities
y = quantity.Symbolic(sqrtm(x.sym()), ...
- 'grid', x(1).grid, 'variable', x(1).variable, ...
+ 'domain', x(1).domain, ...
'name', ['sqrtm(', x(1).name, ')']);
end % sqrtm()
- function result = diff(obj, k, diffGridName)
- % diff applies the kth-derivative for the variable specified with
- % the input gridName to the obj. If no gridName is specified, then diff
- % applies the derivative w.r.t. to all gridNames / variables.
- if nargin == 1 || isempty(k)
- k = 1; % by default, only one derivatve per diffGridName is applied
- end
- if nargin <= 2 % if no diffGridName is specified, then the derivative
- % w.r.t. to all gridNames is applied
- diffGridName = obj(1).gridName;
- end
-
- result = obj.copy();
- if iscell(diffGridName)
- for it = 1 : numel(diffGridName)
- result = result.diff(k, diffGridName{it});
- end
- else
- diffVariable = obj.gridName2variable(diffGridName);
- [result.name] = deal(['(d_{' char(diffVariable) '} ' result(1).name, ')']);
- [result.valueDiscrete] = deal([]);
- for l = 1:numel(obj)
- result(l).valueSymbolic = diff(obj(l).valueSymbolic, diffVariable, k);
- result(l).valueContinuous = obj.setValueContinuous(result(l).valueSymbolic, obj(1).variable);
- end
- end
- end % diff()
-
function b = flipGrid(a, myGridName)
if ~iscell(myGridName)
myGridName = {myGridName};
@@ -538,7 +542,7 @@ classdef Symbolic < quantity.Function
variableNew{it} = 1-variableOld{it};
end
b = quantity.Symbolic(subs(a.sym, variableOld, variableNew), ...
- 'grid', a(1).grid, 'variable', a(1).variable, ...
+ 'domain', a(1).domain, ...
'name', ['flip(', a(1).name, ')']);
end % flipGrid()
@@ -569,35 +573,13 @@ classdef Symbolic < quantity.Function
function mObj = uminus(obj)
% unitary minus: C = -A
mObj = quantity.Symbolic(-obj.sym, 'grid', obj(1).grid, ...
- 'variable', obj(1).variable, 'name', ['-', obj(1).name]);
+ 'domain', obj(1).domain, 'name', ['-', obj(1).name]);
end % uminus()
function mObj = uplus(obj)
% unitary plus: C = +A
mObj = copy(obj);
end % uplus()
-
- function parameters = combineGridGridNameVariable(A, B, parameters)
- assert(isa(A, 'quantity.Symbolic') && isa(B, 'quantity.Symbolic'), ...
- 'A and B must be quantity.Symbolic');
-
- if nargin < 3
- parameters = struct(A(1));
- end
- % combine grids of two quantities. This is often used in mathematical
- % operations to obtain the parameters of the resulting quantity.
- idx = computePermutationVectors(A, B);
- parameters.grid = ...
- [A(1).grid(idx.A.common), ...
- A(1).grid(~idx.A.common), ...
- B(1).grid(~idx.B.common)];
- parameters.gridName = [A(1).gridName(idx.A.common), ...
- A(1).gridName(~idx.A.common), ...
- B(1).gridName(~idx.B.common)];
- parameters.variable = [A(1).variable(idx.A.common), ...
- A(1).variable(~idx.A.common), ...
- B(1).variable(~idx.B.common)];
- end % combineGridGridNameVariable()
-
+
function C = mtimes(A, B)
% if one input is ordinary matrix, this is very simple.
if isempty(B) || isempty(A)
@@ -609,13 +591,13 @@ classdef Symbolic < quantity.Function
end
if isnumeric(B)
C = quantity.Symbolic(A.sym() * B, ...
- 'grid', A(1).grid, 'variable', A(1).variable, ...
+ 'domain', A(1).domain, ...
'name', [A(1).name, ' c']);
return
end
if isnumeric(A)
C = quantity.Symbolic(A * B.sym(), ...
- 'grid', B(1).grid, 'variable', B(1).variable, ...
+ 'domain', B(1).domain, ...
'name', ['c ', B(1).name]);
return
end
@@ -624,7 +606,7 @@ classdef Symbolic < quantity.Function
return;
end
- parameters = combineGridGridNameVariable(A, B);
+ parameters.domain = gridJoin(A(1).domain, B(1).domain);
parameters.name = [A(1).name, ' ', B(1).name];
parameters = misc.struct2namevaluepair(parameters);
@@ -635,16 +617,32 @@ classdef Symbolic < quantity.Function
% minus() see quantity.Discrete.
assert(isequal(size(a), size(b)), 'terms must be of same size');
if ~isa(a, 'quantity.Discrete')
- a = quantity.Symbolic(a, 'name', 'c');
+ if isnumeric(a)
+ a = quantity.Symbolic(a, 'name', 'c', 'domain', quantity.Domain.empty());
+ else
+ error('Not yet implemented')
+ end
parameters = struct(b(1));
else
if ~isa(b, 'quantity.Discrete')
- b = quantity.Symbolic(b, 'name', 'c');
+ if isnumeric(b)
+ b = quantity.Symbolic(b, 'name', 'c', 'domain', quantity.Domain.empty());
+ else
+ error('Not yet implemented')
+ end
end
parameters = struct(a(1));
end
- parameters = combineGridGridNameVariable(a, b, parameters);
+ % remove the not required parameters:
+ for k = {'domain', 'grid', 'gridName', 'variable'}
+ if isfield(parameters, k)
+ parameters = rmfield(parameters, k);
+ end
+ end
+
+ % set new parameter values:
+ parameters.domain = gridJoin(a(1).domain, b(1).domain);
parameters.name = [a(1).name, '+' , b(1).name];
parameters = misc.struct2namevaluepair(parameters);
@@ -654,56 +652,50 @@ classdef Symbolic < quantity.Function
function absQuantity = abs(obj)
% abs returns the absolut value of the quantity as a quantity
absQuantity = quantity.Symbolic(abs(obj.sym()), ...
- 'grid', obj(1).grid, 'variable', obj(1).variable, ...
- 'name', ['|', obj(1).name, '|']);
+ 'name', ['|', obj(1).name, '|'], ...
+ 'domain', obj(1).domain);
end % abs()
function y = real(obj)
% real() returns the real part of the obj.
y = quantity.Symbolic(real(obj.sym()), ...
'name', ['real(', obj(1).name, ')'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % real()
function y = imag(obj)
% imag() returns the real part of the obj.
y = quantity.Symbolic(imag(obj.sym()), ...
'name', ['real(', obj(1).name, ')'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % imag()
function y = inv(obj)
% inv inverts the matrix obj.
y = quantity.Symbolic(inv(obj.sym()), ...
'name', ['(', obj(1).name, ')^{-1}'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % inv()
function y = exp(obj)
% exp() is the exponential function using obj as the exponent.
y = quantity.Symbolic(exp(obj.sym()), ...
'name', ['exp(', obj(1).name, ')'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % exp()
function y = expm(obj)
% exp() is the matrix-exponential function using obj as the exponent.
y = quantity.Symbolic(expm(obj.sym()), ...
'name', ['expm(', obj(1).name, ')'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % expm()
function y = ctranspose(obj)
% ctranspose() or ' is the complex conjugate tranpose
y = quantity.Symbolic(conj(obj.sym().'), ...
'name', ['{', obj(1).name, '}^{H}'], ...
- 'grid', obj(1).grid, ...
- 'variable', obj(1).variable);
+ 'domain', obj(1).domain);
end % expm()
function C = int(obj, z, a, b)
@@ -780,7 +772,7 @@ classdef Symbolic < quantity.Function
% with numeric integral bounds. Therefore, this try-catch block
% is needed.
C = quantity.Symbolic(I, ...
- 'grid', gridI, 'variable', variableI, ...
+ 'grid', gridI, 'gridName', gridNameI, ...
'name', ['int[', obj(1).name, ']']);
catch
C = cumInt(quantity.Discrete(obj), z, a, b);
@@ -872,6 +864,16 @@ classdef Symbolic < quantity.Function
end
end
+ function result = diff_inner(obj, k, diffGridName)
+ result = obj.copy();
+ [result.name] = deal(['(d_{' diffGridName '} ' result(1).name, ')']);
+ [result.valueDiscrete] = deal([]);
+ for l = 1:numel(obj)
+ result(l).valueSymbolic = diff(obj(l).valueSymbolic, diffGridName, k);
+ result(l).valueContinuous = obj.setValueContinuous(result(l).valueSymbolic, obj(1).variable);
+ end
+ end
+
end
methods (Static)
@@ -937,11 +939,7 @@ classdef Symbolic < quantity.Function
end
function f = setValueContinuous(f, symVar)
% converts symVar into a matlabFunction.
- if isa(f, 'sym')
- if nargin == 1
- symVar = quantity.Symbolic.getVariable(f);
- end
-
+ if isa(f, 'sym')
if isempty(symvar(f))
f = @(varargin) double(f) + quantity.Function.zero(varargin{:});
else
diff --git a/+quantity/TimeDelayOperator.m b/+quantity/TimeDelayOperator.m
new file mode 100644
index 0000000000000000000000000000000000000000..4ad46845500091f78be38a0fda10a561284af49e
--- /dev/null
+++ b/+quantity/TimeDelayOperator.m
@@ -0,0 +1,130 @@
+classdef TimeDelayOperator < handle & matlab.mixin.Copyable
+ %TIMEDELAYOPERATOR class to describe an operator of the form
+ % D[h](z,t) = h( t + coefficient(z) )
+ % This comes from the inverse Laplace transfromed value of
+ % exp(-s coefficient(z) ) * H(s) <-> h(t + coefficient(z))
+ %
+
+ properties
+ coefficient {mustBe.quantity};
+ end
+
+ properties (Access = protected)
+ spatialDomain_inner quantity.Domain;
+ timeDomain_inner quantity.Domain;
+ isZero logical = false;
+ end
+
+ methods
+ function obj = TimeDelayOperator(coefficient, timeDomain, optionalArgs)
+ %TIMEDELAYOPERATOR Construct an instance of a time delay
+ %operator
+ % obj = TimeDelayOperator(COEFFICIENT, TIMEDOMAIN, VARARGIN)
+ % will construct a time delay operator of the form
+ % OBJ[h](z,t) = h( t + COEFFICIENT(z) )
+ % with the time delay described by COEFFICIENT on a temporal
+ % domain TIMEDOMAIN. Using the name-value-pairs optional
+ % parameters, a spatial domain can be specified by the name
+ % 'spatialDomain'.
+ arguments
+ coefficient;
+ timeDomain;
+ optionalArgs.spatialDomain = coefficient.domain;
+ optionalArgs.isZero logical = false;
+ end
+
+ obj.coefficient = coefficient;
+ obj.timeDomain_inner = timeDomain;
+ obj.spatialDomain_inner = optionalArgs.spatialDomain;
+ obj.isZero = optionalArgs.isZero;
+ end
+
+ 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.
+ % d = OBJ[h]
+ 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);
+
+ d = quantity.Discrete.zeros([n, l], [obj.spatialDomain, obj.timeDomain]);
+
+ for i = 1 : n
+ 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, 'domain', optionalArgs.domain );
+ end
+ end
+ end
+ end
+ end
+
+ function d = compositionDomain( obj, h)
+
+ % todo assert all coefficients have the same domain
+
+ d = h(1).compositionDomain( obj.t + obj(1).coefficient );
+ end
+
+ function t = t(obj)
+ t = quantity.Symbolic(sym(obj.timeDomain.name), 'domain', obj.timeDomain);
+ end
+
+ function d = diag(obj)
+
+ n = length(obj);
+ idx = 1:n;
+
+ d(1:(n^2)) = obj.zero(obj.timeDomain);
+ d = reshape(d, n, n);
+ d(idx == idx') = obj(:);
+
+ end
+
+ function C = coefficients(obj)
+ C = reshape([obj.coefficient], size(obj));
+ end
+
+ function sd = spatialDomain(obj)
+ sd = join( [obj.spatialDomain_inner] );
+ end
+
+ function td = timeDomain(obj)
+ td = join( [obj.timeDomain_inner] );
+ end
+
+
+ end
+
+ methods (Static)
+ function o = zero(timeDomain, namedArgs)
+
+ arguments
+ timeDomain quantity.Domain;
+ namedArgs.spatialDomain quantity.Domain = quantity.Domain.empty();
+ end
+
+ if isempty( namedArgs.spatialDomain )
+ infValue = inf;
+ else
+ infValue = inf(obj.spatialDomain.gridLength, 1);
+ end
+
+ qinf = quantity.Discrete(infValue, ...
+ 'domain', namedArgs.spatialDomain);
+ o = quantity.TimeDelayOperator(qinf, timeDomain, 'isZero', true);
+ end
+ end
+
+end
+
diff --git a/+signals/+faultmodels/TaylorPolynomial.m b/+signals/+faultmodels/TaylorPolynomial.m
index c5363b2bab45fc6de420e35b07d6b4728923659c..d451e9b9e72eb5b671f1eb77a1955beb8431c005 100644
--- a/+signals/+faultmodels/TaylorPolynomial.m
+++ b/+signals/+faultmodels/TaylorPolynomial.m
@@ -128,12 +128,11 @@ classdef TaylorPolynomial < quantity.Discrete
end
function monomials = getMonomials(degree, t, tStar, myLocalGrid, T)
monomials = quantity.Symbolic.zeros([degree, 1], ...
- 'variable', {t, tStar}, ...
- 'grid', {myLocalGrid, T}, ...
+ 'domain', quantity.Domain.empty(), ...
'name', 'varphi_k');
for k = 1:degree
monomials(k, 1) = quantity.Symbolic((t-tStar)^(k-1), ...
- 'variable', {t, tStar}, ...
+ 'gridName', cellfun(@char, {t, tStar}, 'UniformOutput', false), ...
'grid', {myLocalGrid, T}, ...
'name', 'varphi_k');
end
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index fc20f9e5dab6eb9c45f2b8ae6089972a00c088bc..65fb584378377f8cb3e71efe162818e8a7a73a8b 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -4,6 +4,50 @@ function [tests] = testDiscrete()
tests = functiontests(localfunctions);
end
+function testCompose(testCase)
+
+t = quantity.Domain('t', linspace(0, 2));
+f = quantity.Discrete(sin(t.grid * 2 * pi) , 'domain', t, 'name', 'f' );
+g = quantity.Discrete( 0.5 * t.grid, 'domain', t, 'name', 'g');
+h = quantity.Discrete( 2 * t.grid, 'domain', t, 'name', 'h');
+
+% verify simple evaluation of f( g(t) )
+testCase.verifyEqual( on( f.compose(g) ), sin(t.grid * pi), 'AbsTol', 3e-3 )
+
+% verify that a warning is thrown, if the domain of definition is exceeded
+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('tau', linspace(-1, 3));
+z = quantity.Domain('z', linspace(0, 1, 51));
+
+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('z', linspace(0, 1, 11));
+t = quantity.Domain('t', linspace(0, 1, 21));
+tau = quantity.Domain('tau', linspace(0, 2, 31));
+
+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)
z = linspace(0, 2*pi)';
@@ -24,12 +68,37 @@ tc.verifyEqual(blubQuadraticNorm.on(), 2*ones(11,1));
tc.verifyEqual(blubQuadraticWeightedNorm.on(), 4*ones(11,1));
end % testQuadraticNorm()
+function testSort(tc)
+
+t = linspace(0, pi, 7)';
+domA = quantity.Domain('a', t);
+domB = quantity.Domain('b', t);
+q = quantity.Discrete(sin(t) .* cos(t'), ...
+ 'domain', [domA, domB]);
+
+q1V = q.valueDiscrete();
+
+q.sort('descend');
+
+q2V = q.valueDiscrete();
+
+tc.verifyEqual(q1V , q2V')
+
+
+end % testSort
+
function testBlkdiag(tc)
% init some data
syms z zeta
-A = quantity.Symbolic(...
- [1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
+
+z = quantity.Domain('z', linspace(0, 1, 21));
+zeta = quantity.Domain('zeta', linspace(0, 1, 41));
+
+A = quantity.Discrete(cat(4, ...
+ cat(3, 1 + z.grid .* zeta.grid', - z.ones .* zeta.grid'), ...
+ cat(3, -z.grid .* zeta.ones', z.grid .^2 .* zeta.ones') ), ...
+ 'domain', [z, zeta], 'name', 'q');
+
B = 2*A(:, 1);
C = 3*A(1);
@@ -66,7 +135,7 @@ function testCtranspose(tc)
syms z zeta
qSymbolic = quantity.Symbolic(...
[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
+ 'gridName', {'z', 'zeta'}, 'name', 'q');
qDiscrete = quantity.Discrete(qSymbolic);
qDiscreteCtransp = qDiscrete';
tc.verifyEqual(qDiscrete(1,1).on(), qDiscreteCtransp(1,1).on());
@@ -82,12 +151,14 @@ tc.verifyEqual(qDiscrete2.imag.on(), -qDiscrete2Ctransp.imag.on());
end % testCtranspose
function testTranspose(tc)
+
syms z zeta
qSymbolic = quantity.Symbolic(...
[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
+ 'gridName', {'z', 'zeta'}, 'name', 'q');
qDiscrete = quantity.Discrete(qSymbolic);
qDiscreteTransp = qDiscrete.';
+
tc.verifyEqual(qDiscrete(1,1).on(), qDiscreteTransp(1,1).on());
tc.verifyEqual(qDiscrete(2,2).on(), qDiscreteTransp(2,2).on());
tc.verifyEqual(qDiscrete(1,2).on(), qDiscreteTransp(2,1).on());
@@ -98,16 +169,16 @@ function testFlipGrid(tc)
syms z zeta
myGrid = linspace(0, 1, 11);
f = quantity.Discrete(quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid}));
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
% % flip one grid
fReference = quantity.Discrete(quantity.Symbolic([1+z+(1-zeta); 2*(1-zeta)+sin(z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid}));
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
fFlipped = f.flipGrid('zeta');
tc.verifyEqual(fReference.on(), fFlipped.on(), 'AbsTol', 10*eps);
% flip both grids
fReference2 = quantity.Discrete(quantity.Symbolic([1+(1-z)+(1-zeta); 2*(1-zeta)+sin(1-z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid}));
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
fFlipped2 = f.flipGrid({'z', 'zeta'});
fFlipped3 = f.flipGrid({'zeta', 'z'});
tc.verifyEqual(fReference2.on(), fFlipped2.on(), 'AbsTol', 10*eps);
@@ -115,19 +186,22 @@ tc.verifyEqual(fReference2.on(), fFlipped3.on(), 'AbsTol', 10*eps);
end % testFlipGrid();
function testScalarPlusMinusQuantity(testCase)
-syms z
-myGrid = linspace(0, 1, 7);
-f = quantity.Discrete(quantity.Symbolic([1; 2] + zeros(2, 1)*z, ...
- 'variable', {z}, 'grid', {myGrid}));
+
+z = quantity.Domain('z', linspace(0, 1, 7));
+f = quantity.Discrete( ( [2; 1] + zeros(2, 1) .* z.grid' )', ...
+ 'domain', z);
+
testCase.verifyError(@() 1-f-1, '');
testCase.verifyError(@() 1+f+1, '');
+
end % testScalarPlusMinusQuantity
function testNumericVectorPlusMinusQuantity(testCase)
-syms z
-myGrid = linspace(0, 1, 7);
-f = quantity.Discrete(quantity.Symbolic([1+z; 2+sin(z)] + zeros(2, 1)*z, ...
- 'variable', {z}, 'grid', {myGrid}));
+
+z = quantity.Domain('z', linspace(0, 1, 7));
+f = quantity.Discrete( [1+z.grid, 2+sin(z.grid)] , ...
+ 'domain', z);
+
a = ones(size(f));
testCase.verifyEqual(on(a-f), 1-f.on());
testCase.verifyEqual(on(f-a), f.on()-1);
@@ -136,19 +210,18 @@ testCase.verifyEqual(on(f+a), f.on()+1);
end % testNumericVectorPlusMinusQuantity
function testUnitaryPluasAndMinus(testCase)
-syms z zeta
-qSymbolic = quantity.Symbolic(...
- [1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
-qDiscrete = quantity.Discrete(qSymbolic);
-qDoubleArray = qSymbolic.on();
+z = quantity.Domain('z', linspace(0, 1, 7));
+zeta = quantity.Domain('zeta', linspace(0, 1, 7));
+
+q = quantity.Discrete(cat(4, ...
+ cat(3, 1 + z.grid .* zeta.grid', - z.grid .* zeta.ones'), ...
+ cat(3, - z.ones .* zeta.grid', z.grid.^2 .* zeta.ones') ), ...
+ 'domain', [z, zeta]) ;
-testCase.verifyEqual(on(-qSymbolic), -qDoubleArray);
-testCase.verifyEqual(on(-qDiscrete), -qDoubleArray);
-testCase.verifyEqual(on(+qSymbolic), +qDoubleArray);
-testCase.verifyEqual(on(+qDiscrete), +qDoubleArray);
-testCase.verifyEqual(on(+qDiscrete), on(+qSymbolic));
-testCase.verifyEqual(on(-qDiscrete), on(-qSymbolic));
+qDoubleArray = q.on();
+
+testCase.verifyEqual(on(-q), -qDoubleArray);
+testCase.verifyEqual(on(+q), +qDoubleArray);
end
@@ -194,31 +267,51 @@ testCase.verifyEqual(DE.on(), compareDE.on())
ED = [E, D];
compareED = quantity.Discrete({tan(t*z), sin(t * z); cos(t * z), cos(t*z)}, 'grid', {t, z}, 'gridName', {'t', 'z'});
testCase.verifyEqual(ED.on(), compareED.on())
+
+% concatenation of a constant and a function:
+C = [5; 2];
+const = quantity.Discrete(C, 'domain', quantity.Domain.empty);
+%
+constA = [const, A];
+testCase.verifyTrue(all(constA(1,1).on() == C(1)));
+testCase.verifyTrue(all(constA(2,1).on() == C(2)));
+
+Aconst = [A, const];
+testCase.verifyTrue(all(Aconst(1,2).on() == C(1)));
+testCase.verifyTrue(all(Aconst(2,2).on() == C(2)));
+
+
end
function testExp(testCase)
% 1 spatial variable
-syms z zeta
-s1d = quantity.Discrete(quantity.Symbolic(...
- 1+z*z, 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 's1d'));
+z = quantity.Domain('z', linspace(0, 1, 3));
+s1d = quantity.Discrete(1 + z.grid .* z.grid.', 'domain', z);
testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));
% diagonal matrix
-s2dDiag = quantity.Discrete(quantity.Symbolic(...
- [1+z*zeta, 0; 0, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2dDiag'));
+zeta = quantity.Domain('zeta', linspace(0, 1, 4));
+s2dDiag = quantity.Discrete(cat(4, ...
+ cat(3, 1 + z.grid .* zeta.grid', zeros(z.n, zeta.n)), ....
+ cat(3, zeros(z.n, zeta.n), z.grid.^2 .* zeta.ones')), 'domain', [z, zeta]);
+
testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
end
function testExpm(testCase)
-syms z zeta
-mat2d = quantity.Discrete(quantity.Symbolic(...
- ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2d'));
+
+z = quantity.Domain('z', linspace(0, 1, 3));
+zeta = quantity.Domain('zeta', linspace(0, 1, 4));
+
+mat2d = quantity.Discrete(...
+ {1 + z.grid .* zeta.grid', 3 * z.ones .* zeta.grid'; ...
+ 2 + 5 * z.grid + zeta.grid', z.grid .^2 .* zeta.ones'}, ...
+ 'domain', [z, zeta]);
+
mat2dMat = mat2d.on();
mat2dExpm = 0 * mat2d.on();
-for zIdx = 1 : 21
- for zetaIdx = 1 : 41
+for zIdx = 1 : z.n
+ for zetaIdx = 1 : zeta.n
mat2dExpm(zIdx, zetaIdx, :, :) = expm(reshape(mat2dMat(zIdx, zetaIdx, :, :), [2, 2]));
end
end
@@ -289,7 +382,8 @@ end
function testDiag2Vec(testCase)
% quantity.Symbolic
syms z
-myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2]);
+do = quantity.Domain('z', linspace(0,1,3));
+myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], 'domain', do);
myVectorSymbolic = diag2vec(myMatrixSymbolic);
testCase.verifyEqual(myMatrixSymbolic(1,1).valueContinuous, myVectorSymbolic(1,1).valueContinuous);
testCase.verifyEqual(myMatrixSymbolic(2,2).valueContinuous, myVectorSymbolic(2,1).valueContinuous);
@@ -302,6 +396,7 @@ testCase.verifyEqual(myMatrixDiscrete(1,1).on(), myVectorDiscrete(1,1).on());
testCase.verifyEqual(myMatrixDiscrete(2,2).on(), myVectorDiscrete(2,1).on());
testCase.verifyEqual(numel(myVectorDiscrete), size(myMatrixDiscrete, 1));
end
+
function testVec2Diag(testCase)
% quantity.Discrete
n = 7;
@@ -343,39 +438,38 @@ end
function testSolveDVariableEqualQuantityConstant(testCase)
%% simple constant case
-quan = quantity.Discrete(ones(51, 1), 'grid', linspace(0, 1, 51), ...
+quan = quantity.Discrete(ones(3, 1), 'grid', linspace(0, 1, 3), ...
'gridName', 'z', 'size', [1, 1], 'name', 'a');
solution = quan.solveDVariableEqualQuantity();
-[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 51), linspace(0, 1, 51));
+[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
end
function testSolveDVariableEqualQuantityNegative(testCase)
syms z zeta
assume(z>0 & z<1); assume(zeta>0 & zeta<1);
-myParameterGrid = linspace(0, 1, 51);
+myParameterGrid = linspace(0, 0.1, 5);
Lambda = quantity.Symbolic(-0.1-z^2, ...%, -1.2+z^2]),...1+z*sin(z)
- 'variable', z, ...
'grid', myParameterGrid, 'gridName', 'z', 'name', '\Lambda');
%%
-myGrid = linspace(-2, 2, 101);
-sEval = -0.9;
+myGrid = linspace(-.2, .2, 11);
+sEval = -0.09;
solveThisOdeDiscrete = solveDVariableEqualQuantity(...
diag2vec(quantity.Discrete(Lambda)), 'variableGrid', myGrid);
solveThisOdeSymbolic = solveDVariableEqualQuantity(...
diag2vec(Lambda), 'variableGrid', myGrid);
-testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.5}), solveThisOdeDiscrete.on({sEval, 0.5}), 'RelTol', 5e-4);
+testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.05}), solveThisOdeDiscrete.on({sEval, 0.05}), 'RelTol', 5e-3);
end
function testSolveDVariableEqualQuantityComparedToSym(testCase)
%% compare with symbolic implementation
syms z
assume(z>0 & z<1);
-quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 21)}, ...
- 'gridName', 'z', 'name', 'bSym', 'variable', {z});
-quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 21)}, ...
+quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 5)}, ...
+ 'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
+quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 5)}, ...
'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
solutionBSym = quanBSym.solveDVariableEqualQuantity();
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
@@ -387,10 +481,11 @@ function testSolveDVariableEqualQuantityAbsolut(testCase)
%% compare with symbolic implementation
syms z
assume(z>0 & z<1);
-quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 51)}, ...
- 'gridName', 'z', 'name', 'bSym', 'variable', {z});
-quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 51)}, ...
+quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 11)}, ...
+ 'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
+quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 11)}, ...
'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
+
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
myGrid = solutionBDiscrete.grid{1};
solutionBDiscreteDiff = solutionBDiscrete.diff(1, 's');
@@ -399,9 +494,10 @@ for icIdx = 1 : length(myGrid)
quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
end
%%
-testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-3);
+testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-2);
end
+
function testMtimesDifferentGridLength(testCase)
%%
a = quantity.Discrete(ones(11, 1), 'grid', linspace(0, 1, 11), 'gridName', 'z', ...
@@ -411,10 +507,9 @@ b = quantity.Discrete(ones(5, 1), 'grid', linspace(0, 1, 5), 'gridName', 'z', ..
ab = a*b;
%
syms z
-c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'variable', z, 'name', 'c');
+c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'gridName', 'z', 'name', 'c');
ac = a*c;
-%%
%%
testCase.verifyEqual(ab.on(), ones(11, 1));
testCase.verifyEqual(ac.on(), ones(21, 1));
@@ -472,23 +567,28 @@ testCase.verifyEqual(myQuantityDZzeta(3).on(), zeros(11, 21), 'AbsTol', 10*eps);
end
function testOn(testCase)
-%% init data
-gridVecA = linspace(0, 1, 27);
-gridVecB = linspace(0, 1, 41);
-[value] = createTestData(gridVecA, gridVecB);
+% init data
+domainVecA = quantity.Domain('z', linspace(0, 1, 27));
+domainVecB = quantity.Domain('zeta', linspace(0, 1, 41));
+[value] = createTestData(domainVecA.grid, domainVecB.grid);
a = quantity.Discrete(value, 'size', [2, 3], ...
- 'grid', {gridVecA, gridVecB}, 'gridName', {'z', 'zeta'}, 'name', 'A');
+ 'domain', [domainVecA, domainVecB], 'name', 'A');
-%%
+%
testCase.verifyEqual(value, a.on());
testCase.verifyEqual(permute(value, [2, 1, 3, 4]), ...
- a.on({gridVecB, gridVecA}, {'zeta', 'z'}));
+ a.on([domainVecB, domainVecA]));
testCase.verifyEqual(createTestData(linspace(0, 1, 100), linspace(0, 1, 21)), ...
a.on({linspace(0, 1, 100), linspace(0, 1, 21)}), 'AbsTol', 100*eps);
+
+d24z = quantity.Domain('z', linspace(0, 1, 24));
+d21zeta = quantity.Domain('zeta', linspace(0, 1, 21));
+
testCase.verifyEqual(createTestData(linspace(0, 1, 24), linspace(0, 1, 21)), ...
- a.on({linspace(0, 1, 24), linspace(0, 1, 21)}, {'z', 'zeta'}), 'AbsTol', 24*eps);
+ a.on( [d24z, d21zeta] ), 'AbsTol', 24*eps);
+
testCase.verifyEqual(permute(createTestData(linspace(0, 1, 21), linspace(0, 1, 24)), [2, 1, 3, 4]), ...
- a.on({linspace(0, 1, 24), linspace(0, 1, 21)}, {'zeta', 'z'}), 'AbsTol', 2e-4);
+ a.on( {linspace(0, 1, 24), linspace(0, 1, 21)}, {'zeta', 'z'}), 'AbsTol', 2e-4);
function [value] = createTestData(gridVecA, gridVecB)
[zGridA , zetaGridA] = ndgrid(gridVecA, gridVecB);
@@ -499,6 +599,21 @@ testCase.verifyEqual(permute(createTestData(linspace(0, 1, 21), linspace(0, 1, 2
value(:,:,1,2) = 1+zGridA + zetaGridA.^2;
value(:,:,2,3) = 1+zeros(numel(gridVecA), numel(gridVecB));
end
+
+%% test on on 3-dim domain
+z = quantity.Domain('z', 0:3);
+t = quantity.Domain('t', 0:4);
+x = quantity.Domain('x', 0:5);
+
+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
function testOn2(testCase)
@@ -518,30 +633,37 @@ testCase.verifyEqual([a(1).on({zeta(2), eta(3), z(4)}); a(2).on({zeta(2), eta(3)
end
function testMtimesZZeta2x2(testCase)
-gridVecA = linspace(0, 1, 101);
-[zGridA , ~] = ndgrid(gridVecA, gridVecA);
-a = quantity.Discrete(ones([size(zGridA), 2, 2]), 'size', [2, 2], ...
- 'grid', {gridVecA, gridVecA}, 'gridName', {'z', 'zeta'}, 'name', 'A');
+
+
+%%
+z11_ = linspace(0, 1, 11);
+z11 = quantity.Domain('z', z11_);
+zeta11 = quantity.Domain('zeta', z11_);
+
+a = quantity.Discrete(ones([z11.n, z11.n, 2, 2]), 'size', [2, 2], ...
+ 'domain', [z11, zeta11], 'name', 'A');
syms zeta
-assume(zeta>=0 & zeta>=1)
-gridVecB = linspace(0, 1, 100);
-b = quantity.Symbolic((eye(2, 2)), 'variable', zeta, ...
- 'grid', gridVecB);
+z10_ = linspace(0, 1, 10);
+
+b = quantity.Symbolic((eye(2, 2)), 'gridName', 'zeta', ...
+ 'grid', z10_);
+
c = a*b;
-%%
+%
testCase.verifyEqual(c.on(), a.on());
end
function testMTimesPointWise(tc)
syms z zeta
-Z = linspace(0, 1, 501)';
-ZETA = Z;
-P = quantity.Symbolic([z, z^2; z^3, z^4] * zeta, 'grid', {Z, ZETA});
-B = quantity.Symbolic([sin(zeta); cos(zeta)], 'grid', ZETA);
+Z = linspace(0, pi, 51)';
+ZETA = linspace(0, pi / 2, 71)';
+
+P = quantity.Symbolic( sin(zeta) * z, 'grid', {Z, ZETA}, 'gridName', {'z', 'zeta'});
+B = quantity.Symbolic( sin(zeta), 'grid', ZETA, 'gridName', 'zeta');
PB = P*B;
@@ -549,6 +671,7 @@ p = quantity.Discrete(P);
b = quantity.Discrete(B);
pb = p*b;
+
tc.verifyEqual(MAX(abs(PB-pb)), 0, 'AbsTol', 10*eps);
end
@@ -575,7 +698,10 @@ KLdKR = KL \ KR;
%%
testCase.verifyEqual(vLdVR.on(), 2 .\ (linspace(1, 2, 21).').^2);
-testCase.verifyEqual(KLdKR.on(), permute(repmat(2*eye(2), [1, 1, size(T)]), [4, 3, 1, 2]));
+
+% TODO #discussion: for the following case the operation KL \ KR seems to change the order of the grid in the old version. As this is very unexpected. I would like to change them to not change the order!
+testCase.verifyEqual(KLdKR.on(), permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]));
+
end
function testMrdivide(testCase)
@@ -601,7 +727,9 @@ KLdKR = KL / KR;
%%
testCase.verifyEqual(vLdVR.on(), 2 ./ (linspace(1, 2, 21).').^2);
-testCase.verifyEqual(KLdKR.on(), permute(repmat(0.5*eye(2), [1, 1, size(T)]), [4, 3, 1, 2]));
+
+% TODO #discussion: for the following case the operation KL \ KR seems to change the order of the grid in the old version. As this is very unexpected. I would like to change them to not change the order!
+testCase.verifyEqual(KLdKR.on(), permute(repmat(0.5*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]));
end
function testInv(testCase)
@@ -665,7 +793,7 @@ testCase.verifyEqual(V.on(), f.on(), 'AbsTol', 1e-5);
%% int_s_t a(t,s) * b(s) ds
v = int(a*b, sy, sy, t);
-V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s});
+V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s}, 'gridName', {'t', 's'});
f = cumInt(k * x, 's', 's', 't');
@@ -680,7 +808,7 @@ end
y = quantity.Discrete(C, 'size', size(c), 'grid', {tGrid s}, 'gridName', {'t' 's'});
v = int(a*c, sy, sy, t);
-V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s});
+V = quantity.Symbolic(subs(v, {t, sy}, {'t', 's'}), 'grid', {tGrid, s}, 'gridName', {'t', 's'});
f = cumInt( k * y, 's', 's', 't');
testCase.verifyEqual( f.on(f(1).grid, f(1).gridName), V.on(f(1).grid, f(1).gridName), 'AbsTol', 1e-5 );
@@ -714,12 +842,12 @@ tGrid = linspace(0, pi, 51)';
a = [ 1, sin(t); t, 1];
% compute the numeric version of the volterra integral
-f = quantity.Symbolic(a, 'grid', {tGrid}, 'variable', {'t'});
+f = quantity.Symbolic(a, 'grid', {tGrid}, 'gridName', {'t'});
f = quantity.Discrete(f);
intK = cumInt(f, 't', 0, 't');
-K = quantity.Symbolic( int(a, 0, t), 'grid', {tGrid}, 'variable', {'t'});
+K = quantity.Symbolic( int(a, 0, t), 'grid', {tGrid}, 'gridName', {'t'});
testCase.verifyEqual(intK.on(), K.on(), 'AbsTol', 1e-3);
@@ -753,61 +881,6 @@ testCase.verifyTrue(all(B_z_1(:,:,2) == 1));
end
-function testGridJoin(testCase)
-s = linspace(0, 1);
-z = linspace(0, 1, 71)';
-t = linspace(0, 1, 51);
-
-[Z, T, S] = ndgrid(z, t, s);
-
-a = quantity.Discrete(cat(4, sin(Z.*Z.*S), cos(Z.*T.*S)), ...
- 'size', [2 1], 'grid', {z, t, s}, 'gridName', {'z', 't', 's'});
-b = quantity.Discrete(ones(numel(s), numel(t)), ...
- 'size', [1 1], 'grid', {s, t}, 'gridName', {'p', 's'});
-c = quantity.Discrete(ones(numel(t), 2, 2), ...
- 'size', [2 2], 'grid', {t}, 'gridName', {'p'});
-
-[gridJoinedAB, gridNameJoinedAB] = gridJoin(a, b);
-[gridJoinedCC, gridNameJoinedCC] = gridJoin(c, c);
-[gridJoinedBC, gridNameJoinedBC] = gridJoin(b, c);
-
-testCase.verifyEqual(gridNameJoinedAB, {'p', 's', 't', 'z'});
-testCase.verifyEqual(gridJoinedAB, {s, s, t, z});
-testCase.verifyEqual(gridNameJoinedCC, {'p'});
-testCase.verifyEqual(gridJoinedCC, {t});
-testCase.verifyEqual(gridNameJoinedBC, {'p', 's'});
-testCase.verifyEqual(gridJoinedBC, {s, t});
-end
-
-function testGridIndex(testCase)
-
-
-z = linspace(0, 2*pi, 71)';
-t = linspace(0, 3*pi, 51);
-s = linspace(0, 1);
-
-[Z, T, S] = ndgrid(z, t, s);
-
-a = quantity.Discrete(cat(4, sin(Z.*Z.*S), cos(Z.*T.*S)), ...
- 'size', [2 1], 'grid', {z, t, s}, 'gridName', {'z', 't', 's'});
-
-idx = a.gridIndex('z');
-
-testCase.verifyEqual(idx, 1);
-
-
-idx = a.gridIndex({'z', 't'});
-testCase.verifyEqual(idx, [1 2]);
-
-idx = a.gridIndex({'z', 's'});
-testCase.verifyEqual(idx, [1 3]);
-
-idx = a.gridIndex('t');
-testCase.verifyEqual(idx, 2);
-
-end
-
-
function testMTimes(testCase)
% multiplication of a(z) * a(z)'
% a(z) \in (3, 2), z \in (100)
@@ -836,7 +909,7 @@ testCase.verifyTrue(numeric.near(b1b1', B1B1.on()));
% vector multiplication of b(z, t) * b(v, t)
b2 = cat(3, z .* t, z*0 +1);
b2b2 = misc.multArray(b2, permute(b2, [1 2 4 3]), 4, 3, [1, 2]);
-b2b2 = permute(b2b2, [2 1 3 4]);
+b2b2 = permute(b2b2, [1 2 3 4]);
B2 = quantity.Discrete(b2, 'size', [ 2, 1 ], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
B2B2 = B2*B2';
testCase.verifyTrue(numeric.near(b2b2, B2B2.on()));
@@ -845,12 +918,12 @@ c2 = cat(4, sin(z), sin(t), cos(z));
b2c2 = misc.multArray(b2, c2, 4, 3, [1 2]);
C2 = quantity.Discrete(c2, 'size', [1 3], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
B2C2 = B2 * C2;
-testCase.verifyTrue(numeric.near(permute(b2c2, [2 1 3 4]), B2C2.on()));
+testCase.verifyTrue(numeric.near(permute(b2c2, [1 2 3 4]), B2C2.on()));
% matrix b(z,t) * b(z,t)
b = cat(4, cat(3, z .* t, z*0 +1), cat(3, sin( z ), cos( z .* t )));
bb = misc.multArray(b, b, 4, 3, [1, 2]);
-bb = permute(bb, [2 1 3 4]);
+bb = permute(bb, [1 2 3 4]);
B = quantity.Discrete(b, 'size', [ 2, 2 ], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
BB = B*B;
testCase.verifyTrue(numeric.near(bb, BB.on));
@@ -859,8 +932,8 @@ testCase.verifyTrue(numeric.near(bb, BB.on));
c = cat(4, cat(3, v .* t2, v*0 +1), cat(3, sin( v ), cos( v .* t2 )), cat(3, cos( t2 ), tan( v .* t2 )));
bc = misc.multArray(b, c, 4, 3, 2);
bc = permute(bc, [ 1 2 4 3 5 ]);
-B = quantity.Discrete(b, 'size', [ 2, 2 ], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
-C = quantity.Discrete(c, 'size', [2, 3], 'gridName', {'v', 't'}, 'grid',{v(:,1), t2(1,:)});
+B = quantity.Discrete(b, 'size', [2, 2], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
+C = quantity.Discrete(c, 'size', [2, 3], 'gridName', {'v', 't'}, 'grid', {v(:,1), t2(1,:)});
BC = B*C;
testCase.verifyTrue(numeric.near(bc, BC.on));
@@ -876,7 +949,7 @@ testCase.verifyTrue( numeric.near( aa(2,2).on() , cos(z * pi) .* cos(z * pi)));
%% test multiplicatin with constants:
C = [3 0; 0 5];
-c = quantity.Discrete(C, 'gridName', {}, 'grid', {}, 'name', 'c');
+c = quantity.Discrete(C, 'domain', quantity.Domain.empty(), 'name', 'c');
testCase.verifyTrue( numeric.near( squeeze(double(a*c)), [sin(z * pi) * 3, cos(z * pi) * 5]));
testCase.verifyTrue( numeric.near( squeeze(double(a*[3 0; 0 5])), [sin(z * pi) * 3, cos(z * pi) * 5]));
testCase.verifyTrue( numeric.near( double(c*c), C*C));
@@ -888,16 +961,23 @@ testCase.verifyTrue( numeric.near( squeeze(double(a * 42)), [sin(z * pi) * 42, c
testCase.verifyTrue( numeric.near( squeeze(double(s * a')), [sin(z * pi) * 42, cos(z * pi) * 42]));
testCase.verifyTrue( numeric.near( squeeze(double(a * s)), [sin(z * pi) * 42, cos(z * pi) * 42]));
+end
-%% test
+function testOnConstant(testCase)
+A = [0 1; 2 3];
+const = quantity.Discrete(A, 'domain', quantity.Domain.empty);
+C = const.on([1:5]);
+for k = 1:numel(A)
+ testCase.verifyTrue( all( A(k) == C(:,k) ));
+end
end
function testIsConstant(testCase)
-C = quantity.Discrete(rand(3,7), 'gridName', {});
+C = quantity.Discrete(rand(3,7), 'domain', quantity.Domain.empty());
testCase.verifyTrue(C.isConstant());
z = linspace(0, pi)';
@@ -913,7 +993,7 @@ a = cat(3, [myGrid*0, myGrid, sin(myGrid)], [myGrid.^2, myGrid.^3, cos(myGrid)])
c = rand(2, 3); % dim = (2, 3)
A = quantity.Discrete(a, 'size', [3, 2], 'grid', myGrid, 'gridName', {'z'});
-C = quantity.Discrete(c, 'gridName', {});
+C = quantity.Discrete(c, 'domain', quantity.Domain.empty());
ac = misc.multArray(a, c, 3, 1); % dim = (100, 3, 3)
AC = A*C;
@@ -993,7 +1073,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());
@@ -1009,6 +1088,26 @@ cAB2 = [1 2; 3 4] + AB2;
testCase.verifyEqual(AB2c.on(), tst)
testCase.verifyEqual(cAB2.on(), tst)
+%% test plus on different domains
+z = quantity.Domain('z', 0:3);
+t = quantity.Domain('t', 0:4);
+x = quantity.Domain('x', 0:5);
+
+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;
+testCase.verifyEqual( on( TZX + XTZ ), on( 2 * TZX ), 'AbsTol', eps )
+
end
function testInit(testCase)
@@ -1019,7 +1118,11 @@ v = {sin(z * t * pi); cos(z * t * pi)};
V = cat(3, v{:});
b = quantity.Discrete(v, 'grid', {z, t}, 'gridName', {'z', 't'});
c = quantity.Discrete(V, 'grid', {z, t}, 'gridName', {'z', 't'});
-d = quantity.Discrete(V, 'gridName', {'z', 't'});
+
+Z = quantity.Domain('z', z);
+T = quantity.Domain('t', t);
+
+d = quantity.Discrete(V, 'domain', [Z, T]);
%%
verifyTrue(testCase, misc.alln(b.on() == c.on()));
@@ -1071,23 +1174,12 @@ verifyTrue(testCase, numeric.near(A, Anum, 1e-2));
end
-% function testNDGrid(testCase)
-% %%
-% z = linspace(0,1).';
-% t = linspace(0,1,101);
-% b = quantity.Discrete({sin(z * t * pi); cos(z * t * pi)}, 'grid', {z, t}, 'gridName', {'z', 't'});
-% % #TODO
-% end
-
-function testDefaultGrid(testCase)
-g = quantity.Discrete.defaultGrid([100, 42]);
-testCase.verifyEqual(g{1}, linspace(0, 1, 100).');
-testCase.verifyEqual(g{2}, linspace(0, 1, 42));
-end
-
function testValue2Cell(testCase)
v = rand([100, 42, 2, 3]);
-V = quantity.Discrete( quantity.Discrete.value2cell(v, [100, 42], [2, 3]), 'gridName', {'z', 't'} );
+Z = quantity.Domain('z', linspace(0,1,100));
+T = quantity.Domain('t', linspace(0,1,42));
+V = quantity.Discrete( quantity.Discrete.value2cell(v, [100, 42], [2, 3]), ...
+ 'domain', [Z, T]);
verifyTrue(testCase, misc.alln(v == V.on()));
end
@@ -1139,6 +1231,7 @@ t = linspace(0, 1, 31);
quan = quantity.Discrete({sin(z * t * pi); cos(z * t * pi)}, 'grid', {z, t}, 'gridName', {'z', 't'});
gridSampled = {linspace(0, 1, 11), linspace(0, 1, 21)};
quanCopy = copy(quan);
+
quanSampled = quanCopy.changeGrid(gridSampled, {'z', 't'});
testCase.verifyEqual(quanSampled.on(), quan.on(gridSampled))
@@ -1147,14 +1240,16 @@ quanSampled2 = quanCopy2.changeGrid(gridSampled, {'t', 'z'});
testCase.verifyEqual(quanSampled2.on(), permute(quan.on(gridSampled, {'t', 'z'}), [2, 1, 3]));
end
-
-
function testSubs(testCase)
%%
myTestArray = ones(21, 41, 2, 3);
myTestArray(:,1,1,1) = linspace(0, 1, 21);
myTestArray(:,:,2,:) = 2;
-quan = quantity.Discrete(myTestArray, 'gridName', {'z', 'zeta'}, 'size', [2, 3]);
+z = quantity.Domain('z', linspace(0, 1, 21));
+zeta = quantity.Domain('zeta', linspace(0, 1, 41));
+
+quan = quantity.Discrete(myTestArray, 'domain', [z, zeta]);
+
quanZetaZ = quan.subs({'z', 'zeta'}, {'zeta', 'z'});
quan10 = quan.subs({'z', 'zeta'}, {1, 0});
quan01 = quan.subs({'z', 'zeta'}, {0, 1});
@@ -1162,7 +1257,7 @@ quant1 = quan.subs({'z', 'zeta'}, {'t', 1});
quanZetaZeta = quan.subs({'z'}, 'zeta');
quanPt = quan.subs({'z'}, {1});
quanEta = quan.subs({'z'}, {'eta'});
-%%
+%
testCase.verifyEqual(quan10, shiftdim(quan.on({1, 0})));
testCase.verifyEqual(quan01, shiftdim(quan.on({0, 1})));
testCase.verifyEqual(quanZetaZ.on(), quan.on());
@@ -1175,9 +1270,17 @@ testCase.verifyEqual(quanZetaZeta.on(), quanZetaZetaReference);
testCase.verifyEqual(quanEta(1).gridName, {'eta', 'zeta'})
testCase.verifyEqual(quanPt.on(), shiftdim(quan.on({1, quan(1).grid{2}})));
-%% 4d-case
+testCase.verifyEqual( gridSize( quanZetaZeta.subs( ...
+ quantity.Domain('zeta', linspace(0,1,3))) ), 3);
+
+% 4d-case
myTestArray4d = rand(11, 21, 31, 41, 1, 2);
-quan4d = quantity.Discrete(myTestArray4d, 'gridName', {'z', 'zeta', 'eta', 'beta'}, 'name', 'fun5d');
+x1 = quantity.Domain('z', linspace(0, 1, 11));
+x2 = quantity.Domain('zeta', linspace(0, 1, 21));
+x3 = quantity.Domain('eta', linspace(0, 1, 31));
+x4 = quantity.Domain('beta', linspace(0, 1, 41));
+
+quan4d = quantity.Discrete(myTestArray4d, 'domain', [x1, x2, x3, x4], 'name', 'fun5d');
quan4dAllEta = quan4d.subs({'z', 'zeta', 'eta', 'beta'}, {'eta', 'eta', 'eta', 'eta'});
testCase.verifyEqual(reshape(quan4d.on({1, 1, 1, 1}), size(quan4d)), ...
reshape(quan4dAllEta.on({1}), size(quan4dAllEta)));
@@ -1185,6 +1288,8 @@ testCase.verifyEqual(reshape(quan4d.on({1, 1, 1, 1}), size(quan4d)), ...
quan4dArbitrary = quan4d.subs({'z', 'zeta', 'eta', 'beta'}, {'zeta', 'beta', 'z', 1});
testCase.verifyEqual(reshape(quan4d.on({1, 1, 1, 1}), size(quan4d)), ...
reshape(quan4dArbitrary.on({1, 1, 1}), size(quan4dAllEta)));
+
+
end
function testZTTimes(testCase)
diff --git a/+unittests/+quantity/testDomain.m b/+unittests/+quantity/testDomain.m
new file mode 100644
index 0000000000000000000000000000000000000000..fae2acc0d35e4fdcdf159f778c6dcc56d31de900
--- /dev/null
+++ b/+unittests/+quantity/testDomain.m
@@ -0,0 +1,180 @@
+function [tests] = testDomain()
+ tests = functiontests(localfunctions);
+end
+
+function setupOnce(testCase)
+end
+
+function testEquality(testCase)
+
+d = quantity.Domain('d', 0);
+e = quantity.Domain.empty();
+
+dd = d;
+ee = e;
+
+testCase.verifyNotEqual( d, e );
+testCase.verifyEqual( dd, d );
+testCase.verifyEqual( ee, e );
+testCase.verifyNotEqual( e, d );
+
+end
+
+function testDomainInit(testCase)
+
+Z = linspace(0,pi, 3);
+d = quantity.Domain('z', Z);
+D = [d d];
+
+testCase.verifyEqual( ndims(D), 2);
+testCase.verifyEqual( ndims(d), 1);
+end
+
+function testSort(testCase)
+
+z = linspace(0, pi, 3);
+a = quantity.Domain('a', z);
+b = quantity.Domain('b', z);
+c = quantity.Domain('c', z);
+
+d = [a b c];
+
+[d_, I_] = d.sort('descend');
+[d__, I__] = d.sort('ascend');
+d___ = d.sort();
+
+testCase.verifyEqual({d_.name}, {'c', 'b', 'a'});
+testCase.verifyEqual({d__.name}, {d.name});
+testCase.verifyEqual({d__.name}, {d___.name});
+testCase.verifyEqual( I_, flip(I__) );
+
+
+end
+
+function testDomainGridLength(testCase)
+
+Z = linspace(0,pi, 3);
+d = quantity.Domain('z', Z);
+D = [d d];
+
+testCase.verifyEqual( cellfun(@(v) numel(v), {Z, Z}), D.gridLength)
+
+end
+
+function testDomainNumGridElements(testCase)
+
+Z = linspace(0,pi, 3);
+d = quantity.Domain('z', Z);
+D = [d d];
+
+testCase.verifyEqual( D.numGridElements, prod([length(Z), length(Z)]));
+
+end
+
+function testEmpty(testCase)
+
+d = quantity.Domain();
+e = quantity.Domain.empty();
+o = quantity.Domain('', 1);
+
+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
+
+
+function testDefaultGrid(testCase)
+ g = quantity.Domain.defaultGrid([100, 42]);
+ testCase.verifyEqual(g(1).grid, linspace(0, 1, 100).');
+ testCase.verifyEqual(g(2).grid, linspace(0, 1, 42).');
+end
+
+function testGridJoin(testCase)
+s = linspace(0, 1);
+z = linspace(0, 1, 71)';
+t = linspace(0, 1, 51);
+
+S = quantity.Domain('s', s);
+Z = quantity.Domain('z', z);
+T = quantity.Domain('t', t);
+Ps = quantity.Domain('p', s);
+St = quantity.Domain('s', t);
+Pt = quantity.Domain('p', t);
+
+a = [Z, T, S];
+b = [Ps, St];
+c = Pt;
+%
+%
+% a = quantity.Discrete(cat(4, sin(Z.*Z.*S), cos(Z.*T.*S)), ...
+% 'size', [2 1], 'grid', {z, t, s}, 'gridName', {'z', 't', 's'});
+% b = quantity.Discrete(ones(numel(s), numel(t)), ...
+% 'size', [1 1], 'grid', {s, t}, 'gridName', {'p', 's'});
+% c = quantity.Discrete(ones(numel(t), 2, 2), ...
+% 'size', [2 2], 'grid', {t}, 'gridName', {'p'});
+
+joinedDomainAB = sort( gridJoin(a, b) );
+joinedDomainCC = sort( gridJoin(c, c) );
+joinedDomainBC = sort( gridJoin(b, c) );
+joinedDomain = sort( gridJoin(Z, [Z, T]) );
+
+testCase.verifyEqual({joinedDomainAB.name}, {'p', 's', 't', 'z'});
+testCase.verifyEqual({joinedDomainAB.grid}, {s(:), s(:), t(:), z(:)});
+testCase.verifyEqual({joinedDomainCC.name}, {'p'});
+testCase.verifyEqual({joinedDomainCC.grid}, {t(:)});
+testCase.verifyEqual({joinedDomainBC.name}, {'p', 's'});
+testCase.verifyEqual({joinedDomainBC.grid}, {s(:), t(:)});
+testCase.verifyEqual({joinedDomain.grid}, {t(:), z(:)});
+end
+
+function testGridIndex(testCase)
+z = linspace(0, 2*pi, 71)';
+t = linspace(0, 3*pi, 51);
+s = linspace(0, 1);
+
+S = quantity.Domain('s', s);
+Z = quantity.Domain('z', z);
+T = quantity.Domain('t', t);
+
+a = [Z, T, S];
+
+idx = a.gridIndex('z');
+
+testCase.verifyEqual(idx, 1);
+
+
+idx = a.gridIndex({'z', 't'});
+testCase.verifyEqual(idx, [1 2]);
+
+idx = a.gridIndex({'z', 's'});
+testCase.verifyEqual(idx, [1 3]);
+
+idx = a.gridIndex('t');
+testCase.verifyEqual(idx, 2);
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/+unittests/+quantity/testOperator.m b/+unittests/+quantity/testOperator.m
index 8e853d2226aef03e5e817925914111e0eb4cc749..81e8acfcab8fbeca9fa04ff95391d313cf169a10 100644
--- a/+unittests/+quantity/testOperator.m
+++ b/+unittests/+quantity/testOperator.m
@@ -27,7 +27,7 @@ A = cat(3, [ 0, 0, 0, 0; ...
B = cat(3, [1, 0, 0, 0]', zeros(4, 1), zeros(4, 1));
for k = 1:3
- a{k} = quantity.Symbolic(A(:,:,k), 'grid', Z, 'variable', z);
+ a{k} = quantity.Symbolic(A(:,:,k), 'grid', Z, 'gridName', 'z');
end
A = quantity.Operator(a, 's', s);
@@ -37,6 +37,31 @@ testCase.TestData.a = a;
end
+function testApplyTo(testCase)
+
+ A = [0 1 0; 0 0 1; - (1:3)];
+ B = [0 0 1]';
+ C = [1 0 0];
+
+ sys = ss(A, B, C, []);
+ detsIA = quantity.Operator(num2cell(flip(charpoly(A))));
+ adjB = quantity.Operator( adjoint(sym('s')*eye(size(A)) - A) * B );
+
+ t = sym('t', 'real');
+ sol.y = quantity.Symbolic( t.^3 .* (1 - t).^3, ...
+ 'grid', linspace(0, 1, 1e2), 'gridName', 't');
+
+ sol.u = detsIA.applyTo(sol.y);
+ sol.x = adjB.applyTo(sol.y);
+
+ [sim.y, sim.t, sim.x] = lsim(sys, sol.u.on(), sol.y.domain.grid, 'foh');
+
+ % deviation of final values:
+ dev.y = sim.y(end)
+ dev.x = sim.x(end, :)
+
+end
+
function testCTranspose(testCase)
At = testCase.TestData.A';
testCase.verifyTrue(At(3).coefficient == testCase.TestData.A(3).coefficient');
@@ -47,9 +72,9 @@ function testFundamentalMatrixSpaceDependent(testCase)
%
a = 20;
z = sym('z', 'real');
- Z = linspace(0, 1, 32)';
+ Z = linspace(0, 1, 5)';
A = quantity.Operator(...
- {quantity.Symbolic([1, z; 0, a], 'variable', 'z', 'grid', Z)});
+ {quantity.Symbolic([1, z; 0, a], 'gridName', 'z', 'grid', Z)});
F = A.stateTransitionMatrix();
@@ -60,7 +85,7 @@ function testFundamentalMatrixSpaceDependent(testCase)
f = (exp(z) - exp(a * z) - (1 - a) * z * exp( a * z ) ) / (1 - a)^2;
end
- PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
+ PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'gridName', 'z', 'name', 'Phi_A');
testCase.verifyEqual(double(F), PhiExact.on(), 'RelTol', 1e-6);
@@ -69,10 +94,10 @@ end
function testStateTransitionMatrix(testCase)
N = 2;
Z = linspace(0,1,11)';
-A0 = quantity.Symbolic([0 1; 1 0], 'grid', Z, 'variable', 'z');
-A1 = quantity.Symbolic([0 1; 0 0], 'grid', Z, 'variable', 'z');
+A0 = quantity.Symbolic([0 1; 1 0], 'grid', Z, 'gridName', 'z');
+A1 = quantity.Symbolic([0 1; 0 0], 'grid', Z, 'gridName', 'z');
A = quantity.Operator({A0, A1});
-B = quantity.Operator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'variable', 'z'));
+B = quantity.Operator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'gridName', 'z'));
[Phi1, Psi1] = A.stateTransitionMatrix(N, B);
[Phi2, Psi2] = A.stateTransitionMatrixByOdeSystem(N, B);
@@ -97,7 +122,7 @@ end
function testMTimes(testCase)
- z = linspace(0, pi, 15)';
+ z = linspace(0, pi, 3)';
A0 = quantity.Discrete(reshape((repmat([1 2; 3 4], length(z), 1)), length(z), 2, 2),...
'gridName', 'z', 'grid', z, 'name', 'A');
A1 = quantity.Discrete(reshape((repmat([5 6; 7 8], length(z), 1)), length(z), 2, 2), ...
@@ -115,7 +140,7 @@ function testMTimes(testCase)
As = A0.at(0) + A1.at(0)*s + A2.at(0)*s^2;
Bs = A0.at(0);
- Cs = misc.polynomial2coefficients(As*Bs, s);
+ Cs = polyMatrix.polynomial2coefficients(As*Bs, s);
for k = 1:3
testCase.verifyTrue( numeric.near(Cs(:,:,k), C(k).coefficient.at(0)) );
@@ -133,7 +158,7 @@ function testMTimes(testCase)
end
function testSum(testCase)
-z = linspace(0, pi, 15)';
+z = linspace(0, pi, 5)';
A0 = quantity.Discrete(reshape((repmat([1 2; 3 4], length(z), 1)), length(z), 2, 2),...
'gridName', 'z', 'grid', z, 'name', 'A');
A1 = quantity.Discrete(reshape((repmat([5 6; 7 8], length(z), 1)), length(z), 2, 2), ...
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 5634c445b328a6105b88f0d24e20ed1299677aa2..9d58191833f10fc15862123def7719d015923cfb 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -7,8 +7,11 @@ end
function testVec2Diag(testCase)
% quantity.Symbolic
syms z
-myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2]);
-myVectorSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1; 0.9-z/2]);
+Z = quantity.Domain('z', linspace(0, 1, 3));
+myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], ...
+ 'domain', Z);
+myVectorSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1; 0.9-z/2], ...
+ 'domain', Z);
testCase.verifyTrue( myVectorSymbolic.vec2diag.near(myMatrixSymbolic) );
end
@@ -17,9 +20,9 @@ function testSymbolicEvaluation(tc)
syms z
myGrid = linspace(0, 1, 7);
fS = quantity.Symbolic(z * sinh(z * 1e5) / cosh(z * 1e5), ...
- 'variable', {z}, 'grid', {myGrid}, 'symbolicEvaluation', true);
+ 'gridName', {'z'}, 'grid', {myGrid}, 'symbolicEvaluation', true);
fF = quantity.Symbolic(z * sinh(z * 1e5) / cosh(z * 1e5), ...
- 'variable', {z}, 'grid', {myGrid}, 'symbolicEvaluation', false);
+ 'gridName', {'z'}, 'grid', {myGrid}, 'symbolicEvaluation', false);
tc.verifyTrue(any(isnan(fF.on)))
tc.verifyFalse(any(isnan(fS.on)))
@@ -29,7 +32,7 @@ function testCtranspose(tc)
syms z zeta
qSymbolic = quantity.Symbolic(...
[1+zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
+ 'gridName', {'z', 'zeta'}, 'name', 'q');
qSymbolicCtransp = qSymbolic';
tc.verifyEqual(qSymbolic(1,1).on(), qSymbolicCtransp(1,1).on());
tc.verifyEqual(qSymbolic(2,2).on(), qSymbolicCtransp(2,2).on());
@@ -37,7 +40,7 @@ tc.verifyEqual(qSymbolic(1,2).on(), qSymbolicCtransp(2,1).on());
tc.verifyEqual(qSymbolic(2,1).on(), qSymbolicCtransp(1,2).on());
qSymbolic2 = quantity.Symbolic(sym('z') * 2i + sym('zeta'), ...
- 'grid', {linspace(0, 1, 21), linspace(0, 1, 11)}, 'variable', sym({'z', 'zeta'}), 'name', 'im');
+ 'grid', {linspace(0, 1, 21), linspace(0, 1, 11)}, 'gridName', {'z', 'zeta'}, 'name', 'im');
qSymolic2Ctransp = qSymbolic2';
tc.verifyEqual(qSymbolic2.real.on(), qSymolic2Ctransp.real.on());
tc.verifyEqual(qSymbolic2.imag.on(), -qSymolic2Ctransp.imag.on());
@@ -47,7 +50,7 @@ function testTranspose(tc)
syms z zeta
qSymbolic = quantity.Symbolic(...
[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 'q');
+ 'gridName', {'z', 'zeta'}, 'name', 'q');
qSymbolicTransp = qSymbolic.';
tc.verifyEqual(qSymbolic(1,1).on(), qSymbolicTransp(1,1).on());
tc.verifyEqual(qSymbolic(2,2).on(), qSymbolicTransp(2,2).on());
@@ -59,16 +62,16 @@ function testFlipGrid(tc)
syms z zeta
myGrid = linspace(0, 1, 11);
f = quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
% flip zeta
fReference = quantity.Symbolic([1+z+(1-zeta); 2*(1-zeta)+sin(z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
fFlipped = f.flipGrid('zeta');
tc.verifyEqual(fReference.on, fFlipped.on, 'AbsTol', 10*eps);
% flip z and zeta
fReference2 = quantity.Symbolic([1+(1-z)+(1-zeta); 2*(1-zeta)+sin(1-z)] + zeros(2, 1)*z*zeta, ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
fFlipped2 = f.flipGrid({'z', 'zeta'});
tc.verifyEqual(fReference2.on, fFlipped2.on, 'AbsTol', 10*eps);
end % testFlipGrid();
@@ -76,14 +79,15 @@ end % testFlipGrid();
function testCumInt(testCase)
tGrid = linspace(pi, 1.1*pi, 51)';
sGrid = tGrid;
-syms s t
+s = sym('s');
+t = sym('t');
a = [ 1, s; t, 1];
b = [ s; 2*s];
%% int_0_t a(t,s) * b(s) ds
% compute symbolic version of the volterra integral
-integrandSymbolic = quantity.Symbolic(a*b, 'grid', {tGrid, sGrid}, 'variable', {t, s});
+integrandSymbolic = quantity.Symbolic(a*b, 'grid', {tGrid, sGrid}, 'gridName', {'t', 's'});
integrandDiscrete = quantity.Discrete(integrandSymbolic);
V = cumInt(integrandSymbolic, 's', tGrid(1), 't');
f = cumInt(integrandDiscrete, 's', sGrid(1), 't');
@@ -100,6 +104,7 @@ testCase.verifyEqual( f.on(f(1).grid, f(1).gridName), V.on(f(1).grid, f(1).gridN
fCumInt2Bcs = cumInt(integrandSymbolic, 's', 'zeta', 't');
fCumInt2Cum = cumInt(integrandSymbolic, 's', tGrid(1), 't') ...
- cumInt(integrandSymbolic, 's', tGrid(1), 'zeta');
+
testCase.verifyEqual(fCumInt2Bcs.on(fCumInt2Bcs(1).grid, fCumInt2Bcs(1).gridName), ...
fCumInt2Cum.on(fCumInt2Bcs(1).grid, fCumInt2Bcs(1).gridName), 'AbsTol', 100*eps);
end
@@ -108,7 +113,7 @@ function testScalarPlusMinusQuantity(testCase)
syms z
myGrid = linspace(0, 1, 7);
f = quantity.Symbolic([1; 2] + zeros(2, 1)*z, ...
- 'variable', {z}, 'grid', {myGrid});
+ 'gridName', {'z'}, 'grid', {myGrid});
testCase.verifyError(@() 1-f-1, '');
testCase.verifyError(@() 1+f+1, '');
end
@@ -117,7 +122,7 @@ function testNumericVectorPlusMinusQuantity(testCase)
syms z
myGrid = linspace(0, 1, 7);
f = quantity.Symbolic([1+z; 2+sin(z)] + zeros(2, 1)*z, ...
- 'variable', {z}, 'grid', {myGrid});
+ 'gridName', {'z'}, 'grid', {myGrid});
a = ones(size(f));
testCase.verifyEqual(on(a-f), 1-f.on(), 'RelTol', 10*eps);
testCase.verifyEqual(on(f-a), f.on()-1, 'RelTol', 10*eps);
@@ -129,11 +134,11 @@ function testDiffWith2Variables(testCase)
syms z zeta
myGrid = linspace(0, 1, 7);
f = quantity.Symbolic([z*zeta, z; zeta, 1], ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
fdz = quantity.Symbolic([zeta, 1; 0, 0], ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
fdzeta = quantity.Symbolic([z, 0; 1, 0], ...
- 'variable', {z, zeta}, 'grid', {myGrid, myGrid});
+ 'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
fRefTotal = 0*f.on();
fRefTotal(:,:,1,1) = 1;
@@ -149,7 +154,7 @@ function testGridName2variable(testCase)
syms z zeta eta e a
myGrid = linspace(0, 1, 7);
obj = quantity.Symbolic([z*zeta, eta*z; e*a, a], ...
- 'variable', {z zeta eta e a}, ...
+ 'gridName', {'z' 'zeta' 'eta' 'e' 'a'}, ...
'grid', {myGrid, myGrid, myGrid, myGrid, myGrid});
thisGridName = {'eta', 'e', 'z'};
@@ -168,8 +173,8 @@ end
function testCat(testCase)
syms z zeta
-f1 = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 'f1');
-f2 = quantity.Symbolic(sin(z), 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 'f2');
+f1 = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'gridName', {'z'}, 'name', 'f1');
+f2 = quantity.Symbolic(sin(z), 'grid', {linspace(0, 1, 21)}, 'gridName', {'z'}, 'name', 'f2');
% vertical concatenation
F = [f1; f2];
@@ -189,7 +194,7 @@ testCase.verifyTrue(F(2,1) == f1);
testCase.verifyTrue(F(2,2) == f2);
% concatenation on different grids
-f3 = quantity.Symbolic(cos(z), 'grid', {linspace(0, 1, 13)}, 'variable', {z}, 'name', 'f1');
+f3 = quantity.Symbolic(cos(z), 'grid', {linspace(0, 1, 13)}, 'gridName', {'z'}, 'name', 'f1');
F = [f1, f3];
testCase.verifyEqual(F(2).gridSize, f1.gridSize)
@@ -208,19 +213,19 @@ end
function testExp(testCase)
% 1 spatial variable
syms z zeta
-s1d = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 's1d');
+s1d = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'gridName', {'z'}, 'name', 's1d');
testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));
% diagonal matrix
s2dDiag = quantity.Symbolic([1+z*zeta, 0; 0, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2dDiag');
+ 'gridName', {'z', 'zeta'}, 'name', 's2dDiag');
testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
end
function testExpm(testCase)
syms z zeta
mat2d = quantity.Symbolic(ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2d');
+ 'gridName', {'z', 'zeta'}, 'name', 's2d');
mat2dMat = mat2d.on();
mat2dExpm = 0 * mat2d.on();
for zIdx = 1 : 21
@@ -272,7 +277,7 @@ function testOn(testCase)
syms z zeta
gridZ = linspace(0, 1, 40);
gridZeta = linspace(0, 1, 21);
-A = quantity.Symbolic(2+[z, z^2, z+zeta; -sin(zeta), z*zeta, 1], 'variable', {z, zeta}, ...
+A = quantity.Symbolic(2+[z, z^2, z+zeta; -sin(zeta), z*zeta, 1], 'gridName', {'z', 'zeta'}, ...
'grid', {gridZ, gridZeta}, 'name', 'A');
%%
@@ -305,17 +310,17 @@ end
function testSqrt(testCase)
% 1 spatial variable
syms z zeta
-s1d = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 's1d');
+s1d = quantity.Symbolic(1+z*z, 'grid', {linspace(0, 1, 21)}, 'gridName', {'z'}, 'name', 's1d');
testCase.verifyEqual(s1d.sqrt.on(), sqrt(s1d.on()));
% 2 spatial variables
s2d = quantity.Symbolic(1+z*zeta, 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2d');
+ 'gridName', {'z', 'zeta'}, 'name', 's2d');
testCase.verifyEqual(s2d.sqrt.on(), sqrt(s2d.on()));
% diagonal matrix
s2dDiag = quantity.Symbolic([1+z*zeta, 0; 0, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2dDiag');
+ 'gridName', {'z', 'zeta'}, 'name', 's2dDiag');
testCase.verifyEqual(s2dDiag.sqrt.on(), sqrt(s2dDiag.on()));
end
@@ -323,7 +328,7 @@ end
function testSqrtm(testCase)
syms z zeta
mat2d = quantity.Symbolic(ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
- 'variable', {z, zeta}, 'name', 's2d');
+ 'gridName', {'z', 'zeta'}, 'name', 's2d');
mat2dMat = mat2d.on();
mat2dSqrtm = 0 * mat2d.on();
for zIdx = 1 : 21
@@ -342,7 +347,7 @@ function testInv(testCase)
syms z
zGrid = linspace(0, 1, 21);
v = quantity.Symbolic(1+z*z, ...
- 'grid', {zGrid}, 'variable', {z}, 'name', 's1d');
+ 'grid', {zGrid}, 'gridName', {'z'}, 'name', 's1d');
vInvReference = 1 ./ (1 + zGrid.^2);
vInv = v.inv();
@@ -365,7 +370,7 @@ end
function testZero(testCase)
-x = quantity.Symbolic(0);
+x = quantity.Symbolic(0, 'domain', quantity.Domain.empty());
verifyEqual(testCase, x.on(), 0);
verifyTrue(testCase, misc.alln(on(x) * magic(5) == zeros(5)));
@@ -375,7 +380,10 @@ end
function testCast2Quantity(testCase)
syms z t
-x = quantity.Symbolic(sin(z * t * pi), 'grid', {linspace(0, 1).', linspace(0, 2, 200)});
+d = [quantity.Domain('z', linspace(0, 1, 5)), ...
+ quantity.Domain('t', linspace(0, 2, 2))];
+
+x = quantity.Symbolic(sin(z * t * pi), 'domain', d);
X = quantity.Discrete(x);
@@ -384,7 +392,9 @@ end
function testMrdivide(testCase)
syms z t
-x = quantity.Symbolic(sin(z * t * pi), 'grid', {linspace(0, 1).', linspace(0, 2, 200)});
+d = [quantity.Domain('z', linspace(0, 1, 5)), ...
+ quantity.Domain('t', linspace(0, 2, 2))];
+x = quantity.Symbolic(sin(z * t * pi), 'domain', d);
o = x / x;
verifyTrue(testCase, misc.alln(o.on == 1));
@@ -398,7 +408,9 @@ end
function testMTimesConstants(testCase)
syms z t
-x = quantity.Symbolic(sin(z * t * pi), 'grid', {linspace(0, 1).', linspace(0, 2, 200)});
+d = [quantity.Domain('z', linspace(0, 1, 5)), ...
+ quantity.Domain('t', linspace(0, 2, 2))];
+x = quantity.Symbolic(sin(z * t * pi), 'domain', d);
% test multiplicaiton iwth a constant scalar
x5 = x * 5;
@@ -424,9 +436,9 @@ function testMLRdivide2ConstantSymbolic(testCase)
syms z
assume(z>0 & z<1);
-Lambda = quantity.Symbolic(eye(2), 'variable', z, ...
+Lambda = quantity.Symbolic(eye(2), ...
'grid', linspace(0, 1), 'gridName', 'z', 'name', 'Lambda');
-A = quantity.Symbolic(eye(2), 'variable', z, ...
+A = quantity.Symbolic(eye(2), ...
'grid', linspace(0, 1), 'gridName', 'z', 'name', 'A');
C1 = A / Lambda;
C2 = A \ Lambda;
@@ -451,7 +463,9 @@ end
function testCastToQuantity(testCase)
syms z t
-x = quantity.Symbolic(sin(z * t * pi), 'grid', {linspace(0, 1).', linspace(0, 2, 200)});
+x = quantity.Symbolic(sin(z * t * pi), ...
+ 'grid', {linspace(0, 1).', linspace(0, 2, 200)}, ...
+ 'gridName', {'z', 't'});
X1 = quantity.Discrete(x);
@@ -460,25 +474,14 @@ verifyEqual(testCase, X1.on, x.on);
end
-% function testmTimesFunctionHandleSymbolic(testCase)
-%
-% TODO repair this kind of multiplication!
-%
-% %
-% f = quantity.Function(@(z) z, 'name', 'f');
-% s = quantity.Symbolic(sym('z'), 'name', 'f');
-% verifyTrue(testCase, isa(f*s, 'quantity.Function'));
-%
-% end
-
-
function testDiff(testCase)
%%
z = sym('z', 'real');
f1 = sinh(z * pi);
f2 = cosh(z * pi);
-f = quantity.Symbolic([f1 ; f2 ], 'name', 'f');
+f = quantity.Symbolic([f1 ; f2 ], 'name', 'f', 'domain', ...
+ quantity.Domain('z', linspace(0,1,11)));
d2f = f.diff(2);
d1f = f.diff(1);
@@ -502,23 +505,37 @@ verifyTrue(testCase, numeric.near(double([R3{:}]), f.diff(2).on(), 1e-12));
end
-
function testInit(testCase)
-syms z
-A = [0 1; 1 0];
+z = linspace(0,1).';
+t = linspace(0,1,3);
+
+syms x y
-%% compute comparative solution
-P1 = expm(A * z);
+v = [sin(x * y * pi); cos(x * y * pi)];
+
+try
+ b = quantity.Symbolic(v, 'grid', {z, t}, 'gridName', {'z', 't'});
+ testCase.verifTrue(false);
+end
+
+b = quantity.Symbolic(v, 'grid', {z, t}, 'gridName', {'x', 'y'});
+c = quantity.Symbolic(v, 'grid', {z, t}, 'gridName', {'x', 'y'});
+%%
+verifyTrue(testCase, misc.alln(b.on() == c.on()));
end
function testPlus(testCase)
syms z zeta
+Z = quantity.Domain('z', linspace(0, 1, 3));
+Zeta = quantity.Domain('zeta', linspace(0, 1, 4));
+
%% 1 + 2 variables
-fSymVec3 = quantity.Symbolic([sinh(pi * z)], 'name', 'sinh');
-fSymVec4 = quantity.Symbolic([cosh(zeta * pi * z)], 'name', 'cosh');
+fSymVec3 = quantity.Symbolic([sinh(pi * z)], 'name', 'sinh', 'domain', Z);
+fSymVec4 = quantity.Symbolic([cosh(zeta * pi * z)], 'name', 'cosh', ...
+ 'domain', [Z Zeta]);
result = fSymVec3 - fSymVec4;
resultDiscrete = quantity.Discrete(fSymVec3) - quantity.Discrete(fSymVec4);
@@ -527,7 +544,9 @@ resultDiscrete = quantity.Discrete(fSymVec3) - quantity.Discrete(fSymVec4);
testCase.verifyEqual(result.on(), resultDiscrete.on(), 'AbsTol', 10*eps);
%% quantity + constant
-myQuan = quantity.Symbolic([sinh(pi * z), cosh(zeta * pi * z); 1, z^2], 'name', 'myQuan');
+
+myQuan = quantity.Symbolic([sinh(pi * z), cosh(zeta * pi * z); 1, z^2], ...
+ 'name', 'myQuan', 'domain', [Z Zeta]);
myQuanPlus1 = myQuan + ones(size(myQuan));
my1PlusQuan = ones(size(myQuan)) + myQuan;
@@ -537,7 +556,8 @@ testCase.verifyEqual(myQuan.on()+1, my1PlusQuan.on(), 'AbsTol', 10*eps);
%% 1 variable
-fSymVec = quantity.Symbolic([sinh(z * pi) ; cosh(z * pi)], 'name', 'sinhcosh');
+fSymVec = quantity.Symbolic([sinh(z * pi) ; cosh(z * pi)], ...
+ 'domain', Z, 'name', 'sinhcosh');
F = fSymVec + fSymVec;
val = F.on();
@@ -556,7 +576,8 @@ syms zeta
f1fun = @(z, zeta) sinh(z * pi);
f2fun = @(z, zeta) cosh(zeta * pi) + 0*z;
-fSymVec = quantity.Symbolic([sinh(z * pi); cosh(zeta * pi)], 'name', 'sinhcosh');
+fSymVec = quantity.Symbolic([sinh(z * pi); cosh(zeta * pi)], ...
+ 'domain', [Z Zeta], 'name', 'sinhcosh');
fSymVec2 = [fSymVec(2); fSymVec(1)];
F = fSymVec + fSymVec;
@@ -579,7 +600,8 @@ syms z
f1 = sinh(z * pi);
f2 = cosh(z * pi);
-f = quantity.Symbolic([f1 ; f2 ], 'name', 'sinhcosh');
+f = quantity.Symbolic([f1 ; f2 ], 'name', 'sinhcosh', 'domain', ...
+ quantity.Domain('z', linspace(0,1,3)));
F = f * f.';
val = F.on();
@@ -605,36 +627,42 @@ end
function testVariable(testCase)
syms z
-q = quantity.Symbolic(z);
+q = quantity.Symbolic(z, 'domain', quantity.Domain('z', 1));
verifyEqual(testCase, q.variable, z);
end
-function testUMinus(testCase)
+function testUMinusAndUPlus(testCase)
%%
syms z;
f1 = (sin(z .* pi) );
f2 = z;
f3 = (cos(z .* pi) );
-zz = linspace(0, 1, 42).';
-f = quantity.Symbolic([f1, f2; 0, f3], ...
- 'grid', {zz});
+d = quantity.Domain('z', linspace(0, 1, 3).');
-mf = -f;
+f = quantity.Symbolic([f1, f2; 0, f3], ...
+ 'domain', d);
+fDisc = quantity.Discrete(f);
+values = f.on();
%%
-verifyEqual(testCase, [f.valueDiscrete], -[mf.valueDiscrete]);
+testCase.verifyEqual(on( -f ), -values);
+testCase.verifyEqual(on( + ( - f ) ), -values );
+testCase.verifyEqual(on( + fDisc), on(+f));
+testCase.verifyEqual(on( - fDisc), on(-f));
+
+
end
function testConstantValues(testCase)
%%
syms z
-q = quantity.Symbolic(magic(7), 'grid', {linspace(0, 1)'}, 'variable', z);
+q = quantity.Symbolic(magic(3), 'grid', {linspace(0, 1, 4)'}, 'gridName', 'z');
%%
-magic7 = magic(7);
-magic700 = zeros(100, 7, 7);
+magic7 = magic(3);
+magic700 = zeros(4, 3, 3);
for k = 1:q.gridSize()
magic700(k, :,:) = magic7;
end
@@ -643,7 +671,7 @@ end
testCase.verifyEqual(magic700, q.on());
testCase.verifyTrue(q.isConstant());
-p = quantity.Symbolic([0 0 0 0 z], 'grid', linspace(0,1)', 'variable', z);
+p = quantity.Symbolic([0 0 0 0 z], 'grid', linspace(0, 1, 4)', 'gridName', 'z');
testCase.verifyFalse(p.isConstant());
end
@@ -651,8 +679,8 @@ end
function testIsConstant(testCase)
%%
syms z
-Q_constant = quantity.Symbolic(1);
-Q_var = quantity.Symbolic(z);
+Q_constant = quantity.Symbolic(1, 'domain', quantity.Domain.empty());
+Q_var = quantity.Symbolic(z, 'domain', quantity.Domain.empty());
%%
verifyTrue(testCase, Q_constant.isConstant);
@@ -662,22 +690,23 @@ end
function testEqual(testCase)
%%
-Q = quantity.Symbolic(0, 'name', 'P');
-P = quantity.Symbolic(0, 'name', 'Q');
+Q = quantity.Symbolic(0, 'name', 'P', 'domain', quantity.Domain.empty());
+P = quantity.Symbolic(0, 'name', 'Q', 'domain', quantity.Domain.empty());
syms Z;
f1 = (sin(Z .* pi) );
f2 = Z;
f3 = (cos(Z .* pi) );
-z = linspace(0, 1, 42).';
+z = quantity.Domain('Z', linspace(0, 1, 4));
+
R1 = quantity.Symbolic([f1, f2; 0, f3], ...
- 'grid', {z});
+ 'domain', z);
R2 = quantity.Symbolic([f2, f1; 0, f3], ...
- 'grid', {z});
+ 'domain', z);
R3 = quantity.Symbolic([f1, f2; 0, f3], ...
- 'grid', {z});
+ 'domain', z);
%%
verifyTrue(testCase, Q == P);
@@ -690,16 +719,18 @@ end
function testEvaluateConstant(testCase)
%%
syms z zeta
-Q = quantity.Symbolic(magic(7), 'grid', {1}, 'variable', z);
-Q2D = quantity.Symbolic(magic(7), 'grid', {1, 1}, 'variable', {z, zeta});
+Z = quantity.Domain('z', linspace(0,1,1)');
+Zeta = quantity.Domain('zeta', linspace(0,1,1)');
+Q = quantity.Symbolic(magic(2), 'domain', Z);
+Q2D = quantity.Symbolic(magic(2), 'domain', [Z, Zeta]);
%%
-verifyEqual(testCase, shiftdim(Q.on(), 1), magic(7))
-verifyEqual(testCase, shiftdim(Q2D.on(), 2), magic(7))
+verifyEqual(testCase, shiftdim(Q.on(), 1), magic(2))
+verifyEqual(testCase, shiftdim(Q2D.on(), 2), magic(2))
end
function testSize(testCase)
-Q = quantity.Symbolic(ones(1,2,3,4));
+Q = quantity.Symbolic(ones(1,2,3,4), 'domain', quantity.Domain.empty());
verifyEqual(testCase, size(Q), [1, 2, 3, 4]);
end
@@ -709,9 +740,10 @@ syms Z;
f1 = (sin(Z .* pi) );
f2 = Z;
f3 = (cos(Z .* pi) );
-z = linspace(0, 1, 42).';
+z = linspace(0, 1, 4).';
+zDomain = quantity.Domain('Z', z);
f = quantity.Symbolic([f1, f2; 0, f3], ...
- 'grid', {z});
+ 'domain', zDomain);
value = f.on();
@@ -725,13 +757,15 @@ end
function testSolveAlgebraic(testCase)
syms x
assume(x>0 & x<1);
+
+X = quantity.Domain('x', linspace(0,1,11));
% scalar
-fScalar = quantity.Symbolic([(1+x)^2]);
+fScalar = quantity.Symbolic([(1+x)^2], 'domain', X);
solutionScalar = fScalar.solveAlgebraic(2, fScalar.gridName{1});
testCase.verifyEqual(solutionScalar, sqrt(2)-1, 'AbsTol', 1e-12);
% array
-f = quantity.Symbolic([2*x, 1]);
+f = quantity.Symbolic([2*x, 1], 'domain', X);
solution = f.solveAlgebraic([1, 1], f(1).gridName{1});
testCase.verifyEqual(solution, 0.5);
end
@@ -739,17 +773,29 @@ end
function testSubs(testCase)
%%
% init
-syms x y z bla
-f = quantity.Symbolic([x, y^2; y+x, 1]);
+syms x y z w
+X = quantity.Domain('x', linspace(0, 1, 11));
+Y = quantity.Domain('y', linspace(0, 1, 7));
+Z = quantity.Domain('z', linspace(0, 1, 5));
+W = quantity.Domain('w', linspace(0, 1, 3));
+
+f = quantity.Symbolic([x, y^2; y+x, 1], 'domain', [X Y]);
% subs on variable
Ff1 = [1, 1^2; 1+1, 1];
-Ffz = quantity.Symbolic([z, z^2; z+z, 1]);
-Fyx = quantity.Symbolic([bla, x^2; x+bla, 1]);
+Ffz = quantity.Symbolic([z, z^2; z+z, 1], 'domain', Z);
+Fyx = quantity.Symbolic([w, x^2; x+w, 1], 'domain', [X W]);
+% test 1: replace all domains by numerical values
fOf1 = f.subs({'x', 'y'}, {1, 1});
-fOfz = f.subs({'x', 'y'}, {'z', 'z'});
-fOfyx = f.subs({'x', 'y'}, {'bla', 'x'});
+
+% test 2: replace one domain by a numerical value
+fOf2 = f.subs({'x'}, {1});
+
+Z0 = quantity.Domain('z', 0);
+
+fOfz = f.subs([X, Y], [Z Z]);
+fOfyx = f.subs({'x', 'y'}, {'w', 'x'});
testCase.verifyEqual(Ff1, fOf1);
testCase.verifyEqual(Ffz.sym(), fOfz.sym());
@@ -758,20 +804,20 @@ testCase.verifyEqual(Fyx.sym(), fOfyx.sym());
%%
syms z zeta
assume(z>0 & z<1); assume(zeta>0 & zeta<1);
-F = quantity.Symbolic(zeros(1), 'grid', {linspace(0,1), linspace(0,1)}, 'variable', {z, zeta});
+F = quantity.Symbolic(zeros(1), 'grid', ...
+ {linspace(0, 1, 3), linspace(0, 1, 3)}, 'gridName', {'z', 'zeta'});
Feta = F.subs('z', 'eta');
testCase.verifyEqual(Feta(1).gridName, {'eta', 'zeta'});
%%
-fMessy = quantity.Symbolic(x^1*y^2*z^4*bla^5, 'variables', {x, y, z, bla}, ...
- 'grid', {linspace(0, 1, 5), linspace(0, 1, 11), linspace(0, 1, 15), linspace(0, 1, 21)});
-fMessyA = fMessy.subs({'x', 'y', 'z', 'bla'}, {'a', 'a', 'a', 'a'});
-fMessyYY = fMessy.subs({'bla', 'x', 'z'}, {'bli', 'y', 'y'});
+fMessy = quantity.Symbolic(x^1*y^2*z^4*w^5, 'domain', [X, Y, Z, W]);
+fMessyA = fMessy.subs({'x', 'y', 'z', 'w'}, {'a', 'a', 'a', 'a'});
+fMessyYY = fMessy.subs({'w', 'x', 'z'}, {'xi', 'y', 'y'});
testCase.verifyEqual(fMessyA.gridName, {'a'});
-testCase.verifyEqual(fMessyA.grid, {linspace(0, 1, 21)});
-testCase.verifyEqual(fMessyYY.gridName, {'y', 'bli'});
-testCase.verifyEqual(fMessyYY.grid, {linspace(0, 1, 15), linspace(0, 1, 21)});
+testCase.verifyEqual(fMessyA.grid, {linspace(0, 1, 11)'});
+testCase.verifyEqual(fMessyYY.gridName, {'y', 'xi'});
+testCase.verifyEqual(fMessyYY.grid, {linspace(0, 1, 11)', linspace(0, 1, 3)'});
% % sub multiple numerics -> not implemented yet
% f11 = f.subs('x', [1; 2]);
% f1 = f.subs('x', 1);
diff --git a/+unittests/+quantity/testTimeDelayOperator.m b/+unittests/+quantity/testTimeDelayOperator.m
new file mode 100644
index 0000000000000000000000000000000000000000..0e6dcf9e8ad745aa6b16b1401472ea254ac16a1f
--- /dev/null
+++ b/+unittests/+quantity/testTimeDelayOperator.m
@@ -0,0 +1,117 @@
+function [tests] = testTimeDelayOperator()
+%TESTGEVREY Summary of this function goes here
+% Detailed explanation goes here
+tests = functiontests(localfunctions);
+end
+
+function setupOnce(testCase)
+
+z = quantity.Domain('z', linspace(0,1));
+zeta = quantity.Domain('zeta', linspace(0,1));
+t = quantity.Domain('t', linspace(0, 2, 31));
+v = quantity.Discrete( z.grid, 'domain', z );
+c_zeta = quantity.Discrete( zeta.grid * 0, 'domain', zeta );
+c = quantity.Discrete( 1, 'domain', quantity.Domain.empty);
+
+testCase.TestData.z = z;
+testCase.TestData.t = t;
+testCase.TestData.zeta = zeta;
+testCase.TestData.coefficient.fun_zeta = c_zeta;
+testCase.TestData.coefficient.variable = v;
+testCase.TestData.coefficient.constant = c;
+
+testCase.TestData.delay.variable = quantity.TimeDelayOperator(v, t);
+testCase.TestData.delay.constant = quantity.TimeDelayOperator(c, t);
+testCase.TestData.delay.zero = quantity.TimeDelayOperator.zero(t);
+testCase.TestData.delay.spatialDomain2 = quantity.TimeDelayOperator( ...
+ v - c_zeta, t);
+
+testCase.TestData.delay.diagonal = diag( ...
+ [testCase.TestData.delay.constant, ...
+ testCase.TestData.delay.variable]);
+end
+
+function testInit(testCase)
+
+td = testCase.TestData;
+
+testCase.verifyEqual(td.delay.constant.coefficient, ...
+ td.coefficient.constant);
+testCase.verifyEqual( td.delay.variable.coefficient, td.coefficient.variable);
+
+end
+
+
+function testDiag(testCase)
+
+ D = testCase.TestData.delay.diagonal;
+ v = testCase.TestData.delay.variable;
+ c = testCase.TestData.delay.constant;
+ o = testCase.TestData.delay.zero;
+
+ testCase.verifyEqual( D(1), c );
+ testCase.verifyEqual( D(2), o );
+ testCase.verifyEqual( D(3), o );
+ testCase.verifyEqual( D(4), v );
+
+end
+
+function testApplyTo(testCase)
+
+ z = testCase.TestData.z;
+ t = testCase.TestData.t;
+ tau = quantity.Domain('tau', linspace(-1, 3, 201));
+ v = testCase.TestData.delay.variable;
+
+ h = quantity.Discrete( sin(tau.grid * 2 * pi), 'domain', tau );
+ H = quantity.Discrete( sin(z.grid * 2 * pi + t.grid' * 2 * pi), ...
+ 'domain', [z, t]);
+ H_ = v.applyTo(h);
+
+ %% test a scalar operator
+ testCase.verifyEqual(H.on(), H_.on(), 'AbsTol', 5e-3);
+
+ d2 = testCase.TestData.delay.spatialDomain2;
+ H2 = d2.applyTo(h);
+ H2 = subs(H2, 'zeta', 0);
+ testCase.verifyEqual(H.on(), H2.on(), 'AbsTol', 5e-3);
+
+ %% 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); ...
+ sin(z.grid * 2 * pi + t.grid' * 2 * pi)}, 'domain', [z, t]);
+ D = testCase.TestData.delay.diagonal;
+ H_ = D.applyTo(h);
+
+ 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
+
+
+
+
+
+
+
+
+
+