%% Efficient Object-Oriented Kronecker Product Manipulation
% Sean's pick this week is KronProd by Matt J.
%
%% Background
% This one is old and very dear to my heart.  Some time in the 2009 era, as I was starting grad school, I needed help calculating subpixel resolution displacement vectors between to CT image volumes.  The displacement vector calculation was one of the main underlying things pieces of most of my grad research.  This was before imregcorr existed, though it still does not support 3D which was my requirement.
% I'd managed to get a 2D variant working, but needed 3D.  Matt replied to my newsreader post to discuss the algorithm and its extension in 3D, and pointed to this to be more efficient than using kron directly.  Performance was important as I'd have to compute over four billion displacement vectors over the next year.
% Matt has provided many algorithms and overloaded a bunch of operators to work with the KronProd class and configured it with custom indexing and size.  In 2009, I didn't understand how custom indexing with subsref/subsasgn/size overloading worked, so it seemed like magic.  With the knowledge of how those work now, and the performance implications in that era, I probably would've pulled the core algorithm into a function for speed since my use was stateless, i.e. a new KronProd each time.
% Here are a few simple operations with a KronProd.  We'll create one from a three arrays.

kp = KronProd({magic(3), [randi(10,[3 4]), zeros(3,1)], randn(9, 15)})

%
% Now we can work with this as if the full expansion you'd get from the kron function.

nzeros = nnz(kp)

% Or a matrix operation where I get back KronProds.

[u, s, v] = svd(kp)

% Matt has overloaded subsref so when you index into the KronProd, you get back a KronProd or just matrix.

kkp = kp(1)
km = kp{1}

% He did not overload subsasgn though, so we can't do the opposite.

try
    kp{1} = rand(3)
catch ME
    fprintf(2, ME.message)
end

% He's also overloaded size so the appearance of size is that of the full Kronecker product.  Under the hood, however, the KronProd is actually a scalar, it just appears to be the expanded size.

sz = size(s)

%
% Starting in MATLAB R2021b, we now have an easier and more formal means for controlling indexing and creating scalar objects using modular indexing.
% Right now subsref overloads (), {}, and dot(.).  It overloads dot, just because it has to with the old approach, so we don't need to do that with the newer indexing.  So the class definition line would need to be:
% classdef KronProd < matlab.mixin.indexing.RedefinesParen & matlab.mixin.indexing.RedefinesBrace
% Now, we'd need to implement the abstract methods for parenAssign, parenReference, braceAssign, braceReference.  These methods would then use similar logic to what's used today, in a protected methods block, something along these lines:

function B = braceReference(obj, indexOp)
    inds = M.opinds(indexOp(1).Indices{1});
    B=M.opset(inds);
    if isnum(inds)
        B=B{:};
    end
end

%% Comments
% Matt, if we ever meet in real life, I owe you a beer!
% Give it a try and let us know what you think here or leave a comment for Matt.