accuracy for CPDAG

notears.py

@@ -99,5 +99,5 @@ def _func(w):
  import igraph as ig
  assert ig.Graph.Weighted_Adjacency(W_est.tolist()).is_dag()
  np.savetxt('W_est.csv', W_est, delimiter=',')
- acc = ut.count_accuracy(W_true, W_est)
+ acc = ut.count_accuracy(B_true, W_est != 0)
  print(acc)

utils.py

@@ -119,13 +119,12 @@ def _simulate_single_equation(X, w):
  return X
 
 
-def count_accuracy(W_true, W, W_und=None):
- """Compute FDR, TPR, and FPR for B, or optionally for CPDAG = B + B_und.
+def count_accuracy(B_true, B_est):
+ """Compute various accuracy metrics for B_est.
 
  Args:
- W_true (np.ndarray): [d, d] ground truth graph
- W (np.ndarray): [d, d] predicted graph
- W_und (np.ndarray): [d, d] predicted undirected edges in CPDAG, asymmetric
+ B_true (np.ndarray): [d, d] ground truth graph, {0, 1}
+ B_est (np.ndarray): [d, d] estimate, {0, 1, -1}, -1 is undirected edge in CPDAG
 
  Returns:
  fdr: (reverse + false positive) / prediction positive
@@ -134,46 +133,35 @@ def count_accuracy(W_true, W, W_und=None):
  shd: undirected extra + undirected missing + reverse
  nnz: prediction positive
  """
- d = W_true.shape[0]
- # convert to binary adjacency matrix
- B_true = (W_true != 0)
- B = (W != 0)
- B_und = None if W_und is None else (W_und != 0)
+ if ((B_est == -1) == (B_est.T == -1)).any():
+ raise ValueError('undirected edge should only appear once')
+ d = B_true.shape[0]
  # linear index of nonzeros
- pred_und = None
- if B_und is not None:
- pred_und = np.flatnonzero(B_und)
- pred = np.flatnonzero(B)
+ pred_und = np.flatnonzero(B_est == -1)
+ pred = np.flatnonzero(B_est == 1)
  cond = np.flatnonzero(B_true)
  cond_reversed = np.flatnonzero(B_true.T)
  cond_skeleton = np.concatenate([cond, cond_reversed])
  # true pos
  true_pos = np.intersect1d(pred, cond, assume_unique=True)
- if B_und is not None:
- # treat undirected edge favorably
- true_pos_und = np.intersect1d(pred_und, cond_skeleton, assume_unique=True)
- true_pos = np.concatenate([true_pos, true_pos_und])
+ # treat undirected edge favorably
+ true_pos_und = np.intersect1d(pred_und, cond_skeleton, assume_unique=True)
+ true_pos = np.concatenate([true_pos, true_pos_und])
  # false pos
  false_pos = np.setdiff1d(pred, cond_skeleton, assume_unique=True)
- if B_und is not None:
- false_pos_und = np.setdiff1d(pred_und, cond_skeleton, assume_unique=True)
- false_pos = np.concatenate([false_pos, false_pos_und])
+ false_pos_und = np.setdiff1d(pred_und, cond_skeleton, assume_unique=True)
+ false_pos = np.concatenate([false_pos, false_pos_und])
  # reverse
  extra = np.setdiff1d(pred, cond, assume_unique=True)
  reverse = np.intersect1d(extra, cond_reversed, assume_unique=True)
  # compute ratio
- pred_size = len(pred)
- if B_und is not None:
- pred_size += len(pred_und)
+ pred_size = len(pred) + len(pred_und)
  cond_neg_size = 0.5 * d * (d - 1) - len(cond)
  fdr = float(len(reverse) + len(false_pos)) / max(pred_size, 1)
  tpr = float(len(true_pos)) / max(len(cond), 1)
  fpr = float(len(reverse) + len(false_pos)) / max(cond_neg_size, 1)
  # structural hamming distance
- B_lower = np.tril(B + B.T)
- if B_und is not None:
- B_lower += np.tril(B_und + B_und.T)
- pred_lower = np.flatnonzero(B_lower)
+ pred_lower = np.flatnonzero(np.tril(B_est + B_est.T))
  cond_lower = np.flatnonzero(np.tril(B_true + B_true.T))
  extra_lower = np.setdiff1d(pred_lower, cond_lower, assume_unique=True)
  missing_lower = np.setdiff1d(cond_lower, pred_lower, assume_unique=True)