Commit 3e1b9829 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Fixed minor bugs to increase runtime

changed misc.diag_nd to use a for loop on the domain, which will be collapsed.

Added misc.rat to compute rational approximations of real numbers

Added quantity.EquidistantDomain to simplify the initialization of a domain with a certain step size

Fixed signals/SmoothStep to use spline wise interpolation instead of symbolic computation
parent 24553886
......@@ -29,7 +29,7 @@ function value = diagNd(value, dimensionsToBeDiagonalized)
% diagonalization are at the beginning
value = permute(value, [dimensionsToBeDiagonalized, ...
originalOrder(all(originalOrder ~= dimensionsToBeDiagonalized(:)))]);
% the result contains the diagonalized dimension first, then the
% remaining
tempValueSize = size(value);
......@@ -60,4 +60,4 @@ function result = allEqual(a)
result = (a{1} == a{2});
end
end
\ No newline at end of file
end
......@@ -36,22 +36,10 @@ if valueSize(1) ~= valueSize(2)
'be calculated must be quadratic but size(value, 1) ~= size(value, 2)']);
end
% valueDiagSize is initialised and written in such weired way in the
% following to ensure that it a vector with at least 2 elements, even if
% the variable value is 2D. If this is not done, the final reshape wont
% work.
valueDiagSize = ones(1:max((numel(valueSize)-1), 2));
valueDiagSize(1:numel(valueSize([1, 3:end]))) = valueSize([1, 3:end]);
% The implemenetation with num2cell and diag(blkdiag(...)) is terribly slow
% because of blkdiag. Since it is slower than for-loops and diag, for-loops
% are used instead.
% valueCell = num2cell(value, 1:2); % valueDiagUnshaped = diag(blkdiag(valueCell{:}));
valueUnshaped = reshape(value, valueSize(1), valueSize(2), []);
valueDiagUnshaped = zeros(valueSize(1), prod(valueSize(3:end)));
for it = 1:prod(valueSize(3:end))
valueDiagUnshaped(:,it) = diag(valueUnshaped(:, :, it));
% just do an iteration over the dimensions for the diagonal evaluation.
valueDiag = zeros([valueSize(2:end), 1]);
for i = 1:valueSize(1)
valueDiag(i,:) = value(i,i,:);
end
valueDiag = reshape(valueDiagUnshaped, valueDiagSize);
end
function [approx] = rat(X, tol)
%RAT_K Rational approximation.
% approx = RAT(X,tol) returns a cell with the rational approximations of
% real valued X.
% X is approximated by
%
% 1
% d1 + ----------------------------------
% 1
% d2 + -----------------------------
% 1
% d3 + ------------------------
% 1
% d4 + -------------------
% 1
% d5 + --------------
% 1
% d6 + ---------
% 1
% d7 + ----
% d8
%
% in order to chose an approximation of an required tolerance, all
% approximation steps will be returned. So APPROX is a cell of the same
% dimension as X. Each entry in X is a struct-array with the fields
% N for the nominator, D for the denominator, j for the index in X and
% error for the absolute error of this approximation.
% For further details, see matlab function rat.
%
% see rat
if nargin < 2
tol = 1.e-6*norm(X(isfinite(X)),1);
end
assert( all( isreal(X(:)) & isfinite(X(:))));
for j = 1:numel(X)
k = 0;
C = [1 0; 0 1];
x = X(j);
while 1
k = k+1;
d = round(x);
x = x - d;
C = [C*[d;1] C(:,1)];
xError = abs(C(1,1)/C(2,1) - X(j));
o.N = C(1,1)/sign(C(2,1));
o.D = abs(C(2,1));
o.j = j;
o.error = xError;
rationals(k) = o;
if (x==0) || (xError <= max(tol,eps(X(j))))
break
end
x = 1/x;
end
approx{j} = rationals(:);
end
end
......@@ -91,21 +91,21 @@ classdef interpolant
end
function value = evaluate(obj, varargin)
% Example:
% test_2p5_1_1 = testInterpolant.evaluate(2.5, 1:5, 1)
%
% INPUT PARAMETERS:
% NUMERIC.INTERPOLANT obj : interpolant object
% VARARGIN varargin : sample points for which the values are
% interpolated. Expected to be a list of
% arrays
% OUTPUT:
% ARRAY value : values at the points specified by varargin
% Example:
% test_2p5_1_1 = testInterpolant.evaluate(2.5, 1:5, 1)
%
% INPUT PARAMETERS:
% obj : interpolant object VARARGIN
% varargin : sample points for which the values are
% interpolated. Expected to be a list of arrays
% OUTPUT:
% value : values at the points specified by varargin
% evaluate is used to get the value of the data at certain
% points specified by varargin. It is assumed that
% numel(varargin) == ndims(valueOriginal)
reducedGrid = varargin(obj.reducedDimension);
value = permute(obj.reducedGriddedInterpolant(reducedGrid), obj.permutationVector);
interpolant = obj.reducedGriddedInterpolant(reducedGrid);
value = permute(interpolant, obj.permutationVector);
end
function value = eval(obj, varargin)
......
......@@ -234,7 +234,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
methods (Access = public)
function [d, I, d_size] = compositionDomain(obj, domainName)
function d = compositionDomain(obj, domainName)
assert(isscalar(obj));
......@@ -243,17 +243,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% the evaluation of obj.on( compositionDomain ) is done by:
d_size = size(d);
% 1) vectorization of the n-d-grid: compositionDomain
d = d(:);
% 2) then it is sorted in ascending order
[d, I] = sort(d);
% verify the domain to be monotonical increasing
deltaCOD = diff(d);
assert(misc.alln(deltaCOD >= 0), 'The domain for the composition f(g(.)) must be monotonically increasing');
d = quantity.Domain(domainName, d);
% vectorization of the n-d-grid: compositionDomain
d = quantity.Domain(domainName, d(:));
end
function obj_hat = compose(obj, g, optionalArgs)
......@@ -276,11 +267,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
'Composition of functions: The domains for the composition must be equal. A grid join is not implemented yet.');
assert( any( logOfDomain ) )
% 2) get the composition domain:
% get the composition domain:
% For the argument y of a function f(y) which should be
% composed by y = g(z,t) a new dommain will be created on the
% basis of evaluation of g(z,t).
[composeOnDomain, I, domainSize] = ...
composeOnDomain = ...
g.compositionDomain(myCompositionDomain.name);
% check if the composition domain is in the range of definition
......@@ -290,14 +281,20 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
'The composition domain is not a subset of obj.domain! The missing values will be extrapolated.');
end
% 3) evaluation on the new grid:
% evaluation on the new grid:
% the order of the domains is important. At first, the
% domains which will not be replaced are taken. The last
% domain must be the composed domain. For example: a function
% f(x, y, z, t), where y should be composed with g(z, t) will
% be resorted to f_(x, z, t, y) and then evaluated with y =
% g(z,t)
% #TODO: optimize the memory consumption of this function.
% 1) only consider the unqiue grid points in evaluationDomain
% 2) do the conversion of the evaluationDomain directly to
% the target domain.
evaluationDomain = [originalDomain( ~logOfDomain ), composeOnDomain ];
newValues = obj.on( evaluationDomain );
% reshape the new values into a 2-d array so that the first
......@@ -306,20 +303,16 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
sizeOldDomain = prod( [originalDomain( ~logOfDomain ).n] );
sizeComposeDomain = composeOnDomain.gridLength;
newValues = reshape(newValues, [sizeOldDomain, sizeComposeDomain]);
% 4) reorder the computed values, dependent on the sort
% position fo the new domain
newValues(:,I) = newValues(:,:);
% 5) rearrange the computed values, to have the same dimension
%rearrange the computed values, to have the same dimension
% as the required domain
% *) consider the domain
% consider the domain
% f( z, g(z,t) ) = f(z, g(zeta,t) )|_{zeta = z}
tmpDomain = [originalDomain( ~logOfDomain ), g(1).domain ];
% newValues will be reshaped into the form
% f(z, t, zeta)
newValues = reshape( newValues, [tmpDomain.gridLength, 1] );
% *) now the common domains, i.e., zeta = z must be merged:
% now the common domains, i.e., zeta = z must be merged:
% for this, find the index of the common domain in list of
% temporary combined domain
......@@ -329,10 +322,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if ~isempty(intersectDomain)
idx = 1:length(tmpDomain);
idxCommon = idx(strcmp({tmpDomain.name}, intersectDomain));
% take the diagonal values of the common domain, i.e., z = zeta
newValues = misc.diagNd( newValues, idxCommon );
logCommon = strcmp({tmpDomain.name}, intersectDomain);
% take the diagonal values of the common domain, i.e., z = zeta
% use the diag_nd function because it seems to be faster
% then the diagNd function, although the values must be
% sorted.
% #TODO: Rewrite the diagNd function, using for loops, in order to be as fast as diag_nd
newValues = permute( newValues, [idx(logCommon), idx(~logCommon)]);
newValues = misc.diag_nd(newValues);
end
% *) build a new valueDiscrete on the correct grid.
......@@ -366,10 +364,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
value = obj.obj2value(obj(1).domain);
elseif nargin == 2 && (iscell(myDomain) || isnumeric(myDomain))
% case 1: a domain is specified by myDomain == grid(-cell)
% case 1: a domain is specified by myDomain as agrid
myDomain = misc.ensureIsCell(myDomain);
newGrid = myDomain;
myDomain = quantity.Domain.empty();
if obj(1).isConstant()
gridNames = repmat({''}, length(newGrid));
......@@ -377,11 +374,14 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
gridNames = {obj(1).domain.name};
end
% initialize the new domain
clear('myDomain');
myDomain(1:length(newGrid)) = quantity.Domain();
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
value = reshape(obj.obj2value(myDomain), ...
[myDomain.gridLength, size(obj)]);
[myDomain.gridLength, size(obj)]);
else
% Since in the remaining cases the order of the domains is not
% neccessarily equal to the order in obj(1).domain, this is
......@@ -903,8 +903,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
gridName2Replace = {gridName2Replace.name};
elseif nargin == 3
gridName2Replace = misc.ensureIsCell(gridName2Replace);
for k = 1:numel( gridName2Replace )
if isa(gridName2Replace{k}, 'quantity.Domain')
gridName2Replace{k} = gridName2Replace{k}.name;
end
end
values = misc.ensureIsCell(values);
end
......@@ -1341,28 +1345,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
for it = 1 : numel(obj)
newObj(it).valueDiscrete = obj(it).on(newDomain);
end
end
function [lowerBounds, upperBounds] = getBoundaryValues(obj, varargin)
% GETBOUNDARYVALUES returns the boundary values of the grid
% [lowerBounds, upperBounds] = getBoundaryValues(obj, gridName)
% returns the boundary values of the grid, specified by
% gridName
%
% [lowerBounds, upperBounds] = getBoundaryValues(obj) returns
% the boudary values of all grids defined for the object.
grids = obj(1).grid(obj(1).domain.gridIndex(varargin{:}));
lowerBounds = zeros(1, numel(grids));
upperBounds = zeros(1, numel(grids));
for k = 1:numel(grids)
lowerBounds(k) = grids{k}(1);
upperBounds(k) = grids{k}(end);
end
end
end
end
%% math
......@@ -1681,8 +1664,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% xNorm = sqrt(int_0^1 x.' * weight * x dz).
% The integral domain is specified by integralGrid.
if obj.nargin > 1
assert(ischar(integralGridName), 'integralGrid must specify a gridName as a char');
if nargin > 1
integralGridName = misc.ensureIsCell(integralGridName);
assert(ischar([integralGridName{:}]), 'integralGrid must specify a gridName as a char');
else
integralGridName = obj(1).gridName;
end
......@@ -1693,8 +1677,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if obj.nargin == 1 && all(strcmp(obj(1).gridName, integralGridName))
xNorm = sqrtm(on(int(obj.' * myParser.Results.weight * obj), integralGridName));
else
xNorm = sqrtm(int(obj.' * myParser.Results.weight * obj, integralGridName));
[xNorm.name] = deal(['||', obj(1).name, '||_{L2}']);
xNorm = sqrtm( int(obj.' * myParser.Results.weight * obj, integralGridName) );
if isa(xNorm, 'quantity.Discrete')
[xNorm.name] = deal(['||', obj(1).name, '||_{L2}']);
end
end
end % l2norm()
......@@ -1926,69 +1912,40 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
return
end
intGridName = obj(1).gridName;
[a, b] = obj.getBoundaryValues;
if nargin == 1 || nargin == 3
% obj.int() -> integrate over all dimensions of this
% quantity.
intDomain = obj(1).domain;
elseif nargin == 2
% obj.int(z) OR obj.int(z, a, b)
% integrate over the domain
if isa(varargin{1}, 'quantity.Domain')
intDomain = varargin{1};
else
intDomain = obj(1).domain.find( varargin{1} );
end
elseif nargin == 4
% obj.int(z, <domain>, lowerBound, upperBOund)
%
I = cumInt( obj, varargin{:});
return;
end
if nargin == 1
% see default settings.
elseif nargin == 2 || nargin == 4
% (obj, z) OR (obj, z, a, b)
intGridName = varargin{1};
elseif nargin == 3
if nargin == 3
% (obj, a, b)
a = varargin{1};
b = varargin{2};
elseif nargin == 1 || nargin == 2
a = [intDomain.lower];
b = [intDomain.upper];
end
if iscell(intGridName) && numel(intGridName) > 1
obj = obj.int(intGridName{1}, a(1), b(1));
I = obj.int(intGridName{2:end}, a(2:end), b(2:end));
return
end
[idxGrid, isGrid] = obj(1).domain.gridIndex(intGridName);
if nargin == 4
I = cumInt( obj, varargin{:});
return;
%
% assert(all((varargin{2} == a(idxGrid)) & (varargin{3} == b(idxGrid))), ...
% 'only integration from beginning to end of domain is implemented');
end
%% do the integration over one dimension
I = obj.copy();
% get index of the variables that should be considered for
% integration
I = obj.cumInt(intDomain(1), a(1), b(1));
[domainIdx{1:obj(1).nargin}] = deal(':');
currentGrid = obj(1).grid{ idxGrid };
domainIdx{idxGrid} = find(currentGrid >= a(idxGrid) & currentGrid <= b(idxGrid));
myGrid = currentGrid(domainIdx{idxGrid});
for kObj = 1:numel(obj)
J = numeric.trapz_fast_nDim(myGrid, obj(kObj).atIndex(domainIdx{:}), idxGrid);
% If the quantity is integrated, so that the first
% dimension collapses, this dimension has to be shifted
% away, so that the valueDiscrete can be set correctly.
if size(J,1) == 1
J = shiftdim(J);
end
% create result:
I(kObj).valueDiscrete = J;
I(kObj).name = ['int(' obj(kObj).name ')'];
end
if all(isGrid)
newDomain = quantity.Domain();
else
newDomain = quantity.Domain(obj(kObj).gridName{~isGrid}, ...
obj(kObj).grid{~isGrid});
if numel(intDomain) > 1
I = I.int(intDomain(2:end), a(2:end), b(2:end));
return
end
[I.domain] = deal(newDomain);
end
function result = cumInt(obj, domain, lowerBound, upperBound)
......@@ -2023,7 +1980,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% int_lowerBound^upperBound f(.) =
% F(upperBound) - F(lowerBound)
result = result.subs({domain}, upperBound) ...
result = result.subs( {domain}, upperBound) ...
- result.subs({domain}, lowerBound);
if isa(result, 'quantity.Discrete')
result.setName(deal(['int(', obj(1).name, ')']));
......@@ -2055,8 +2012,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
% combine both domains with finest grid
joinedDomain = gridJoin(A(1).domain, B(1).domain);
joinedDomain = join(A(1).domain, B(1).domain);
[aDiscrete] = A.expandValueDiscrete(joinedDomain);
[bDiscrete] = B.expandValueDiscrete(joinedDomain);
......
classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
%DOMAIN class to describes a range of values on which a function can be defined.
% todo:
......@@ -16,10 +16,11 @@ classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
% a speaking name for this domain; Should be unique, so that the
% domain can be identified by the name.
name char;
n; % number of discretization points == gridLength
end
properties (Dependent)
n; % number of discretization points == gridLength
lower; % lower bound of the domain == grid(1)
upper; % upper bound of the domain == grid(end)
end
......@@ -35,6 +36,7 @@ classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
if nargin >= 1
obj.grid = grid(:);
obj.name = name;
obj.n = length(grid);
end
end % Domain()
......@@ -42,14 +44,6 @@ classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
i = ones(size(obj.grid));
end
function n = get.n(obj)
if isempty(obj.grid)
n = [];
else
n = length(obj.grid);
end
end
function lower = get.lower(obj)
% minimum value of grid = grid(1), since grid must be ascending
lower = obj.grid(1);
......@@ -139,7 +133,7 @@ classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
[joinedGrid, index] = unique({domain1.name, domain2.name}, 'stable');
joinedDomain = quantity.Domain.empty;
joinedDomain(1 : numel(joinedGrid)) = quantity.Domain();
% #todo@domain: the gird comparison seems to be very
% complicated
......@@ -164,8 +158,8 @@ classdef Domain < handle & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
tempDomain1 = domain1(index1);
tempDomain2 = domain2(index2);
assert(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
assert( numeric.near( tempDomain1.lower, tempDomain2.lower), 'Grids must have same domain for gridJoin')
assert( numeric.near( tempDomain1.upper, tempDomain2.upper), 'Grids must have same domain for gridJoin')
if tempDomain1.gridLength > tempDomain2.gridLength
joinedDomain(i) = tempDomain1;
......
......@@ -3,21 +3,34 @@ classdef EquidistantDomain < quantity.Domain
%definition of a function. The discretization points are equally
%distributed over the domain.
properties
Property1
properties (Access = protected)
stepSize;
end
methods
function obj = untitled(inputArg1,inputArg2)
%UNTITLED Construct an instance of this class
% Detailed explanation goes here
obj.Property1 = inputArg1 + inputArg2;
function obj = EquidistantDomain(name, lower, upper, optArgs)
arguments
name,
lower,
upper,
optArgs.N = [];
optArgs.stepSize = [];
end
if ~isempty(optArgs.stepSize)
grid = lower : optArgs.stepSize : upper;
elseif ~isempty(optArgs.N)
grid = linspace(lower, upper, optArgs.N);
else
error('Either number of discretization points or the step size must be defined!')
end
obj@quantity.Domain(name, grid);
obj.stepSize = optArgs.stepSize;
end
function outputArg = method1(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
function d = quantity.Domain(obj)
d = arrayfun(@(o) quantity.Domain(o.name, o.grid), obj);
end
end
end
......
......@@ -87,8 +87,8 @@ classdef Function < quantity.Discrete
% otherwise the function has to be evaluated on the new
% domain
value = cell(size(obj));
ndGrd = myDomain.ndgrid;
for k = 1:numel(obj)
ndGrd = myDomain.ndgrid;
tmp = obj(k).evaluateFunction( ndGrd{:} );
value{k} = tmp(:);
end
......
......@@ -87,7 +87,6 @@ classdef Symbolic < quantity.Function
% if ~all(strcmp(sort({myDomain.name}), sort(variableName)))
% error('If the gridName(s) are set explicitly, they have to be the same as the variable names!')
% end
parentVarargin = [{fun}, varargin(:)', {'domain'}, {myDomain}];
end
% call parent class
......@@ -99,21 +98,27 @@ classdef Symbolic < quantity.Function
obj = quantity.Symbolic.empty(size(obj));
end
mySymbolicVariable = ...
quantity.Symbolic.getVariable(valueContinuous);
% special properties
for k = 1:numel(valueContinuous)
obj(k).variable = ...
quantity.Symbolic.getVariable(valueContinuous); %#ok<AGROW>
obj(k).variable = mySymbolicVariable; %#ok<AGROW>
obj(k).valueSymbolic = symb{k}; %#ok<AGROW>
obj(k).symbolicEvaluation = ...
variableParser.Results.symbolicEvaluation; %#ok<AGROW>
end
assert( ~( ~isempty( mySymbolicVariable ) && isempty( myDomain ) ), ...
'quantity:Symbolic:Initialisation', ...
'If a quantity.Symbolic is initialized with a symbolic variable, a domain is required.');
end
end
function itIs = isConstant(obj)
itIs = true;
for k = 1:numel(obj)
itIs = itIs && isempty(symvar(obj(k).sym));
itIs = itIs && isempty( obj(k).variable );
if ~itIs
break