%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-12-07 09:13:42
 % @LastEditTime: 2022-01-13 21:12:16

 % sensitivity on align-CVT of MMC
 % dcdp = dcdH * dHdp
%--------------------------
function [c_dH, c_interEles] = sensitivity_mmc_align(nodes_cell, faces_cell, seeds0, step, T, ...
    bpoly, epsilon, conn_list, c_rho, c_phi, c_interPolySeed, maxx, minx, maxy, miny)

    npnts = size(seeds0, 1);    
    c_dH        = cell(npnts, 3);
    c_interEles = cell(npnts, 3);

    seeds = [seeds0, T];
    
    parfor i = 1 : npnts
        
        neighbor = conn_list{i};
        local_idx = int32([neighbor; i]);

        for j = 1:3
%             tic
            % minus
            seeds_minus = seeds;
            seeds_minus(i, j) = seeds_minus(i, j) - step(j);
            
            if j == 1 && seeds_minus(i, j) < minx
                seeds_minus(i, j) = minx;
            elseif j == 2 && seeds_minus(i, j) < miny
                seeds_minus(i, j) = miny;
            end

            % only update the neighbor cells
            [~, edges, c_nnbEdgeIDX, nnbSeedIDX] = intersection_info(seeds_minus(:,1:2), ...
                local_idx, bpoly);
            [MMCs, ~] = update_MMCs(edges, seeds_minus(:,3), nnbSeedIDX, seeds_minus(:,1:2));

            [rhoMinus, interEles] = update_projection(nodes_cell, faces_cell, ...
                MMCs, epsilon, c_rho, c_phi, c_nnbEdgeIDX, c_interPolySeed, local_idx);
            
            % plus
            seeds_plus = seeds;
            seeds_plus(i, j) = seeds_plus(i, j) + step(j);
            
            if j == 1 && seeds_plus(i, j) > maxx
                seeds_plus(i, j) = maxx;
            elseif j == 2 && seeds_plus(i, j) > maxy
                seeds_plus(i, j) = maxy;
            end

            [~, edges, c_nnbEdgeIDX, nnbSeedIDX] = intersection_info(seeds_plus(:,1:2), ...
                local_idx, bpoly);
            [MMCs, ~] = update_MMCs(edges, seeds_plus(:,3), nnbSeedIDX, seeds_plus(:,1:2));

            [rhoPlus, interEles] = update_projection(nodes_cell, faces_cell, ...
                MMCs, epsilon, c_rho, c_phi, c_nnbEdgeIDX, c_interPolySeed, local_idx);

            nnb = length(interEles);
            dH_ele = cell(nnb, 1); 
            for k = 1 : nnb
                ele = interEles(k);
                dH_ele{k} = (rhoPlus{ele} - rhoMinus{ele}) / 2 / step(j);
%                 max(dH_ele{k}(:))
            end
            c_dH{i,j} = dH_ele;
            c_interEles{i,j} = interEles;

%             mp = [];
%             figure; 
%             for testi = 1:size(MMCs_poly,1)
%                 mp = [mp; polyshape(reshape(MMCs_poly(testi,:), 2, [])')];
%             end
%             plot(mp);
%             bp = [];
%             for testi = 1:nnb
%                 bp = [bp; polyshape(c_polys{interEles(testi)})];
%             end
%             hold on; plot(bp);
%             toc
        end
    end
end