Vote for neighbors not being outliers. The key method of SOS.

@param ignore Object to ignore
@param di Neighbor object IDs.
@param p Probabilities
@param norm Normalization factor (1/sum)
@param scores Output score storage
771918
public static void nominateNeighbors(DBIDIter ignore, DBIDArrayIter di, double[] p, double norm, WritableDoubleDataStore scores) {
    for(di.seek(0); di.valid(); di.advance()) {
      if(DBIDUtil.equal(ignore, di)) {
        continue;
      }
      double v = p[di.getOffset()] * norm; // Normalize
      if(!(v > 0)) {
        break;
      }
      scores.increment(di, FastMath.log1p(-v));
    }
  }
347014
protected void computeScoreFive( int top[] , int middle[] , int bottom[] , int score[] , int width ) {

		// disparity as the outer loop to maximize common elements in inner loops, reducing redundant calculations
		for( int d = minDisparity; d < maxDisparity; d++ ) {

			// take in account the different in image border between the sub-regions and the effective region
			int indexSrc = (d-minDisparity)*width + (d-minDisparity) + radiusX;
			int indexDst = (d-minDisparity)*width + (d-minDisparity);
			int end = indexSrc + (width-d-4*radiusX);
			while( indexSrc < end ) {
				int s = 0;

				// sample four outer regions at the corners around the center region
				int val0 = top[indexSrc-radiusX];
				int val1 = top[indexSrc+radiusX];
				int val2 = bottom[indexSrc-radiusX];
				int val3 = bottom[indexSrc+radiusX];

				// select the two best scores from outer for regions
				if( val1 < val0 ) {
					int temp = val0;
					val0 = val1;
					val1 = temp;
				}

				if( val3 < val2 ) {
					int temp = val2;
					val2 = val3;
					val3 = temp;
				}

				if( val3 < val0 ) {
					s += val2;
					s += val3;
				} else if( val2 < val1 ) {
					s += val2;
					s += val0;
				} else {
					s += val0;
					s += val1;
				}

				score[indexDst++] = s + middle[indexSrc++];
			}
		}
	}
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);
  }
18144
public OutlierResult run(Database database, Relation<O> relation) {
    DistanceQuery<O> distFunc = database.getDistanceQuery(relation, getDistanceFunction());
    KNNQuery<O> knnQuery = database.getKNNQuery(distFunc, k);

    // track the maximum value for normalization
    DoubleMinMax ldofminmax = new DoubleMinMax();
    // compute the ldof values
    WritableDoubleDataStore ldofs = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_TEMP);

    // compute LOF_SCORE of each db object
    if(LOG.isVerbose()) {
      LOG.verbose("Computing LDOFs");
    }
    FiniteProgress progressLDOFs = LOG.isVerbose() ? new FiniteProgress("LDOF for objects", relation.size(), LOG) : null;

    Mean dxp = new Mean(), Dxp = new Mean();
    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      KNNList neighbors = knnQuery.getKNNForDBID(iditer, k);
      dxp.reset();
      Dxp.reset();
      DoubleDBIDListIter neighbor1 = neighbors.iter(),
          neighbor2 = neighbors.iter();
      for(; neighbor1.valid(); neighbor1.advance()) {
        // skip the point itself
        if(DBIDUtil.equal(neighbor1, iditer)) {
          continue;
        }
        dxp.put(neighbor1.doubleValue());
        for(neighbor2.seek(neighbor1.getOffset() + 1); neighbor2.valid(); neighbor2.advance()) {
          // skip the point itself
          if(DBIDUtil.equal(neighbor2, iditer)) {
            continue;
          }
          Dxp.put(distFunc.distance(neighbor1, neighbor2));
        }
      }
      double ldof = dxp.getMean() / Dxp.getMean();
      if(Double.isNaN(ldof) || Double.isInfinite(ldof)) {
        ldof = 1.0;
      }
      ldofs.putDouble(iditer, ldof);
      // update maximum
      ldofminmax.put(ldof);

      LOG.incrementProcessed(progressLDOFs);
    }
    LOG.ensureCompleted(progressLDOFs);

    // Build result representation.
    DoubleRelation scoreResult = new MaterializedDoubleRelation("LDOF Outlier Score", "ldof-outlier", ldofs, relation.getDBIDs());
    OutlierScoreMeta scoreMeta = new QuotientOutlierScoreMeta(ldofminmax.getMin(), ldofminmax.getMax(), 0.0, Double.POSITIVE_INFINITY, LDOF_BASELINE);
    return new OutlierResult(scoreMeta, scoreResult);
  }
