Skip to content

Revise misplaced columns in dynotears estimation, and fix a minor bug in dataset generation. #119

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ThinkNaive wants to merge 9 commits into
base: develop
Choose a base branch
from
Open

Revise misplaced columns in dynotears estimation, and fix a minor bug in dataset generation. #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
1ed7976
Revise misplaced columns in dynotears estimation, and fix a minor bug…
ThinkNaive Jun 19, 2021
b9b5e4d
Merge branch 'develop' into develop
oentaryorj Nov 3, 2021
1c56978
Replace w_threshold with tau_w
oentaryorj Nov 3, 2021
eb61fff
Update test_dynotears.py
oentaryorj Nov 3, 2021
32372c1
Update .gitignore
oentaryorj Nov 3, 2021
e734905
Merge branch 'develop' into develop
oentaryorj Aug 26, 2022
d128d82
Merge branch 'develop' into develop
oentaryorj Aug 26, 2022
0e39493
Merge branch 'develop' into develop
oentaryorj Aug 26, 2022
36727b5
Merge branch 'develop' into develop
oentaryorj Aug 29, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions causalnex/structure/data_generators/wrappers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -533,10 +533,12 @@ def generate_dataframe_dynamic( # pylint: disable=R0914
ValueError: if sem_type isn't linear-gauss/linear_exp/linear-gumbel
"""
s_types = ("linear-gauss", "linear-exp", "linear-gumbel")

if sem_type not in s_types:
raise ValueError(f"unknown sem type {sem_type}. Available types are: {s_types}")
intra_nodes = sorted(el for el in g.nodes if "_lag0" in el)
inter_nodes = sorted(el for el in g.nodes if "_lag0" not in el)

intra_nodes = sorted([el for el in g.nodes if "_lag0" in el], key=lambda t: t.split('_lag')[1])

@mzjp2 mzjp2 Jun 21, 2021

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorted accepts a generator expression, doesn't need to (waste memory!) coercing it into a list first :)

@ThinkNaive ThinkNaive Jun 21, 2021

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorted(el for el in g.nodes if "_lag0" not in el) will lead to misplaced columns, e.g. [0_lag1, 0_lag2, 0_lag3, 1_lag1, 1_lag2, 1_lag3, 2_lag1, 2_lag2, 2_lag3] (node first, then lag, case1). But the dynotears model accepts columns like [0_lag1, 1_lag1, 2_lag1, 0_lag2, 1_lag2, 2_lag2, 0_lag3, 1_lag3, 2_lag3] (lag first, then node, case2). The default sorted expression takes a case1 style. I've test the revised code using the code below.

__test__.py

# __test__.py
from causalnex.structure.transformers import DynamicDataTransformer
import os
import random
from matplotlib import ticker
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from causalnex.structure.data_generators import gen_stationary_dyn_net_and_df
from netgraph import InteractiveGraph
from causalnex.structure.dynotears import from_pandas_dynamic
import scipy.linalg as slin

# metrics
def count_accuracy(B_true, B_est) -> tuple:
    B_true = B_true != 0
    B_est = B_est != 0
    d = min(B_est.shape)

    # fdr
    fp = np.sum(~B_true & B_est)
    pp = np.sum(B_est)
    fdr = fp / pp

    # tpr
    tp = np.sum(B_true & B_est)
    tt = np.sum(B_true)
    tpr = tp / tt

    # fpr
    tf = d * (d - 1) / 2 - np.sum(B_true) + \
        B_est.shape[0] * B_est.shape[1] - d * d
    fpr = fp / tf

    # shd
    shd = np.sum(B_true != B_est)

    # nnz
    nnz = pp

    return fdr, tpr, fpr, shd, nnz

# reproduce
def set_seed(seed):
    """
    Referred from:
    - https://stackoverflow.com/questions/38469632/tensorflow-non-repeatable-results
    """
    # Reproducibility
    random.seed(seed)
    np.random.seed(seed)
    try:
        os.environ["PYTHONHASHSEED"] = str(seed)
    except:
        pass

# plot weighted directed graph
def plot_graph(g, node_size=8, font_size=8):
    # define layout
    pos = {}
    for node in g.nodes():
        # get node vertex index and slice index
        index = node.split('_lag')
        vertex = int(index[0]) + 1
        slice = int(index[1])

        # generate position with each slice as round shape
        x = np.random.uniform(low=-0.0, high=0.0) + 0.3 * \
            np.cos(2 * np.pi * vertex / d) + 1.0 * (p - slice)
        y = np.random.uniform(low=-0.1, high=0.1) + \
            np.sin(2 * np.pi * vertex / d)
        pos[node] = (x, y)

    g = nx.DiGraph(g)
    fig, ax = plt.subplots()
    plot_instance = InteractiveGraph(
        g,
        node_size=node_size,
        node_labels=True,
        node_label_fontdict=dict(size=font_size),
        node_layout=pos,
        arrows=True,
        ax=ax
    )

# transform graph to matrix
def to_matrix(g):
    a = nx.to_numpy_array(g)
    # a /= np.abs(a).max()

    # permute array by order from max lag to instant.
    sorted_nodes = sorted(g.nodes(), key=lambda k: d * int(k[-1]) + int(k[0]))
    sorted_index = [list(g.nodes()).index(i) for i in sorted_nodes]
    p = np.zeros(a.shape)
    for y, x in enumerate(sorted_index):
        p[y, x] = 1.0

    a = p.dot(a).dot(p.T)
    return a[:, :d]

# plot matrix
def plot_matrix(gt, ge=None, names=None, rng=2.0):
    # transform graph to matrix and compute layout
    bt = to_matrix(gt)
    d = bt.shape[1]
    n_col = bt.shape[0] // d
    if ge:
        be = to_matrix(ge)
        n_row = 2
        b = [bt, be]
    else:
        n_row = 1
        b = [bt]

    # plot matrix
    fig, ax = plt.subplots(n_row, n_col, figsize=(16, 6))
    ax = ax.flatten()
    for row in range(n_row):
        for col in range(n_col):
            # split matrix for intra and inter-p
            mat = b[row][d * col: d * (col + 1), :]
            # plot matrix
            im = ax[row * n_col + col].imshow(mat, cmap="seismic",
                                            interpolation="none", vmin=-rng, vmax=rng)
            # add value labels
            for i in range(d):
                for j in range(d):
                    color = 'white' if abs(mat[i, j]) > 0.3 * rng else 'black'
                    value = '{:g}'.format(round(mat[i, j], 3))
                    ax[row * n_col + col].text(j, i, value,
                                            ha="center", va="center", color=color, fontsize=9)
            # hide axis
            ax[row * n_col + col].set_xticks([])
            ax[row * n_col + col].set_yticks([])
            # show graph type labels
            if col == 0:
                ax[row * n_col + col].set_ylabel(names[row], fontsize=10)
            # show instant and lagged labels
            if row == 0:
                if col == 0:
                    label = '$W$ (Intra-slice)'
                else:
                    label = '$A_{%d}$ (Inter-slice)' % col
                ax[row * n_col + col].set_xlabel(label, fontsize=10)
                ax[row * n_col + col].xaxis.set_label_position('top')

    # add colorbar
    cb = fig.colorbar(im, ax=ax, shrink=0.7)
    tick_locator = ticker.MaxNLocator(nbins=5)
    cb.locator = tick_locator
    cb.set_ticks([t - rng for t in range(int(rng) * 2 + 1)])
    cb.update_ticks()

    return fig

# compute loss
def compute_loss(dataset, g, p):
    X, Xlags = DynamicDataTransformer(
        p=p).fit_transform(dataset, return_df=False)
    n, d_vars = X.shape

    wa = to_matrix(g)
    w_mat, a_mat = wa[: d_vars, :], wa[d_vars:, :]

    loss = (
        0.5
        / n
        * np.square(
            np.linalg.norm(
                X.dot(np.eye(d_vars, d_vars) - w_mat)
                - Xlags.dot(a_mat), "fro"
            )
        )
    )

    h = np.trace(slin.expm(w_mat * w_mat)) - d_vars

    return loss, h

# main
set_seed(0)
d = 5
p = 3

# generate dataset
g, df, intra, inter = gen_stationary_dyn_net_and_df(
    n_samples=500,
    num_nodes=5,
    p=3,
    w_min_intra=0.5,
    w_max_intra=2.0,
    w_min_inter=0.3,
    w_max_inter=0.5,
    degree_intra=4,
    degree_inter=1,
    w_decay=1.3
)

# convert dataset for learning
dataset = df.iloc[:, :d]
dataset.columns = [t.split('_lag')[0] for t in dataset.columns]

# estimate graph
gg = from_pandas_dynamic(
    time_series=[dataset], p=3, tau_w=0.01, tau_a=0.01, lambda_w=0.05, lambda_a=0.05)

# compute metrics
fdr, tpr, fpr, shd, nnz = count_accuracy(
    to_matrix(g), to_matrix(gg))
print("FDR:{:.2%} TPR:{:.2%} FPR:{:.2%} SHD:{:2d} NNZ:{:2d}".format(
    fdr, tpr, fpr, shd, nnz
))

# compute loss for both true graph and estimated graph
print('TrueGraph: loss={} h={}'.format(*compute_loss(dataset, g, p)))
print('Estimated: loss={} h={}'.format(*compute_loss(dataset, gg, p)))

# plot graph and matrix
# plot_graph(g)
# plot_graph(gg)
fig = plot_matrix(g, gg, names=('True', 'Estimated'))

# save image
imgfile = "graph.svg"
fig.savefig(imgfile, transparent=True)

plt.show(block=True)

inter_nodes = sorted([el for el in g.nodes if "_lag0" not in el], key=lambda t: t.split('_lag')[1])
w_mat = nx.to_numpy_array(g, nodelist=intra_nodes)
a_mat = nx.to_numpy_array(g, nodelist=intra_nodes + inter_nodes)[
len(intra_nodes) :, : len(intra_nodes)
Expand Down Expand Up @@ -612,7 +614,7 @@ def gen_stationary_dyn_net_and_df( # pylint: disable=R0913, R0914
- full: constructs a fully-connected graph - degree has no effect
graph_type_inter:
- erdos-renyi: constructs a graph such that the probability of any given edge is degree / (num_nodes - 1)
- full: connect all past nodes to all present nodesw_min_intra:
- full: connect all past nodes to all present nodes
w_min_intra: minimum weight on intra-slice adjacency matrix
w_max_intra: maximum weight on intra-slice adjacency matrix
w_min_inter: minimum weight on inter-slice adjacency matrix
Expand Down
13 changes: 8 additions & 5 deletions causalnex/structure/dynotears.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,8 @@ def from_pandas_dynamic( # pylint: disable=too-many-arguments
lambda_a: float = 0.1,
max_iter: int = 100,
h_tol: float = 1e-8,
w_threshold: float = 0.0,
tau_w: float = 0.0,
tau_a: float = 0.0,
tabu_edges: List[Tuple[int, int, int]] = None,
tabu_parent_nodes: List[int] = None,
tabu_child_nodes: List[int] = None,
Expand Down Expand Up @@ -116,7 +117,8 @@ def from_pandas_dynamic( # pylint: disable=too-many-arguments
lambda_a,
max_iter,
h_tol,
w_threshold,
tau_w,
tau_a,
tabu_edges,
tabu_parent_nodes,
tabu_child_nodes,
Expand Down Expand Up @@ -162,7 +164,8 @@ def from_numpy_dynamic( # pylint: disable=too-many-arguments
lambda_a: float = 0.1,
max_iter: int = 100,
h_tol: float = 1e-8,
w_threshold: float = 0.0,
tau_w: float = 0.0,
tau_a: float = 0.0,
tabu_edges: List[Tuple[int, int, int]] = None,
tabu_parent_nodes: List[int] = None,
tabu_child_nodes: List[int] = None,
Expand Down Expand Up @@ -246,8 +249,8 @@ def from_numpy_dynamic( # pylint: disable=too-many-arguments
X, Xlags, bnds, lambda_w, lambda_a, max_iter, h_tol
)

w_est[np.abs(w_est) < w_threshold] = 0
a_est[np.abs(a_est) < w_threshold] = 0
w_est[np.abs(w_est) < tau_w] = 0
a_est[np.abs(a_est) < tau_a] = 0
sm = _matrices_to_structure_model(w_est, a_est)
return sm

Expand Down
44 changes: 22 additions & 22 deletions tests/structure/test_dynotears.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ def test_expected_structure_learned_p1(self, data_dynotears_p1):
"""

sm = from_numpy_dynamic(
data_dynotears_p1["X"], data_dynotears_p1["Y"], w_threshold=0.2
data_dynotears_p1["X"], data_dynotears_p1["Y"], tau_w=0.2
)
w_edges = [
(f"{i}_lag0", f"{j}_lag0")
Expand Down Expand Up @@ -145,7 +145,7 @@ def test_expected_structure_learned_p2(self, data_dynotears_p2):
"""

sm = from_numpy_dynamic(
data_dynotears_p2["X"], data_dynotears_p2["Y"], w_threshold=0.25
data_dynotears_p2["X"], data_dynotears_p2["Y"], tau_w=0.25
)
w_edges = [
(f"{i}_lag0", f"{j}_lag0")
Expand Down Expand Up @@ -250,7 +250,7 @@ def test_all_columns_in_structure(self, data_dynotears_p2):
def test_isolated_nodes_exist(self, data_dynotears_p2):
"""Isolated nodes should still be in the learned structure"""
sm = from_numpy_dynamic(
data_dynotears_p2["X"], data_dynotears_p2["Y"], w_threshold=1
data_dynotears_p2["X"], data_dynotears_p2["Y"], tau_w=1
)
assert len(sm.edges) == 2
assert len(sm.nodes) == 15
Expand All @@ -273,7 +273,7 @@ def test_certain_relationships_get_near_certain_weight(self):
[[np.sqrt(el), np.sqrt(el)] for el in np.random.choice(100, size=500)],
columns=["a", "b"],
)
sm = from_numpy_dynamic(data.values[1:], data.values[:-1], w_threshold=0.1)
sm = from_numpy_dynamic(data.values[1:], data.values[:-1], tau_w=0.1)
edge = (
sm.get_edge_data("1_lag0", "0_lag0") or sm.get_edge_data("0_lag0", "1_lag0")
)["weight"]
Expand All @@ -287,7 +287,7 @@ def test_inverse_relationships_get_negative_weight(self):
data = pd.DataFrame(
[[el, -el] for el in np.random.choice(100, size=500)], columns=["a", "b"]
)
sm = from_numpy_dynamic(data.values[1:], data.values[:-1], w_threshold=0.1)
sm = from_numpy_dynamic(data.values[1:], data.values[:-1], tau_w=0.1)
edge = (
sm.get_edge_data("1_lag0", "0_lag0") or sm.get_edge_data("0_lag0", "1_lag0")
)["weight"]
Expand All @@ -299,20 +299,20 @@ def test_no_cycles(self, data_dynotears_p2):
"""

sm = from_numpy_dynamic(
data_dynotears_p2["X"], data_dynotears_p2["Y"], w_threshold=0.05
data_dynotears_p2["X"], data_dynotears_p2["Y"], tau_w=0.05
)
assert nx.algorithms.is_directed_acyclic_graph(sm)

def test_tabu_edges_on_non_existing_edges_do_nothing(self, data_dynotears_p2):
"""If tabu edges do not exist in the original unconstrained network then nothing changes"""
sm = from_numpy_dynamic(
data_dynotears_p2["X"], data_dynotears_p2["Y"], w_threshold=0.2
data_dynotears_p2["X"], data_dynotears_p2["Y"], tau_w=0.2
)

sm_2 = from_numpy_dynamic(
data_dynotears_p2["X"],
data_dynotears_p2["Y"],
w_threshold=0.2,
tau_w=0.2,
tabu_edges=[(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3)],
)
assert set(sm_2.edges) == set(sm.edges)
Expand Down Expand Up @@ -391,7 +391,7 @@ def test_expected_structure_learned_p1(self, data_dynotears_p1):
sm = from_pandas_dynamic(
df,
p=1,
w_threshold=0.2,
tau_w=0.2,
)
map_ = dict(zip(range(5), ["a", "b", "c", "d", "e"]))
w_edges = [
Expand Down Expand Up @@ -430,7 +430,7 @@ def test_expected_structure_learned_p2(self, data_dynotears_p2):
sm = from_pandas_dynamic(
df,
p=2,
w_threshold=0.25,
tau_w=0.25,
)
map_ = dict(zip(range(5), ["a", "b", "c", "d", "e"]))
w_edges = [
Expand Down Expand Up @@ -532,7 +532,7 @@ def test_all_columns_in_structure(self, data_dynotears_p2):
sm = from_pandas_dynamic(
pd.DataFrame(data_dynotears_p2["X"], columns=["a", "b", "c", "d", "e"]),
p=2,
w_threshold=0.4,
tau_w=0.4,
)
assert sorted(sm.nodes) == [
f"{var}_lag{l_val}"
Expand All @@ -547,7 +547,7 @@ def test_isolated_nodes_exist(self, data_dynotears_p2):
df.loc[-2, :] = data_dynotears_p2["Y"][0, 5:10]
df = df.sort_index()

sm = from_pandas_dynamic(df, p=2, w_threshold=1)
sm = from_pandas_dynamic(df, p=2, tau_w=1)
assert len(sm.edges) == 2
assert len(sm.nodes) == 15

Expand All @@ -572,7 +572,7 @@ def test_certain_relationships_get_near_certain_weight(self):
[[np.sqrt(el), np.sqrt(el)] for el in np.random.choice(100, size=500)],
columns=["a", "b"],
)
sm = from_pandas_dynamic(data, p=1, w_threshold=0.1)
sm = from_pandas_dynamic(data, p=1, tau_w=0.1)
edge = (
sm.get_edge_data("b_lag0", "a_lag0") or sm.get_edge_data("a_lag0", "b_lag0")
)["weight"]
Expand All @@ -586,7 +586,7 @@ def test_inverse_relationships_get_negative_weight(self):
data = pd.DataFrame(
[[el, -el] for el in np.random.choice(100, size=500)], columns=["a", "b"]
)
sm = from_pandas_dynamic(data, p=1, w_threshold=0.1)
sm = from_pandas_dynamic(data, p=1, tau_w=0.1)
edge = (
sm.get_edge_data("b_lag0", "a_lag0") or sm.get_edge_data("a_lag0", "b_lag0")
)["weight"]
Expand All @@ -599,7 +599,7 @@ def test_no_cycles(self, data_dynotears_p2):
sm = from_pandas_dynamic(
pd.DataFrame(data_dynotears_p2["X"], columns=["a", "b", "c", "d", "e"]),
p=2,
w_threshold=0.05,
tau_w=0.05,
)
assert nx.algorithms.is_directed_acyclic_graph(sm)

Expand All @@ -613,12 +613,12 @@ def test_tabu_edges_on_non_existing_edges_do_nothing(self, data_dynotears_p2):
sm = from_pandas_dynamic(
df,
p=2,
w_threshold=0.2,
tau_w=0.2,
)
sm_2 = from_pandas_dynamic(
df,
p=2,
w_threshold=0.2,
tau_w=0.2,
tabu_edges=[(0, "a", "a"), (0, "a", "b"), (0, "a", "c"), (0, "a", "d")],
)
assert set(sm_2.edges) == set(sm.edges)
Expand All @@ -636,9 +636,9 @@ def test_list_of_dfs_as_input(self, data_dynotears_p2):
df_ = df.copy()
df_.index = range(100, 152)
df = pd.concat([df, df_])
sm = from_pandas_dynamic(df, p=2, w_threshold=0.05)
sm_1 = from_pandas_dynamic([df], p=2, w_threshold=0.05)
sm_2 = from_pandas_dynamic([df, df], p=2, w_threshold=0.05)
sm = from_pandas_dynamic(df, p=2, tau_w=0.05)
sm_1 = from_pandas_dynamic([df], p=2, tau_w=0.05)
sm_2 = from_pandas_dynamic([df, df], p=2, tau_w=0.05)

assert list(sm_2.edges) == list(sm_1.edges)
assert list(sm.edges) == list(sm_1.edges)
Expand All @@ -664,8 +664,8 @@ def test_discondinuity(self):
index=np.arange(200, 300),
)

sm = from_pandas_dynamic(pd.concat([df, df_2], axis=0), p=2, w_threshold=0.05)
sm_1 = from_pandas_dynamic([df, df_2], p=2, w_threshold=0.05)
sm = from_pandas_dynamic(pd.concat([df, df_2], axis=0), p=2, tau_w=0.05)
sm_1 = from_pandas_dynamic([df, df_2], p=2, tau_w=0.05)

assert [(u, v, round(w, 3)) for u, v, w in sm_1.edges(data="weight")] == [
(u, v, round(w, 3)) for u, v, w in sm.edges(data="weight")
Expand Down