A Matlab class for Commutative Hypercomplex number: How to overload the basic Matlab function?

Question

I am programming commutative hypercomplex number class.

First of all: If that already exist, please tell. I have not found. PS: Quaternion exists but are not commutative.

I am trying to reproduce all that Matlab can do with real and complex, but with hypercomplex. The solution that I am working on is to define a class. If a more smart method may be used, please tell.

Many codes with complex numbers exist, the idea is to overload smartly to make working all existing program with Hypercomplex. This work might be used by many users.

please see at the end of the post the main maths principle found in Ablamowicz1996 _ Clifford Algebras with Numeric and Symbolic Computations

My Hypercomplex class works for basic operation: I have overloaded the mtimes(a,b) = a*b operation. For instance, I can do:

A = MyClass(1);
B = MyClass(2);
A*B; % ->that works

As a main problem, for instance in 2D, I would like to be able to use

[A B C].*[D ; E]

while using the element-wise behavior of Matlab, and for any dimension , and with my definition of *. For instance, while dealing with matrices of compatible size but not identical:

A = [1 2 3]; B = [1;2]

Matlab gives

A.*B    
ans =

 1     2     3
 2     4     6

That works also for N dimension when compatible. How can I overload times(a,b) = a.*b , power, rdivide.. in order to have the element-wise behavior of Matlab?

Here is my Class

classdef Hypercomplex < handle
    properties (SetAccess = public, GetAccess = public)
        A
    end
    properties (SetAccess = private, GetAccess = private)
        eps,nu,e1,e2
    end
    
    methods
        function s = Hypercomplex(Z)             %% Conctruc
            % Z might be real, complex or hypercomplex
            if isa(Z,'Hypercomplex')            % is HyperComplex
                s = Z;
            elseif ~isreal(Z) && length(Z) == 1 % is Complex
                s = Hypercomplex([real(Z) imag(Z) 0 0]);
            elseif isreal(Z) && length(Z) < 5
                s.A = zeros(1,4);
                for ii = 1:length(Z)            % is vecteur
                    s.A(ii) = Z(ii);
                end
            end
        end
        function s = epsnu(s)
            s.eps = (s.A(1) - s.A(4)) + 1i*(s.A(2)+s.A(3));
            s.nu  = (s.A(1) + s.A(4)) + 1i*(s.A(2)-s.A(3));
        end
        function s = e1e2(s)
            s.e1 = Hypercomplex([0.5 0 0 -0.5]);
            s.e2 = Hypercomplex([0.5 0 0  0.5]);
        end
        %% Overload
        function out = cos(a)
            a.epsnu.e1e2;
            out = cos(a.eps)*a.e1 + cos(a.nu)*a.e2;
            % I might also do for  sin exp sqrt ...
        end
        
        function out = mtimes(a,b) % '*'
            a = Hypercomplex(a).epsnu;
            b = Hypercomplex(b).epsnu;
            eps12 = 0.5*(a.eps * b.eps);% produit complex
            nu12  = 0.5*(a.nu  * b.nu );% produit complex
            
            h1 = Hypercomplex([real(eps12) imag(eps12) imag(eps12) -real(eps12)]);
            h2 = Hypercomplex([real(nu12)  imag(nu12) -imag(nu12)   real(nu12) ]);
            out = h1+h2;
        end
        function out = mrdivide(a,b)
            a = Hypercomplex(a);
            b = Hypercomplex(b);
            out = a * b.inverse;
        end
        function s = inverse(s)
            s = Hypercomplex(s).epsnu.e1e2;
            s = 1/(s.eps)*s.e1 + 1/(s.nu)*s.e2;
        end
        function out = plus(a,b)
            a = Hypercomplex(a);
            b = Hypercomplex(b);
            out = Hypercomplex(a.A + b.A);
        end
        function out = minus(a,b)
            a = Hypercomplex(a);
            b = Hypercomplex(b);
            out = Hypercomplex(a.A - b.A);
        end
    end
end

Example i^2 = -1 :

>> Ni =  Hypercomplex([0 1 0 0]) ;
>> Ni*Ni
  Hypercomplex with properties:

    A: [-1 0 0 0]

user_7328782 · Answer

I would think you need to loop over your arrays, and have a special case where, if for one dimension one of the inputs has size 1 and the other doesn't, then the same element for that array is used. A simple way to accomplish this is to index with (i-1)*step + 1, where step is 0 if the array has size 1, and 1 otherwise.

One thing to note is that, in your overloaded method, indexing doesn't use your custom indexing definition (subsref and subsasgn), but the built-in one.

For example, assuming you only allow 2D arrays of your class:

function C = times(A, B)
szA = size(A);
szB = size(B);
if numel(szA)~=2 || numel(szB)~=2
   error('Only 2D arrays are supported')
end
szC = max(szA, szB);
stepA = [matchSizes(szA(1), szC(1)), matchSizes(szA(2), szC(2))];
stepB = [matchSizes(szB(1), szC(1)), matchSizes(szB(2), szC(2))];
C = MyClass.empty(szC);
for jj = 1:szC(2)
   for ii = 1:szC(1)
      C(ii,jj) = A((ii-1)*stepA(1)+1,(jj-1)*stepA(2)+1) * B((ii-1)*stepB(1)+1,(jj-1)*stepB(2)+1);
   end
end

% private helper function, keep in the same file
function step = matchSizes(szIn, szOut)
step = 1;
if szIn ~= szOut
   if szIn == 1
      step = 0;
   else
      error('Sizes not compatible')
   end
end

NOTE: Code not tested.

NOTE2: This is one of the cases where 0-based indexing is more convenient!

user_6579744 · Answer

If MyClass has a property val to store numeric values I think (hope) that this should be possible (not tested):

function obj = times(A, B)
  r = reshape([A.val], size(A)) .* reshape([B.val], size(B));
  obj = arrayfun(@MyClass, r);
end

EDIT (because question is edited):

You can follow two ways to implement array of those numbers:

Array of objects
In this method each object just contains a hypercomplex number. It is the way that you have tried. As @CrisLuengo pointed out you can use loops to create the desired result or you can use a combination of vectorization and use of functions like arrayfun and cellfun. You need to know how to work with cell arrays and comma separated list and be familiar with verctorization operations like implicit expansion. Here I implemented plus and epsnu using vectorization and cellfun. They can work with both scalar and object arrays:
```
 function out = plus(a,b)
     a = Hypercomplex(a);
     b = Hypercomplex(b);
     AA = reshape([a.A], [4, size(a)]);
     BB = reshape([b.A], [4, size(b)]])
     CC = AA + BB;
     A1 = reshape()
     out = cellfun(@Hypercomplex, num2cell(reshape(CC, 4, []).', 2));
 end 

 function s = epsnu(s)
   AA = reshape([s.A], 4, []);
   eps = num2cell(AA(1, :) - AA(4, :) + 1i.*(AA(2, :) + AA(3, :)));
   nu = num2cell(AA(1, :) + AA(4, :) + 1i.*(AA(2, :) - AA(3, :)));
   [s.eps] = eps{:};
   [s.nu] = nu{:};
end
```
Object containing array of numbers
Instead of defining A as a [4 x 1] array you can define each coordinate (x,y,z,et) as a separate property of class each can contain array of numbers. you need to overload the cat operator to concatenate objects. As a result no object array is created but the underlying verctors are concatenated. You need tooverload subsref and subsasgn if you want to simulate element access. I think this method is simpler to implement than the first method.

A Matlab class for Commutative Hypercomplex number: How to overload the basic Matlab function?

2 Answers