From e37e9a541d195d2a3bc05d55a8f0c201f5aa63f7 Mon Sep 17 00:00:00 2001
From: Jakob Gabriel <jakob.gabriel@uni-ulm.de>
Date: Tue, 11 Aug 2020 16:20:28 +0200
Subject: [PATCH] quantity.Discrete: simplified implementation of power() +
wrote unittests (also for quantity.Symbolic)
---
+quantity/Discrete.m | 22 ++++------------------
+unittests/+quantity/testDiscrete.m | 23 +++++++++++++++++++++++
+unittests/+quantity/testSymbolic.m | 25 ++++++++++++++++++++++++-
3 files changed, 51 insertions(+), 19 deletions(-)
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index b1ae25f..96490e0 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -1609,24 +1609,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end % sqrtm()
function P = power(obj, p)
- % a.^p elementwise power implemented by multiplication
- % #todo unittest required
- P = quantity.Discrete.zeros(size(obj), obj(1).domain);
-
- if round(p) == p
- % power of integer order:
- for k = 1:numel(obj)
- P(k) = obj(k)^p;
- end
- else
- % power of not integer order
- for k = 1:numel(obj)
- P(k) = obj(k).copy();
- P(k).valueDiscrete = obj(k).valueDiscrete.^p;
- P(k).name = obj(k).name + ".^" + num2str(p);
- end
- end
- end % power
+ % a.^p elementwise power
+ P = quantity.Discrete(obj.on().^(p), obj(1).domain, ...
+ "name", obj(1).name + ".^{" + num2str(p) + "}");
+ end % power()
function P = mpower(a, p)
% Matrix power a^p is matrix or scalar a to the power p.
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index 1a1014f..1cf32c1 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -4,6 +4,29 @@ function [tests] = testDiscrete()
tests = functiontests(localfunctions);
end
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Discrete(sin(z.grid * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Discrete + zeta.Discrete, A + 0*zeta.Discrete; ...
+ z.Discrete+zeta.Discrete, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+ 'AbsTol', 1e-15);
+end % testPower()
+
function setupOnce(testCase)
syms z zeta t sy
testCase.TestData.sym.z = z;
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 5206643..25f1b37 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -13,12 +13,35 @@ B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
B2 = [zeta.Symbolic + 0*z.Symbolic, A + 0*zeta.Symbolic; ...
z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
tc.verifyEqual(B.on([z, zeta]), B2.on([z, zeta]));
-% TODO: fix this error!
+% TODO: fix this error! See issue #41
% this problem is maybe caused as input arguments of B2(1).valueContinuous is somehow ordered
% differently compared to B2(2:4).valueContinuous or all B(:).valueContinuous,
% i.e. B2(1).valueContinuous is (zeta, z) instead of (z, zeta).
end % testValueContinuousBug()
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Symbolic(sin(sym("z") * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
+ z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+ 'AbsTol', 1e-15);
+end % testPower()
+
function testMPower(tc)
%% scalar case
z = quantity.Domain("z", linspace(0, 1, 11));
--
GitLab