480845
public boolean refine(List<AssociatedTriple> listObs , DMatrixRMaj P2 , DMatrixRMaj P3 )
	{
		CommonOps_DDRM.setIdentity(P1);
		initializeStructure(listObs, P2, P3);

		if( scale ) {
			scaler.applyScale(structure,observations);
		}

//		sba.setVerbose(System.out,0);
		sba.setParameters(structure,observations);
		sba.configure(converge.ftol,converge.gtol,converge.maxIterations);

		if( !sba.optimize(structure) ) {
			return false;
		}

		// save the results
		P2.set(structure.views[1].worldToView);
		P3.set(structure.views[2].worldToView);

		if( scale ) {
			// don't use built in unscaling function because it undoes scaling on points. Those are disposable
			scaler.pixelScaling.get(1).remove(P2,P2);
			scaler.pixelScaling.get(2).remove(P3,P3);
		}

		return true;
	}
304644
def _outlier_rejection(self, params, model, signal, ii):
        
        # Z score across repetitions:
        z_score = (params - np.mean(params, 0))/np.std(params, 0)
        # Silence warnings: 
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            outlier_idx = np.where(np.abs(z_score)>3.0)[0]
            nan_idx = np.where(np.isnan(params))[0]
            outlier_idx = np.unique(np.hstack([nan_idx, outlier_idx]))
            ii[outlier_idx] = 0
            model[outlier_idx] = np.nan
            signal[outlier_idx] = np.nan
            params[outlier_idx] = np.nan

        return model, signal, params, ii
552215
public double outresScore(final int s, long[] subspace, DBIDRef id, KernelDensityEstimator kernel, DBIDs cands) {
    double score = 1.0; // Initial score is 1.0
    final SubspaceEuclideanDistanceFunction df = new SubspaceEuclideanDistanceFunction(subspace);
    MeanVariance meanv = new MeanVariance();
    ModifiableDoubleDBIDList neighcand = DBIDUtil.newDistanceDBIDList(cands.size());
    ModifiableDoubleDBIDList nn = DBIDUtil.newDistanceDBIDList(cands.size());

    for(int i = s; i < kernel.dim; i++) {
      assert !BitsUtil.get(subspace, i);
      BitsUtil.setI(subspace, i);
      df.setSelectedDimensions(subspace);
      final double adjustedEps = kernel.adjustedEps(kernel.dim);
      DoubleDBIDList neigh = initialRange(id, cands, df, adjustedEps * 2, kernel, neighcand);
      // Relevance test
      if(neigh.size() > 2) {
        if(relevantSubspace(subspace, neigh, kernel)) {
          final double density = kernel.subspaceDensity(subspace, neigh);
          // Compute mean and standard deviation for densities of neighbors.
          meanv.reset();
          for(DoubleDBIDListIter neighbor = neigh.iter(); neighbor.valid(); neighbor.advance()) {
            subsetNeighborhoodQuery(neighcand, neighbor, df, adjustedEps, kernel, nn);
            meanv.put(kernel.subspaceDensity(subspace, nn));
          }
          final double deviation = (meanv.getMean() - density) / (2. * meanv.getSampleStddev());
          // High deviation:
          if(deviation >= 1) {
            score *= density / deviation;
          }
          // Recursion
          score *= outresScore(i + 1, subspace, id, kernel, neighcand);
        }
      }
      BitsUtil.clearI(subspace, i);
    }
    return score;
  }
