function [pert, C21, all_T12, all_C21] = HOPE(M, N, options)
%% get number of nodes
n1 = M.n;
n2 = N.n;
num = n1;
if n1 > n2
    num = n2;
end

if nargout > 2, all_T12 = {}; all_C21 = {}; end

adj1 = sparse(M.adj);
adj2 = sparse(N.adj);

% G1 = build_witnesses(adj1, n1, 2);
% G2 = build_witnesses(adj2, n2, 2);

%% initialize map with shot
T12 = knnsearch(N.shots, M.shots,'NSMethod','kdtree');
pert = sparse(T12, 1:length(T12), 1, n2, n1);
if nargout > 2, all_T12{end+1} = pert; all_C21{end+1} = 0;end

%% refining params
% cc = ceil(n1*0.002);
% nums = [1:n1]';
% size_L2 = 0;
% count = 0;

[T12, col] = find(pert);
% num = length(T12);
% nn = length(col);
MA = full(diag(M.Phi.A));
e = zeros(num,1);
R_max = max(max(N.distances));
th = linspace(100,0.2, 10);%[100, 50,  25, 12.5, 6.25, 3.125, 1.5625 0.2];
max_T_o = 1;
max_T = 1;
%% smoothness, continuity and global accuracy
for kk = 1:60   
    D1 = build_witnesses(adj1, n1, max_T);
    D2 = build_witnesses(adj2, n2, max_T);
    max_T = max_T + 1;
    if max_T > 8
        max_T = max_T_o;
    end
    % begining = 1;
    % vec = double(int32(linspace(nn/cc,nn,cc)));
    % for k=1:length(vec)
    %     ending = vec(k);
    %     idx = col(begining:ending);
    %     begining = ending + 1;                
    %     pertF = greedy_match(G2*(pert*G1(:, idx))); 
    %     [idx2, ~] = find(pertF);                
    %     idx = idx(1:size(idx2, 1));
    %     T12(idx) = idx2;
    %     pert = sparse(T12, col(1:length(T12)), 1, n2, n1); 
    % end 
    
    % pert = greedy_match(D2*(pert*D1(:, 1:n1))); 
    % [idx2, idx] = find(pert); 
    % T12(col(idx)) = idx2;
    % pert = sparse(T12, col(1:length(T12)), 1, n2, n1); 

    DD1 = cell(n1,1);
    DD2 = cell(n2,1);
    idx = cell(num,1);
    ee = zeros(num,1);
    good = 1:num;
    D1T = M.distances(good,good); 
    D2T = N.distances(good,T12(good));
    for i = 1:length(good)
        idx{i} = find(D1T(:,i) ~= 0);
        DD1{i} = D1T(idx{i},i);
        DD2{i} = D2T(idx{i},i);
        DD2{i}(DD2{i}==0) = R_max;
        r = max(DD1{i});
        ee(i) = sum(((abs(DD1{i}-DD2{i}))/r).*MA(idx{i}))/sum(MA(idx{i}));        
    end
    e(good)=ee; 
    if kk<=length(th)
        landmarks = find(e<th(kk));
    else
        landmarks = find(e<th(end));
    end

    size_L = length(landmarks);
    bad = setdiff(1:n1, landmarks);
    size_b = length(bad);

    % patience = size_L - size_L2;
    % if ((patience < 1) && (kk > 2))
    %     count = count + 1;
    % end
    % if count > 10
    %     break
    % end 
    % size_L2 = size_L;

    % fprintf("[GLOBAL] iter: %i, num_lmks: %i, num_bad: %i\n", kk, size_L, size_b);
    if size_b > 0  
        pertF = greedy_match(D2*(pert*D1(:, bad)));
        [idx2, idx] = find(pertF);
        bad = bad(idx);
        T12(col(bad)) = idx2;           
        pert = sparse(T12, col(1:length(T12)), 1, n2, n1); 
    end  

    C21 = 0;
    if nargout > 2, all_T12{end+1} = pert; all_C21{end+1} = C21;end
end
