Description:
            Takes in a set of rankings and computes the model
            parameters for a mixture of Plackett-Luce models.
        Parameters:
            rankings:   set of rankings to aggregate
            K:          number of mixture components to compute
            epsilon:    convergence condition threshold value for overall EM algorithm
            epsilon_mm: convergence condition threshold value for MM algorithm
            iters:      dict, iterations configuration for EM and MM algorithms

def aggregate(self, rankings, K, epsilon, epsilon_mm, iters):
        """
        Description:
            Takes in a set of rankings and computes the model
            parameters for a mixture of Plackett-Luce models.
        Parameters:
            rankings:   set of rankings to aggregate
            K:          number of mixture components to compute
            epsilon:    convergence condition threshold value for overall EM algorithm
            epsilon_mm: convergence condition threshold value for MM algorithm
            iters:      dict, iterations configuration for EM and MM algorithms
        """
        x = rankings # shorter pseudonym for voting data
        self.n = len(rankings) # number of votes

        # "fixed" iterations type variables
        outer_iters = None
        inner_iters = None
        inner_range = None

        # Additional "scaling" iterations type variables
        inner_iters_base = None
        scaling_divisor = None

        # Additional "total" iterations type variables
        total_iters = None
        isIncremented = False

        if "type" not in iters:
            raise ValueError("iters dict must contain key \"type\"")
        iters_type = iters["type"]
        if iters_type == "fixed":
            outer_iters = iters["em_iters"]
            inner_iters = iters["mm_iters"]
        elif iters_type == "scaling":
            outer_iters = iters["em_iters"]
            inner_iters_base = iters["mm_iters_base"]
            scaling_divisor = iters["scaling_divisor"]
        elif iters_type == "total":
            total_iters = iters["total_iters"]
            outer_iters = iters["em_iters"]
            inner_iters = total_iters // outer_iters
        else:
            raise ValueError("iters dict value for key \"type\" is invalid: " + str(iters_type))

        # pre-compute the delta values
        delta_i_j_s = np.empty((self.n, self.m, self.m + 1))
        for i in range(self.n):
            for j in range(self.m):
                for s in range(self.m + 1):
                    delta_i_j_s[i][j][s] = EMMMixPLAggregator.delta(x[i], j, s, self.m)

        # generate initial values for p and pi:
        p_h0 = np.random.rand(K, self.m)
        p_h0 /= np.sum(p_h0, axis=1, keepdims=True)

        pi_h0 =  np.random.rand(K)
        pi_h0 /= np.sum(pi_h0)

        p_h = np.copy(p_h0)
        pi_h = np.copy(pi_h0)

        for g in range(outer_iters):

            p_h1 = np.empty((K, self.m))
            pi_h1 = np.empty(K)
            z_h1 = np.empty((self.n, K))

            # E-Step:
            self._EStep(K, x, z_h1, pi_h, p_h)

            # M-Step:
            if iters_type == "fixed":
                inner_range = range(inner_iters)
            elif iters_type == "total" and not isIncremented:
                test = (g + 1) * inner_iters + (outer_iters - g - 1) * (inner_iters + 1)
                if test < total_iters:
                    inner_iters += 1
                    isIncremented = True
                inner_range = range(inner_iters)
            elif iters_type == "scaling":
                inner_range = range(int(g/scaling_divisor) + inner_iters_base)

            for l in inner_range:
                self._MStep(l, K, x, z_h1, pi_h1, p_h1, p_h, delta_i_j_s)

                if (epsilon_mm != None and
                    np.all(np.absolute(p_h1 - p_h) < epsilon_mm)):
                        break

                p_h = np.copy(p_h1) # deep copy p for next MM iteration
                # pi does not change across MM iterations, no copy needed

            if (epsilon != None and
                np.all(np.absolute(p_h1 - p_h) < epsilon) and
                np.all(np.absolute(pi_h1 - pi_h) < epsilon)):
                    break

            # remember that assignments below are references only, not copies
            p_h = p_h1
            pi_h = pi_h1

        return (pi_h1, p_h1, pi_h0, p_h0)