1145796
public OutlierResult run(Database database, Relation<NumberVector> relation) {
    final int dbdim = RelationUtil.dimensionality(relation);
    final int mindim = dbdim >> 1;
    final int maxdim = dbdim - 1;
    final Random rand = rnd.getSingleThreadedRandom();

    ArrayList<OutlierResult> results = new ArrayList<>(num);
    {
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("LOF iterations", num, LOG) : null;
      for(int i = 0; i < num; i++) {
        long[] dimset = randomSubspace(dbdim, mindim, maxdim, rand);
        SubspaceEuclideanDistanceFunction df = new SubspaceEuclideanDistanceFunction(dimset);
        LOF<NumberVector> lof = new LOF<>(k, df);

        // run LOF and collect the result
        OutlierResult result = lof.run(database, relation);
        results.add(result);
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
    }

    WritableDoubleDataStore scores = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_STATIC);
    DoubleMinMax minmax = new DoubleMinMax();
    if(breadth) {
      FiniteProgress cprog = LOG.isVerbose() ? new FiniteProgress("Combining results", relation.size(), LOG) : null;
      @SuppressWarnings("unchecked")
      Pair<DBIDIter, DoubleRelation>[] IDVectorOntoScoreVector = (Pair<DBIDIter, DoubleRelation>[]) new Pair[results.size()];

      // Mapping score-sorted DBID-Iterators onto their corresponding scores.
      // We need to initialize them now be able to iterate them "in parallel".
      {
        int i = 0;
        for(OutlierResult r : results) {
          IDVectorOntoScoreVector[i] = new Pair<DBIDIter, DoubleRelation>(r.getOrdering().order(relation.getDBIDs()).iter(), r.getScores());
          i++;
        }
      }

      // Iterating over the *lines* of the AS_t(i)-matrix.
      for(int i = 0; i < relation.size(); i++) {
        // Iterating over the elements of a line (breadth-first).
        for(Pair<DBIDIter, DoubleRelation> pair : IDVectorOntoScoreVector) {
          DBIDIter iter = pair.first;
          // Always true if every algorithm returns a complete result (one score
          // for every DBID).
          if(iter.valid()) {
            double score = pair.second.doubleValue(iter);
            if(Double.isNaN(scores.doubleValue(iter))) {
              scores.putDouble(iter, score);
              minmax.put(score);
            }
            iter.advance();
          }
          else {
            LOG.warning("Incomplete result: Iterator does not contain |DB| DBIDs");
          }
        }
        // Progress does not take the initial mapping into account.
        LOG.incrementProcessed(cprog);
      }
      LOG.ensureCompleted(cprog);
    }
    else {
      FiniteProgress cprog = LOG.isVerbose() ? new FiniteProgress("Combining results", relation.size(), LOG) : null;
      for(DBIDIter iter = relation.iterDBIDs(); iter.valid(); iter.advance()) {
        double sum = 0.0;
        for(OutlierResult r : results) {
          final double s = r.getScores().doubleValue(iter);
          if(!Double.isNaN(s)) {
            sum += s;
          }
        }
        scores.putDouble(iter, sum);
        minmax.put(sum);
        LOG.incrementProcessed(cprog);
      }
      LOG.ensureCompleted(cprog);
    }
    OutlierScoreMeta meta = new BasicOutlierScoreMeta(minmax.getMin(), minmax.getMax());
    DoubleRelation scoreres = new MaterializedDoubleRelation("Feature bagging", "fb-outlier", scores, relation.getDBIDs());
    return new OutlierResult(meta, scoreres);
  }
