Commit 9fe59c24 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

cumInt of Function & Symbolic: only supports string, no char -> updated also unittests.

parent 6c5018c9
......@@ -1132,158 +1132,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
[obj.domain] = deal(newDomainComplete);
end
end % subsDomain()
% 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 or an object-array with
% % quantity.Domain objects 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;
% elseif isempty(obj)
% % if the object is empty, nothing must be done.
% solution = obj;
% else
% % input checks
% 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);
% for k = 1:numel( gridName2Replace )
% if isa(gridName2Replace{k}, 'quantity.Domain')
% gridName2Replace{k} = gridName2Replace{k}.name;
% end
% end
% values = misc.ensureIsCell(values);
% end
%
% assert(numel(values) == numel(gridName2Replace), ...
% 'gridName2Replace and values must be of same size');
%
% % set the newDomain once. If obj(1).domain is a quantity.Equidistant domain, it can
% % not be mixed with other quantity.Domains in an array. Hence, it must be casted to
% % a quantity.Domain. The following strange form of concatenation an empty Domain
% % with the required domain, ensures that the result is an array of quantity.Domain
% % objects.
% newDomain = [quantity.Domain.empty(), obj(1).domain];
%
% % here substitution starts:
% % The first (gridName2Replace{1}, values{1})-pair is
% % replaced. If there are more cell-elements in those inputs
% % then subs() is called again for the remaining pairs
% % (gridName2Replace{2:end}, values{2:end}).
% if ischar(values{1}) || isstring(values{1})
% % if values{1} is a char-array, then the gridName is
% % replaced
% if any(strcmp(values{1}, gridName2Replace(2:end)))
% % in the case if a quantity f(z, zeta) should be
% % substituted like subs(f, {z, zeta}, {zeta, z})
% % this would cause an error, since after the first
% % substituion subs(f, z, zeta) the result would be
% % f(zeta, zeta) -> the 2nd subs(f, zeta, z) will
% % result in f(z, z) and not in f(zeta, z) as
% % intended. This is solved, by an additonal
% % substitution:
% % f.subs(z,zetabackUp).subs(zeta,z).subs(zetabackUp,zeta)
% values{end+1} = values{1};
% gridName2Replace{end+1} = [gridName2Replace{1}, 'backUp'];
% values{1} = gridName2Replace{end};
% end
% if isequal(values{1}, gridName2Replace{1})
% % replace with same variable... everything stays the
% % same.
% % Do not use "return", since, later subs might need to be
% % called recursively!
% 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 and zeta are
% % merged.
% domainIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
% obj(1).domain.gridIndex(values{1})];
% newDomainForOn = obj(1).domain;
% if obj(1).domain(domainIndices(1)).n > obj(1).domain(domainIndices(2)).n
% newDomainForOn(domainIndices(2)) = quantity.Domain(...
% newDomainForOn(domainIndices(2)).name, ...
% newDomainForOn(domainIndices(1)).grid);
% elseif obj(1).domain(domainIndices(1)).n < obj(1).domain(domainIndices(2)).n
% newDomainForOn(domainIndices(1)) = quantity.Domain(...
% newDomainForOn(domainIndices(1)).name, ...
% newDomainForOn(domainIndices(2)).grid);
% end
% newValue = misc.diagNd(obj.on(newDomainForOn), domainIndices);
% newDomain = [newDomainForOn(domainIndices(2)), ...
% newDomainForOn(all(1:1:numel(newDomainForOn) ~= domainIndices(:)))];
% else
% % this is the default case. just grid name is changed.
% %newDomain = obj(1).domain;
% newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
% quantity.Domain(values{1}, ...
% obj(1).domain(obj(1).domain.gridIndex(gridName2Replace{1})).grid);
% newValue = obj.on();
% end
%
% elseif isnumeric(values{1}) && numel(values{1}) == 1
% % if values{1} is a scalar, then obj is evaluated and
% % the resulting quantity looses that spatial grid and
% % gridName
% newDomain = newDomain(~strcmp(gridName2Replace{1}, [newDomain.name]));
% % newGrid is the similar to the original grid, but the
% % grid of gridName2Replace is removed.
% newGridSize = newDomain.gridLength();
% % newGridForOn is the similar to the original grid, but
% % the grid of gridName2Replace is set to values{1} for
% % evaluation of obj.on().
% newGridForOn = obj(1).grid;
% newGridForOn{obj(1).domain.gridIndex(gridName2Replace{1})} = values{1};
% newValue = reshape(obj.on(newGridForOn), [newGridSize, size(obj)]);
%
% elseif isnumeric(values{1}) && numel(values{1}) > 1
% % if values{1} is a double vector, then the grid is
% % replaced.
% newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
% quantity.Domain(gridName2Replace{1}, values{1});
% newValue = obj.on(newDomain);
% else
% error('value must specify a gridName or a gridPoint');
% end % if ischar(values{1}) || isstring(values{1})
% if isempty(newDomain)
% % TODO@Jakob: Is it possible that this case occurrs?
% % If I see it right, the only case where the newDomain is not set is the case
% % where a error is thrown.
% solution = newValue;
% else
% solution = quantity.Discrete(newValue, newDomain, ...
% 'name', obj(1).name);
% end
% if numel(gridName2Replace) > 1
% solution = solution.subs(gridName2Replace(2:end), values(2:end));
% end
% end
% end
function [idx, logicalIdx] = gridIndex(obj, varargin)
warning('DEPRICATED: use quantity.Domain.gridIndex method instead')
......@@ -2409,7 +2257,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
b = [intDomain.upper];
end
I = obj.cumInt(intDomain(1), a(1), b(1));
I = obj.cumInt(intDomain(1).name, a(1), b(1));
if numel(intDomain) > 1
I = I.int(intDomain(2:end), a(2:end), b(2:end));
......@@ -2419,13 +2267,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function result = cumInt(obj, domainName, lowerBound, upperBound)
% CUMINT cumulative integration
% result = cumInt(obj, domain, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and
% the specified 'bounds'. 'domain' must be a gridName of the
% object. 'lowerBound' and 'upperBound' define the boundaries
% of the integration domain. These can be either doubles or
% gridNames. If it is double, the bound is constant. If it is
% a gridName, the bound is variable and the result is a
% result = cumInt(obj, domainName, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and the specified 'bounds'.
% domainName must be a name of one domain of the object.
% 'lowerBound' and 'upperBound' define the boundaries of the integration domain.
% These can be either doubles or names of a domain as a string. If they are doubles,
% the bound is constant. If it is a string, the bound is variable and the result is a
% function dependent of this variable.
arguments
......
......@@ -186,34 +186,40 @@ classdef Function < quantity.Discrete
end
function I = cumInt(obj, domain, lowerBound, upperBound)
function I = cumInt(obj, domainName, lowerBound, upperBound)
% CUMINT cumulative integration
% result = cumInt(obj, domainName, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and the specified 'bounds'.
% domainName must be a name of one domain of the object.
% 'lowerBound' and 'upperBound' define the boundaries of the integration domain.
% These can be either doubles or names of a domain as a string. If they are doubles,
% the bound is constant. If it is a string, the bound is variable and the result is a
% function dependent of this variable.
% input parser since some default values depend on intGridName.
myParser = misc.Parser;
myParser.addRequired('domain', @(d) obj(1).domain.gridIndex(d) ~= 0);
myParser.addRequired('lowerBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l));
myParser.addRequired('upperBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l));
myParser.parse(domain, lowerBound, upperBound)
arguments
obj;
domainName string;
lowerBound;
upperBound;
end
assert(isnumeric(lowerBound) || isstring(lowerBound), ...
"integral bounds must be string or numeric");
assert(isnumeric(upperBound) || isstring(upperBound), ...
"integral bounds must be string or numeric");
% get the desired domain. This step is necessary, as the input
% argument DOMAIN can be either a string or a quantity.Domain
% object.
domain = obj(1).domain.find(domain);
domain = obj(1).domain.find(domainName);
for k = 1:numel(obj)
I_val = numeric.cumIntegral( ...
@(varargin)obj(k).evaluateFunction(varargin{:}), ...
domain.grid);
@(varargin)obj(k).evaluateFunction(varargin{:}), domain.grid);
result = quantity.Discrete(I_val, obj(1).domain);
% int_lowerBound^upperBound f(.) =
% F(upperBound) - F(lowerBound)
I(k) = result.subs(domain, upperBound) ...
- result.subs(domain, lowerBound);
% int_lowerBound^upperBound f(.) = F(upperBound) - F(lowerBound)
I(k) = result.subs(domainName, upperBound) - result.subs(domainName, lowerBound);
if isa(result, 'quantity.Discrete')
result.setName("int(" + obj(1).name + ")");
end
......@@ -221,9 +227,9 @@ classdef Function < quantity.Discrete
I = reshape(I, size(obj));
end
end % cumInt()
end
end % methods (Access = public)
methods (Access = protected)
function v = evaluateFunction(obj, varargin)
......
......@@ -685,65 +685,64 @@ classdef Symbolic < quantity.Function
'name', "{" + obj(1).name + "}^{H}");
end % expm()
function C = cumInt(obj, domain, lowerBound, upperBound)
function C = cumInt(obj, domainName, lowerBound, upperBound)
% CUMINT cumulative integration
% result = cumInt(obj, domain, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and
% the specified 'bounds'. 'domain' must be a gridName of the
% object. 'lowerBound' and 'upperBound' define the boundaries
% of the integration domain. These can be either doubles or
% gridNames. If it is double, the bound is constant. If it is
% a gridName, the bound is variable and the result is a
% function dependent of this variable. Since cumInt can be
% done with the usual Symbolic int(), this method just checks
% the input and calls int().
% result = cumInt(obj, domainName, lowerBound, upperBound)
% performes the integration over 'obj' for the 'domain' and the specified 'bounds'.
% domainName must be a name of one domain of the object.
% 'lowerBound' and 'upperBound' define the boundaries of the integration domain.
% These can be either doubles or names of a domain as a string. If they are doubles,
% the bound is constant. If it is a string, the bound is variable and the result is a
% function dependent of this variable.
%input parser since some default values depend on intGridName.
myParser = misc.Parser;
myParser.addRequired('domain', @(d) obj(1).domain.gridIndex(d) ~= 0);
myParser.addRequired('lowerBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l) );
myParser.addRequired('upperBound', ...
@(l) isnumeric(l) || ischar(l) || isstring(l) );
myParser.parse(domain, lowerBound, upperBound)
arguments
obj;
domainName string;
lowerBound;
upperBound;
end
assert(isnumeric(lowerBound) || isstring(lowerBound), ...
"integral bounds must be string or numeric");
assert(isnumeric(upperBound) || isstring(upperBound), ...
"integral bounds must be string or numeric");
% get the desired domain. This step is necessary, as the input
% argument DOMAIN can be either a string or a quantity.Domain
% object.
domain = obj(1).domain.find(domain);
domainName = obj(1).domain.find(domainName);
% do the integration element-wise
if obj.nargin == numel(domain) && isnumeric(lowerBound) && isnumeric(upperBound)
if obj.nargin == numel(domainName) && isnumeric(lowerBound) && isnumeric(upperBound)
% if the integration has fixed upper and lower bounds, and
% should be applied for all independent variables, use the
% implementation based on function handles, because it is
% much faster:
C = cumInt@quantity.Function(obj, domain, lowerBound, upperBound);
C = cumInt@quantity.Function(obj, domainName, lowerBound, upperBound);
else
try
% try to compute the symbolic integration
I = quantity.Symbolic(int(obj.sym, domain.name), obj(1).domain);
I = quantity.Symbolic(int(obj.sym, domainName.name), obj(1).domain);
C = subs(I, domain.name, upperBound) - ...
subs(I, domain.name, lowerBound);
C = subs(I, domainName.name, upperBound) - ...
subs(I, domainName.name, lowerBound);
catch ex
try
% try to compute the integration by function
% handles
C = cumInt@quantity.Function(obj, domain, lowerBound, upperBound);
C = cumInt@quantity.Function(obj, domainName, lowerBound, upperBound);
catch ex
% last fall back to compute the integration by
% discrete values. This should work always
C = cumInt(quantity.Discrete(obj), domain, lowerBound, upperBound);
C = cumInt(quantity.Discrete(obj), domainName, lowerBound, upperBound);
end
end
end % try-catch
end
end
end
end % if-else
end % cumInt
end % methods (Access = public)
methods (Access = protected)
......
......@@ -57,14 +57,14 @@ f = quantity.Function( @(z) sin( sin( z ) ), z);
F_fn = int(f, z, 0, "z");
F_dk = int(quantity.Discrete(f), z, 0, "z");
testCase.verifyEqual( F_fn.on, F_dk.on, 'AbsTol', 6e-6)
testCase.verifyEqual( F_fn.on, F_dk.on, "AbsTol", 6e-6)
f2 = quantity.Function({@(z)sin(z); @(z)cos(z)}, z);
F2_fn = int(f2, z, "zeta", "z");
F2_dk = int( quantity.Discrete(f2), z, "zeta", "z");
testCase.verifyEqual( F2_fn.on, F2_dk.on, 'AbsTol', 7.5e-6);
testCase.verifyEqual( F2_fn.on, F2_dk.on, "AbsTol", 7.5e-6);
end
......@@ -116,10 +116,10 @@ f1 = @(z) sinh(z * pi);
f2 = @(z) cosh(z * pi);
f3 = @(z) sin(z * pi);
f = quantity.Function({f1, f3 ; f2, f2; f1, f3 }, z, 'name', 'sinhcosh');
f = quantity.Function({f1, f3 ; f2, f2; f1, f3 }, z, "name", "sinhcosh");
F = -f;
%
z = sym('z', 'real');
z = sym("z", "real");
fsym = sym({f1, f3 ; f2, f2; f1, f3});
fsym = symfun(- fsym, z);
......@@ -153,12 +153,12 @@ f1 = @(z) sinh(z * pi);
f2 = @(z) cosh(z * pi);
f3 = @(z) sin(z * pi);
z = sym('z', 'real');
z = sym("z", "real");
fsym = sym({f1, f3 ; f2, f2; f1, f3});
fsym = symfun(fsym * fsym.', z);
f = quantity.Function({f1, f3 ; f2, f2; f1, f3 }, quantity.Domain.defaultDomain(10, "z"),...
'name', 'sinhcosh');
"name", "sinhcosh");
F = f * f.';
......@@ -174,7 +174,7 @@ end
function testInit(testCase)
fun = @(z) z.^2;
q = quantity.Function(fun, quantity.Domain.defaultDomain(19, "z"), 'name', "test");
q = quantity.Function(fun, quantity.Domain.defaultDomain(19, "z"), "name", "test");
verifyEqual(testCase, q.valueContinuous, fun);
verifyEqual(testCase, q.name, "test");
......@@ -185,8 +185,8 @@ a = quantity.Function(@(z,t) z + t, [z, t]);
b = quantity.Function(@(t,z) z + t, [z, t]);
testCase.verifyError(@() quantity.Function(@(z,t,x) z + t, [z,t]), 'quantity:Function:domain')
testCase.verifyError(@() quantity.Function({ @(z,t) z + t, @(z) z, @(z,t,x) z + t + x}, [z,t]), 'quantity:Function:domain')
testCase.verifyError(@() quantity.Function(@(z,t,x) z + t, [z,t]), "quantity:Function:domain")
testCase.verifyError(@() quantity.Function({ @(z,t) z + t, @(z) z, @(z,t,x) z + t + x}, [z,t]), "quantity:Function:domain")
end
......
This diff is collapsed.
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