function [A, G, repu_obj, quad_const] = compute_parts(Y,P,weights,Pnz,I,J,attr,theta,method,constant,varargin)
% wrapper of computing parts for specific optimizers (e.g. 'mm' or 'fphssne')
% 
%   [A, G, repu_obj, quad_const] = compute_parts(Y,P,weights,Pnz,I,J,attr,theta,method,constant,varargin)
% 
% Input:
%   Y : N x 2, 2-D coordinates
%   P : N x N, the pairwise similarities or network (weighted) adjacency
%   weights: weights of the data points (or graph nodes)
%   Pnz : = nonzeros(P)
%   I,J : the indices of non-zero entries in P
%   attr : the over-attraction scaling factor (default=1)
%   theta : the Barnes-Hut tree farness-factor parameter
%   method : the visualization method
%   constant : the constant part of the objective function
%   varargin : visualization method specific parameters
%              for LinLog, lambda = varargin{1}
%              for EE, lambda = varargin{1}
%              for NeRV, lambda = varargin{1}, epsilon = varargin{2}
%
% Output:
%   A : the coefficients of the quadratic part, i.e. W+W' in the paper
%   G : the gradient of the Lipschitz part (usually the repulsive part)
%   repu_obj: the repulsive objective
%   quad_const: the constant induced in the quadratification step
%
% Copyright (c) 2014, Zhirong Yang (Aalto University)
% All rights reserved.

n = size(P,1);
switch method
    case 'linlog'
        lambda = varargin{1};
%         myeps = 1e-50;
        myeps = 0;
        qnz = sqrt(sum((Y(I,:)-Y(J,:)).^2,2))+myeps;
        A = attr * sparse(I,J,0.5*Pnz./qnz,n,n);
        [repu_obj, G] = compute_linlog_obj_grad_repulsive_barneshut(Y, theta, lambda, 2);
        quad_const = attr * 0.5*sum(Pnz.*qnz);
    case 'ee'
        lambda = varargin{1};
        A = attr * P;
        [repu_obj, G] = compute_ee_obj_grad_repulsive_barneshut(Y, theta, lambda, 2);
        quad_const = 0;
    case 'sne'
        A = attr * 0.5 * (P+P');
        [repu_obj, G] = compute_sne_obj_grad_repulsive_barneshut(Y, theta, 2);
        quad_const = 0;
    case 'ssne'
        A = attr * P;
        [repu_obj, G] = compute_ssne_obj_grad_repulsive_barneshut(Y, theta, 2);
        quad_const = 0;
    case 'tsne'
        qnz = 1./(1+sum((Y(I,:)-Y(J,:)).^2,2));
        A = attr * sparse(I,J,Pnz.*qnz,n,n);
        [repu_obj, G] = compute_tsne_obj_grad_repulsive_barneshut(Y, theta, 2);
        quad_const = attr * (sum(Pnz.*(qnz-log(qnz)-1)));
    case 'wtsne'
        qnz = 1./(1+sum((Y(I,:)-Y(J,:)).^2,2));
        A = attr * sparse(I,J,Pnz.*qnz,n,n);
        [repu_obj, G] = compute_wtsne_obj_grad_repulsive_barneshut(Y, weights, theta, 2);
        quad_const = attr * (sum(Pnz.*(qnz-log(qnz)-1)));
    case 'nerv'
        lambda = varargin{1};
        epsilon = varargin{2};
        A = attr * 0.5 * (P+P');
        [obj, grad] = compute_nerv_obj_grad_barneshut(Y, P', lambda, epsilon, theta, 2);
        repu_obj = obj - trace(Y'*GraphLaplacian(P+P')*Y);
        G = reshape(grad, n,2) - 4 * GraphLaplacian(A) * Y / attr;
        quad_const = 0;
    case 'mdsks'
        dist2nz = sum((Y(I,:)-Y(J,:)).^2,2);
        qnz = exp(-dist2nz);
        pqg = Pnz.*qnz;
        A = attr * sparse(I,J,pqg/sum(pqg),n,n);
        [repu_obj, G] = compute_mdsks_obj_grad_repulsive_barneshut(Y, theta, 2);
        qpq = qnz./sum(pqg);
        quad_const = attr * sum(Pnz.*(qpq.*log(qpq)));
    otherwise
        error('unknown method');
end