Matlab implementation of the "deleted pivot" operation described by Greg Kuperberg in "Kasteleyn cokernels"
1; % script file | |
function res = deleted_pivot(M, i, j) | |
% Calculates the deleted pivot of matrix M at row i, column j | |
% as described by Greg Kuperberg in "Kasteleyn cokernels" | |
% X = the jth column of M with element i removed | |
X = [M(:,j)(1:i-1); M(:,j)(i+1:end)]; | |
% Yt = the ith row of M with element j removed | |
Yt = [M(i,:)(1:j-1), M(i,:)(j+1:end)]; | |
% Mxy is what we're calling M', which is the matrix M with row i | |
% and column j removed. | |
Mx = [M(:, 1:j-1), M(:, j+1:end)]; | |
Mxy = [Mx(1:i-1, :); Mx(i+1:end, :)]; | |
% The result is given as M' - X*Yt, where * is the outer product | |
res = Mxy - X*Yt; | |
endfunction | |
function res = deleted_pivots(M, varargin) | |
% Calculates the deleted pivot of matrix M at rows i_k, and columns | |
% j_k given by arguments [i_k, j_k]. | |
% Example: | |
% M = eye(5) | |
% Mdp = deleted_pivots(M, [2, 4], [1, 3], [5, 5]) | |
% Stop recursion if we're out of varargin | |
if (numel(varargin) == 0) | |
res = M; | |
return; | |
endif | |
% Validate and then do the actual computation | |
ij = varargin{1}; | |
assert(numel(ij) == 2, "Arguments to deleted_pivots must be two-element arrays.") | |
M_new = deleted_pivot(M, ij(1), ij(2)); | |
% Compute varargin for the next recursion | |
next_varargin = varargin(2:end); | |
reduce_indices_curried = @(i) @(j) @(indices) reduce_indices(i, j, indices); | |
reduce_indices_ij = reduce_indices_curried(ij(1))(ij(2)); | |
next_varargin = cellfun(reduce_indices_ij, next_varargin, "UniformOutput", false); | |
% Recurse! | |
res = deleted_pivots(M_new, next_varargin{:}); | |
endfunction | |
function res = reduce_indices(i, j, indices) | |
% Reduce indices [i_k, j_k] by one if they are greater than i or j, | |
% respectively. This is useful for updating indices when a row or | |
% column of a matrix has been removed. | |
res = indices; | |
if (indices(1) > i) | |
res(1) -= 1; | |
endif | |
if (indices(2) > j) | |
res(2) -= 1; | |
endif | |
endfunction |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment