function indx = get_target_index(obj)
% Choose which subfunction to update in certain learning iteration
    gd = find((obj.eval_count < 2) & obj.active);  % If an active subfunction has less than two observations, then 
% evaluate it again. Two evaluations per subfunction are needed to estimate its Hessian; 
    if ~isempty(gd)
        indx = gd(randperm(length(gd), 1));
        return
    end
    % choose the subfunction evaluated farthest using either total weighted Hessian, or by the subfunction Hessian
    if randn() < 0
        max_dist = -1;
        indx = -1;
        for i = 1:obj.N
            dtheta = obj.theta_proj - obj.last_theta(:,i);
            bdtheta = obj.b(:,:,i).' * dtheta;
            dist = sum(bdtheta.^2) + sum(dtheta.^2)*obj.min_eig_sub(i);
            if (dist > max_dist) && obj.active(i)
                max_dist = dist;
                indx = i;
            end
        end
    else
        dtheta = bsxfun(@plus, obj.theta_proj, -obj.last_theta); % difference between current theta & all subfunctions recent evaluation;
        full_H_combined = obj.get_full_H_with_diagonal(); % full Hessian;
        distance = sum(dtheta.*(full_H_combined * dtheta), 1); % squared distance;
        [~, dist_ord] = sort(-distance); % sort the distances from largest to smallest;
        dist_ord = dist_ord(obj.active(dist_ord)); % keep indices of active subfunctions;
        indx = dist_ord(1); % choose the active subfunction that is farthest from current location;;
    end
end
