diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 123bbc1b80fa5ed2770b1d9c521de24e669c2b4e..650e062a842267a307a26218ebbdd11e7ad18495 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -90,23 +90,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
domainParser.parse(varargin{:});
if domainParser.isDefault('domain') && ...
- domainParser.isDefault('grid') && ...
- domainParser.isDefault('gridName')
+ ( domainParser.isDefault('grid') || ...
+ domainParser.isDefault('gridName') )
% case 1: nothing about the grid is defined
% -> use default grid
- assert(iscell(valueOriginal), 'If no grid is specified, valueOriginal must be a cell-array')
- valueOriginalGridSize = size(valueOriginal{1});
- myDomain = quantity.Domain.defaultGrid(valueOriginalGridSize);
-
- elseif domainParser.isDefault('domain') && ...
- domainParser.isDefault('grid')
- % case 2: the gridNames are specified
- % -> use default grid with specified names
- assert(iscell(valueOriginal), 'If no grid is specified, valueOriginal must be a cell-array')
- valueOriginalGridSize = size(valueOriginal{1});
- myDomain = quantity.Domain.defaultGrid(valueOriginalGridSize, ...
- misc.ensureIsCell(domainParser.Results.gridName));
-
+ 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
@@ -224,7 +212,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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
@@ -287,46 +275,61 @@ 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 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();
+ for k = 1:length(newGrid)
+ myDomain(k) = quantity.Domain('grid', newGrid{k}, 'name', obj(1).domain(k).name);
+ 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);
+ elseif nargin == 3
+ assert(iscell(myDomain), 'If the domain is specified by cell-array pairs, the value myDomain must be a cell-array with grid entries')
+ 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('grid', newGrid{k}, 'name', gridNames{k});
+ end
+ else
+ myDomain = obj(1).domain;
end
+
+ % verify the domain
+ assert(numel(myDomain) == numel(obj(1).domain), ['Wrong grid for the evaluation of the object']);
+ [myDomain, gridPermuteIdx] = obj(1).domain.permute(myDomain);
+
+ % get the valueDiscrete data for this object. Apply the
+ % permuted myDomain. Then the obj2value will be evaluated
+ % in the order of the original domain. The permuatation to
+ % the new order will be done in the next step.
+ value = obj.obj2value(myDomain(gridPermuteIdx));
- value = 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);
- end
- value = permute(reshape(value, [cellfun(@(v) numel(v), myGrid), size(obj)]), ...
+ value = permute(reshape(value, [cellfun(@(v) numel(v), {myDomain(gridPermuteIdx).grid}), size(obj)]), ...
[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
end
end
@@ -422,37 +425,6 @@ 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.
@@ -814,6 +786,26 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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 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;
@@ -822,14 +814,19 @@ 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');
@@ -865,8 +862,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)};
@@ -886,7 +883,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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
@@ -900,7 +897,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
@@ -930,9 +927,12 @@ 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
% with the same size as size(obj).
@@ -1017,20 +1017,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();
@@ -1222,28 +1208,37 @@ 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.
+ %
+ % 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};
+
+ gridNew = misc.ensureIsCell(gridNew);
+ gridNameNew = misc.ensureIsCell(gridNameNew);
+
+ [gridIndexNew, logIdx] = obj(1).domain.gridIndex(gridNameNew);
+ newDomain = obj(1).domain;
+
+ for i = 1 : length(gridIndexNew)
+ newDomain(gridIndexNew(i)) = quantity.Domain(...
+ 'grid', gridNew{i}, 'name', gridNameNew{i});
end
- assert(isequal(myGridName(:), obj(1).gridName(:)), 'rearranging grids failed');
+ assert(isequal({newDomain.name}, 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)
@@ -1438,6 +1433,9 @@ 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
@@ -1469,38 +1467,48 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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)];
+ parameters.figureID = a(1).figureID;
+
+
+ domainA = a(1).domain;
+ domainB = b(1).domain;
% select finest grid from a and b
- parameters.grid = cell(1, numel(parameters.gridName));
- [gridJoined, gridNameJoined] = gridJoin(a, b);
+ joinedDomain = gridJoin(domainA, domainB);
+ newDomainName = {joinedDomain.name};
+
% 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};
+
+ finalGridOrder = [domainA(idx.A.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(newDomainName)
+
+ parameters.domain(i) = finalGridOrder(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(finalGridOrder(i).name);
+ idxB = domainB.gridIndex(finalGridOrder(i).name);
+
+ if idxA
+ newDomainA = [newDomainA, finalGridOrder(i)];
end
- if b.gridIndex(parameters.gridName{it})
- bGrid{b.gridIndex(parameters.gridName{it})} = ...
- parameters.grid{it};
+
+ if idxB
+ newDomainB = [newDomainB, finalGridOrder(i)];
end
end
parameters = misc.struct2namevaluepair(parameters);
- valueA = permute(a.on(aGrid), idx.A.permute);
- valueB = permute(b.on(bGrid), idx.B.permute);
+ %%
+ valueA = a.on(newDomainA);
+ valueB = b.on(newDomainB);
C = misc.multArray(valueA, valueB, idx.A.value(end), idx.B.value(1), idx.common);
C = permute(C, permuteGrid);
@@ -1712,7 +1720,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
@@ -1721,7 +1729,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
@@ -1731,37 +1739,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
@@ -1888,8 +1865,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
myGrid = currentGrid(domainIdx{idxGrid});
-
-
for kObj = 1:numel(obj)
J = numeric.trapz_fast_nDim(myGrid, obj(kObj).atIndex(domainIdx{:}), idxGrid);
@@ -1900,21 +1875,22 @@ 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)
+ grdName = '';
+ newGrid = {};
+ else
+ grdName = obj(kObj).gridName{~isGrid};
+ newGrid = obj(kObj).grid{~isGrid};
+ end
+
+ newDomain = quantity.Domain('grid', newGrid, 'name', grdName);
+ [I.domain] = deal(newDomain);
+
end
function result = cumInt(obj, domain, lowerBound, upperBound)
@@ -1930,7 +1906,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', ...
@@ -1939,7 +1915,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}, ...
@@ -1968,54 +1944,57 @@ 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);
+ [aDiscrete] = A.expandValueDiscrete(joinedDomain);
+ [bDiscrete] = B.expandValueDiscrete(joinedDomain);
% create result object
C = quantity.Discrete(aDiscrete + bDiscrete, ...
- 'grid', gridJoined, 'gridName', gridNameJoined, ...
- 'size', size(A), 'name', [A(1).name, '+', B(1).name]);
+ 'domain', joinedDomain, ...
+ 'name', [A(1).name, '+', B(1).name]);
end
- function [valDiscrete] = expandValueDiscrete(obj, gridJoined, gridNameJoined, gridJoinedLength)
+ function [valDiscrete] = expandValueDiscrete(obj, newDomain)
% EXPANDVALUEDISCRETE
% [valDiscrete, gridNameSorted] = ...
% expandValueDiscrete(obj, gridIndex, valDiscrete)
% Expanses the value of obj, so that
- % todo
+
+
+ gridNameJoined = {newDomain.name};
+ gridJoinedLength = newDomain.gridLength;
% get the index of obj.grid in the joined grid
- gridLogical = logical( obj.gridIndex(gridNameJoined) );
+ [~, logicalIdx] = newDomain.gridIndex({obj(1).domain.name});
- valDiscrete = obj.on(gridJoined(gridLogical > 0), gridNameJoined(gridLogical > 0));
+ valDiscrete = obj.on( newDomain(logicalIdx) );
oldDim = ndims(valDiscrete);
- valDiscrete = permute(valDiscrete, [(1:sum(~gridLogical)) + oldDim, 1:oldDim] );
- valDiscrete = repmat(valDiscrete, [gridJoinedLength(~gridLogical), ones(1, ndims(valDiscrete))]);
+ valDiscrete = permute(valDiscrete, [(1:sum(~logicalIdx)) + oldDim, 1:oldDim] );
+ valDiscrete = repmat(valDiscrete, [gridJoinedLength(~logicalIdx), ones(1, ndims(valDiscrete))]);
%
valDiscrete = reshape(valDiscrete, ...
- [gridJoinedLength(~gridLogical), gridJoinedLength(gridLogical), size(obj)]);
+ [gridJoinedLength(~logicalIdx), gridJoinedLength(logicalIdx), 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))]);
+ % permute valDiscrete such that grids are in the order specified
+ % by gridNameJoined.
+ gridIndex = 1:numel(logicalIdx);
+ gridOrder = [gridIndex(~logicalIdx), gridIndex(logicalIdx)];
+ valDiscrete = permute(valDiscrete, [gridOrder, numel(logicalIdx)+(1:ndims(obj))]);
end
@@ -2155,15 +2134,32 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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 myDomain ~= obj(1).domain
+ % if a new domain is specified for the evaluation of
+ % the quantity, 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{:});
+ elseif obj.isConstant
+ value = repmat(value, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
+ end
end
@@ -2225,7 +2221,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
myGridSize = num2cell( fullSize(1:length(gridSize)) );
if nargin == 2
- valueSize = [fullSize(length(myGridSize)+1:length(fullSize)), 1];
+ valueSize = [fullSize(length(myGridSize)+1:length(fullSize)), 1, 1];
end
myValueSize = arrayfun(@(n) ones(n, 1), valueSize, 'UniformOutput', false);
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index 1484f5c5e9645fe3aaf11cfbb27766423c64598f..fb8647eed11b52e69942e617b2ff9e9b2af21e63 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -1,7 +1,7 @@
classdef Domain
%DOMAIN class to describes a range of values on which a function can be defined.
- % todo:
+ % todo:
% * EquidistantDomain
% * multi dimensional
@@ -9,16 +9,18 @@ classdef Domain
% 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};
+ % generated by X_grid = linspace(0, 1).
+ grid double {mustBeReal};
% a speaking name for this domain; Should be unique, so that the
% domain can be identified by the name.
name char;
+
+ isConstant = false ;
end
properties (Dependent)
- n; % number of discretization points
+ n; % number of discretization points
lower; % lower bound of the domain
upper; % upper bound of the domain
end
@@ -29,13 +31,13 @@ classdef Domain
%
if nargin >= 1
parser = misc.Parser();
- parser.addParameter('grid', [], @isvector);
+ parser.addParameter('grid', [], @(g) isvector(g) || isempty(g));
parser.addParameter('name', '', @ischar);
parser.parse(varargin{:});
-
+
% todo: assertions
% * ascending ?
-
+
obj.grid = parser.Results.grid(:);
obj.name = parser.Results.name;
else
@@ -72,6 +74,124 @@ classdef Domain
s = [obj.n];
end
+ function s = eq(obj, obj2)
+ 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
+
+ function l = ne(obj, obj2)
+ l = ~(obj == obj2);
+ end
+
+ function [joinedDomain] = 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.
+
+ [joinedGrid] = union({obj1.name}, {obj2.name}, 'sorted');
+
+ joinedDomain = quantity.Domain.empty;
+
+ % 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] = obj1.gridIndex(currentGridName);
+ [index2, logicalIndex2] = obj2.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 = obj1(index1);
+ tempDomain2 = obj2(index2);
+
+ if ~tempDomain1.isConstant && ~tempDomain2.isConstant
+ assert(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
+ assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
+ end
+ 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) = obj1(index1);
+ elseif any(logicalIndex2)
+ joinedDomain(i) = obj2(index2);
+ end
+
+ end
+ end % gridJoin()
+
+ 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.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)
+ idx(k) = nArgIdx(log);
+ else
+ idx(k) = 0;
+ end
+ end
+ end % gridIndex
+
+ 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 [newDomain, idx] = permute(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
end
methods (Static)
@@ -97,7 +217,7 @@ classdef Domain
g(k) = quantity.Domain('grid', O, 'name', name{k});
end
- end
+ end
end
end
diff --git a/+quantity/Function.m b/+quantity/Function.m
index 4b090acdc57b53388e2b9c58fb142af9540ba511..c8eba8e971b35257a8e7ce2d83c70211470beb8e 100644
--- a/+quantity/Function.m
+++ b/+quantity/Function.m
@@ -74,57 +74,25 @@ 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:
+ function value = obj2value(obj, myDomain, recalculate)
- gridSize = cellfun('length', myGrid);
- if isempty(gridSize), gridSize = 1; end
-
- 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 && (myDomain == obj(1).domain)
+ value = obj2value@quantity.Discrete(obj, myDomain);
+ else
+ % otherwise the function has to be evaluated on the new
+ % domain
+ value = nan([numel(obj), myDomain.gridLength]);
+ for k = 1:numel(obj)
+ ndGrd = myDomain.ndgrid;
+ tmp = obj(k).evaluateFunction( ndGrd{:} );
+ value(k,:) = tmp(:);
+ end
end
- value = reshape(permute(value, [1+(gridPermuteIdx), 1]), ...
- [gridSize(gridPermuteIdx), size(obj)]);
end
function n = nargin(obj)
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index fc20f9e5dab6eb9c45f2b8ae6089972a00c088bc..e18b0fc973d26609401aaa282130e7fbd8fa5031 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -472,23 +472,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('grid', linspace(0, 1, 27), 'name', 'z');
+domainVecB = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
+[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('grid', linspace(0, 1, 24), 'name', 'z');
+d21zeta = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'zeta');
+
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);
@@ -531,7 +536,7 @@ b = quantity.Symbolic((eye(2, 2)), 'variable', zeta, ...
'grid', gridVecB);
c = a*b;
-%%
+%
testCase.verifyEqual(c.on(), a.on());
end
@@ -575,7 +580,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 +609,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)
@@ -753,61 +763,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 +791,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 +800,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 +814,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 +831,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));
@@ -897,7 +852,7 @@ 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 +868,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;
@@ -1019,7 +974,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('grid', z, 'name', 'z');
+T = quantity.Domain('grid', t, 'name', 't');
+
+d = quantity.Discrete(V, 'domain', [Z, T]);
%%
verifyTrue(testCase, misc.alln(b.on() == c.on()));
@@ -1071,23 +1030,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('grid', linspace(0,1,100), 'name', 'z');
+T = quantity.Domain('grid', linspace(0,1,42), 'name', 't');
+V = quantity.Discrete( quantity.Discrete.value2cell(v, [100, 42], [2, 3]), ...
+ 'domain', [Z, T]);
verifyTrue(testCase, misc.alln(v == V.on()));
end
@@ -1139,6 +1087,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({linspace(0,1,11)}, {'t'});
quanSampled = quanCopy.changeGrid(gridSampled, {'z', 't'});
testCase.verifyEqual(quanSampled.on(), quan.on(gridSampled))
@@ -1147,14 +1096,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]);
+d1 = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'z');
+d2 = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
+
+quan = quantity.Discrete(myTestArray, 'domain', [d1, d2]);
+
quanZetaZ = quan.subs({'z', 'zeta'}, {'zeta', 'z'});
quan10 = quan.subs({'z', 'zeta'}, {1, 0});
quan01 = quan.subs({'z', 'zeta'}, {0, 1});
@@ -1162,7 +1113,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 +1126,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('grid', linspace(0,1,3), 'name', 'zeta')) ), 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('grid', linspace(0, 1, 11), 'name', 'z');
+x2 = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'zeta');
+x3 = quantity.Domain('grid', linspace(0, 1, 31), 'name', 'eta');
+x4 = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'beta');
+
+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 +1144,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
index 3f14eb8a4ca10a2f3192b285e2841e6a45a6d9f1..da316cad0b01cbea9e3daebea28a35ec9cf86b63 100644
--- a/+unittests/+quantity/testDomain.m
+++ b/+unittests/+quantity/testDomain.m
@@ -42,4 +42,94 @@ d = quantity.Domain();
testCase.verifyTrue( isempty(d) )
-end
\ No newline at end of file
+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('grid', s, 'name', 's');
+Z = quantity.Domain('grid', z, 'name', 'z');
+T = quantity.Domain('grid', t, 'name', 't');
+Ps = quantity.Domain('grid', s, 'name', 'p');
+St = quantity.Domain('grid', t, 'name', 's');
+Pt = quantity.Domain('grid', t, 'name', 'p');
+
+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 = gridJoin(a, b);
+joinedDomainCC = gridJoin(c, c);
+joinedDomainBC = gridJoin(b, c);
+
+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(:)});
+end
+
+function testGridIndex(testCase)
+z = linspace(0, 2*pi, 71)';
+t = linspace(0, 3*pi, 51);
+s = linspace(0, 1);
+
+S = quantity.Domain('grid', s, 'name', 's');
+Z = quantity.Domain('grid', z, 'name', 'z');
+T = quantity.Domain('grid', t, 'name', '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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+