How do I implement a numeric data type?

Hello MatLab experts,

I'm implementing a numeric data type. It must (eventually) implement all of the mathematical operators and functions +, *, -, /, log, exp, etc...

The problem I'm running into is that objects in matlab can implicitly be stored as vectors or matrices, and I'm having trouble figuring out how best to deal with that. This project is currently in-progress and not very well organized, but the full file is currently available here https://github.com/jeffersoncarpenter/RandomVariable/blob/master/RandomVariable.m

The numeric data type is for "random variables" or specifically for observed values with known uncertainties that I'm collecting as physics lab data.

There is a formula for error propagation involving the squares of partial derivatives, and it amounts to applying the chain rule at each step of the computation in which the numbers are pairs (observational value, standard deviation) and add, multiply, and so on as random variables. Ergo the "random variable" numeric data type.

Here is a sample class method. It's pretty long, but the thing to notice is the top-level if statement checking whether the argument is a Datum (random variable) or a standard number, followed by more if statements to determine whether to loop over the left-hand argument, the right-hand argument, or both. Ultimately I will have to support multiplying on the left by standard numbers, multiplying on the right by standard numbers, and multiplying two random variables together, together with every combination of vector or matrix rank.

Is there a better way to do this? I expect to have to branch to check whether I'm operating on the left-hand side or the right-hand side with a standard number instead of a Datum (random variable), but I would like to not have to further check - in the body of each operator, that would be a lot of code - whether the left-hand side, right-hand side, or both are vectors or matrices. It would be great if that could be written more compactly.

        function r = mtimes(a, b)
            if isa(a, 'Datum')
                if 1 == numel(b)
                    r = a;
                    for i = 1:numel(a)
                        r(i) = Datum(a(i).Value * b.Value, a(i).Value^2 * b.Variance + a(i).Variance * b.Value^2);
                    end
                elseif (1 == numel(a))
                    r = b;
                    for i = 1:numel(b)
                        r(i) = Datum(a.Value * b(i).Value, a.Value^2 * b(i).Variance + a.Variance * b(i).Value^2);
                    end
                else
                    r = a;
                    for i = 1:numel(a)
                        r(i) = Datum(a(i).Value * b(i).Value, a(i).Value^2 * b(i).Variance + a(i).Variance * b(i).Value^2);
                    end
                end
            else
                if 1 == numel(b)
                    if 1 == numel(a)
                        r = Datum(a * b.Value, a^2 * b.Variance);
                    else
                        for i = 1:numel(a)
                            r(i) = Datum(a(i) * b.Value, a(i)^2 * b.Variance);
                        end
                    end
                else
                    r = b;
                    for i = 1:numel(b)
                        r(i) = Datum(a * b(i).Value, a^2 * b(i).Variance);
                    end
                end
            end
        end

Answer 1

[Ameer Hamza]

See the following definition of your class. This uses vectorized operations, which make the code simple, faster, and more readable. Also, there is no need to use if-condition for each case, as did in your code.

% WIP
classdef Datum
    properties
        Value
        Variance
    end
    methods
        function obj = Datum(Val_ini, Var_ini)
            % Following allow faster initialization. See here: https://blogs.mathworks.com/loren/2012/03/26/considering-performance-in-object-oriented-matlab-code/
            if nargin > 0
                for i = numel(Val_ini):-1:1
                    obj(i).Value = Val_ini(i);
                    obj(i).Variance = Var_ini(i);
                end
            end
        end
        
        function r = plus(a, b)
            if isa(a, 'Datum') && isa(b, 'Datum')
                r = Datum(a.value + b.value, a.variance + b.variance);
            elseif isa(a, 'Datum')
                r = Datum(a.value + b, a.variance);
            else % only when b is Datum
                r = Datum(a + b.value, b.variance);
            end
        end
        
        function r = minus(a, b)
            r = plus(a, -1 * b);
        end
        
        function r = mtimes(a, b)
            % matlab's automatically making everything matrices all the
            % time is a disaster
            if isa(a, 'Datum') && isa(b, 'Datum')
                r = Datum(a.value .* b.value, a.value.^2 .* b.variance + a.variance .* b.value.^2);
            elseif isa(a, 'Datum')
                r = Datum(a.value * b, b^2 * a.variance);
            else
                r = Datum(a * b.value, a^2 * b.variance);
            end
        end
        function r = mpower(a, b)
            if isa(a, 'Datum')
                r =  Datum(a.value.^b, b^2 * a.value.^(2*b-2) .* a.variance);
            elseif isa(b, 'Datum')
                r = Datum(a.^b.value, log(a)^2 * a.^(2 * v.value) .* b.variance);
            end
        end
        function r = mrdivide(a, b)
            r = mtimes(a, b^-1);
        end
        function r = rdivide(a, b)
            r = arrayfun(@mrdivide, a, b);
        end
        function r = log(a)
            r = Datum(log(a.value), a.variance ./ a.value);
        end
        function r = exp(a)
            r = Datum(exp(a.value), exp(2 * a.value .* a.variance));
        end
        function r = variance(a)
            r = [a.Variance];
        end
        function r = stdev(a)
            r = sqrt([a.Variance]);
        end
        function r = value(a)
            r = [a.Value];
        end
    end
end

Few test cases

x1 = Datum(1:10, zeros(1, 10)); % 1 x 10 Datum
x2 = Datum(2*ones(1, 10), 2*ones(1,10)); % 1 x 10 Datum
x3 = Datum(1, 3); % 1 x 1 Datum

y1 = 2 + x1;
y2 = x2 - 3;
y3 = x1 + x2;
y4 = x2 - x1;
y5 = x1 - x3;
y6 = x1 * x2;
y7 = x1 * x3;
y8 = 2 * x1;
y9 = x1 * 5;

[Jefferson Carpenter]

Ah, that's good. I hadn't thought of making sure to put the arrays inside the datum, rather than having an array of datums.

[Ameer Hamza]

No, this does not put an array inside a Datum class. It creates an array of Datum objects, and each element still holds a single value. If you check the values of x1 and x2, you will see that these are Datum arrays.

[Jefferson Carpenter]

It seems to be working.