Skip to content
GitLab
Commit 605efca2 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

tmp

parent 6da3dcdf
Hide whitespace changes
Inline Side-by-side
+quantity/Discrete.m
View file @ 605efca2
......@@ -175,11 +175,19 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% --- getter and setters ---
%---------------------------
function gridName = get.gridName(obj)
gridName = {obj.domain.name};
if isempty(obj.domain)
gridName = {};
else
gridName = {obj.domain.name};
end
end
function grid = get.grid(obj)
grid = {obj.domain.grid};
if isempty(obj.domain)
grid = {};
else
grid = {obj.domain.grid};
end
end
function itIs = isConstant(obj)
......@@ -190,23 +198,6 @@ 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.
......@@ -430,39 +421,19 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% 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));
[sortedDomain, I] = obj(1).domain.sort(varargin{:});
[obj.domain] = deal(sortedDomain);
% assign the new grid names
[obj.gridName] = deal(sortedNames);
% permute the value discrete
for k = 1:numel(obj)
obj(k).valueDiscrete = permute(obj(k).valueDiscrete, I);
end
end
end% sort()
function c = horzcat(a, varargin)
%HORZCAT Horizontal concatenation.
% [A B] is the horizontal concatenation of objects A and B
......@@ -733,15 +704,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
myParser.addParameter('AbsTol', 1e-6);
myParser.parse(varargin{:});
variableGrid = myParser.Results.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);
......@@ -1006,7 +977,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
......@@ -1449,7 +1420,7 @@ 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
......@@ -1470,46 +1441,48 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
parameters.name = [a(1).name, ' ', b(1).name];
parameters.figureID = a(1).figureID;
domainA = a(1).domain;
domainB = b(1).domain;
% select finest grid from a and b
joinedDomain = gridJoin(domainA, domainB);
newDomainName = {joinedDomain.name};
% gridJoin combines the grid of a and b while using the finer
% of both.
finalGridOrder = [domainA(idx.A.common), domainA(~idx.A.common), domainB(~idx.B.common)];
% 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(newDomainName)
for i = 1 : numel(joinedDomain)
parameters.domain(i) = finalGridOrder(i);
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(finalGridOrder(i).name);
idxB = domainB.gridIndex(finalGridOrder(i).name);
idxA = domainA.gridIndex(joinedDomain(i).name);
idxB = domainB.gridIndex(joinedDomain(i).name);
if idxA
newDomainA = [newDomainA, finalGridOrder(i)];
newDomainA = [newDomainA, joinedDomain(i)]; %#ok<AGROW>
end
if idxB
newDomainB = [newDomainB, finalGridOrder(i)];
newDomainB = [newDomainB, joinedDomain(i)]; %#ok<AGROW>
end
end
parameters = misc.struct2namevaluepair(parameters);
%%
% 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{:});
......@@ -1674,11 +1647,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()
......@@ -2148,17 +2123,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
value = reshape(v, [obj(1).gridSize(), size(obj)]);
if myDomain ~= obj(1).domain
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.
% the quantity, ...
if obj.isConstant
% ... duplicate the constant value on the desired
% domain
value = repmat(value, [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{:});
elseif obj.isConstant
value = repmat(value, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
end
end
end
......@@ -2349,10 +2328,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
......@@ -2367,7 +2342,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
s.PropertyList.valueDiscrete = ...
[sprintf('%ix', gridSize(obj), size(obj)) sprintf('\b')];
end
end
end % getPropertyGroups
end % methods (Access = protected)
methods ( Static, Access = protected )
......
+quantity/Domain.m
View file @ 605efca2
classdef Domain
classdef Domain < matlab.mixin.CustomDisplay
%DOMAIN class to describes a range of values on which a function can be defined.
% todo:
......@@ -15,8 +15,6 @@ classdef Domain
% 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)
......@@ -24,6 +22,10 @@ classdef Domain
lower; % lower bound of the domain
upper; % upper bound of the domain
end
properties (Dependent, Hidden)
ones;
end
methods
function obj = Domain(varargin)
......@@ -45,6 +47,10 @@ classdef Domain
end
end
function i = get.ones(obj)
i = ones(size(obj.grid));
end
function n = get.n(obj)
n = length(obj.grid);
end
......@@ -86,12 +92,12 @@ classdef Domain
l = ~(obj == obj2);
end
function [joinedDomain] = gridJoin(obj1, obj2)
function [joinedDomain, index] = 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');
[joinedGrid, index] = unique({obj1.name, obj2.name}, 'stable');
joinedDomain = quantity.Domain.empty;
......@@ -115,10 +121,9 @@ classdef Domain
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
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
......@@ -192,8 +197,83 @@ classdef Domain
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);
for k = 1:length(myGrids)
newObj(k) = quantity.Domain('grid', myGrids{k}, ...
'name', sortedNames{k}); %#ok<AGROW>
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 g = defaultGrid(gridSize, name)
......
+quantity/Symbolic.m
View file @ 605efca2
......@@ -576,28 +576,28 @@ classdef Symbolic < quantity.Function
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 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)
......@@ -624,7 +624,7 @@ classdef Symbolic < quantity.Function
return;
end
parameters = combineGridGridNameVariable(A, B);
parameters.domain = gridJoin(A.domain, B.domain);
parameters.name = [A(1).name, ' ', B(1).name];
parameters = misc.struct2namevaluepair(parameters);
......@@ -644,7 +644,7 @@ classdef Symbolic < quantity.Function
parameters = struct(a(1));
end
parameters = combineGridGridNameVariable(a, b, parameters);
parameters.domain = gridJoin(a.domain, b.domain);
parameters.name = [a(1).name, '+' , b(1).name];
parameters = misc.struct2namevaluepair(parameters);
......
+unittests/+quantity/testDiscrete.m
View file @ 605efca2
......@@ -24,12 +24,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('grid', t, 'name', 'a');
domB = quantity.Domain('grid', t, 'name', 'b');
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('grid', linspace(0, 1, 21), 'name', 'z');
zeta = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
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);
......@@ -82,12 +107,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');
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());
......@@ -115,19 +142,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('grid', linspace(0, 1, 7), 'name', 'z');
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('grid', linspace(0, 1, 7), 'name', 'z');
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 +166,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('grid', linspace(0, 1, 7), 'name', 'z');
zeta = quantity.Domain('grid', linspace(0, 1, 7), 'name', 'zeta');
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
......@@ -198,27 +227,33 @@ 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('grid', linspace(0, 1, 3), 'name', 'z');
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('grid', linspace(0, 1, 4), 'name', 'zeta');
s2dDiag = quantity.Discrete(cat(4, ...
cat(3, 1 + z.grid .* zeta.grid