600728
private static void outsideAlgorithm(final double[][] scores, final ProjTreeChart inChart, final ProjTreeChart outChart, 
            boolean singleRoot) {             
        final int n = scores.length;
        final int startIdx = singleRoot ? 1 : 0;         

        if (singleRoot) {
            // Initialize.
            double goalScore = 0.0;
            outChart.updateCell(0, n-1, RIGHT, COMPLETE, goalScore, -1);
            
            // The inside algorithm is effectively doing this...
            //
            // wallScore[r] log+=  inChart.scores[0][r][LEFT][COMPLETE] +
            //                     inChart.scores[r][n-1][RIGHT][COMPLETE] + 
            //                     fracRoot[r];
            //
            // goalScore log+= wallScore[r];            
            for (int r=1; r<n; r++) {
                outChart.updateCell(0, r, RIGHT, INCOMPLETE, goalScore, -1);
            }
            
            // Un-build goal constituents by combining left and right complete
            // constituents, on the left and right respectively. This corresponds to
            // left and right triangles. (Note: this is the opposite of how we
            // build an incomplete constituent.)
            for (int r=1; r<n; r++) {
                // Left child.
                double leftScore = outChart.getScore(0, r, RIGHT, INCOMPLETE) + 
                                   inChart.getScore(r, n - 1, RIGHT, COMPLETE) + 
                                   scores[0][r];
                outChart.updateCell(1, r, LEFT, COMPLETE, leftScore, -1);
                // Right child.
                double rightScore = outChart.getScore(0, r, RIGHT, INCOMPLETE) + 
                                    inChart.getScore(1, r, LEFT, COMPLETE) + 
                                    scores[0][r];
                outChart.updateCell(r, n - 1, RIGHT, COMPLETE, rightScore, -1);
            }
        } else  {
            // Base case.
            for (int d=0; d<2; d++) {
                outChart.updateCell(0, n-1, d, COMPLETE, 0.0, -1);
            }
        }
        
        // Parse.
        for (int width = n - 1; width >= 1; width--) {
            for (int s = startIdx; s < n - width; s++) {
                int t = s + width;
                
                // First create complete items (opposite of inside).
                // -- Left side.
                for (int r=s; r<t; r++) {
                    final int d = LEFT;
                    // Left child.
                    double leftScore = outChart.getScore(s, t, d, COMPLETE) + 
                            inChart.getScore(r, t, d, INCOMPLETE);
                    outChart.updateCell(s, r, d, COMPLETE, leftScore, -1);
                    // Right child.
                    double rightScore = outChart.getScore(s, t, d, COMPLETE) + 
                            inChart.getScore(s, r, d, COMPLETE);
                    outChart.updateCell(r, t, d, INCOMPLETE, rightScore, -1);
                }
                // -- Right side.
                for (int r=s+1; r<=t; r++) {
                    final int d = RIGHT;
                    // Left child.
                    double leftScore = outChart.getScore(s, t, d, COMPLETE) + 
                            inChart.getScore(r, t, d, COMPLETE);
                    outChart.updateCell(s, r, d, INCOMPLETE, leftScore, -1);
                    // Right child.
                    double rightScore = outChart.getScore(s, t, d, COMPLETE) + 
                            inChart.getScore(s, r, d, INCOMPLETE);
                    outChart.updateCell(r, t, d, COMPLETE, rightScore, -1);
                }
                
                // Second create incomplete items (opposite of inside).
                for (int r=s; r<t; r++) {
                    for (int d=0; d<2; d++) {
                        double edgeScore = (d == LEFT) ? scores[t][s] : scores[s][t];
                        // Left child.
                        double leftScore = outChart.getScore(s, t, d, INCOMPLETE)
                                + inChart.getScore(r + 1, t, LEFT, COMPLETE) + edgeScore;
                        outChart.updateCell(s, r, RIGHT, COMPLETE, leftScore, -1);
                        // Right child.
                        double rightScore = outChart.getScore(s, t, d, INCOMPLETE)
                                + inChart.getScore(s, r, RIGHT, COMPLETE) + edgeScore;
                        outChart.updateCell(r + 1, t, LEFT, COMPLETE, rightScore, -1);
                    }
                }
            }
        }
    }
