Get the marginal distribution for a node.
332798
def marginal(repertoire, node_index):
    """"""
    index = tuple(i for i in range(repertoire.ndim) if i != node_index)

    return repertoire.sum(index, keepdims=True)
987219
def _call_marginalizevperp(self,o,integrate_method='dopr54_c',**kwargs):
        """"""
        #Get d, l, vlos
        l= o.ll(obs=[1.,0.,0.],ro=1.)*_DEGTORAD
        vlos= o.vlos(ro=1.,vo=1.,obs=[1.,0.,0.,0.,0.,0.])
        R= o.R(use_physical=False)
        phi= o.phi(use_physical=False)
        #Get local circular velocity, projected onto the los
        if isinstance(self._pot,list):
            vcirc= calcRotcurve([p for p in self._pot if not p.isNonAxi],R)[0]
        else:
            vcirc= calcRotcurve(self._pot,R)[0]
        vcirclos= vcirc*math.sin(phi+l)
        #Marginalize
        alphalos= phi+l
        if not 'nsigma' in kwargs or ('nsigma' in kwargs and \
                                          kwargs['nsigma'] is None):
            nsigma= _NSIGMA
        else:
            nsigma= kwargs['nsigma']
        kwargs.pop('nsigma',None)
        #BOVY: add asymmetric drift here?
        if math.fabs(math.sin(alphalos)) < math.sqrt(1./2.):
            sigmaR1= nu.sqrt(self._initdf.sigmaT2(R,phi=phi,
                                                  use_physical=False)) #Slight abuse
            cosalphalos= math.cos(alphalos)
            tanalphalos= math.tan(alphalos)
            return integrate.quad(_marginalizeVperpIntegrandSinAlphaSmall,
                                  -nsigma,nsigma,
                                  args=(self,R,cosalphalos,tanalphalos,
                                        vlos-vcirclos,vcirc,
                                        sigmaR1,phi),
                                  **kwargs)[0]/math.fabs(cosalphalos)*sigmaR1
        else:
            sigmaR1= nu.sqrt(self._initdf.sigmaR2(R,phi=phi,
                                                  use_physical=False))
            sinalphalos= math.sin(alphalos)
            cotalphalos= 1./math.tan(alphalos)
            return integrate.quad(_marginalizeVperpIntegrandSinAlphaLarge,
                                  -nsigma,nsigma,
                                  args=(self,R,sinalphalos,cotalphalos,
                                        vlos-vcirclos,vcirc,sigmaR1,phi),
                                  **kwargs)[0]/math.fabs(sinalphalos)*sigmaR1
987402
private final double getNodeParentsAlpha(Node<?> node)
    {
        double alpha = 1;

        node = node.getParent();

        while (null != node)
        {
            alpha = alpha * node.getAttributes().getAlpha();

            node = node.getParent();

            if ((null != node) && (node.getNodeType() == NodeType.LAYER))
            {
                node = null;
            }
        }
        return alpha;
    }
1163795
def visit_Constant(self, node):
        """"""
        nnode = self.dnode(node)
        node.first_line, node.first_col = ast_pos(nnode, self.bytes_pos_to_utf8)
        node.last_line = node.first_line
        node.last_col = node.first_col + len(repr(node.value))
        node.uid = (node.last_line, node.last_col)
1272287
private CfgParseChart marginal(CfgParseChart chart, List<?> terminals,
      Factor rootDist) {
    initializeChart(chart, terminals);
    upwardChartPass(chart);
    // Set the initial outside probabilities
    Factor root = rootDistribution.product(rootDist);
    chart.updateOutsideEntry(0, chart.chartSize() - 1, 
        root.coerceToDiscrete().getWeights().getValues(), root, parentVar);
    downwardChartPass(chart);
    return chart;
  }
