gt:         if (existsKnowledge()) {
context: ///////////////////////////////////////////////////////////////////////////////
  // For information as to what this class does, see the Javadoc, below.       //
  // Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006,       //
  // 2007, 2008, 2009, 2010, 2014, 2015 by Peter Spirtes, Richard Scheines, Joseph   //
  // Ramsey, and Clark Glymour.                                                //
  //                                                                           //
  // This program is free software; you can redistribute it and/or modify      //
  // it under the terms of the GNU General Public License as published by      //
  // the Free Software Foundation; either version 2 of the License, or         //
  // (at your option) any later version.                                       //
  //                                                                           //
  // This program is distributed in the hope that it will be useful,           //
  // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
  // GNU General Public License for more details.                              //
  //                                                                           //
  // You should have received a copy of the GNU General Public License         //
  // along with this program; if not, write to the Free Software               //
  // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA //
  ///////////////////////////////////////////////////////////////////////////////
  
  package edu.cmu.tetrad.search;
  
  import edu.cmu.tetrad.data.*;
  import edu.cmu.tetrad.graph.*;
  import edu.cmu.tetrad.util.ChoiceGenerator;
  import edu.cmu.tetrad.util.ForkJoinPoolInstance;
  import edu.cmu.tetrad.util.TaskManager;
  import edu.cmu.tetrad.util.TetradLogger;
  
  import java.io.PrintStream;
  import java.text.DecimalFormat;
  import java.text.NumberFormat;
  import java.util.*;
  import java.util.concurrent.*;
  
  
  /**
   * GesSearch is an implementation of the GES algorithm, as specified in Chickering (2002) "Optimal structure
   * identification with greedy search" Journal of Machine Learning Research. It works for both BayesNets and SEMs.
   * <p>
   * Some code optimization could be done for the scoring part of the graph for discrete models (method scoreGraphChange).
   * Some of Andrew Moore's approaches for caching sufficient statistics, for instance.
   * <p>
   * To speed things up, it has been assumed that variables X and Y with zero correlation do not correspond to edges in
   * the graph. This is a restricted form of the heuristicSpeedup assumption, something GES does not assume. This
   * the graph. This is a restricted form of the heuristicSpeedup assumption, something GES does not assume. This
   * heuristicSpeedup assumption needs to be explicitly turned on using setHeuristicSpeedup(true).
   * <p>
   * A number of other optimizations were added 5/2015. See code for details.
   *
   * @author Ricardo Silva, Summer 2003
   * @author Joseph Ramsey, Revisions 5/2015
   */
  public final class Fgs2 implements GraphSearch, GraphScorer {
  
  
      /**
       * Internal.
       */
      private enum Mode {
          allowUnfaithfulness, heuristicSpeedup, coverNoncolliders
      }
  
      /**
       * Specification of forbidden and required edges.
       */
      private IKnowledge knowledge = new Knowledge2();
  
      /**
       * List of variables in the data set, in order.
       */
      private List<Node> variables;
  
      /**
       * The true graph, if known. If this is provided, asterisks will be printed out next to false positive added edges
       * (that is, edges added that aren't adjacencies in the true graph).
       */
      private Graph trueGraph;
  
      /**
       * An initial graph to start from.
       */
      private Graph initialGraph;
  
      /**
       * If non-null, edges not adjacent in this graph will not be added.
       */
      private Graph boundGraph = null;
  
      /**
       * Elapsed time of the most recent search.
       */
      private long elapsedTime;
  
      /**
       * A bound on cycle length.
       */
      private int cycleBound = -1;
  
      /**
       * The totalScore for discrete searches.
       */
      private Score score;
  
      /**
       * The logger for this class. The config needs to be set.
       */
      private TetradLogger logger = TetradLogger.getInstance();
  
      /**
       * The top n graphs found by the algorithm, where n is numPatternsToStore.
       */
      private LinkedList<ScoredGraph> topGraphs = new LinkedList<>();
  
      /**
       * The number of top patterns to store.
       */
      private int numPatternsToStore = 0;
  
      /**
       * True if verbose output should be printed.
       */
      private boolean verbose = false;
  
      // Potential arrows sorted by bump high to low. The first one is a candidate for adding to the graph.
      private SortedSet<Arrow> sortedArrows = null;
  
      // Arrows added to sortedArrows for each <i, j>.
      private Map<OrderedPair<Node>, Set<Arrow>> lookupArrows = null;
  
      // A utility map to help with orientation.
      private Map<Node, Set<Node>> neighbors = null;
  
      // Map from variables to their column indices in the data set.
      private ConcurrentMap<Node, Integer> hashIndices;
  
      // The static ForkJoinPool instance.
      private ForkJoinPool pool = ForkJoinPoolInstance.getInstance().getPool();
  
      // A running tally of the total BIC totalScore.
      private double totalScore;
  
      // A graph where X--Y means that X and Y have non-zero total effect on one another.
      private Graph effectEdgesGraph;
  
      // The minimum number of operations to do before parallelizing.
      private final int minChunk = 100;
  
      // Where printed output is sent.
      private PrintStream out = System.out;
  
      // A initial adjacencies graph.
      private Graph adjacencies = null;
  
      // The graph being constructed.
      private Graph graph;
  
      // Arrows with the same totalScore are stored in this list to distinguish their order in sortedArrows.
      // The ordering doesn't matter; it just have to be transitive.
      int arrowIndex = 0;
  
      // The final totalScore after search.
      private double modelScore;
  
      // Internal.
      private Mode mode = Mode.heuristicSpeedup;
  
      /**
       * True if one-edge faithfulness is assumed. Speeds the algorithm up.
       */
      private boolean faithfulnessAssumed = true;
  
      // Bounds the indegree of the graph.
      private int maxIndegree = -1;
  
      final int maxThreads = ForkJoinPoolInstance.getInstance().getPool().getParallelism();
  
      //===========================CONSTRUCTORS=============================//
  
      /**
       * Construct a Score and pass it in here. The totalScore should return a
       * positive value in case of conditional dependence and a negative
       * values in case of conditional independence. See Chickering (2002),
       * locally consistent scoring criterion.
       */
      public Fgs2(Score score) {
          if (score == null) throw new NullPointerException();
          setScore(score);
          this.graph = new EdgeListGraphSingleConnections(getVariables());
      }
  
      //==========================PUBLIC METHODS==========================//
  
      /**
       * Set to true if it is assumed that all path pairs with one length 1 path do not cancel.
       */
      public void setFaithfulnessAssumed(boolean faithfulnessAssumed) {
          this.faithfulnessAssumed = true;
      }
  
      /**
       * @return true if it is assumed that all path pairs with one length 1 path do not cancel.
       */
      public boolean isFaithfulnessAssumed() {
          return faithfulnessAssumed;
      }
  
      /**
       * Greedy equivalence search: Start from the empty graph, add edges till model is significant. Then start deleting
       * edges till a minimum is achieved.
       *
       * @return the resulting Pattern.
       */
      public Graph search() {
          topGraphs.clear();
  
          lookupArrows = new ConcurrentHashMap<>();
          final List<Node> nodes = new ArrayList<>(variables);
          graph = new EdgeListGraphSingleConnections(nodes);
  
          if (adjacencies != null) {
              adjacencies = GraphUtils.replaceNodes(adjacencies, nodes);
          }
  
          if (initialGraph != null) {
              graph = new EdgeListGraphSingleConnections(initialGraph);
              graph = GraphUtils.replaceNodes(graph, nodes);
          }
  
          addRequiredEdges(graph);
  
          if (faithfulnessAssumed) {
              initializeForwardEdgesFromEmptyGraph(getVariables());
  
              // Do forward search.
              this.mode = Mode.heuristicSpeedup;
              fes();
              bes();
  
              this.mode = Mode.coverNoncolliders;
              initializeTwoStepEdges(getVariables());
              fes();
              bes();
          } else {
              initializeForwardEdgesFromEmptyGraph(getVariables());
  
              // Do forward search.
              this.mode = Mode.heuristicSpeedup;
              fes();
              bes();
  
              this.mode = Mode.allowUnfaithfulness;
              initializeForwardEdgesFromExistingGraph(getVariables());
              fes();
              bes();
          }
  
          long start = System.currentTimeMillis();
          totalScore = 0.0;
  
          long endTime = System.currentTimeMillis();
          this.elapsedTime = endTime - start;
          this.logger.log("graph", "\nReturning this graph: " + graph);
  
          this.logger.log("info", "Elapsed time = " + (elapsedTime) / 1000. + " s");
          this.logger.flush();
  
          this.modelScore = totalScore;
  
          return graph;
      }
  
      /**
       * @return the background knowledge.
       */
  
      public IKnowledge getKnowledge() {
          return knowledge;
      }
  
      /**
       * Sets the background knowledge.
       *
       * @param knowledge the knowledge object, specifying forbidden and required edges.
       */
      public void setKnowledge(IKnowledge knowledge) {
          if (knowledge == null) throw new NullPointerException();
          this.knowledge = knowledge;
      }
  
      public long getElapsedTime() {
          return elapsedTime;
      }
  
      /**
       * If the true graph is set, askterisks will be printed in log output for the true edges.
       */
      public void setTrueGraph(Graph trueGraph) {
          this.trueGraph = trueGraph;
      }
  
      /**
       * @return the totalScore of the given DAG, up to a constant.
       */
      public double getScore(Graph dag) {
          return scoreDag(dag);
      }
  
      /**
       * @return the list of top scoring graphs.
       */
      public LinkedList<ScoredGraph> getTopGraphs() {
          return topGraphs;
      }
  
      /**
       * @return the number of patterns to store.
       */
      public int getNumPatternsToStore() {
          return numPatternsToStore;
      }
  
      /**
       * Sets the number of patterns to store. This should be set to zero for fast search.
       */
      public void setNumPatternsToStore(int numPatternsToStore) {
          if (numPatternsToStore < 0) {
              throw new IllegalArgumentException("# graphs to store must at least 0: " + numPatternsToStore);
          }
  
          this.numPatternsToStore = numPatternsToStore;
      }
  
      /**
       * @return the initial graph for the search. The search is initialized to this graph and
       * proceeds from there.
       */
      public Graph getInitialGraph() {
          return initialGraph;
      }
  
      /**
       * Sets the initial graph.
       */
      public void setInitialGraph(Graph initialGraph) {
          if (initialGraph != null) {
              initialGraph = GraphUtils.replaceNodes(initialGraph, variables);
  
              if (verbose) {
                  out.println("Initial graph variables: " + initialGraph.getNodes());
                  out.println("Data set variables: " + variables);
              }
  
              if (!new HashSet<>(initialGraph.getNodes()).equals(new HashSet<>(variables))) {
                  throw new IllegalArgumentException("Variables aren't the same.");
              }
          }
  
          this.initialGraph = initialGraph;
      }
  
      /**
       * Sets whether verbose output should be produced.
       */
      public void setVerbose(boolean verbose) {
          this.verbose = verbose;
      }
  
      /**
       * Sets the output stream that output (except for log output) should be sent to.
       * By detault System.out.
       */
      public void setOut(PrintStream out) {
          this.out = out;
      }
  
      /**
       * @return the output stream that output (except for log output) should be sent to.
       */
      public PrintStream getOut() {
          return out;
      }
  
      /**
       * @return the set of preset adjacenies for the algorithm; edges not in this adjacencies graph
       * will not be added.
       */
      public Graph getAdjacencies() {
          return adjacencies;
      }
  
      /**
       * Sets the set of preset adjacenies for the algorithm; edges not in this adjacencies graph
       * will not be added.
       */
      public void setAdjacencies(Graph adjacencies) {
          this.adjacencies = adjacencies;
      }
  
      /**
       * A bound on cycle length.
       */
      public int getCycleBound() {
          return cycleBound;
      }
  
      /**
       * A bound on cycle length.
       *
       * @param cycleBound The bound, >= 1, or -1 for unlimited.
       */
      public void setCycleBound(int cycleBound) {
          if (!(cycleBound == -1 || cycleBound >= 1))
              throw new IllegalArgumentException("Cycle bound needs to be -1 or >= 1: " + cycleBound);
          this.cycleBound = cycleBound;
      }
  
      /**
       * Creates a new processors pool with the specified number of threads.
       */
      public void setParallelism(int numProcessors) {
          this.pool = new ForkJoinPool(numProcessors);
      }
  
      /**
       * If non-null, edges not adjacent in this graph will not be added.
       */
      public void setBoundGraph(Graph boundGraph) {
          this.boundGraph = GraphUtils.replaceNodes(boundGraph, getVariables());
      }
  
      /**
       * For BIC totalScore, a multiplier on the penalty term. For continuous searches.
       *
       * @deprecated Use the getters on the individual scores instead.
       */
      public double getPenaltyDiscount() {
          if (score instanceof ISemBicScore) {
              return ((ISemBicScore) score).getPenaltyDiscount();
          } else {
              return 2.0;
          }
      }
  
      /**
       * @deprecated Use the setters on the individual scores instead.
       */
      public void setSamplePrior(double samplePrior) {
          if (score instanceof LocalDiscreteScore) {
              ((LocalDiscreteScore) score).setSamplePrior(samplePrior);
          }
      }
  
      /**
       * @deprecated Use the setters on the individual scores instead.
       */
      public void setStructurePrior(double expectedNumParents) {
          if (score instanceof LocalDiscreteScore) {
              ((LocalDiscreteScore) score).setStructurePrior(expectedNumParents);
          }
      }
  
      /**
       * For BIC totalScore, a multiplier on the penalty term. For continuous searches.
       *
       * @deprecated Use the setters on the individual scores instead.
       */
      public void setPenaltyDiscount(double penaltyDiscount) {
          if (score instanceof ISemBicScore) {
              ((ISemBicScore) score).setPenaltyDiscount(penaltyDiscount);
          }
      }
  
      /**
       * The maximum of parents any nodes can have in output pattern.
       * @return -1 for unlimited.
       */
      public int getMaxIndegree() {
          return maxIndegree;
      }
  
      /**
       * The maximum of parents any nodes can have in output pattern.
       * @param maxIndegree -1 for unlimited.
       */
      public void setMaxIndegree(int maxIndegree) {
          if (maxIndegree < -1) throw new IllegalArgumentException();
          this.maxIndegree = maxIndegree;
      }
  
      //===========================PRIVATE METHODS========================//
  
      //Sets the discrete scoring function to use.
      private void setScore(Score totalScore) {
          this.score = totalScore;
  
          this.variables = new ArrayList<>();
  
          for (Node node : totalScore.getVariables()) {
              if (node.getNodeType() == NodeType.MEASURED) {
                  this.variables.add(node);
              }
          }
  
          buildIndexing(totalScore.getVariables());
  
          this.maxIndegree = score.getMaxIndegree();
      }
  
      final int[] count = new int[1];
  
      public int getMinChunk(int n) {
          return Math.max(n / maxThreads, minChunk);
      }
  
      class NodeTaskEmptyGraph extends RecursiveTask<Boolean> {
          private final int from;
          private final int to;
          private final List<Node> nodes;
          private final Set<Node> emptySet;
  
          public NodeTaskEmptyGraph(int from, int to, List<Node> nodes, Set<Node> emptySet) {
              this.from = from;
              this.to = to;
              this.nodes = nodes;
              this.emptySet = emptySet;
          }
  
          @Override
          protected Boolean compute() {
              for (int i = from; i < to; i++) {
                  if ((i + 1) % 1000 == 0) {
                      count[0] += 1000;
                      out.println("Initializing effect edges: " + (count[0]));
                  }
  
                  Node y = nodes.get(i);
                  neighbors.put(y, emptySet);
  
                  for (int j = i + 1; j < nodes.size(); j++) {
                      Node x = nodes.get(j);
  
                      if (existsKnowledge()) {
                          if (getKnowledge().isForbidden(x.getName(), y.getName()) && getKnowledge().isForbidden(y.getName(), x.getName())) {
                              continue;
                          }
  
                          if (!validSetByKnowledge(y, emptySet)) {
                              continue;
                          }
                      }
  
                      if (adjacencies != null && !adjacencies.isAdjacentTo(x, y)) {
                          continue;
                      }
  
                      int child = hashIndices.get(y);
                      int parent = hashIndices.get(x);
                      double bump = score.localScoreDiff(parent, child);
  
                      if (boundGraph != null && !boundGraph.isAdjacentTo(x, y)) continue;
  
                      if (bump > 0) {
                          final Edge edge = Edges.undirectedEdge(x, y);
                          effectEdgesGraph.addEdge(edge);
                      }
  
                      if (bump > 0.0) {
                          addArrow(x, y, emptySet, emptySet, bump);
                          addArrow(y, x, emptySet, emptySet, bump);
                      }
                  }
              }
  
              return true;
          }
      }
  
      private void initializeForwardEdgesFromEmptyGraph(final List<Node> nodes) {
  //        if (verbose) {
  //            System.out.println("heuristicSpeedup = true");
  //        }
  
          sortedArrows = new ConcurrentSkipListSet<>();
          lookupArrows = new ConcurrentHashMap<>();
          neighbors = new ConcurrentHashMap<>();
          final Set<Node> emptySet = new HashSet<>();
  
          long start = System.currentTimeMillis();
          this.effectEdgesGraph = new EdgeListGraphSingleConnections(nodes);
  
          class InitializeFromEmptyGraphTask extends RecursiveTask<Boolean> {
  
              public InitializeFromEmptyGraphTask() {
              }
  
              @Override
              protected Boolean compute() {
                  Queue<NodeTaskEmptyGraph> tasks = new ArrayDeque<>();
  
                  int numNodesPerTask = Math.max(100, nodes.size() / maxThreads);
  
                  for (int i = 0; i < nodes.size(); i += numNodesPerTask) {
                      NodeTaskEmptyGraph task = new NodeTaskEmptyGraph(i, Math.min(nodes.size(), i + numNodesPerTask),
                              nodes, emptySet);
                      tasks.add(task);
                      task.fork();
  
                      for (NodeTaskEmptyGraph _task : new ArrayList<>(tasks)) {
                          if (_task.isDone()) {
                              _task.join();
                              tasks.remove(_task);
                          }
                      }
  
                      while (tasks.size() > maxThreads) {
                          NodeTaskEmptyGraph _task = tasks.poll();
                          _task.join();
                      }
                  }
  
                  for (NodeTaskEmptyGraph task : tasks) {
                      task.join();
                  }
  
                  return true;
              }
          }
  
          pool.invoke(new InitializeFromEmptyGraphTask());
  
          long stop = System.currentTimeMillis();
  
          if (verbose) {
              out.println("Elapsed initializeForwardEdgesFromEmptyGraph = " + (stop - start) + " ms");
          }
      }
  
      private void initializeTwoStepEdges(final List<Node> nodes) {
  //        if (verbose) {
  //            System.out.println("heuristicSpeedup = false");
  //        }
  
          count[0] = 0;
  
          sortedArrows = new ConcurrentSkipListSet<>();
          lookupArrows = new ConcurrentHashMap<>();
          neighbors = new ConcurrentHashMap<>();
  
          if (this.effectEdgesGraph == null) {
              this.effectEdgesGraph = new EdgeListGraph(nodes);
          }
  
          if (initialGraph != null) {
              for (Edge edge : initialGraph.getEdges()) {
                  if (!effectEdgesGraph.isAdjacentTo(edge.getNode1(), edge.getNode2())) {
                      effectEdgesGraph.addUndirectedEdge(edge.getNode1(), edge.getNode2());
                  }
              }
          }
  
          final Set<Node> emptySet = new HashSet<>(0);
  
          class InitializeFromExistingGraphTask extends RecursiveTask<Boolean> {
              private int chunk;
              private int from;
              private int to;
  
              public InitializeFromExistingGraphTask(int chunk, int from, int to) {
                  this.chunk = chunk;
                  this.from = from;
                  this.to = to;
              }
  
              @Override
              protected Boolean compute() {
                  if (TaskManager.getInstance().isCanceled()) return false;
  
                  if (to - from <= chunk) {
                      for (int i = from; i < to; i++) {
                          if ((i + 1) % 1000 == 0) {
                              count[0] += 1000;
                              out.println("Initializing effect edges: " + (count[0]));
                          }
  
                          Node y = nodes.get(i);
  
                          Set<Node> g = new HashSet<>();
  
                          for (Node n : graph.getAdjacentNodes(y)) {
                              for (Node m : graph.getAdjacentNodes(n)) {
                                  if (graph.isAdjacentTo(y, m)) {
                                      continue;
                                  }
  
                                  if (graph.isDefCollider(m, n, y)) {
                                      continue;
                                  }
  
                                  g.add(m);
                              }
                          }
  
                          for (Node x : g) {
                              if (existsKnowledge()) {
                                  if (getKnowledge().isForbidden(x.getName(), y.getName()) && getKnowledge().isForbidden(y.getName(), x.getName())) {
                                      continue;
                                  }
  
                                  if (!validSetByKnowledge(y, emptySet)) {
                                      continue;
                                  }
                              }
  
                              if (adjacencies != null && !adjacencies.isAdjacentTo(x, y)) {
                                  continue;
                              }
  
                              if (removedEdges.contains(Edges.undirectedEdge(x, y))) {
                                  continue;
                              }
  
                              calculateArrowsForward(x, y);
                          }
                      }
  
                      return true;
                  } else {
                      int mid = (to + from) / 2;
  
                      InitializeFromExistingGraphTask left = new InitializeFromExistingGraphTask(chunk, from, mid);
                      InitializeFromExistingGraphTask right = new InitializeFromExistingGraphTask(chunk, mid, to);
  
                      left.fork();
                      right.compute();
                      left.join();
  
                      return true;
                  }
              }
          }
  
          pool.invoke(new InitializeFromExistingGraphTask(getMinChunk(nodes.size()), 0, nodes.size()));
      }
  
      private void initializeForwardEdgesFromExistingGraph(final List<Node> nodes) {
  //        if (verbose) {
  //            System.out.println("heuristicSpeedup = false");
  //        }
  
          count[0] = 0;
  
          sortedArrows = new ConcurrentSkipListSet<>();
          lookupArrows = new ConcurrentHashMap<>();
          neighbors = new ConcurrentHashMap<>();
  
          if (this.effectEdgesGraph == null) {
              this.effectEdgesGraph = new EdgeListGraph(nodes);
          }
  
          if (initialGraph != null) {
              for (Edge edge : initialGraph.getEdges()) {
                  if (!effectEdgesGraph.isAdjacentTo(edge.getNode1(), edge.getNode2())) {
                      effectEdgesGraph.addUndirectedEdge(edge.getNode1(), edge.getNode2());
                  }
              }
          }
  
          final Set<Node> emptySet = new HashSet<>(0);
  
          class InitializeFromExistingGraphTask extends RecursiveTask<Boolean> {
              private int chunk;
              private int from;
              private int to;
  
              public InitializeFromExistingGraphTask(int chunk, int from, int to) {
                  this.chunk = chunk;
                  this.from = from;
                  this.to = to;
              }
  
              @Override
              protected Boolean compute() {
                  if (TaskManager.getInstance().isCanceled()) return false;
  
                  if (to - from <= chunk) {
                      for (int i = from; i < to; i++) {
                          if ((i + 1) % 1000 == 0) {
                              count[0] += 1000;
                              out.println("Initializing effect edges: " + (count[0]));
                          }
  
                          Node y = nodes.get(i);
                          Set<Node> D = new HashSet<>();
                          List<Node> cond = new ArrayList<>();
                          D.addAll(GraphUtils.getDconnectedVars(y, cond, graph));
                          D.remove(y);
                          D.removeAll(effectEdgesGraph.getAdjacentNodes(y));
  
                          for (Node x : D) {
                              if (existsKnowledge()) {
                                  if (getKnowledge().isForbidden(x.getName(), y.getName()) && getKnowledge().isForbidden(y.getName(), x.getName())) {
                                      continue;
                                  }
  
                                  if (!validSetByKnowledge(y, emptySet)) {
                                      continue;
                                  }
                              }
  
                              if (adjacencies != null && !adjacencies.isAdjacentTo(x, y)) {
                                  continue;
                              }
  
                              calculateArrowsForward(x, y);
                          }
                      }
  
                      return true;
                  } else {
                      int mid = (to + from) / 2;
  
                      InitializeFromExistingGraphTask left = new InitializeFromExistingGraphTask(chunk, from, mid);
                      InitializeFromExistingGraphTask right = new InitializeFromExistingGraphTask(chunk, mid, to);
  
                      left.fork();
                      right.compute();
                      left.join();
  
                      return true;
                  }
              }
          }
  
          pool.invoke(new InitializeFromExistingGraphTask(getMinChunk(nodes.size()), 0, nodes.size()));
      }
  
      private void fes() {
          TetradLogger.getInstance().log("info", "** FORWARD EQUIVALENCE SEARCH");
  
          while (!sortedArrows.isEmpty()) {
              Arrow arrow = sortedArrows.first();
              sortedArrows.remove(arrow);
  
              Node x = arrow.getA();
              Node y = arrow.getB();
  
              if (graph.isAdjacentTo(x, y)) {
                  continue;
              }
  
              if (!arrow.getNaYX().equals(getNaYX(x, y))) {
                  continue;
              }
  
              if (!getTNeighbors(x, y).containsAll(arrow.getHOrT())) {
                  continue;
              }
  
              if (!validInsert(x, y, arrow.getHOrT(), getNaYX(x, y))) {
                  continue;
              }
  
              Set<Node> T = arrow.getHOrT();
              double bump = arrow.getBump();
  
              boolean inserted = insert(x, y, T, bump);
              if (!inserted) continue;
  
              totalScore += bump;
  
              Set<Node> visited = reapplyOrientation(x, y, null);
              Set<Node> toProcess = new HashSet<>();
  
              for (Node node : visited) {
                  final Set<Node> neighbors1 = getNeighbors(node);
                  final Set<Node> storedNeighbors = this.neighbors.get(node);
  
                  if (!(neighbors1.equals(storedNeighbors))) {
                      toProcess.add(node);
                  }
              }
  
              toProcess.add(x);
              toProcess.add(y);
  
              storeGraph();
              reevaluateForward(toProcess, arrow);
          }
      }
  
      private void bes() {
          TetradLogger.getInstance().log("info", "** BACKWARD EQUIVALENCE SEARCH");
  
          sortedArrows = new ConcurrentSkipListSet<>();
          lookupArrows = new ConcurrentHashMap<>();
          neighbors = new ConcurrentHashMap<>();
  
          initializeArrowsBackward();
  
          while (!sortedArrows.isEmpty()) {
              Arrow arrow = sortedArrows.first();
              sortedArrows.remove(arrow);
  
              Node x = arrow.getA();
              Node y = arrow.getB();
  
              if (!arrow.getNaYX().equals(getNaYX(x, y))) {
                  continue;
              }
  
              if (!graph.isAdjacentTo(x, y)) continue;
  
              Edge edge = graph.getEdge(x, y);
              if (edge.pointsTowards(x)) continue;
  
              HashSet<Node> diff = new HashSet<>(arrow.getNaYX());
              diff.removeAll(arrow.getHOrT());
  
              if (!validDelete(x, y, arrow.getHOrT(), arrow.getNaYX())) continue;
  
              Set<Node> H = arrow.getHOrT();
              double bump = arrow.getBump();
  
              boolean deleted = delete(x, y, H, bump, arrow.getNaYX());
              if (!deleted) continue;
  
              totalScore += bump;
  
              clearArrow(x, y);
  
              Set<Node> visited = reapplyOrientation(x, y, H);
  
              Set<Node> toProcess = new HashSet<>();
  
              for (Node node : visited) {
                  final Set<Node> neighbors1 = getNeighbors(node);
                  final Set<Node> storedNeighbors = this.neighbors.get(node);
  
                  if (!(neighbors1.equals(storedNeighbors))) {
                      toProcess.add(node);
                  }
              }
  
              toProcess.add(x);
              toProcess.add(y);
              toProcess.addAll(getCommonAdjacents(x, y));
  
              storeGraph();
              reevaluateBackward(toProcess);
          }
  
          meekOrientRestricted(getVariables(), getKnowledge());
      }
  
      private Set<Node> getCommonAdjacents(Node x, Node y) {
          Set<Node> commonChildren = new HashSet<>(graph.getAdjacentNodes(x));
          commonChildren.retainAll(graph.getAdjacentNodes(y));
          return commonChildren;
      }
  
      private Set<Node> reapplyOrientation(Node x, Node y, Set<Node> newArrows) {
          Set<Node> toProcess = new HashSet<>();
          toProcess.add(x);
          toProcess.add(y);
  
          if (newArrows != null) {
              toProcess.addAll(newArrows);
          }
  
          return meekOrientRestricted(new ArrayList<>(toProcess), getKnowledge());
      }
  
      // Returns true if knowledge is not empty.
      private boolean existsKnowledge() {
          return !knowledge.isEmpty();
      }
  
  
      // Initiaizes the sorted arrows lists for the backward search.
      private void initializeArrowsBackward() {
          for (Edge edge : graph.getEdges()) {
              Node x = edge.getNode1();
              Node y = edge.getNode2();
  
              if (existsKnowledge()) {
                  if (!getKnowledge().noEdgeRequired(x.getName(), y.getName())) {
                      continue;
                  }
              }
  
              clearArrow(x, y);
              clearArrow(y, x);
  
              if (edge.pointsTowards(y)) {
                  calculateArrowsBackward(x, y);
              } else if (edge.pointsTowards(x)) {
                  calculateArrowsBackward(y, x);
              } else {
                  calculateArrowsBackward(x, y);
                  calculateArrowsBackward(y, x);
              }
  
              this.neighbors.put(x, getNeighbors(x));
              this.neighbors.put(y, getNeighbors(y));
          }
      }
  
      // Calcuates new arrows based on changes in the graph for the forward search.
      private void reevaluateForward(final Set<Node> nodes, final Arrow arrow) {
          class AdjTask extends RecursiveTask<Boolean> {
              private final List<Node> nodes;
              private int from;
              private int to;
              private int chunk;
  
              public AdjTask(int chunk, List<Node> nodes, int from, int to) {
                  this.nodes = nodes;
                  this.from = from;
                  this.to = to;
                  this.chunk = chunk;
              }
  
              @Override
              protected Boolean compute() {
                  if (to - from <= chunk) {
                      for (int _w = from; _w < to; _w++) {
                          Node x = nodes.get(_w);
  
                          List<Node> adj;
  
                          if (mode == Mode.heuristicSpeedup) {
                              adj = effectEdgesGraph.getAdjacentNodes(x);
                          } else if (mode == Mode.coverNoncolliders) {
                              Set<Node> g = new HashSet<>();
  
                              for (Node n : graph.getAdjacentNodes(x)) {
                                  for (Node m : graph.getAdjacentNodes(n)) {
                                      if (graph.isAdjacentTo(x, m)) {
                                          continue;
                                      }
  
                                      if (graph.isDefCollider(m, n, x)) {
                                          continue;
                                      }
  
                                      g.add(m);
                                  }
                              }
  
                              adj = new ArrayList<>(g);
                          } else if (mode == Mode.allowUnfaithfulness) {
                              HashSet<Node> D = new HashSet<>();
                              D.addAll(GraphUtils.getDconnectedVars(x, new ArrayList<Node>(), graph));
                              D.remove(x);
                              adj = new ArrayList<>(D);
                          } else {
                              throw new IllegalStateException();
                          }
  
                          for (Node w : adj) {
                              if (adjacencies != null && !(adjacencies.isAdjacentTo(w, x))) {
                                  continue;
                              }
  
                              if (w == x) continue;
  
                              if (!graph.isAdjacentTo(w, x)) {
                                  clearArrow(w, x);
                                  calculateArrowsForward(w, x);
                              }
                          }
                      }
  
                      return true;
                  } else {
                      int mid = (to - from) / 2;
  
                      List<AdjTask> tasks = new ArrayList<>();
  
                      tasks.add(new AdjTask(chunk, nodes, from, from + mid));
                      tasks.add(new AdjTask(chunk, nodes, from + mid, to));
  
                      invokeAll(tasks);
  
                      return true;
                  }
              }
          }
  
          final AdjTask task = new AdjTask(getMinChunk(nodes.size()), new ArrayList<>(nodes), 0, nodes.size());
          pool.invoke(task);
      }
  
      // Calculates the new arrows for an a->b edge.
      private void calculateArrowsForward(Node a, Node b) {
          if (mode == Mode.heuristicSpeedup && !effectEdgesGraph.isAdjacentTo(a, b)) return;
          if (adjacencies != null && !adjacencies.isAdjacentTo(a, b)) return;
          this.neighbors.put(b, getNeighbors(b));
  
          if (existsKnowledge()) {
              if (getKnowledge().isForbidden(a.getName(), b.getName())) {
                  return;
              }
          }
  
          Set<Node> naYX = getNaYX(a, b);
          if (!isClique(naYX)) return;
  
          List<Node> TNeighbors = getTNeighbors(a, b);
          int _maxIndegree = maxIndegree == -1 ? 1000 : maxIndegree;
  
          final int _max = Math.min(TNeighbors.size(), _maxIndegree - graph.getIndegree(b));
  
          Set<Set<Node>> previousCliques = new HashSet<>();
          previousCliques.add(new HashSet<Node>());
          Set<Set<Node>> newCliques = new HashSet<>();
  
          FOR:
          for (int i = 0; i <= _max; i++) {
              final ChoiceGenerator gen = new ChoiceGenerator(TNeighbors.size(), i);
              int[] choice;
  
              while ((choice = gen.next()) != null) {
                  Set<Node> T = GraphUtils.asSet(choice, TNeighbors);
  
                  Set<Node> union = new HashSet<>(naYX);
                  union.addAll(T);
  
                  boolean foundAPreviousClique = false;
  
                  for (Set<Node> clique : previousCliques) {
                      if (union.containsAll(clique)) {
                          foundAPreviousClique = true;
                          break;
                      }
                  }
  
                  if (!foundAPreviousClique) {
                      break FOR;
                  }
  
                  if (!isClique(union)) continue;
                  newCliques.add(union);
  
                  double bump = insertEval(a, b, T, naYX, hashIndices);
  
                  if (bump > 0.0) {
                      addArrow(a, b, naYX, T, bump);
                  }
  
  //                if (mode == Mode.heuristicSpeedup && union.isEmpty() && score.isEffectEdge(bump) &&
  //                        !effectEdgesGraph.isAdjacentTo(a, b) && graph.getParents(b).isEmpty()) {
  //                    effectEdgesGraph.addUndirectedEdge(a, b);
  //                }
              }
  
              previousCliques = newCliques;
              newCliques = new HashSet<>();
          }
      }
  
      private void addArrow(Node a, Node b, Set<Node> naYX, Set<Node> hOrT, double bump) {
          Arrow arrow = new Arrow(bump, a, b, hOrT, naYX, arrowIndex++);
          sortedArrows.add(arrow);
          addLookupArrow(a, b, arrow);
      }
  
      // Reevaluates arrows after removing an edge from the graph.
      private void reevaluateBackward(Set<Node> toProcess) {
          class BackwardTask extends RecursiveTask<Boolean> {
              private final Node r;
              private List<Node> adj;
              private Map<Node, Integer> hashIndices;
              private int chunk;
              private int from;
              private int to;
  
              public BackwardTask(Node r, List<Node> adj, int chunk, int from, int to,
                                  Map<Node, Integer> hashIndices) {
                  this.adj = adj;
                  this.hashIndices = hashIndices;
                  this.chunk = chunk;
                  this.from = from;
                  this.to = to;
                  this.r = r;
              }
  
              @Override
              protected Boolean compute() {
                  if (to - from <= chunk) {
                      for (int _w = from; _w < to; _w++) {
                          final Node w = adj.get(_w);
                          Edge e = graph.getEdge(w, r);
  
                          if (e != null) {
                              if (e.pointsTowards(r)) {
                                  clearArrow(w, r);
                                  clearArrow(r, w);
  
                                  calculateArrowsBackward(w, r);
                              } else if (Edges.isUndirectedEdge(graph.getEdge(w, r))) {
                                  clearArrow(w, r);
                                  clearArrow(r, w);
  
                                  calculateArrowsBackward(w, r);
                                  calculateArrowsBackward(r, w);
                              }
                          }
                      }
  
                      return true;
                  } else {
                      int mid = (to - from) / 2;
  
                      List<BackwardTask> tasks = new ArrayList<>();
  
                      tasks.add(new BackwardTask(r, adj, chunk, from, from + mid, hashIndices));
                      tasks.add(new BackwardTask(r, adj, chunk, from + mid, to, hashIndices));
  
                      invokeAll(tasks);
  
                      return true;
                  }
              }
          }
  
          for (Node r : toProcess) {
              this.neighbors.put(r, getNeighbors(r));
              List<Node> adjacentNodes = graph.getAdjacentNodes(r);
              pool.invoke(new BackwardTask(r, adjacentNodes, getMinChunk(adjacentNodes.size()), 0,
                      adjacentNodes.size(), hashIndices));
          }
      }
  
      // Calculates the arrows for the removal in the backward direction.
      private void calculateArrowsBackward(Node a, Node b) {
          if (existsKnowledge()) {
              if (!getKnowledge().noEdgeRequired(a.getName(), b.getName())) {
                  return;
              }
          }
  
          Set<Node> naYX = getNaYX(a, b);
  
          List<Node> _naYX = new ArrayList<>(naYX);
  
          final int _depth = _naYX.size();
  
          for (int i = 0; i <= _depth; i++) {
              final ChoiceGenerator gen = new ChoiceGenerator(_naYX.size(), i);
              int[] choice;
  
              while ((choice = gen.next()) != null) {
                  Set<Node> diff = GraphUtils.asSet(choice, _naYX);
  
                  Set<Node> h = new HashSet<>(_naYX);
                  h.removeAll(diff);
  
                  if (existsKnowledge()) {
                      if (!validSetByKnowledge(b, h)) {
                          continue;
                      }
                  }
  
                  double bump = deleteEval(a, b, diff, naYX, hashIndices);
  
                  if (bump > 0.0) {
                      addArrow(a, b, naYX, h, bump);
                  }
              }
          }
      }
  
      public double getModelScore() {
          return modelScore;
      }
  
      // Basic data structure for an arrow a->b considered for additiom or removal from the graph, together with
      // associated sets needed to make this determination. For both forward and backward direction, NaYX is needed.
      // For the forward direction, T neighbors are needed; for the backward direction, H neighbors are needed.
      // See Chickering (2002). The totalScore difference resulting from added in the edge (hypothetically) is recorded
      // as the "bump".
      private static class Arrow implements Comparable<Arrow> {
          private double bump;
          private Node a;
          private Node b;
          private Set<Node> hOrT;
          private Set<Node> naYX;
          private int index = 0;
  
          public Arrow(double bump, Node a, Node b, Set<Node> hOrT, Set<Node> naYX, int index) {
              this.bump = bump;
              this.a = a;
              this.b = b;
              this.hOrT = hOrT;
              this.naYX = naYX;
              this.index = index;
          }
  
          public double getBump() {
              return bump;
          }
  
          public Node getA() {
              return a;
          }
  
          public Node getB() {
              return b;
          }
  
          public Set<Node> getHOrT() {
              return hOrT;
          }
  
          public Set<Node> getNaYX() {
              return naYX;
          }
  
          // Sorting by bump, high to low. The problem is the SortedSet contains won't add a new element if it compares
          // to zero with an existing element, so for the cases where the comparison is to zero (i.e. have the same
          // bump, we need to determine as quickly as possible a determinate ordering (fixed) ordering for two variables.
          // The fastest way to do this is using a hash code, though it's still possible for two Arrows to have the
          // same hash code but not be equal. If we're paranoid, in this case we calculate a determinate comparison
          // not equal to zero by keeping a list. This last part is commened out by default.
          public int compareTo(Arrow arrow) {
              if (arrow == null) throw new NullPointerException();
  
              final int compare = Double.compare(arrow.getBump(), getBump());
  
              if (compare == 0) {
                  return Integer.compare(getIndex(), arrow.getIndex());
              }
  
              return compare;
          }
  
          public String toString() {
              return "Arrow<" + a + "->" + b + " bump = " + bump + " t/h = " + hOrT + " naYX = " + naYX + ">";
          }
  
          public int getIndex() {
              return index;
          }
      }
  
      // Get all adj that are connected to Y by an undirected edge and not adjacent to X.
      private List<Node> getTNeighbors(Node x, Node y) {
          List<Edge> yEdges = graph.getEdges(y);
          List<Node> tNeighbors = new ArrayList<>();
  
          for (Edge edge : yEdges) {
              if (!Edges.isUndirectedEdge(edge)) {
                  continue;
              }
  
              Node z = edge.getDistalNode(y);
  
              if (graph.isAdjacentTo(z, x)) {
                  continue;
              }
  
              tNeighbors.add(z);
          }
  
          return tNeighbors;
      }
  
      // Get all adj that are connected to Y.
      private Set<Node> getNeighbors(Node y) {
          List<Edge> yEdges = graph.getEdges(y);
          Set<Node> neighbors = new HashSet<>();
  
          for (Edge edge : yEdges) {
              if (!Edges.isUndirectedEdge(edge)) {
                  continue;
              }
  
              Node z = edge.getDistalNode(y);
  
              neighbors.add(z);
          }
  
          return neighbors;
      }
  
      // Evaluate the Insert(X, Y, T) operator (Definition 12 from Chickering, 2002).
      private double insertEval(Node x, Node y, Set<Node> t, Set<Node> naYX,
                                Map<Node, Integer> hashIndices) {
          Set<Node> set = new HashSet<>(naYX);
          set.addAll(t);
          set.addAll(graph.getParents(y));
          return scoreGraphChange(y, set, x, hashIndices);
      }
  
      // Evaluate the Delete(X, Y, T) operator (Definition 12 from Chickering, 2002).
      private double deleteEval(Node x, Node y, Set<Node> diff, Set<Node> naYX,
                                Map<Node, Integer> hashIndices) {
          Set<Node> set = new HashSet<>(diff);
          set.addAll(graph.getParents(y));
          set.remove(x);
          return -scoreGraphChange(y, set, x, hashIndices);
      }
  
      // Do an actual insertion. (Definition 12 from Chickering, 2002).
      private boolean insert(Node x, Node y, Set<Node> T, double bump) {
          if (graph.isAdjacentTo(x, y)) {
              return false; // The initial graph may already have put this edge in the graph.
          }
  
          Edge trueEdge = null;
  
          if (trueGraph != null) {
              Node _x = trueGraph.getNode(x.getName());
              Node _y = trueGraph.getNode(y.getName());
              trueEdge = trueGraph.getEdge(_x, _y);
          }
  
          if (boundGraph != null && !boundGraph.isAdjacentTo(x, y)) return false;
  
          graph.addDirectedEdge(x, y);
  
          if (verbose) {
              String label = trueGraph != null && trueEdge != null ? "*" : "";
              TetradLogger.getInstance().log("insertedEdges", graph.getNumEdges() + ". INSERT " + graph.getEdge(x, y) +
                      " " + T + " " + bump + " " + label);
          }
  
          int numEdges = graph.getNumEdges();
  
  //        if (verbose) {
          if (numEdges % 1000 == 0) out.println("Num edges added: " + numEdges);
  //        }
  
          if (verbose) {
              String label = trueGraph != null && trueEdge != null ? "*" : "";
              out.println(graph.getNumEdges() + ". INSERT " + graph.getEdge(x, y) +
                      " " + T + " " + bump + " " + label
                      + " degree = " + GraphUtils.getDegree(graph)
                      + " indegree = " + GraphUtils.getIndegree(graph));
          }
  
          for (Node _t : T) {
              graph.removeEdge(_t, y);
              if (boundGraph != null && !boundGraph.isAdjacentTo(_t, y)) continue;
  
              graph.addDirectedEdge(_t, y);
  
              if (verbose) {
                  String message = "--- Directing " + graph.getEdge(_t, y);
                  TetradLogger.getInstance().log("directedEdges", message);
                  out.println(message);
              }
          }
  
          return true;
      }
  
      Set<Edge> removedEdges = new HashSet<>();
  
      // Do an actual deletion (Definition 13 from Chickering, 2002).
      private boolean delete(Node x, Node y, Set<Node> H, double bump, Set<Node> naYX) {
          Edge trueEdge = null;
  
          if (trueGraph != null) {
              Node _x = trueGraph.getNode(x.getName());
              Node _y = trueGraph.getNode(y.getName());
              trueEdge = trueGraph.getEdge(_x, _y);
          }
  
          Edge oldxy = graph.getEdge(x, y);
  
          Set<Node> diff = new HashSet<>(naYX);
          diff.removeAll(H);
  
          graph.removeEdge(oldxy);
          removedEdges.add(Edges.undirectedEdge(x, y));
  
  //        if (verbose) {
          int numEdges = graph.getNumEdges();
          if (numEdges % 1000 == 0) out.println("Num edges (backwards) = " + numEdges);
  //        }
  
          if (verbose) {
              String label = trueGraph != null && trueEdge != null ? "*" : "";
              String message = (graph.getNumEdges()) + ". DELETE " + x + "-->" + y +
                      " H = " + H + " NaYX = " + naYX + " diff = " + diff + " (" + bump + ") " + label;
              TetradLogger.getInstance().log("deletedEdges", message);
              out.println(message);
          }
  
          for (Node h : H) {
              if (graph.isParentOf(h, y) || graph.isParentOf(h, x)) continue;
  
              Edge oldyh = graph.getEdge(y, h);
  
              graph.removeEdge(oldyh);
  
              graph.addEdge(Edges.directedEdge(y, h));
  
              if (verbose) {
                  TetradLogger.getInstance().log("directedEdges", "--- Directing " + oldyh + " to " +
                          graph.getEdge(y, h));
                  out.println("--- Directing " + oldyh + " to " + graph.getEdge(y, h));
              }
  
              Edge oldxh = graph.getEdge(x, h);
  
              if (Edges.isUndirectedEdge(oldxh)) {
                  graph.removeEdge(oldxh);
  
                  graph.addEdge(Edges.directedEdge(x, h));
  
                  if (verbose) {
                      TetradLogger.getInstance().log("directedEdges", "--- Directing " + oldxh + " to " +
                              graph.getEdge(x, h));
                      out.println("--- Directing " + oldxh + " to " + graph.getEdge(x, h));
                  }
              }
          }
  
          return true;
      }
  
      // Test if the candidate insertion is a valid operation
      // (Theorem 15 from Chickering, 2002).
      private boolean validInsert(Node x, Node y, Set<Node> T, Set<Node> naYX) {
          boolean violatesKnowledge = false;
  
  