function [A, Bx, By] = operator2matrix(Mx, x, y, F, Bx, By, My)
% Forms a matrix representing a linear operator between two tangent spaces
%
% function [A, Bx, By] = operator2matrix(M,  x, y, F)
% function [A, Bx, By] = operator2matrix(Mx, x, y, F, [], [], My)
% function [A, Bx, By] = operator2matrix(M,  x, y, F, Bx, By)
% function [A, Bx, By] = operator2matrix(Mx, x, y, F, Bx, By, My)
%
% Given a manifold structure M, two points x and y on that manifold, and a
% function F encoding a linear operator from the tangent space T_x M to the
% tangent space T_y M, this tool generates two random orthonormal bases
% (one for T_x M, and one for T_y M), and forms the matrix A which
% represents the operator F in those bases. In particular, the singular
% values of A are equal to the singular values of F. If two manifold
% structures are passed, then x is a point on Mx and y is a point on My.
%
% If Bx and By are supplied (as cells containing orthonormal vectors in
% T_x M and T_y M respectively), then these bases are used, and the matrix
% A represents the linear operator F restricted to the span of Bx, composed
% with orthogonal projection to the span of By. Of course, if Bx and By are
% orthonormal bases of T_x M and T_y M, then this is simply a
% representation of F. Same comment if two manifolds are passed.
%
% See also: tangentorthobasis hessianmatrix

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sep. 13, 2019.
% Contributors: 
% Change log: 

    % By default, the two points are on the same manifold
    if ~exist('My', 'var') || isempty(My)
        My = Mx;
    end

    if ~exist('Bx', 'var') || isempty(Bx)
        Bx = tangentorthobasis(Mx, x);
    end
    if ~exist('By', 'var') || isempty(By)
        By = tangentorthobasis(My, y);
    end

    assert(iscell(Bx) && iscell(By), 'Bx and By should be cells.');

    n_in = numel(Bx);
    n_out = numel(By);
    
    A = zeros(n_out, n_in);
    
    for k = 1 : n_in
        A(:, k) = tangent2vec(My, y, By, F(Bx{k}));
    end

end