973071
public static int cudnnDivisiveNormalizationForward(
        cudnnHandle handle, 
        cudnnLRNDescriptor normDesc, 
        int mode, 
        Pointer alpha, 
        cudnnTensorDescriptor xDesc, /** same desc for means, temp, temp2 */
        Pointer x, 
        Pointer means, /** if NULL, means are assumed to be zero */
        Pointer temp, 
        Pointer temp2, 
        Pointer beta, 
        cudnnTensorDescriptor yDesc, 
        Pointer y)
    {
        return checkResult(cudnnDivisiveNormalizationForwardNative(handle, normDesc, mode, alpha, xDesc, x, means, temp, temp2, beta, yDesc, y));
    }
1102814
private void accumulateNormals(final RandomElementReader vertices,
            final Element face) throws IOException {
        int[] indices = face.getIntList("vertex_index");
        for (int i = 0; i < indices.length; i++) {
            int pre;
            int post;
            if (counterClockwise) {
                pre = (i + indices.length - 1) % indices.length;
                post = (i + 1) % indices.length;
            } else {
                pre = (i + 1) % indices.length;
                post = (i + indices.length - 1) % indices.length;
            }

            Element centerVertex;
            Element preVertex;
            Element postVertex;
            try {
                centerVertex = vertices.readElement(indices[i]);
                preVertex = vertices.readElement(indices[pre]);
                postVertex = vertices.readElement(indices[post]);
                accumulateNormal(centerVertex, preVertex, postVertex);
            } catch (IndexOutOfBoundsException e) {
                // we ignore defects in the normals.
            }
        }
    }
1020521
protected void checkNeighbors( int cx , int cy , ImageFlow.D flow , ImageFlow image , float score ) {
		for( int i = -regionRadius; i <= regionRadius; i++ ) {
			int index = image.width*(cy+i) + (cx-regionRadius);
			for( int j = -regionRadius; j <= regionRadius; j++ , index++ ) {
				float s = scores[ index ];
				ImageFlow.D f = image.data[index];
				if( s > score ) {
					f.set(flow);
					scores[index] = score;
				} else if( s == score ) {
					// Pick solution with the least motion when ambiguous
					float m0 = f.x*f.x + f.y*f.y;
					float m1 = flow.x*flow.x + flow.y*flow.y;
					if( m1 < m0 ) {
						f.set(flow);
						scores[index] = score;
					}
				}
			}
		}
	}
802836
private void normalizeL2(double[] scores) {
		double square_sum = 0.0;
		for (int i = 0; i < scores.length; i++) {
			square_sum += (scores[i] * scores[i]);
		}
		double norm = Math.sqrt(square_sum);
		if (norm != 0)
			for (int i = 0; i < scores.length; i++) {
				scores[i] = scores[i] / norm;
			}
	}
1348709
private static long specialNStar(final long numSamplesI, final double p, final double delta) {
    final double q, numSamplesF;
    double tot, curTerm;
    long m;
    assertTrue(numSamplesI >= 1);
    assertTrue((0.0 < p) && (p < 1.0));
    assertTrue((0.0 < delta) && (delta < 1.0));
    q = 1.0 - p;
    numSamplesF = numSamplesI;
    // Use a different algorithm if the following isn't true; this one will be too slow, or worse.
    assertTrue((numSamplesF / p) < 500.0);
    curTerm = Math.pow(p, numSamplesF);  // curTerm = posteriorProbability (k, k, p)
    assertTrue(curTerm > 1e-100); // sanity check for non-use of logarithms
    tot = curTerm;
    m = numSamplesI;
    while (tot <= delta) { // this test can fail even the first time
      curTerm = (curTerm * q * (m)) / ((m + 1) - numSamplesI);
      tot += curTerm;
      m += 1;
    }
    // we have reached a state where tot > delta, so back up one
    return (m - 1);
  }