913519
def get_node (self, url_data):
        """"""
        if not url_data.url:
            return None
        elif url_data.url in self.nodes:
            return None
        node = {
            "url": url_data.url,
            "parent_url": url_data.parent_url,
            "id": self.nodeid,
            "label": quote(url_data.title if url_data.title else url_data.name),
            "extern": 1 if url_data.extern else 0,
            "checktime": url_data.checktime,
            "size": url_data.size,
            "dltime": url_data.dltime,
            "edge": quote(url_data.name),
            "valid": 1 if url_data.valid else 0,
        }
        self.nodes[node["url"]] = node
        self.nodeid += 1
        return node
492282
private int findOptimalNodeOrder(final Node<K, V> node) {
		int low = MIN_B_ORDER; // minimal b-order
		int high = (blockSize / node.getStructEstimateSize(1)) << 2; // estimate high b-order

		while (low <= high) {
			int mid = ((low + high) >>> 1);
			mid += (1 - (mid % 2));
			int nodeSize = node.getStructEstimateSize(mid);

			if (log.isDebugEnabled()) {
				log.debug(this.getClass().getName() + "::findOptimalNodeOrder(" + node.getClass().getName()
						+ ") blockSize=" + blockSize + " nodeSize=" + nodeSize + " b_low=" + low
						+ " b_order=" + mid + " b_high=" + high);
			}

			if (nodeSize < blockSize) {
				low = mid + 2;
			} else if (nodeSize > blockSize) {
				high = mid - 2;
			} else {
				return mid;
			}
		}
		return low - 2;
	}
812455
def _logprior(self):
        """"""
        logj = self.logjacobian
        logp = self.prior_distribution(**self.current_params) + logj
        if numpy.isnan(logp):
            logp = -numpy.inf
        return logp
700205
def get_node(self, node_name):
        """"""
        for node in self.nodes:
            if node.__name__ == node_name:
                return node
414321
def Marginal(self, i, name=''):
        
        pmf = Pmf(name=name)
        for vs, prob in self.Items():
            pmf.Incr(vs[i], prob)
        return pmf
803996
def get_lower_bound(self):
        
        lower_bounds = []

        for distribution in self.distribs.values():
            lower_bound = distribution.percent_point(distribution.mean / 10000)
            if not pd.isnull(lower_bound):
                lower_bounds.append(lower_bound)

        return min(lower_bounds)
359462
private Node<T> min(Node<T> x) {
    	do {
    		if (x.left == null) return x;
    		x = x.left;
    	} while (true);
    }
357605
def probability_density(self, X):
        """"""
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 1:
            return np.multiply(U, V)

        else:
            a = np.power(np.multiply(U, V), -1)
            tmp = np.power(-np.log(U), self.theta) + np.power(-np.log(V), self.theta)
            b = np.power(tmp, -2 + 2.0 / self.theta)
            c = np.power(np.multiply(np.log(U), np.log(V)), self.theta - 1)
            d = 1 + (self.theta - 1) * np.power(tmp, -1.0 / self.theta)
            return self.cumulative_distribution(X) * a * b * c * d
150372
def collapse_nodes(self, min_dist=1e-6, min_support=0):
        
        nself = self.copy()
        for node in nself.treenode.traverse():
            if not node.is_leaf():
                if (node.dist <= min_dist) | (node.support < min_support):
                    node.delete()
        nself._coords.update()
        return nself
969104
def _posterior(self, x):
        
        yedge = self._nuis_pdf.marginalization_bins()
        yc = 0.5 * (yedge[1:] + yedge[:-1])
        yw = yedge[1:] - yedge[:-1]

        like_array = self.like(x[:, np.newaxis], yc[np.newaxis, :]) * yw
        like_array /= like_array.sum()

        self._post = like_array.sum(1)
        self._post_interp = castro.Interpolator(x, self._post)
        return self._post
