Commit 7393369c authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

quantity.Discrete: rearranged on() such that the easy cases are faster, since...

quantity.Discrete: rearranged on() such that the easy cases are faster, since the permutation of the grid / domain can be skipped
parent e1553a64
......@@ -157,9 +157,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
function itIs = isConstant(obj)
% the quantity is interpreted as constant if it has no grid or
% it has a grid that is only one point.
% it has a grid that is only defined at one point.
itIs = isempty(obj(1).domain);
end
end % isConstant()
function doNotCopy = get.doNotCopy(obj)
doNotCopy = obj.doNotCopyPropertiesName();
end
......@@ -357,65 +358,83 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
if isempty(obj)
value = zeros(size(obj));
else
if nargin == 2
% case 1: a domain is specified by myDomain or by
% myDomain as a cell-array with grid entries
if iscell(myDomain) || isnumeric(myDomain)
if nargin == 1
% case 0: no domain was specified, hence the value is requested
% on the default grid defined by obj(1).domain.
value = obj.obj2value(obj(1).domain);
elseif nargin == 2 && (iscell(myDomain) || isnumeric(myDomain))
% case 1: a domain is specified by myDomain == grid(-cell)
myDomain = misc.ensureIsCell(myDomain);
newGrid = myDomain;
myDomain = quantity.Domain.empty();
if obj(1).isConstant()
gridNames = repmat({''}, length(newGrid));
else
gridNames = {obj(1).domain.name};
end
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
value = reshape(obj.obj2value(myDomain), ...
[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
% more involved:
if nargin == 2
% case 2: a domain is specified by a myDomain = domain-array
% nothing has to be done to obtain the domain.
elseif nargin == 3
% case 3: a domain is specified by a grid and a grid
% name. Then, the first input parameter is the grid,
% i.e., myGrid = myDomain and the second is the grid
% name.
% Since the order of the domains is not neccessarily equal to the
% order in obj(1).domain, this is more involved:
myDomain = misc.ensureIsCell(myDomain);
% assert(all(cellfun(@(v)isvector(v), myDomain)), 'The cell entries for a new grid have to be vectors')
gridNames = misc.ensureIsCell(gridNames);
assert(all(cellfun(@(v)isvector(v), myDomain)), ...
'The cell entries for a new grid have to be vectors')
assert(iscell(gridNames), ...
'The gridNames parameter must be cell array')
assert(all(cellfun(@ischar, gridNames)), ...
'The gridNames must be strings')
newGrid = myDomain;
myDomain = quantity.Domain.empty();
if obj(1).isConstant()
gridNames = repmat({''}, length(newGrid));
else
gridNames = {obj(1).domain.name};
end
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
else
error('wrong number of input arguments')
end
elseif nargin == 3
% case 2: a domain is specified by a grid and a grid
% name. Then, the first input parameter is the grid,
% i.e., myGrid = myDomain and the second is the grid
% name.
myDomain = misc.ensureIsCell(myDomain);
gridNames = misc.ensureIsCell(gridNames);
assert(all(cellfun(@(v)isvector(v), myDomain)), 'The cell entries for a new grid have to be vectors')
assert(iscell(gridNames), 'The gridNames parameter must be cell array')
assert(all(cellfun(@ischar, gridNames)), 'The gridNames must be strings')
newGrid = myDomain;
myDomain = quantity.Domain.empty();
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
else
myDomain = obj(1).domain;
% verify the domain
if obj(1).isConstant
gridPermuteIdx = 1:length(myDomain);
else
assert(numel(myDomain) == numel(obj(1).domain), ...
'Wrong grid for the evaluation of the object');
% compute the permutation index, in order to bring the
% new domain in the same order as the original one.
gridPermuteIdx = obj(1).domain.getPermutationIdx(myDomain);
end
% get the valueDiscrete data for this object. Apply the
% permuted myDomain. Then the obj2value will be evaluated
% in the order of the original domain. The permutation to
% the new order will be done in the next step.
originalOrderedDomain(gridPermuteIdx) = myDomain;
value = obj.obj2value(originalOrderedDomain);
value = permute(reshape(value, [originalOrderedDomain.gridLength, size(obj)]), ...
[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
end
% verify the domain
if obj(1).isConstant
gridPermuteIdx = 1:length(myDomain);
else
assert(numel(myDomain) == numel(obj(1).domain), ['Wrong grid for the evaluation of the object']);
% compute the permutation index, in order to bring the
% new domain in the same order as the original one.
gridPermuteIdx = obj(1).domain.getPermutationIdx(myDomain);
end
% get the valueDiscrete data for this object. Apply the
% permuted myDomain. Then the obj2value will be evaluated
% in the order of the original domain. The permuatation to
% the new order will be done in the next step.
originalOrderedDomain(gridPermuteIdx) = myDomain;
value = obj.obj2value(originalOrderedDomain);
value = permute(reshape(value, [originalOrderedDomain.gridLength, size(obj)]), ...
[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
end
end
end % if isempty(obj)
end % on()
function interpolant = interpolant(obj)
% get the interpolant of the obj;
......@@ -779,21 +798,19 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
function solution = solveDVariableEqualQuantity(obj, varargin)
%% solves the first order ODE
% solves the first order ODE
% dvar / ds = obj(var(s))
% var(s=0) = ic
% for var(s, ic). Herein, var is the (only) continuous variale
% obj.variable. The initial condition of the IVP is a variable
% of the result var(s, ic).
assert(numel(obj(1).gridName) == 1, ...
'this method is only implemented for quanitities with one gridName');
% var(0) = ic
% to obtain var(s, ic) depending on both the argument s and the initial
% condition ic. Herein, obj may only depend on one variable / gridName / ...
% domain.
assert(numel(obj(1).domain) == 1, ...
'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('newGridName', 's');
myParser.addParameter('RelTol', 1e-6);
myParser.addParameter('AbsTol', 1e-6);
myParser.parse(varargin{:});
variableGrid = myParser.Results.variableGrid(:);
......@@ -807,7 +824,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
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);
odeSolution = zeros([myGridSize, numel(obj)]);
for it = 1:numel(obj)
for icIdx = 1:numel(myParser.Results.initialValueGrid)
......@@ -819,13 +835,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
[resultGridPositive, odeSolutionPositive] = ...
ode45(@(y, z) obj(it).on(z), ...
positiveVariableGrid, ...
myParser.Results.initialValueGrid(icIdx), options);
myParser.Results.initialValueGrid(icIdx));
end
if numel(negativeVariableGrid) >1
[resultGridNegative, odeSolutionNegative] = ...
ode45(@(y, z) obj(it).on(z), ...
negativeVariableGrid, ...
myParser.Results.initialValueGrid(icIdx), options);
myParser.Results.initialValueGrid(icIdx));
end
if any(variableGrid == 0)
resultGrid = [flip(resultGridNegative(2:end)); 0 ; resultGridPositive(2:end)];
......@@ -843,10 +859,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% return result as quantity-object
solution = quantity.Discrete(...
reshape(odeSolution, [myGridSize, size(obj)]), ...
'gridName', {myParser.Results.newGridName, 'ic'}, 'grid', ...
{variableGrid, myParser.Results.initialValueGrid}, ...
'domain', [quantity.Domain(myParser.Results.newGridName, variableGrid), ...
quantity.Domain('ic', myParser.Results.initialValueGrid)], ...
'size', size(obj), 'name', ['solve(', obj(1).name, ')']);
end
end % solveDVariableEqualQuantity()
function solution = subs(obj, gridName2Replace, values)
% SUBS substitute variables of a quantity
......@@ -2200,11 +2216,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
value = repmat(v, [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{:});
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{:});
end
end
......
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