Merge remote-tracking branch 'remotes/origin/Branch_25807034' into...

Merge remote-tracking branch 'remotes/origin/Branch_25807034' into 28-domain-argument-as-required-argument
-> new methods, no conflicts

+misc/Gains.m

@@ -425,11 +425,16 @@ classdef Gains < handle & matlab.mixin.Copyable
 			
 			arguments
 				obj;
-				outputType;
+				outputType = "";
 				inputName = obj.inputType;
 				leftHandSide = "none";
 			end % arguments
 			
+			if strlength(outputType) == 0
+				texString = outputType;
+				return;
+			end
+			
 			% put string together
 			if isnumeric(leftHandSide)
 				if isscalar(leftHandSide) && obj.lengthOutput > 1

+misc/plotRope.m

+function f = plotRope(myOptions, rope)
+%PLOTSTRING Plots the strings defined by string
+arguments 
+	myOptions = [];
+end
+arguments (Repeating)
+	rope quantity.Discrete;
+end
+
+[myOptions, rope] = plotRopeParser(myOptions, rope);
+
+f = figure();
+for it = 1 : numel(rope)
+	ropeTemp = rope{it};
+	if myOptions.t.n <= 1
+		ropeDiscrete = ropeTemp.on(myOptions.s);
+		ropeDiscrete = reshape(ropeDiscrete, [1, size(ropeDiscrete)]);
+	else
+		ropeDiscrete = ropeTemp.on([myOptions.t, myOptions.s]);
+	end
+	
+	% plot
+	for tIdx = 1 : size(ropeDiscrete, 1)
+		if numel(ropeTemp) == 2
+			ax = plot(ropeDiscrete(tIdx, :, 1), ropeDiscrete(tIdx, :, 2));
+			
+		elseif numel(ropeTemp) == 3
+			ax = plot3(ropeDiscrete(tIdx, :, 1), ropeDiscrete(tIdx, :, 2), ropeDiscrete(tIdx, :, 3), ...
+				"LineWidth", myOptions.LineWidth, "Color", myOptions.Color);
+			
+		else
+			error("rope must be a quantity.Discrete with numel(rope) == 2 or == 3");
+		end
+	end % tIdx = 1 : size(ropeDiscrete, 1)
+	
+	if numel(ropeTemp) == 2
+		horLim = [-MAX(-ropeTemp(1)), MAX(ropeTemp(1))];
+		xlim(horLim);
+	else
+		horLim = [-MAX(-ropeTemp(1:2)), MAX(ropeTemp(1:2))];
+		xlim(horLim); ylim(horLim);
+	end
+end
+
+end % plotRope()
+
+function [myOptions, rope] = plotRopeParser(myOptions, rope)
+if isa(myOptions, "quantity.Discrete")
+	rope = {myOptions; rope{:}};
+	myOptions = [];
+end
+if isempty(myOptions)
+	if any(rope{1}(1).domain.gridIndex("t"))
+		myOptions = plotRopeOptionParser(rope{1}, "t", rope{1}(1).domain.find("t"));
+	else
+		myOptions = plotRopeOptionParser(rope{1});
+	end
+elseif isstruct(myOptions)
+	if ~isfield(myOptions, "t") && any(rope(1).domain.contains("t"))
+		myOptions.t = rope{1}(1).domain.find("t");
+	end
+	myOptions = misc.struct2namevaluepair(myOptions);
+	myOptions = plotRopeOptionParser(rope{1}, myOptions{:});
+else
+	error("myOptions is invalid data-type")
+end
+end % plotRopeParser()
+
+function myOptions = plotRopeOptionParser(rope, NameValueArgs)
+% misc.plotRopeOptionParser is a helper function for misc.plotRope -> see this function for
+% documentation.
+arguments
+	rope quantity.Discrete
+	NameValueArgs.shade = [0.4, 1];
+	NameValueArgs.LineWidth = 1;
+	NameValueArgs.Color = [0 0.4470 0.7410];
+	NameValueArgs.t = quantity.Domain("t", 0);
+	NameValueArgs.s = rope(1).domain.find("s");
+end % argumentss
+
+myOptions = NameValueArgs;
+
+end % plotRopeOptionParser()
\ No newline at end of file

