Commit 50f6d678 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

quantity.Discrete & Symbolic:

- bugfix in get.grid()
- horzcat, vertcat, cat now ensure, that grid and gridName of all quantities to be concatenated are compatible, by changing the grid to the finest grid.
parent b2978d30
......@@ -27,8 +27,8 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% ID of the figure handle in which the handle is plotted
figureID double = 1;
% TODO@ff vermutlich ist es schöner einen converter auf dieses
% Objekt zu schreiben, als es hier als Eigenschaft dran zu hängen.
% TODO@ff vermutlich ist es schner einen converter auf dieses
% Objekt zu schreiben, als es hier als Eigenschaft dran zu hngen.
exportData export.Data;
% Name of this object
......@@ -184,7 +184,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
grid = {grid};
end
isV = cellfun(@(v) isvector(v), grid);
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);
......@@ -206,7 +206,6 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
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!
......@@ -216,6 +215,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
else
if nargin == 1
myGrid = obj(1).grid;
myGridName = obj(1).gridName;
elseif nargin >= 2 && ~iscell(myGrid)
myGrid = {myGrid};
end
......@@ -247,6 +247,139 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
end
function assertSameGrid(a, varargin)
% check if all elements of a have same grid and gridName. If
% further quantites are inputs via varargin, it is verified if
% that quantity has same grid and gridName as quantity a
% as well.
referenceGridName = a(1).gridName;
referenceGrid= a(1).grid;
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');
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');
end
end
if nargin > 2
% if more then 1 quantity is in varargin, they are checked
% iteratively by calling assertSameGrid() again.
assertSameGrid(varargin{:});
end
end
function [referenceGrid, referenceGridName] = getFinestGrid(a, varargin)
% find the finest grid of all input quantities by comparing
% gridSize for each iteratively.
referenceGridName = a(1).gridName;
referenceGrid = a(1).grid;
referenceGridSize = a(1).gridSize(referenceGridName);
for it = 1 : numel(varargin)
assert(numel(referenceGridName) == numel(varargin{it}(1).gridName), ...
['For getFinestGrid, the gridName of all objects must be equal', ...
'. Maybe gridJoin() does what you want?']);
comparisonGridSize = varargin{it}.gridSize(referenceGridName);
for jt = 1 : numel(referenceGridName)
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(jt) > referenceGridSize(jt)
referenceGrid{jt} = comparisonGrid;
referenceGridSize(jt) = comparisonGridSize(jt);
end
end
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};
assert(tempGrid1(1) == tempGrid2(1), 'Grids must have same domain for gridJoin')
assert(tempGrid1(end) == tempGrid2(end), 'Grids must have same domain for gridJoin')
if numel(tempGrid1) > numel(tempGrid2)
gridJoined{it} = tempGrid1;
else
gridJoined{it} = tempGrid2;
end
end
end
end
function c = horzcat(a, varargin)
% Before Concatenate, check if all quantites have same grid.
% Furthermore, copies of the input quantites are used for
% concatenation.
if nargin == 1
objCell = {a};
else
objCell = {a, varargin{:}};
end
[fineGrid, fineGridName] = getFinestGrid(objCell{:});
for it = 1 : nargin
objCopyTemp = copy(objCell{it});
objCell{it} = objCopyTemp.changeGrid(fineGrid, fineGridName);
end
assertSameGrid(objCell{:});
c = builtin('horzcat', objCell{:});
end
function c = vertcat(a, varargin)
% Before Concatenate, check if all quantites have same grid.
% Furthermore, copies of the input quantites are used for
% concatenation.
if nargin == 1
objCell = {a};
else
objCell = {a, varargin{:}};
end
[fineGrid, fineGridName] = getFinestGrid(objCell{:});
for it = 1 : nargin
objCopyTemp = copy(objCell{it});
objCell{it} = objCopyTemp.changeGrid(fineGrid, fineGridName);
end
assertSameGrid(objCell{:});
c = builtin('vertcat', objCell{:});
end
function c = cat(a, varargin)
% Before Concatenate, check if all quantites have same grid.
% Furthermore, copies of the input quantites are used for
% concatenation.
if nargin == 1
objCell = {a};
else
objCell = {a, varargin{:}};
end
[fineGrid, fineGridName] = getFinestGrid(objCell{:});
for it = 1 : nargin
objCopyTemp = copy(objCell{it});
objCell{it} = objCopyTemp.changeGrid(fineGrid, fineGridName);
end
assertSameGrid(objCell{:});
c = builtin('cat', objCell{:});
end
function solution = solveAlgebraic(obj, rhs, gridName, objLimit)
%% this method solves
% obj(gridName) == rhs
......@@ -395,9 +528,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
if ~iscell(values)
values = {values};
assert(numel(values) == numel(gridName2Replace), ...
'gridName2Replace and values must be of same size');
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
......@@ -420,7 +553,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
gridName2Replace{end+1} = [gridName2Replace{1}, 'backUp'];
values{1} = [gridName2Replace{1}, 'backUp'];
end
if any(strcmp(values{1}, obj(1).gridName))
if isequal(values{1}, gridName2Replace{1})
% replace with same variable... everything stay the
% same.
newGrid = obj(1).grid;
newGridName = obj(1).gridName;
newValue = obj.on();
elseif any(strcmp(values{1}, obj(1).gridName))
% if for a quantity f(z, zeta) this method is
% called with subs(f, zeta, z), then
% g(z) = f(z, z) results, hence the dimensions z
......@@ -435,11 +574,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
newValue = misc.diagNd(obj.on(newGridForOn), gridIndices);
newGrid = {newGridForOn{gridIndices(1)}, ...
newGridForOn{any(1:1:numel(newGridForOn)) ~= gridIndices(1) ...
& any(1:1:numel(newGridForOn)) ~= gridIndices(2)}};
newGridForOn{1:1:numel(newGridForOn) ~= gridIndices(1) ...
& 1:1:numel(newGridForOn) ~= gridIndices(2)}};
newGridName = {values{1}, ...
obj(1).gridName{any(1:1:numel(obj(1).gridName)) ~= gridIndices(1) ...
& any(1:1:numel(obj(1).gridName)) ~= gridIndices(2)}};
obj(1).gridName{1:1:numel(obj(1).gridName) ~= gridIndices(1) ...
& 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
else
% this is the default case. just grid name is
......@@ -491,101 +630,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
end
% function solution = subs(obj, gridName, values)
% % 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 gridNames or or numerics.
% % In contrast to on() or at(), only some, but not necessarily
% % all variables are evaluated.
% if ~iscell(gridName)
% gridName = {gridName};
% end
% if ~iscell(values)
% values = {values};
% end
% isNumericValue = cellfun(@isnumeric, values);
% if any((cellfun(@(v) numel(v(:)), values)>1) & isNumericValue)
% error('only implemented for one value per grid');
% end
% numericValues = values(isNumericValue);
% if numel(obj(1).gridName) == numel(numericValues)
% % if all grids are evaluated, solution is a double array.
% solution = reshape(obj.on(numericValues), size(obj));
% else
% % evaluate numeric values
% subsGrid = obj(1).grid;
% selectRemainingGrid = false(1, numel(obj(1).grid));
% for currentGridName = gridName(isNumericValue)
% selectGrid = strcmp(obj(1).gridName, currentGridName);
% subsGrid{selectGrid} = values{strcmp(gridName, currentGridName)};
% selectRemainingGrid = selectRemainingGrid | selectGrid;
% end
% newGrid = obj(1).grid(~selectRemainingGrid);
% newGridName = obj(1).gridName(~selectRemainingGrid);
% for it = 1 : numel(values)
% if ~isNumericValue(it) && ~isempty(obj(1).gridName(~selectRemainingGrid))
% newGridName{strcmp(obj(1).gridName(~selectRemainingGrid), gridName{it})} ...
% = values{it};
% end
% end
% % before creating a new quantity, it is checked that
% % newGridName is unique. If there are non-unique
% % gridNames, multiple ones are removed and the finest grid
% % from newGrid is taken.
% [uniqueGridName] = unique(newGridName);
% for kt = 1 : numel(uniqueGridName)
% % select finest grid
% sameGridSelector = strcmp(newGridName, uniqueGridName{kt});
% gridCandidates = subsGrid(sameGridSelector);
% gridCandidate = gridCandidates{1};
% for asdf = 2 : numel(gridCandidates)
% if numel(gridCandidate) < numel(gridCandidates{asdf})
% gridCandidate = gridCandidates{asdf};
% end
% subsGrid(sameGridSelector) = {gridCandidate};
% end
% [newGrid{strcmp(newGridName, uniqueGridName{kt})}] = deal(gridCandidate);
% end
%
% newValue = reshape(obj.on(subsGrid), [cellfun(@numel, newGrid), size(obj)]);
% k = 1;
% while numel(newGridName) > numel(uniqueGridName)
% % elemenate non-diagonal elements
% sameGridSelector = strcmp(newGridName, uniqueGridName{k});
% if sum(sameGridSelector) > 1
% sameGridIndex = 1:numel(newGridName);
% sameGridIndex = sameGridIndex(sameGridSelector);
% newValue = misc.diagNd(newValue, sameGridIndex);
% newGrid = [newGrid(sameGridIndex(1)), ...
% newGrid(all((1:numel(newGrid) ~= sameGridIndex(:) )))];
% newGridName = {newGridName{sameGridIndex(1)}, ...
% newGridName{all((1:numel(newGridName) ~= sameGridIndex(:) ))}};
% end
% k = k+1;
% end
% solution = quantity.Discrete(newValue, 'size', size(obj), ...
% 'gridName', newGridName, 'grid', newGrid, 'name', obj(1).name);
%
% end
% end
% at evaluates the given function at the given grid points.
% value = at(obj, grid, index) returns the evaluated
% function as an array with the dimensions:
% dim(value) = [length(z), size(obj)].
% For example, the quantity describes an array valued function
% f(z) = [f1(z), f2(z), f3(z); ...
% f4(z), f5(z), f6(z)],
% i.e., it has the dimensions (2, 3) and a vector Z for the N =
% 11 discretization points Z = [0; 0.1; ...; 1]. By calling the
% evaluate function for this quantity an array of dimensions (11,
% 2, 3) will be returned.
function value = at(obj, point)
value = shiftdim(obj.on(point), 1);
end
......@@ -625,12 +670,16 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
d = d(1);
end
function s = gridSize(obj)
function s = gridSize(obj, myGridName)
% GRIDSIZE returns the size of all grid entries.
if isempty(obj(1).grid)
s = [];
else
s = cellfun('length', obj(1).grid);
if nargin == 1
s = cellfun('length', obj(1).grid);
elseif nargin == 2
s = cellfun('length', obj(1).gridOf(myGridName));
end
end
end
......@@ -667,7 +716,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
set(h, 'WindowStyle', 'Docked');
end
assert(obj.nargin() <= 2);
assert(obj.nargin() <= 2, 'plot only supports quantities with 2 gridNames');
subplotRowIdx = 1:size(obj, 1);
subpotColumnIdx = 1:size(obj, 2);
......@@ -682,7 +731,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if obj.nargin() == 1
plot(...
obj(rowIdx, columnIdx, figureIdx).grid{1}, ...
obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
obj(rowIdx, columnIdx, figureIdx).valueDiscrete );
elseif obj.nargin() == 2
misc.isurf(obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
......@@ -779,6 +828,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
[obj.derivatives] = deal({});
[obj.grid] = deal(myGrid);
end
end
%% math
......@@ -1064,32 +1114,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
P = quantity.Discrete(quantity.Discrete.value2cell(P, [p, q], [n, m]), argin{:});
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};
assert(tempGrid1(1) == tempGrid2(1), 'Grids must have same domain for gridJoin')
assert(tempGrid1(end) == tempGrid2(end), 'Grids must have same domain for gridJoin')
if numel(tempGrid1) > numel(tempGrid2)
gridJoined{it} = tempGrid1;
else
gridJoined{it} = tempGrid2;
end
function myGrid = gridOf(obj, myGridName)
if ~iscell(myGridName)
myGrid = obj(1).grid{obj(1).gridIndex(myGridName)};
else
myGrid = cell(size(myGridName));
for it = 1 : numel(myGrid)
myGrid{it} = obj(1).grid{obj(1).gridIndex(myGridName{it})};
end
end
end
......@@ -1233,8 +1264,8 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if numel(boundaryGridName) == 1
result = result.subs(intGridName, boundaryGridName{1});
else
result = result.subs(intGridName, boundaryGridName{2}) ...
- result.subs(intGridName, boundaryGridName{1});
result = result.subs({intGridName}, boundaryGridName(2)) ...
- result.subs({intGridName}, boundaryGridName(1));
end
end
......@@ -1352,7 +1383,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if ip.isDefault('variable') && x.nargin > 1
error(['For x quantities wrt. two independent' ...
' variables, the name of the variable has to be' ...
' parametrized explicitly']);
'parametrized explicitly']);
end
......@@ -1636,10 +1667,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% #TODO insert code here if some properties should not be copied.
[cpObj.exportData] = deal(export.Data.empty);
end
end
end % methods (Access = protected)
methods ( Static, Access = protected )
function doNotCopy = doNotCopyPropertiesName()
% #DEPRECATED
doNotCopy = {'valueContinuous', 'valueDiscrete', 'exportData', 'derivatives'};
......
......@@ -404,7 +404,8 @@ classdef Symbolic < quantity.Function
assert(~any( diff(cellfun(@(c) length(c), [a.grid, varargin{k}.grid])) ), 'Grids must have the same size');
end
c = builtin('vertcat', a, varargin{:});
c = vertcat@quantity.Discrete(a, varargin{:});
%c = builtin('vertcat', a, varargin{:});
end
......
......@@ -946,6 +946,17 @@ quanZetaZetaReference = misc.diagNd(quan.on({linspace(0, 1, 41), linspace(0, 1,
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
myTestArray4d = rand(11, 21, 31, 41, 1, 2);
quan4d = quantity.Discrete(myTestArray4d, 'gridName', {'z', 'zeta', 'eta', 'beta'}, '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)));
%
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)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment