Commit 68484bab authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

issue almost solved... some unittests in quantity.Function left

parent eb0f7d62
......@@ -17,20 +17,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
domain;
end
properties ( Dependent )
% Name of the domains that generate the grid.
gridName {mustBe.unique};
% Grid for the evaluation of the continuous quantity. For the
% example with the function f(x,t), the grid would be
% {[<spatial domain>], [<temporal domain>]}
% whereas <spatial domain> is the discrete description of the
% spatial domain x and <temporal domain> the discrete description of
% the temporal domain t.
grid; % in set.grid it is ensured that, grid is a (1,:)-cell-array
end
methods
%--------------------
% --- Constructor ---
......@@ -139,25 +125,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
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 = isNumber(obj)
% the quantity is interpreted as constant if it has no grid.
itIs = isempty(obj(1).domain);
......@@ -181,12 +148,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% make the object names:
if obj.nargin == 1
headers = cell(1, numel(obj) + 1);
headers{1} = obj(1).gridName{1};
headers{1} = obj(1).domain(1).name;
for i= 1:numel(obj) %TODO use easier to read headers
headers{i+1} = obj(i).name + "" + num2str(i);
end
thisGrid = {obj(1).domain.grid};
exportData = export.dd(...
'M', [obj(1).grid{:}, obj.valueDiscrete], ...
'M', [thisGrid{:}, obj.valueDiscrete], ...
'header', headers, varargin{:});
elseif obj.nargin == 2
error('Not yet implemented')
......@@ -513,7 +481,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
interpolant = numeric.interpolant(...
[indexGrid{:}], value);
else
myGrid = obj(1).grid;
myGrid = {obj(1).domain.grid};
value = obj.obj2value();
indexGrid = misc.indexGrid(size(obj));
interpolant = numeric.interpolant(...
......@@ -533,22 +501,22 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
return;
else
referenceGridName = a(1).gridName;
referenceGrid= a(1).grid;
referenceDomainName = [a(1).domain.name];
referenceGrid= {a(1).domain.grid};
end
for it = 1 : numel(a)
assert(isequal(referenceGridName, a(it).gridName), ...
'All elements of a quantity must have same gridNames');
assert(isequal(referenceGrid, a(it).grid), ...
'All elements of a quantity must have same grid');
for it = 2 : numel(a)
assert(isequal(referenceDomainName, [a(it).domain.name]), ...
'All elements of a quantity must have same domain.name');
assert(isequal(referenceGrid, {a(it).domain.grid}), ...
'All elements of a quantity must have same domain.grid');
end
if nargin > 1
b = varargin{1};
for it = 1 : numel(b)
assert(isequal(referenceGridName, b(it).gridName), ...
'All elements of a quantity must have same gridNames');
assert(isequal(referenceGrid, b(it).grid), ...
'All elements of a quantity must have same grid');
assert(isequal(referenceDomainName, [b(it).domain.name]), ...
'All elements of a quantity must have same domain.name');
assert(isequal(referenceGrid, {b(it).domain.grid}), ...
'All elements of a quantity must have same domain.grid');
end
end
if nargin > 2
......@@ -556,46 +524,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% iteratively by calling assertSameGrid() again.
assertSameGrid(varargin{:});
end
end
function [referenceGrid, referenceGridName] = getFinestGrid(a, varargin)
warning("DEPRICATED: Use quantity.Domain/join instead")
% find the finest grid of all input quantities by comparing
% gridSize for each iteratively.
if isempty(a) || isempty(a(1).grid)
if nargin > 1
[referenceGrid, referenceGridName] = varargin{1}.getFinestGrid(varargin{2:end});
else
referenceGrid = {};
referenceGridName = '';
end
return;
else
referenceGridName = a(1).gridName;
referenceGrid = a(1).grid;
referenceGridSize = [a(1).domain.n];
end
for it = 1 : numel(varargin)
if isempty(varargin{it}) || isempty(varargin{it}(1).domain)
continue;
end
assert(numel(referenceGridName) == numel(varargin{it}(1).gridName), ...
['For getFinestGrid, the gridName of all objects must be equal', ...
'. Maybe join() does what you want?']);
for jt = 1 : numel(referenceGridName)
comparisonGridSize = varargin{it}(1).domain.find(referenceGridName{jt}).n;
comparisonGrid = varargin{it}.gridOf(referenceGridName{jt});
assert(referenceGrid{jt}(1) == comparisonGrid(1), 'Grids must have same domain for combining them')
assert(referenceGrid{jt}(end) == comparisonGrid(end), 'Grids must have same domain for combining them')
if comparisonGridSize > referenceGridSize(jt)
referenceGrid{jt} = comparisonGrid;
referenceGridSize(jt) = comparisonGridSize;
end
end
end
end
end % assertSameGrid()
function [B, I] = sort(obj, direction)
% sort Sort in ascending or descending order.
......@@ -869,7 +798,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
if isempty(B)
Y = A;
else
assert(isequal(A(1).gridName, B(1).gridName), 'only implemented for same gridName');
assert(isequal([A(1).domain.name], [B(1).domain.name]), ...
"only implemented for same gridName");
% otherwise cat() throws an error
Y = [A, zeros(size(A, 1), size(B, 2)); ...
zeros(size(B, 1), size(A, 2)), B];
......@@ -880,7 +810,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
end % blkdiag()
function solution = solveAlgebraic(obj, rhs, gridName, objLimit)
function solution = solveAlgebraic(obj, rhs, domainName, objLimit)
%% this method solves
% obj(gridName) == rhs
% for the variable specified by gridName.
......@@ -892,44 +822,33 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% (grid) (see quantity.invert-method).
% The input objLimit specifies minimum and maximum of the
% values of obj, between which the solution should be searched.
assert(numel(obj(1).gridName) == 1);
assert(isequal(size(obj), [1, 1]));
if ~isequal(size(rhs), size(obj))
error('rhs has not the same size as quantity');
end
if ~iscell(gridName)
gridName = {gridName};
end
if numel(gridName) ~= 1
error('this function can only solve for one variable');
end
if isempty(strcmp(obj(1).gridName, gridName{1}))
error('quantity does not depend on variable');
arguments
obj (1, 1);
rhs (1, 1) double;
domainName (1, 1) string = obj.domain(1).name;
objLimit = [];
end
assert(obj(1).nargin == 1, "only implemented for quantites on 1 domain");
if nargin == 4
assert(numel(objLimit)==2, 'a lower and upper limit must be specified (or neither)');
assert(all(strcmp([obj(1).domain.name], domainName)), ...
"quantity does not depend on variable");
if ~isempty(objLimit)
assert(numel(objLimit)==2, "a lower and upper limit must be specified (or neither)");
objValueTemp = obj.on();
gridSelector = (objValueTemp >= objLimit(1)) & (objValueTemp <= objLimit(2));
gridSelector([max(1, find(gridSelector, 1, 'first')-1), ...
min(find(gridSelector, 1, 'last')+1, numel(gridSelector))]) = 1;
limitedGrid = obj(1).grid{1}(gridSelector);
limitedGrid = obj(1).domain.grid(gridSelector);
objCopy = obj.copy();
objCopy = objCopy.changeDomain(quantity.Domain(gridName, limitedGrid));
objInverseTemp = objCopy.invert(gridName);
objCopy = objCopy.changeDomain(quantity.Domain(domainName, limitedGrid));
objInverseTemp = objCopy.invert(domainName);
else
objInverseTemp = obj.invert(gridName);
objInverseTemp = obj.invert(domainName);
end
solution = objInverseTemp.on(rhs);
% solution = zeros(numel(obj), 1);
% for it = 1 : numel(obj)
% objInverseTemp = obj(it).invert(gridName);
% solution(it) = objInverseTemp.on(rhs(it));
% end
% solution = reshape(solution, size(obj));
end % solveAlgebraic()
function inverse = invert(obj, domainName)
......@@ -947,7 +866,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
assert(isscalar(obj), "invert is only implemented for scalars");
assert(numel(obj.nargin) == 1, "invert is only implemented for quantities on 1 domain");
inverse = quantity.Discrete(...
repmat(obj(1).grid{obj(1).domain.index(domainName)}(:), [1, size(obj)]), ...
repmat(obj(1).domain(obj(1).domain.index(domainName)).grid(:), [1, size(obj)]), ...
quantity.Domain([obj(1).name], obj.on()), ...
'name', domainName);
end % invert()
......@@ -963,8 +882,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
'this method is only implemented for quanitities with one domain');
myParser = misc.Parser();
myParser.addParameter('initialValueGrid', obj(1).grid{1});
myParser.addParameter('variableGrid', obj(1).grid{1});
myParser.addParameter('initialValueGrid', obj(1).domain(1).grid);
myParser.addParameter('variableGrid', obj(1).domain(1).grid);
myParser.addParameter('newGridName', 's');
myParser.parse(varargin{:});
......@@ -1180,25 +1099,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end % atIndex()
function n = nargin(obj)
% nargin is the number of domains defined for the object.
% FIXME: check if all funtions in this object have the same
% number of input values.
% FIXME: for some combinations of constant objects, it seems to be
% possible, that the quantity has a gridName but no grid.
% Actually this should not be allowed. This is quick and dirty
% work around.
n = min(numel(obj(1).gridName), numel(obj(1).grid));
n = numel(obj(1).domain);
end
function d = gridDiff(obj)
% #FIXME:
% 1) test for multidimensional grids
% 2) check that the grid is equally spaced
d = diff(obj(1).grid{:});
d = d(1);
end
% function d = gridDiff(obj)
%
% % #FIXME:
% % 1) test for multidimensional grids
% % 2) check that the grid is equally spaced
%
% d = diff(obj(1).grid{:});
% d = d(1);
% end
function H = plot(obj, varargin)
H = [];
......@@ -1346,7 +1261,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
%
assert(numel(obj) == 1, "argList must NOT be called for an array object");
s = struct(obj);
s = rmfield(s, ["domain", "grid", "gridName"]);
s = rmfield(s, "domain");
if nargin == 2
s = rmfield(s, exclude);
end
......@@ -1440,8 +1355,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for it = 1 : length( indexNew )
newDomain(indexNew(it)) = domain(it);
end
assert(isequal([newDomain.name], obj(1).gridName), ...
'rearranging grids failed');
assert(all(strcmp([newDomain.name], [obj(1).domain.name])), ...
"rearranging grids failed");
end
[newObj.domain] = deal(newDomain);
......@@ -1461,13 +1376,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% newObj = REPLACEGRID(obj, domain) changes the domain of the
% object specified by the name of DOMAIN into DOMAIN.
%
% newObj = REPLACEGRID(obj, domain, 'gridName', NEWNAME) changes the domain of the
% newObj = REPLACEGRID(obj, domain, 'domainName', NEWNAME) changes the domain of the
% object specified by NEWNAME into DOMAIN.
arguments
obj
myNewDomain quantity.Domain
optArgs.gridName = [myNewDomain.name];
optArgs.domainName = [myNewDomain.name];
end
if isempty(obj)
......@@ -1475,12 +1390,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
return;
end
assert( intersect(obj(1).gridName, optArgs.gridName) == optArgs.gridName );
assert( intersect([obj(1).domain.name], optArgs.domainName) == optArgs.domainName );
if obj(1).isNumber()
error("Not yet implemented")
else
indexNew = obj(1).domain.index(optArgs.gridName);
indexNew = obj(1).domain.index(optArgs.domainName);
% initialization of the newDomain array as quantity.Domain
% array. This is required in order to handle also
% quantity.EquidistantDomains:
......@@ -1534,17 +1449,17 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
gridName2 = ip.Results.gridName2;
if ip.isDefault('grid')
assert(numel(obj(1).grid) == 1, 'No grid is defined for the computation of the state-transition matrix. Use the name-value pair option "grid" to define the grid.');
myGrid = obj(1).grid{1};
assert(obj(1).nargin == 1, 'No grid is defined for the computation of the state-transition matrix. Use the name-value pair option "grid" to define the grid.');
myGrid = obj(1).domain(1).grid;
end
if ip.isDefault('gridName1')
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};
assert(obj(1).nargin == 1, 'No gridName is defined. Use name-value-pairs property "gridName1" to define a grid name');
gridName1 = [obj(1).domain.name];
end
if ip.isDefault('gridName2')
gridName2 = [gridName1 '0'];
gridName2 = [char(gridName1) '0'];
end
......@@ -2016,7 +1931,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
integralGridName = misc.ensureString(integralGridName);
if obj.nargin == 1 && all(strcmp(obj(1).gridName, integralGridName))
if obj.nargin == 1 && all(strcmp([obj(1).domain.name], integralGridName))
integratedValue = int(obj.' * optArg.weight * obj);
......@@ -2245,7 +2160,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
if ~iscell(myDomainName)
gridIdx = obj(1).domain.index(myDomainName);
if gridIdx > 0
myGrid = obj(1).grid{gridIdx};
myGrid = obj(1).domain(gridIdx).grid;
else
myGrid = [];
end
......@@ -2254,7 +2169,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for it = 1 : numel(myGrid)
gridIdx = obj(1).domain.index(myDomainName{it});
if gridIdx > 0
myGrid{it} = obj(1).grid{gridIdx};
myGrid{it} = obj(1).domain(gridIdx).grid;
else
myGrid{it} = [];
end
......@@ -2270,11 +2185,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% then diff applies the derivative w.r.t. all gridNames.
arguments
obj;
diffGridName (1,1) string = obj(1).gridName;
diffGridName (1,1) string = [obj(1).domain.name];
k uint64 = 1;
end
if obj.isNumber() && isempty(obj(1).gridName)
if obj.isNumber() && isempty(obj(1).domain)
result = quantity.Discrete(zeros(size(obj)), obj(1).domain, ...
"name", "(d_{.}" + obj(1).name + ")");
return
......@@ -2297,10 +2212,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function I = int(obj, varargin)
% INT integrate
% INT(obj) is the definite integral of obj with respect to
% all of its independent variables over the grid obj(1).grid.
% INT(obj, gridName) is the definite integral of obj with
% all of its independent variables over the grid obj(1).domain.grid.
% INT(obj, domainName) is the definite integral of obj with
% respect to the spatial coordinate grid name, that has to
% be one of obj(1).gridName. The names must be a string or a
% be one of obj(1).domain.name. The names must be a string or a
% array of strings.
% INT(obj, a, b) (question: on which domain? what if numel(gridName)>1?)
% INT(obj, gridName, a, b)
......@@ -2367,11 +2282,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
"integral bounds must be string or numeric");
% get grid
myGrid = obj(1).grid;
intGridIdx = obj(1).domain.index(domainName);
% integrate
F = numeric.cumtrapz_fast_nDim(myGrid{intGridIdx}, ...
F = numeric.cumtrapz_fast_nDim(obj(1).domain(intGridIdx).grid, ...
obj.on(), intGridIdx);
result = quantity.Discrete(F, obj(1).domain);
......@@ -2487,11 +2401,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
elseif isnumeric(B)
b = B;
else
b = B.on(obj(1).grid, obj(1).gridName);
b = B.on(obj(1).domain);
end
[i, m, s] = numeric.near(...
obj.on(obj(1).grid, obj(1).gridName), ...
obj.on(), ...
b, varargin{:});
if nargout == 0
i = s;
......@@ -2638,7 +2552,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
value = reshape(cat(numel(obj(1).domain)+1, obj(:).valueDiscrete), ...
[obj(1).domain.gridLength(), size(obj)]);
indexGrid = misc.indexGrid(size(obj));
myInterpolant = numeric.interpolant([obj(1).grid, indexGrid{:}], value);
myInterpolant = numeric.interpolant([{obj(1).domain.grid}, indexGrid{:}], value);
end % getInterpolant()
function result = diag2vec(obj)
......@@ -2893,7 +2807,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
gridJoinedLength = newDomain.gridLength;
% get the index of obj.grid in the joined grid
% get the index of obj.domain in the joined grid
[~, logicalIdx] = newDomain.index([obj(1).domain.name]);
% evaluate the
valDiscrete = obj.on( newDomain(logicalIdx) );
......@@ -2917,7 +2831,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
index = obj(1).domain.index(diffGridName);
permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
permutationVector = 1 : (obj(1).nargin+ndims(obj));
objDiscrete = permute(obj.on(), ...
[permutationVector(index), ...
......@@ -2971,17 +2885,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% Computes the required permutation vectors to use
% misc.multArray for multiplication of matrices
% #todo this assertion should not be required!
% uA = unique(a(1).gridName, 'stable');
% assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
% uB = unique(b(1).gridName, 'stable');
% assert(numel(uB) == numel(b(1).gridName), 'Gridnames have to be unique!');
% 1) find common entries
if isempty( b(1).gridName ) || isempty( a(1).gridName )
if isempty(b(1).domain) || isempty(a(1).domain)
common = [];
else
common = intersect(a(1).gridName, b(1).gridName);
common = intersect([a(1).domain.name], [b(1).domain.name]);
end
commonA = false(1, a.nargin());
......@@ -2998,13 +2906,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
c = common{k};
% 2) find logical indices for the common entries
cA = strcmp(c, a(1).gridName);
cA = strcmp(c, [a(1).domain.name]);
commonA = commonA | cA;
% 3) exchange the order of the logical indices
gridA(k) = idxA0(cA);
cB = strcmp(c, b(1).gridName);
cB = strcmp(c, [b(1).domain.name]);
commonB = commonB | cB;
gridB(k) = idxB0(cB);
......
......@@ -155,16 +155,16 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
end
function i = isequal(obj, a)
function isEq = isequal(obj, a)
% isequal checks if two domains are equal.
i = numel(obj) == numel(a);
isEq = numel(obj) == numel(a);
for k = 1:numel(obj)
i = i && all(obj(k).name == a(k).name);
i = i && obj(k).n == a(k).n;
i = i && all( obj(k).grid == a(k).grid );
isEq = isEq && all(obj(k).name == a(k).name);
isEq = isEq && isequal(obj(k).grid, a(k).grid);
end
end
end % isequal()
function d = rename(obj, newName, oldName)
% RENAME renames the name of the domain OBJ.
......
......@@ -147,11 +147,6 @@ classdef Function < quantity.Discrete
end
end % at()
function n = nargin(obj)
% FIXME: check if all funtions in this object have the same
% number of input values.
n = numel(obj(1).grid);
end % nargin()
end % methods
% -------------
......@@ -167,7 +162,7 @@ classdef Function < quantity.Discrete
end
mObj.setName("-" + obj(1).name);
end
end % uminus
function p = inner(A, B)
......@@ -177,14 +172,14 @@ classdef Function < quantity.Discrete
m = size(B, 2);
p = zeros(n, m);
z0 = A(1).grid{1}(1);
z1 = A(1).grid{1}(end);
z0 = A(1).domain(1).lower;
z1 = A(1).domain(1).upper;
for k = 1 : n*m
p(k) = AB(k).int(z0, z1);
end
end
end % inner()
function I = cumInt(obj, domainName, lowerBound, upperBound)
% CUMINT cumulative integration
......
......@@ -2,7 +2,7 @@ classdef Symbolic < quantity.Function
properties (SetAccess = protected)
valueSymbolic sym;
variable sym;
variable (:, 1) sym;
end
properties
......@@ -102,7 +102,7 @@ classdef Symbolic < quantity.Function
end
end % Symbolic-Constructor
function i = eq(A, B)
function isEq = eq(A, B)
%== Equal.
% A == B does element by element comparisons between A and B and returns
% an array with elements set to logical 1 (TRUE) where the relation is
......@@ -118,33 +118,33 @@ classdef Symbolic < quantity.Function
if isempty(A)
if isempty(B)
if isequal(size(A),size(B))
i=1;
isEq = true;
return
else
i=0;
isEq = false;
return
end
else
i=0;
isEq = false;
return;
end
else
if isempty(B)
i=0;
isEq = false;
return
end
i = all(size(A) == size(B));
isEq = all(size(A) == size(B));
if ~i
if ~isEq