+quantity/Discrete.m

@@ -1844,21 +1844,91 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 				for it = 1 : obj(1).domain.n
 					VDisc(it, :) = eig(objDisc(:, :, it));
 				end % for it = 1 : obj(1).domain.n
-			else	
+				V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
+			else
 				% [V,D] = eig(A) or [V,D,W] = eig(A)
 				VDisc = zeros([obj(1).domain.n, size(obj)]);
 				DDisc = zeros([obj(1).domain.n, size(obj)]);
 				WDisc = zeros([obj(1).domain.n, size(obj)]);
-				for it = 1 : obj(1).domain.n
-					[VDisc(it, :, :), DDisc(it, :, :), WDisc(it, :, :)] = ...
-						eig(objDisc(:, :, it));
+				VDiscOld = zeros(size(obj));
+				WDiscOld = zeros(size(obj));
+				for zt = 1 : obj(1).domain.n
+					[VDiscTemp, DDisc(zt, :, :), WDiscTemp] = eig(objDisc(:, :, zt));
+					if zt > 1
+						VDiscOld(:, :) = VDisc(zt-1, :, :);
+						WDiscOld(:, :) = WDisc(zt-1, :, :);
+						for jt = 1 : size(obj, 1)
+							% as the builtin eig() scales the eigenvectors, as slight change of the
+							% values can cause a change of the sign of the eigenvector. As eig is
+							% called pointwise, often a non-smooth eigenvector-function results.
+							% Due to this, every eigenvector is compared with the previous one and
+							% if there it makes the result more smooth, the sign is changed.
+							if norm(VDiscTemp(:, jt) - VDiscOld(:, jt)) ...
+									> norm(VDiscTemp(:, jt) + VDiscOld(:, jt))
+								VDiscTemp(:, jt) = -VDiscTemp(:, jt);
+							end
+							if norm(WDiscTemp(:, jt) - WDiscOld(:, jt)) ...
+									> norm(WDiscTemp(:, jt) + WDiscOld(:, jt))
+								WDiscTemp(:, jt) = -WDiscTemp(:, jt);
+							end
+						end % for jt = 1 : size(obj, 1)
+					end
+					VDisc(zt, :, :) = VDiscTemp;
+					WDisc(zt, :, :) = WDiscTemp;
 				end % for it = 1 : obj(1).domain.n
 				D = quantity.Discrete(DDisc, obj(1).domain, "name", "D");
 				W = quantity.Discrete(WDisc, obj(1).domain, "name", "W");
+				V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
 			end
-			V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
 		end % eig()
 		
+		function [l, u] = luDecomposition(obj)
+			% luDecomposition split matrix OBJ into lower triangular matrix l and upper triangular
+			% matrix u, such that l * u = obj.
+			n = size(obj, 1);
+			assert(numel(obj) == n*n, "obj must be a quadratic matrix");
+			u = obj;
+			l = obj.eye([n, n], obj(1).domain);
+
+			% simple algorithm from wikipedia:
+			for it = 1 : n-1
+				for k = it+1 : n
+					l(k, it) = u(k, it) / u(it, it);
+					for jt = it : n
+						u(k, jt) = u(k, jt) - l(k, it) * u(it, jt);
+					end % for jt
+				end % for k
+			end % for it
+		end % luDecomposition()
+		
+		function [V, J, l] = jordanReal(obj)
+			% jordanReal Real Jordan form of the matrix OBJ pointwise over the (scalar) domain.
+			% [V, J] = jordanReal(OBJ) computes the transformation matrix V and the jordan
+			% block matrix J, so that OBJ*V = V*J. The columns of OBJ are the generalised 
+			% eigenvectors.
+			%
+			% [V, J, l] = jordanReal(OBJ) also returns the lengths of the jordan
+			% blocks. l is a vector containing the lengths. So length(l) is the
+			% number of jordan blocks.
+			assert(numel(obj(1).domain)==1, "quantity.Discrete/jordanReal is only implemented for one domain");
+			assert((ndims(obj) <= 2) && (size(obj, 1) == size(obj, 2)), ...
+				"quantity.Discrete/jordanReal is only implemented for quadratic quantities");
+			
+			objDisc = permute(obj.on(), [2, 3, 1]);
+			
+			VDisc = zeros([obj(1).domain.n, size(obj)]);
+			JDisc = zeros([obj(1).domain.n, size(obj)]);
+			lDisc = zeros([obj(1).domain.n, size(obj, 1)]);
+			for it = 1 : obj(1).domain.n
+				[VDisc(it, :, :), JDisc(it, :, :), lTemp] = ...
+					misc.jordanReal(objDisc(:, :, it));
+				lDisc(it, 1:numel(lTemp)) = lTemp;
+			end % for it = 1 : obj(1).domain.n
+			V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
+			J = quantity.Discrete(JDisc, obj(1).domain, "name", "J");
+			l = quantity.Discrete(lDisc, obj(1).domain, "name", "l");
+		end % jordanReal()
+		
 		function y = log(obj)
 		%  log    Natural logarithm.
 		%     log(X) is the natural logarithm of the elements of X.
