r"""
        Makes a counter for central moments (n_counter) and a counter for raw moment (k_counter).
        Each is a list of :class:`~means.approximation.ode_problem.Moment`s.
        Therefore, each :class:`~means.approximation.ode_problem.Moments` is represented by both
        a vector of integer and a symbol.

        :param max_order: the maximal order of moment to be computer (will generate a list of moments up to `max_order + 1`)
        :param species: the name of the species
        :return: a pair of lists of :class:`~means.core.descriptors.Moment`s corresponding to central,
        and raw moments, respectively.
        :rtype: (list[:class:`~mea ns.core.descriptors.Moment`],list[:class:`~mea ns.core.descriptors.Moment`])

def generate_n_and_k_counters(max_order, species, central_symbols_prefix="M_", raw_symbols_prefix="x_"):
        r"""
        Makes a counter for central moments (n_counter) and a counter for raw moment (k_counter).
        Each is a list of :class:`~means.approximation.ode_problem.Moment`s.
        Therefore, each :class:`~means.approximation.ode_problem.Moments` is represented by both
        a vector of integer and a symbol.

        :param max_order: the maximal order of moment to be computer (will generate a list of moments up to `max_order + 1`)
        :param species: the name of the species
        :return: a pair of lists of :class:`~means.core.descriptors.Moment`s corresponding to central,
        and raw moments, respectively.
        :rtype: (list[:class:`~mea ns.core.descriptors.Moment`],list[:class:`~mea ns.core.descriptors.Moment`])
        """
        n_moments = max_order + 1
        # first order moments are always 1
        k_counter = [Moment([0] * len(species), sp.Integer(1))]
        n_counter = [Moment([0] * len(species), sp.Integer(1))]

        # build descriptors for first order raw moments aka expectations (e.g. [1, 0, 0], [0, 1, 0] and [0, 0, 1])
        descriptors = []
        for i in range(len(species)):
            row = [0]*len(species)
            row[i] = 1
            descriptors.append(row)

        # We use species name as symbols for first order raw moment
        k_counter += [Moment(d, s) for d,s in zip(descriptors, species)]

        # Higher order raw moment descriptors
        k_counter_descriptors = [i for i in itertools.product(range(n_moments + 1), repeat=len(species))
                                 if 1 < sum(i) <= n_moments]

        #this mimics the order in the original code
        k_counter_descriptors = sorted(k_counter_descriptors,lambda x,y: sum(x) - sum(y))
        #k_counter_descriptors = [[r for r in reversed(k)] for k in k_counter_descriptors]

        k_counter_symbols = [sp.Symbol(raw_symbols_prefix + "_".join([str(s) for s in count]))
                             for count in k_counter_descriptors]
        k_counter += [Moment(d, s) for d,s in zip(k_counter_descriptors, k_counter_symbols)]

        #  central moments
        n_counter_descriptors = [m for m in k_counter_descriptors if sum(m) > 1]
        # arbitrary symbols
        n_counter_symbols = [sp.Symbol(central_symbols_prefix + "_".join([str(s) for s in count]))
                             for count in n_counter_descriptors]

        n_counter += [Moment(c, s) for c,s in zip(n_counter_descriptors, n_counter_symbols)]

        return n_counter, k_counter