154361
public OutlierResult run(Database database, Relation<O> relation) {
    StepProgress stepprog = LOG.isVerbose() ? new StepProgress("LOF", 3) : null;
    DBIDs ids = relation.getDBIDs();

    LOG.beginStep(stepprog, 1, "Materializing nearest-neighbor sets.");
    KNNQuery<O> knnq = DatabaseUtil.precomputedKNNQuery(database, relation, getDistanceFunction(), k);

    // Compute LRDs
    LOG.beginStep(stepprog, 2, "Computing Local Reachability Densities (LRD).");
    WritableDoubleDataStore lrds = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_TEMP);
    computeLRDs(knnq, ids, lrds);

    // compute LOF_SCORE of each db object
    LOG.beginStep(stepprog, 3, "Computing Local Outlier Factors (LOF).");
    WritableDoubleDataStore lofs = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_DB);
    // track the maximum value for normalization.
    DoubleMinMax lofminmax = new DoubleMinMax();
    computeLOFScores(knnq, ids, lrds, lofs, lofminmax);

    LOG.setCompleted(stepprog);

    // Build result representation.
    DoubleRelation scoreResult = new MaterializedDoubleRelation("Local Outlier Factor", "lof-outlier", lofs, ids);
    OutlierScoreMeta scoreMeta = new QuotientOutlierScoreMeta(lofminmax.getMin(), lofminmax.getMax(), 0.0, Double.POSITIVE_INFINITY, 1.0);
    return new OutlierResult(scoreMeta, scoreResult);
  }
844819
@Override
	public void findNeighbor(P target, int searchN, FastQueue<KdTreeResult> results) {
		if( searchN <= 0 )
			throw new IllegalArgumentException("I'm sorry, but I refuse to search for less than or equal to 0 neighbors.");

		if( tree.root == null )
			return;

		this.searchN = searchN;
		this.target = target;
		this.mostDistantNeighborSq = maxDistanceSq;

		stepClosest(tree.root,results);
	}
234467
public OutlierResult run(Database database, Relation<O> relation) {
    StepProgress stepprog = LOG.isVerbose() ? new StepProgress(5) : null;

    Pair<KNNQuery<O>, KNNQuery<O>> pair = getKNNQueries(database, relation, stepprog);
    KNNQuery<O> knnComp = pair.getFirst();
    KNNQuery<O> knnReach = pair.getSecond();

    // Assert we got something
    if(knnComp == null) {
      throw new AbortException("No kNN queries supported by database for comparison distance function.");
    }
    if(knnReach == null) {
      throw new AbortException("No kNN queries supported by database for density estimation distance function.");
    }

    // FIXME: tie handling!

    // Probabilistic distances
    WritableDoubleDataStore pdists = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_DB);
    LOG.beginStep(stepprog, 3, "Computing pdists");
    computePDists(relation, knnReach, pdists);
    // Compute PLOF values.
    WritableDoubleDataStore plofs = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_TEMP);
    LOG.beginStep(stepprog, 4, "Computing PLOF");
    double nplof = computePLOFs(relation, knnComp, pdists, plofs);

    // Normalize the outlier scores.
    DoubleMinMax mm = new DoubleMinMax();
    {// compute LOOP_SCORE of each db object
      LOG.beginStep(stepprog, 5, "Computing LoOP scores");

      FiniteProgress progressLOOPs = LOG.isVerbose() ? new FiniteProgress("LoOP for objects", relation.size(), LOG) : null;
      final double norm = 1. / (nplof * MathUtil.SQRT2);
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        double loop = NormalDistribution.erf((plofs.doubleValue(iditer) - 1.) * norm);
        plofs.putDouble(iditer, loop);
        mm.put(loop);
        LOG.incrementProcessed(progressLOOPs);
      }
      LOG.ensureCompleted(progressLOOPs);
    }

    LOG.setCompleted(stepprog);

    // Build result representation.
    DoubleRelation scoreResult = new MaterializedDoubleRelation("Local Outlier Probabilities", "loop-outlier", plofs, relation.getDBIDs());
    OutlierScoreMeta scoreMeta = new ProbabilisticOutlierScore(mm.getMin(), mm.getMax(), 0.);
    return new OutlierResult(scoreMeta, scoreResult);
  }