113198
private MarginalResult messagePassing(MarginalizationMethod marginalize, boolean includeJointMarginalsAndPartition) {

    // Using the behavior of brute force factor multiplication as ground truth, the desired
    // outcome of marginal calculation with an impossible factor is a uniform probability dist.,
    // since we have a resulting factor of all 0s. That is of course assuming that normalizing
    // all 0s gives you uniform, which is not real math, but that's a useful tolerance to include, so we do.

    boolean impossibleObservationMade = false;

    // Message passing will look at fully observed cliques as non-entities, but their
    // log-likelihood (the log-likelihood of the single observed value) is still relevant for the
    // partition function.

    double partitionFunction = 1.0;

    if (includeJointMarginalsAndPartition) {

      // Find only factors which are fully observed
      outer:
      for (GraphicalModel.Factor f : model.factors) {
        for (int n : f.neigborIndices) {
          if (!model.getVariableMetaDataByReference(n).containsKey(VARIABLE_OBSERVED_VALUE)) continue outer;
        }

        int[] assignment = new int[f.neigborIndices.length];
        for (int i = 0; i < f.neigborIndices.length; i++) {
          assignment[i] = Integer.parseInt(model.getVariableMetaDataByReference(f.neigborIndices[i]).get(VARIABLE_OBSERVED_VALUE));
        }

        double assignmentValue = f.getAssignmentValue(assignment, weights);
        if (Double.isInfinite(assignmentValue)) {
          impossibleObservationMade = true;
        } else {
          assert !Double.isNaN(assignmentValue);
          partitionFunction *= Math.exp(assignmentValue);
          assert !Double.isNaN(partitionFunction);
        }
      }
    }

    // Create the cliques by multiplying out table factors
    // TODO:OPT This could be made more efficient by observing first, then dot product

    List<TableFactor> cliquesList = new ArrayList<>();
    Map<Integer, GraphicalModel.Factor> cliqueToFactor = new HashMap<>();

    for (GraphicalModel.Factor f : model.factors) {
      boolean allObserved = true;
      int maxVar = 0;
      for (int n : f.neigborIndices) {
        if (!model.getVariableMetaDataByReference(n).containsKey(VARIABLE_OBSERVED_VALUE)) allObserved = false;
        if (n > maxVar) maxVar = n;
      }
      if (allObserved) continue;

      TableFactor clique = null;

      // Retrieve cache if exists and none of the observations have changed

      if (cachedFactors.containsKey(f)) {
        CachedFactorWithObservations obs = cachedFactors.get(f);
        boolean allConsistent = true;
        for (int i = 0; i < f.neigborIndices.length; i++) {
          int n = f.neigborIndices[i];
          if (model.getVariableMetaDataByReference(n).containsKey(VARIABLE_OBSERVED_VALUE) &&
              (obs.observations[i] == -1 ||
                  Integer.parseInt(model.getVariableMetaDataByReference(n).get(VARIABLE_OBSERVED_VALUE)) != obs.observations[i])) {
            allConsistent = false;
            break;
          }
          // NOTE: This disqualifies lots of stuff for some reason...
          if (!model.getVariableMetaDataByReference(n).containsKey(VARIABLE_OBSERVED_VALUE) && (obs.observations[i] != -1)) {
            allConsistent = false;
            break;
          }
        }
        if (allConsistent) {
          clique = obs.cachedFactor;
          if (obs.impossibleObservation) {
            impossibleObservationMade = true;
          }
        }
      }

      // Otherwise make a new cache

      if (clique == null) {
        int[] observations = new int[f.neigborIndices.length];
        for (int i = 0; i < observations.length; i++) {
          Map<String, String> metadata = model.getVariableMetaDataByReference(f.neigborIndices[i]);
          if (metadata.containsKey(VARIABLE_OBSERVED_VALUE)) {
            int value = Integer.parseInt(metadata.get(VARIABLE_OBSERVED_VALUE));
            observations[i] = value;
          } else {
            observations[i] = -1;
          }
        }

        clique = new TableFactor(weights, f, observations);

        CachedFactorWithObservations cache = new CachedFactorWithObservations();
        cache.cachedFactor = clique;
        cache.observations = observations;

        // Check for an impossible observation
        boolean nonZeroValue = false;
        for (int[] assignment : clique) {
          double assignmentValue = clique.getAssignmentValue(assignment);
          assert assignmentValue >= 0.0;
          if (assignmentValue > 0) {
            nonZeroValue = true;
            break;
          }
        }
        if (!nonZeroValue) {
          impossibleObservationMade = true;
          cache.impossibleObservation = true;
        }

        cachedFactors.put(f, cache);
      }

      cliqueToFactor.put(cliquesList.size(), f);
      cliquesList.add(clique);
    }

    TableFactor[] cliques = cliquesList.toArray(new TableFactor[cliquesList.size()]);

    // If we made any impossible observations, we can just return a uniform distribution for all the variables that
    // weren't observed, since that's the semantically correct thing to do (our 'probability' is broken at this
    // point).

    if (impossibleObservationMade) {
      int maxVar = 0;
      for (TableFactor c : cliques) {
        for (int i : c.neighborIndices) if (i > maxVar) maxVar = i;
      }

      double[][] result = new double[maxVar + 1][];

      for (TableFactor c : cliques) {
        for (int i = 0; i < c.neighborIndices.length; i++) {
          result[c.neighborIndices[i]] = new double[c.getDimensions()[i]];
          for (int j = 0; j < result[c.neighborIndices[i]].length; j++) {
            result[c.neighborIndices[i]][j] = 1.0 / result[c.neighborIndices[i]].length;
            assert !Double.isNaN(result[c.neighborIndices[i]][j]);
          }
        }
      }

      // Create a bunch of uniform joint marginals, constrained by observations, and fill up the joint marginals
      // with them

      Map<GraphicalModel.Factor, TableFactor> jointMarginals = new IdentityHashMap<>();
      if (includeJointMarginalsAndPartition) {
        for (GraphicalModel.Factor f : model.factors) {
          TableFactor uniformZero = new TableFactor(f.neigborIndices, f.getDimensions());

          for (int[] assignment : uniformZero) {
            uniformZero.setAssignmentValue(assignment, 0.0);
          }

          jointMarginals.put(f, uniformZero);
        }
      }

      return new MarginalResult(result, 1.0, jointMarginals, new int[result.length]);
    }

    // Find the number of variables in the model, so that we can size arrays appropriately

    int maxVar = 0;
    for (GraphicalModel.Factor fac : model.factors) {
      for (int i : fac.neigborIndices) if (i > maxVar) maxVar = i;
    }

    // Indexed by (start-clique, end-clique), this array will remain mostly null in most graphs

    TableFactor[][] messages = new TableFactor[cliques.length][cliques.length];
    boolean[][] backwardPassedMessages = new boolean[cliques.length][cliques.length];

    // Create the data structures to hold the tree pattern

    boolean[] visited = new boolean[cliques.length];
    int numVisited = 0;
    int[] visitedOrder = new int[cliques.length];

    int[] parent = new int[cliques.length];
    for (int i = 0; i < parent.length; i++) parent[i] = -1;
    // Figure out which cliques are connected to which trees. This is important for calculating the partition
    // function later, since each tree will converge to its own partition function by multiplication, and we will
    // need to multiply the partition function of each of the trees to get the global one.
    int[] trees = new int[cliques.length];

    // Forward pass, record a BFS forest pattern that we can use for message passing

    int treeIndex = -1;
    boolean[] seenVariable = new boolean[maxVar + 1];
    while (numVisited < cliques.length) {
      treeIndex++;

      // Pick the largest connected graph remaining as the root for message passing

      int root = -1;
      for (int i = 0; i < cliques.length; i++) {
        if (!visited[i] &&
            (root == -1 || cliques[i].neighborIndices.length > cliques[root].neighborIndices.length)) {
          root = i;
        }
      }
      assert (root != -1);

      Queue<Integer> toVisit = new ArrayDeque<>();
      toVisit.add(root);
      boolean[] toVisitArray = new boolean[cliques.length];
      toVisitArray[root] = true;

      while (toVisit.size() > 0) {
        int cursor = toVisit.poll();
        // toVisitArray[cursor] = false;
        trees[cursor] = treeIndex;
        if (visited[cursor]) {
          log.info("Visited contains: " + cursor);
          log.info("Visited: " + Arrays.toString(visited));
          log.info("To visit: " + toVisit);
        }
        assert (!visited[cursor]);
        visited[cursor] = true;
        visitedOrder[numVisited] = cursor;
        for (int i : cliques[cursor].neighborIndices) seenVariable[i] = true;
        numVisited++;

        childLoop:
        for (int i = 0; i < cliques.length; i++) {
          if (i == cursor) continue;
          if (i == parent[cursor]) continue;
          if (domainsOverlap(cliques[cursor], cliques[i])) {

            // Make sure that for every variable that we've already seen somewhere in the graph, if it's
            // in the child, it's in the parent. Otherwise we'll break the property of continuous
            // transmission of information about variables through messages.

            childNeighborLoop:
            for (int child : cliques[i].neighborIndices) {
              if (seenVariable[child]) {
                for (int j : cliques[cursor].neighborIndices) {
                  if (j == child) {
                    continue childNeighborLoop;
                  }
                }
                // If we get here it means that this clique is not good as a child, since we can't pass
                // it all the information it needs from other elements of the tree
                continue childLoop;
              }
            }

            if (parent[i] == -1 && !visited[i]) {
              if (!toVisitArray[i]) {
                toVisit.add(i);
                toVisitArray[i] = true;
                for (int j : cliques[i].neighborIndices) seenVariable[j] = true;
              }
              parent[i] = cursor;
            }
          }
        }
      }
      // No cycles in the tree
      assert (parent[root] == -1);
    }

    assert (numVisited == cliques.length);

    // Backward pass, run the visited list in reverse

    for (int i = numVisited - 1; i >= 0; i--) {
      int cursor = visitedOrder[i];
      if (parent[cursor] == -1) continue;

      backwardPassedMessages[cursor][parent[cursor]] = true;

      // Calculate the message to the clique's parent, given all incoming messages so far

      TableFactor message = cliques[cursor];
      for (int k = 0; k < cliques.length; k++) {
        if (k == parent[cursor]) continue;
        if (messages[k][cursor] != null) {
          message = message.multiply(messages[k][cursor]);
        }
      }

      // Marginalize all irrelevant variables, and we have our message

      messages[cursor][parent[cursor]] = marginalizeMessage(message, cliques[parent[cursor]].neighborIndices, marginalize);
    }

    // Setup for MAP calculation

    int[] bestGlobalAssignment = new int[maxVar+1];
    for (int i = 0; i < maxVar+1; i++) {
      Map<String, String> metadata = model.getVariableMetaDataByReference(i);
      if (metadata.containsKey(VARIABLE_OBSERVED_VALUE)) {
        int value = Integer.parseInt(metadata.get(VARIABLE_OBSERVED_VALUE));
        bestGlobalAssignment[i] = value;
      } else {
        bestGlobalAssignment[i] = -1;
      }
    }

    // Forward pass, run the visited list forward

    for (int i = 0; i < numVisited; i++) {
      int cursor = visitedOrder[i];

      // For each child of the cursor
      for (int j = 0; j < cliques.length; j++) {
        if (parent[j] != cursor) continue;

        // For every incoming message through the cursor that isn't coming from j, which is a child, multiply the
        // incoming message
        TableFactor message = cliques[cursor];
        for (int k = 0; k < cliques.length; k++) {
          if (k == j) continue;
          if (messages[k][cursor] != null) {
            message = message.multiply(messages[k][cursor]);
          }
        }

        messages[cursor][j] = marginalizeMessage(message, cliques[j].neighborIndices, marginalize);
      }

      // If we're also looking for the MAP estimate, then we need to also do a version of this for the cursor, where we
      // multiply in any messages that are coming from everyone except the parent

      if (marginalize == MarginalizationMethod.MAX) {
        TableFactor mapFactor = cliques[cursor];
        for (int j = 0; j < cliques.length; j++) {
          if (parent[cursor] == j) continue;
          if (messages[j][cursor] != null) {
            mapFactor = mapFactor.multiply(messages[j][cursor]);
          }
        }

        // OPTIMIZATION:
        // Rather than use the standard iterator, which creates lots of int[] arrays on the heap, which need to be GC'd,
        // we use the fast version that just mutates one array. Since this is read once for us here, this is ideal.

        int[] bestAssignment = null;
        double bestAssignmentValue = Double.NEGATIVE_INFINITY;

        Iterator<int[]> fastPassByReferenceIterator = mapFactor.fastPassByReferenceIterator();
        do {
          int[] assignment = fastPassByReferenceIterator.next();
          double v = mapFactor.getAssignmentValue(assignment);
          if (v > bestAssignmentValue) {
            // Check if this is consistent with our global state
            boolean isConsistent = true;
            for (int k = 0; k < assignment.length; k++) {
              if (bestGlobalAssignment[mapFactor.neighborIndices[k]] != -1 &&
                  bestGlobalAssignment[mapFactor.neighborIndices[k]] != assignment[k]) {
                isConsistent = false;
                break;
              }
            }
            if (isConsistent) {
              bestAssignmentValue = v;
              bestAssignment = assignment.clone();
            }
          }
        } while (fastPassByReferenceIterator.hasNext());

        assert(bestAssignment != null);
        for (int k = 0; k < bestAssignment.length; k++) {
          assert(bestGlobalAssignment[mapFactor.neighborIndices[k]] == -1 ||
              bestGlobalAssignment[mapFactor.neighborIndices[k]] == bestAssignment[k]);
          bestGlobalAssignment[mapFactor.neighborIndices[k]] = bestAssignment[k];
        }
      }
    }

    // Calculate final marginals for each variable

    double[][] marginals = new double[maxVar + 1][];

    // Include observed variables as deterministic

    for (GraphicalModel.Factor fac : model.factors) {
      for (int i = 0; i < fac.neigborIndices.length; i++) {
        int n = fac.neigborIndices[i];
        if (model.getVariableMetaDataByReference(n).containsKey(VARIABLE_OBSERVED_VALUE)) {
          double[] deterministic = new double[fac.getDimensions()[i]];
          int assignment = Integer.parseInt(model.getVariableMetaDataByReference(n).get(VARIABLE_OBSERVED_VALUE));
          if (assignment > deterministic.length) {
            throw new IllegalStateException("Variable " + n + ": Can't have as assignment (" + assignment + ") that is out of bounds for dimension size (" + deterministic.length + ")");
          }
          deterministic[assignment] = 1.0;
          marginals[n] = deterministic;
        }
      }
    }

    Map<GraphicalModel.Factor, TableFactor> jointMarginals = new IdentityHashMap<>();

    if (includeJointMarginalsAndPartition) {
      boolean[] partitionIncludesTrees = new boolean[treeIndex + 1];
      double[] treePartitionFunctions = new double[treeIndex + 1];

      for (int i = 0; i < cliques.length; i++) {
        TableFactor convergedClique = cliques[i];

        for (int j = 0; j < cliques.length; j++) {
          if (i == j) continue;
          if (messages[j][i] == null) continue;
          convergedClique = convergedClique.multiply(messages[j][i]);
        }

        // It doesn't make sense to calculate a partition function for MAP message passing

        if (marginalize == MarginalizationMethod.SUM) {
          // Calculate the partition function when we're calculating marginals
          // We need one contribution per tree in our forest graph

          if (!partitionIncludesTrees[trees[i]]) {
            partitionIncludesTrees[trees[i]] = true;
            treePartitionFunctions[trees[i]] = convergedClique.valueSum();
            assert !Double.isNaN(treePartitionFunctions[trees[i]]);
            partitionFunction *= treePartitionFunctions[trees[i]];
          } else {

            // This is all just an elaborate assert
            // Check that our partition function is the same as the trees we're attached to, or with %.1, for numerical reasons.
            // Sometimes the partition function will explode in value, which can make a non-%-based assert worthless here

            if (assertsEnabled() && !TableFactor.USE_EXP_APPROX) {
              double valueSum = convergedClique.valueSum();
              if (Double.isFinite(valueSum) && Double.isFinite(treePartitionFunctions[trees[i]])) {
                if (Math.abs(treePartitionFunctions[trees[i]] - valueSum) >= 1.0e-3 * treePartitionFunctions[trees[i]]) {
                  log.info("Different partition functions for tree " + trees[i] + ": ");
                  log.info("Pre-existing for tree: " + treePartitionFunctions[trees[i]]);
                  log.info("This clique for tree: " + valueSum);
                }
                assert (Math.abs(treePartitionFunctions[trees[i]] - valueSum) < 1.0e-3 * treePartitionFunctions[trees[i]]);
              }
            }
          }
        }

        // Calculate the factor this clique corresponds to, and put in an entry for joint marginals

        GraphicalModel.Factor f = cliqueToFactor.get(i);
        assert (f != null);
        if (!jointMarginals.containsKey(f)) {
          int[] observedAssignments = getObservedAssignments(f);

          // Collect back pointers and check if this factor matches the clique we're using

          int[] backPointers = new int[observedAssignments.length];
          int cursor = 0;
          for (int j = 0; j < observedAssignments.length; j++) {
            if (observedAssignments[j] == -1) {
              backPointers[j] = cursor;
              cursor++;
            }
            // This is not strictly necessary but will trigger array OOB exception if things go wrong, so is nice
            else backPointers[j] = -1;
          }

          double sum = convergedClique.valueSum();
          assert(sum > 0);

          TableFactor jointMarginal = new TableFactor(f.neigborIndices, f.getDimensions());

          // OPTIMIZATION:
          // Rather than use the standard iterator, which creates lots of int[] arrays on the heap, which need to be GC'd,
          // we use the fast version that just mutates one array. Since this is read once for us here, this is ideal.
          Iterator<int[]> fastPassByReferenceIterator = convergedClique.fastPassByReferenceIterator();
          do {
            int[] assignment = fastPassByReferenceIterator.next();
            if (backPointers.length == assignment.length) {
              // Normalize if we're doing SUM (inference) message passing
              if (marginalize == MarginalizationMethod.SUM) {
                jointMarginal.setAssignmentValue(assignment, convergedClique.getAssignmentValue(assignment) / sum);
              }
              // Don't normalize if we're doing MAX (MAP) message passing
              else if (marginalize == MarginalizationMethod.MAX) {
                jointMarginal.setAssignmentValue(assignment, convergedClique.getAssignmentValue(assignment));
              }
            } else {
              int[] jointAssignment = new int[backPointers.length];
              for (int j = 0; j < jointAssignment.length; j++) {
                if (observedAssignments[j] != -1) jointAssignment[j] = observedAssignments[j];
                else jointAssignment[j] = assignment[backPointers[j]];
              }
              // Normalize if we're doing SUM (inference) message passing
              if (marginalize == MarginalizationMethod.SUM) {
                jointMarginal.setAssignmentValue(jointAssignment, convergedClique.getAssignmentValue(assignment) / sum);
              }
              // Don't normalize if we're doing MAX (MAP) message passing
              else if (marginalize == MarginalizationMethod.MAX) {
                jointMarginal.setAssignmentValue(jointAssignment, convergedClique.getAssignmentValue(assignment));
              }
            }

            // Set the assignment arrays correctly
          } while (fastPassByReferenceIterator.hasNext());

          jointMarginals.put(f, jointMarginal);
        }

        boolean anyNull = false;
        for (int j = 0; j < convergedClique.neighborIndices.length; j++) {
          int k = convergedClique.neighborIndices[j];
          if (marginals[k] == null) {
            anyNull = true;
          }
        }

        if (anyNull) {
          double[][] cliqueMarginals = null;
          switch (marginalize) {
            case SUM:
              cliqueMarginals = convergedClique.getSummedMarginals();
              break;
            case MAX:
              cliqueMarginals = convergedClique.getMaxedMarginals();
              break;
          }
          for (int j = 0; j < convergedClique.neighborIndices.length; j++) {
            int k = convergedClique.neighborIndices[j];
            if (marginals[k] == null) {
              marginals[k] = cliqueMarginals[j];
            }
          }
        }
      }
    }
    // If we don't care about joint marginals, we can be careful about not calculating more cliques than we need to,
    // by explicitly sorting by which cliques are most profitable to calculate over. In this way we can avoid, in
    // the case of a chain CRF, calculating almost half the joint factors.
    else {
      // First do a pass where we only calculate all-null neighbors
      for (int i = 0; i < cliques.length; i++) {
        boolean allNull = true;
        for (int k : cliques[i].neighborIndices) {
          if (marginals[k] != null) allNull = false;
        }
        if (allNull) {
          TableFactor convergedClique = cliques[i];

          for (int j = 0; j < cliques.length; j++) {
            if (i == j) continue;
            if (messages[j][i] == null) continue;
            convergedClique = convergedClique.multiply(messages[j][i]);
          }

          double[][] cliqueMarginals = null;
          switch (marginalize) {
            case SUM:
              cliqueMarginals = convergedClique.getSummedMarginals();
              break;
            case MAX:
              cliqueMarginals = convergedClique.getMaxedMarginals();
              break;
          }
          for (int j = 0; j < convergedClique.neighborIndices.length; j++) {
            int k = convergedClique.neighborIndices[j];
            if (marginals[k] == null) {
              marginals[k] = cliqueMarginals[j];
            }
          }
        }
      }
      // Now we calculate any remaining cliques with any non-null variables
      for (int i = 0; i < cliques.length; i++) {
        boolean anyNull = false;
        for (int j = 0; j < cliques[i].neighborIndices.length; j++) {
          int k = cliques[i].neighborIndices[j];
          if (marginals[k] == null) {
            anyNull = true;
          }
        }

        if (anyNull) {
          TableFactor convergedClique = cliques[i];

          for (int j = 0; j < cliques.length; j++) {
            if (i == j) continue;
            if (messages[j][i] == null) continue;
            convergedClique = convergedClique.multiply(messages[j][i]);
          }

          double[][] cliqueMarginals = null;
          switch (marginalize) {
            case SUM:
              cliqueMarginals = convergedClique.getSummedMarginals();
              break;
            case MAX:
              cliqueMarginals = convergedClique.getMaxedMarginals();
              break;
          }
          for (int j = 0; j < convergedClique.neighborIndices.length; j++) {
            int k = convergedClique.neighborIndices[j];
            if (marginals[k] == null) {
              marginals[k] = cliqueMarginals[j];
            }
          }
        }
      }
    }

    // Add any factors to the joint marginal map that were fully observed and so didn't get cliques
    if (marginalize == MarginalizationMethod.SUM && includeJointMarginalsAndPartition) {
      for (GraphicalModel.Factor f : model.factors) {
        if (!jointMarginals.containsKey(f)) {
          // This implies that every variable in the factor is observed. If that's the case, we need to construct
          // a one hot TableFactor representing the deterministic distribution.
          TableFactor deterministicJointMarginal = new TableFactor(f.neigborIndices, f.getDimensions());
          int[] observedAssignment = getObservedAssignments(f);
          for (int i : observedAssignment) assert (i != -1);
          deterministicJointMarginal.setAssignmentValue(observedAssignment, 1.0);

          jointMarginals.put(f, deterministicJointMarginal);
        }
      }
    }

    assert !Double.isNaN(partitionFunction);
    return new MarginalResult(marginals, partitionFunction, jointMarginals, bestGlobalAssignment);
  }
553588
def average(self, rho):
        """
        def marginal(f, rho):
            remaining = len(f.shape)
            if remaining == 0:
                return rho
            rho = sum([f[i]*rho[i] for i in range(rho.shape[0])])
            f = np.sum(f, 0)
            return marginal(f, rho)

        return marginal(self.distribution, rho)