@@ -2502,7 +2572,31 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 				O = ones([domain.gridLength, valueSize(:)']);
 				P = quantity.Discrete(O, domain, 'size', valueSize, varargin{:});
 			end
-		end
+		end % ones()
+		
+		function I = eye(N, domain, varargin)
+			%  eye Identity matrix.
+			%     eye(N, domain) is the N-by-N identity matrix on the domain domain.
+			%  
+			%     eye([M,N], domain) is an M-by-N matrix with 1's on the diagonal and zeros 
+			%		elsewhere.
+			%     eye([M,N], domain, varargin) passes further inputs to the quantity.Discrete
+			%     constructor.
+			
+			if isscalar(N)
+				N = [N, N];
+			elseif ~ismatrix(N)
+				error("eye only supports ndims(N) == 1 or 2");
+			end
+			
+			if any( N == 0)
+				P = quantity.Discrete.empty(N);
+			else
+				IvalueDiscrete = reshape(eye(N), [ones(1, numel(domain)), N]);
+				IvalueDiscrete = repmat(IvalueDiscrete, [domain.n, 1]);
+				I = quantity.Discrete(IvalueDiscrete, domain, varargin{:});
+			end
+		end % eye()
 		
 		function P = zeros(valueSize, domain, varargin)
 			%ZEROS initializes an zero quantity.Discrete object

+quantity/Symbolic.m

@@ -828,7 +828,25 @@ classdef Symbolic < quantity.Function
 			%ZEROS initializes an zero quantity.Discrete object of the size
 			%specified by the input valueSize.
 			P = quantity.Symbolic(zeros(valueSize), varargin{:});
-		end
+		end % zeros()
+		
+		function P = ones(valueSize, varargin)
+			%ONES initializes an ones-quantity.Discrete object
+			%	P = ones(VALUESIZE, DOMAIN) creates a matrix of size
+			%	VALUESIZE on the DOMAIN with ones as entries.
+			P = quantity.Symbolic(ones(valueSize), varargin{:});
+		end % ones()
+		
+		function P = eye(N, varargin)
+			%  eye Identity matrix.
+			%     eye(N, domain) is the N-by-N identity matrix on the domain domain.
+			%  
+			%     eye([M,N], domain) is an M-by-N matrix with 1's on the diagonal and zeros 
+			%		elsewhere.
+			%     eye([M,N], domain, varargin) passes further inputs to the quantity.Discrete
+			%     constructor.
+			P = quantity.Symbolic(eye(N), varargin{:});
+		end % eye()
 
 	end % methods (Static)
 	

+unittests/+quantity/testDiscrete.m

@@ -4,8 +4,36 @@ function [tests] = testDiscrete()
 tests = functiontests(localfunctions);
 end
 
-function testCatBug(tc)
 
+
+function testEye(tc)
+myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
+tc.verifyEqual(MAX(abs(quantity.Discrete.eye(1, myDomain) - 1)), 0);
+tc.verifyEqual(MAX(abs(quantity.Discrete.eye(2, myDomain) - eye(2))), 0);
+tc.verifyEqual(MAX(abs(quantity.Discrete.eye([2, 1], myDomain) - eye([2, 1]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Discrete.eye([2, 2], myDomain(2)) - eye([2, 2]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Discrete.eye([4, 4], myDomain) - eye(4))), 0);
+end % testEye()
+
+function testLuDecomposition(tc)
+myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
+syms z zeta
+quanSymbolic = quantity.Symbolic([z*zeta, z+1, 2; z^2, zeta^2, -2; 1, 2, 3], myDomain);
+quan = quantity.Discrete(quanSymbolic);
+% discrete case
+[l, u] = quan.luDecomposition();
+tc.verifyEqual(MAX(abs(l*u - quan)), 0, 'AbsTol', 1e3*eps);
+tc.verifyEqual(MAX(abs(l.*triu(ones(size(l)), 1))), 0, 'AbsTol', 1e3*eps);
+tc.verifyEqual(MAX(abs(u.*tril(ones(size(u)), -1))), 0, 'AbsTol', 1e3*eps);
+% symbolic case
+[l, u] = quanSymbolic.luDecomposition();
+tc.verifyEqual(MAX(abs(l*u - quanSymbolic)), 0);
+tc.verifyEqual(MAX(abs(l.*triu(ones(size(l)), 1))), 0);
+tc.verifyEqual(MAX(abs(u.*tril(ones(size(u)), -1))), 0);
+tc.verifyTrue(isa(quanSymbolic, "quantity.Symbolic"));
+end % testLuDecomposition()
+
+function testCatBug(tc)
 z = quantity.Domain("z", linspace(0, 1, 7));
 zeta = quantity.Domain("zeta", linspace(0, 1, 5));
 A = quantity.Discrete(sin(z.grid * pi), z);
@@ -19,10 +47,7 @@ x = A * subs(A, "z", "zeta");
 C1 = [0*z.Symbolic + zeta.Symbolic; x];
 C2 = [zeta.Symbolic + 0*z.Symbolic; x];
 tc.verifyEqual(C1.on([z, zeta]), C2.on([z, zeta]));
-
 end
-
-
 function testEigScalar(tc)
 z = quantity.Domain("z", linspace(-1, 1, 11));
 [Vz, Dz, Wz] = z.Discrete.eig();
@@ -50,6 +75,19 @@ tc.verifyEqual(eig(quan.at(0.4)), Eq.at(0.4), "AbsTol", 1e-15);
 tc.verifyEqual(on(Eq), on(Dq.diag2vec()), "AbsTol", 1e-15);
 end % testEigMat()
 
+function testJordanRealMat(tc)
+z = quantity.Domain("z", linspace(-1, 1, 11));
+zSym = sym("z");
+quan = quantity.Discrete(quantity.Symbolic([1+zSym^2, 2; -zSym, zSym], z));
+% all output arguments
+[V, J, l] = quan.jordanReal();
+[V0p4, J0p4, l0p4] = misc.jordanReal(quan.at(0.4));
+tc.verifyEqual(J0p4, J.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual(V0p4, V.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual([l0p4; 0], l.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual(on(quan * V), on(V * J), "AbsTol", 1e-15);
+end % testJordanRealMat()
+
 function testPower(tc)
 %% scalar case
 z = quantity.Domain("z", linspace(0, 1, 11));

+unittests/+quantity/testSymbolic.m

@@ -4,6 +4,25 @@ function [tests ] = testSymbolic()
 tests = functiontests(localfunctions());
 end
 
+
+function testOnes(tc)
+myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
+tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(1, myDomain) - 1)), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.ones(2, myDomain) - ones(2))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 1], myDomain) - ones([2, 1]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([2, 2], myDomain(2)) - ones([2, 2]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.ones([4, 4], myDomain) - ones(4))), 0);
+end % testOnes()
+
+function testEye(tc)
+myDomain = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 5))];
+tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(1, myDomain) - 1)), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.eye(2, myDomain) - eye(2))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 1], myDomain) - eye([2, 1]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([2, 2], myDomain(2)) - eye([2, 2]))), 0);
+tc.verifyEqual(MAX(abs(quantity.Symbolic.eye([4, 4], myDomain) - eye(4))), 0);
+end % testEye()
+
 function testValueContinuousBug(tc)
 z = quantity.Domain("z", linspace(0, 1, 7));
 zeta = quantity.Domain("zeta", linspace(0, 1, 5));