Skip to content

Commit

Permalink
Merge pull request #177 from jakobrunge/developer
Browse files Browse the repository at this point in the history 
fixed bug in causal_effects.py
  • Loading branch information
@jakobrunge
jakobrunge authored Feb 24, 2022
2 parents 0c6a023 + 24e671a commit b833174
Show file tree
Hide file tree
Showing 5 changed files with 131 additions and 220 deletions.
11 changes: 8 additions & 3 deletions release_on_pip.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,16 +33,21 @@ pip install --use-feature=2020-resolver --upgrade setuptools wheel twine auditwh
pip install -e .['dev']

# Rebuild the .c files from the .pyc files
python setup.py develop
# python setup.py develop

# Build the distribution
python setup.py sdist bdist_wheel


# Linux:
auditwheel repair --plat manylinux2014_x86_64 ./dist/*linux_x86_64.whl -w ./dist/
# auditwheel repair --plat manylinux2014_x86_64 ./dist/*linux_x86_64.whl -w ./dist/
# auditwheel repair --plat manylinux2014_x86_64 ./dist/*.whl -w ./dist/


# Upload the distro
twine upload dist/*manylinux*.whl dist/*tar.gz
# twine upload dist/*manylinux*.whl dist/*tar.gz
twine upload dist/*.whl dist/*tar.gz


# Windows (EDIT names to upload all things .whl and tat.gz in dist/ !!!):
# twine upload dist/*manylinux*.whl dist/*tar.gz
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ def run(self):
# Run the setup
setup(
name="tigramite",
version="5.0.0.0",
version="5.0.0.1",
packages=["tigramite", "tigramite.independence_tests", "tigramite.toymodels"],
license="GNU General Public License v3.0",
description="Tigramite causal discovery for time series",
Expand Down
84 changes: 60 additions & 24 deletions tigramite/causal_effects.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,10 +122,10 @@ def __init__(self,

# Construct internal graph from input graph depending on graph type
# and hidden variables
(self.graph, self.graph_type,
self.tau_max, self.hidden_variables) = self._construct_graph(
graph=graph, graph_type=graph_type,
hidden_variables=hidden_variables)
# (self.graph, self.graph_type,
# self.tau_max, self.hidden_variables) =
self._construct_graph(graph=graph, graph_type=graph_type,
hidden_variables=hidden_variables)

# print(self.graph.shape)
self._check_graph(self.graph)
Expand Down Expand Up @@ -212,20 +212,21 @@ def _construct_graph(self, graph, graph_type, hidden_variables):


if graph_type in ['dag', 'admg']:
tau_max = 0
if graph.ndim != 2:
raise ValueError("graph_type in ['dag', 'admg'] assumes graph.shape=(N, N).")
# Convert to shape [N, N, 1, 1] with dummy dimension
# to process as tsg_dag or tsg_admg with potential hidden variables
self.graph = np.expand_dims(graph, axis=(2, 3))

# tau_max needed in _get_latent_projection_graph
self.tau_max = 0

if len(hidden_variables) > 0:
graph = self._get_latent_projection_graph() # stationary=False)
graph_type = "tsg_admg"
self.graph = self._get_latent_projection_graph() # stationary=False)
self.graph_type = "tsg_admg"
else:
graph = self.graph
graph_type = 'tsg_' + graph_type
# graph = self.graph
self.graph_type = 'tsg_' + graph_type

elif graph_type in ['tsg_dag', 'tsg_admg']:
if graph.ndim != 4:
Expand All @@ -237,10 +238,10 @@ def _construct_graph(self, graph, graph_type, hidden_variables):
self.tau_max = graph.shape[2] - 1

if len(hidden_variables) > 0:
graph = self._get_latent_projection_graph() #, stationary=False)
graph_type = "tsg_admg"
self.graph = self._get_latent_projection_graph() #, stationary=False)
self.graph_type = "tsg_admg"
else:
graph_type = graph_type
self.graph_type = graph_type

elif graph_type in ['stationary_dag']:
# Currently only stationary_dag without hidden variables is supported
Expand All @@ -260,16 +261,16 @@ def _construct_graph(self, graph, graph_type, hidden_variables):
for varlag in self.X.union(self.Y).union(self.S):
maxlag_XYS = max(maxlag_XYS, abs(varlag[1]))

tau_max = maxlag_XYS + statgraph_tau_max
self.tau_max = maxlag_XYS + statgraph_tau_max

stat_graph = deepcopy(graph)

# Construct tsg_graph
graph = np.zeros((self.N, self.N, tau_max + 1, tau_max + 1), dtype='<U3')
graph = np.zeros((self.N, self.N, self.tau_max + 1, self.tau_max + 1), dtype='<U3')
graph[:] = ""
for (i, j) in itertools.product(range(self.N), range(self.N)):
for jt, tauj in enumerate(range(0, tau_max + 1)):
for it, taui in enumerate(range(tauj, tau_max + 1)):
for jt, tauj in enumerate(range(0, self.tau_max + 1)):
for it, taui in enumerate(range(tauj, self.tau_max + 1)):
tau = abs(taui - tauj)
if tau == 0 and j == i:
continue
Expand All @@ -293,9 +294,11 @@ def _construct_graph(self, graph, graph_type, hidden_variables):
# graph[i, j, taui, tauj] = "<--"
# graph[j, i, tauj, taui] = "-->"

graph_type = 'tsg_dag'
self.graph_type = 'tsg_dag'
self.graph = graph


return (graph, graph_type, tau_max, hidden_variables)
# return (graph, graph_type, self.tau_max, hidden_variables)

# max_lag = self._get_maximum_possible_lag(XYZ=list(X.union(Y).union(S)), graph=graph)

Expand Down Expand Up @@ -1980,7 +1983,7 @@ def fit_wright_effect(self,
coeffs = {}
for medy in [med for med in mediators] + [y for y in self.listY]:
coeffs[medy] = {}
mediator_parents = self._get_all_parents([medy]).intersection(mediators.union(self.X)) - set([medy])
mediator_parents = self._get_all_parents([medy]).intersection(mediators.union(self.X).union(self.Y)) - set([medy])
all_parents = self._get_all_parents([medy]) - set([medy])
for par in mediator_parents:
Sprime = set(all_parents) - set([par, medy])
Expand All @@ -1998,6 +2001,7 @@ def fit_wright_effect(self,
cut_off='max_lag_or_tau_max',
return_data=False)
coeffs[medy][par] = fit_res[medy]['model'].coef_[0]
# print(mediators, par, medy, coeffs[medy][par])

elif method == 'parents':
if 'dag' not in self.graph_type:
Expand Down Expand Up @@ -2030,6 +2034,7 @@ def fit_wright_effect(self,
effect[(x, y)] = 0.
for causal_path in causal_paths[x][y]:
effect_here = 1.
# print(x, y, causal_path)
for index, node in enumerate(causal_path[:-1]):
i, taui = node
j, tauj = causal_path[index + 1]
Expand Down Expand Up @@ -2089,9 +2094,40 @@ def predict_wright_effect(self,
print("No causal path from X to Y exists.")
return np.zeros((len(intervention_data), len(self.Y)))

effect = self.model.get_general_prediction(
intervention_data=intervention_data,
conditions_data=None,
pred_params=pred_params)
intervention_T, lenX = intervention_data.shape

lenY = len(self.Y)

predicted_array = np.zeros((intervention_T, lenY))
pred_dict = {}
for iy, y in enumerate(self.Y):
# Print message
if self.verbosity > 1:
print("\n## Predicting target %s" % str(y))
if pred_params is not None:
for key in list(pred_params):
print("%s = %s" % (key, pred_params[key]))
# Default value for pred_params
if pred_params is None:
pred_params = {}
# Check this is a valid target
if y not in self.model.fit_results:
raise ValueError("y = %s not yet fitted" % str(y))

# Transform the data if needed
a_transform = self.model.fit_results[y]['data_transform']
if a_transform is not None:
intervention_data = a_transform.transform(X=intervention_data)


# Now iterate through interventions (and potentially S)
for index, dox_vals in enumerate(intervention_data):
# Construct XZS-array
intervention_array = dox_vals.reshape(1, lenX)
predictor_array = intervention_array

predicted_vals = self.model.fit_results[y]['model'].predict(
X=predictor_array, **pred_params)
predicted_array[index, iy] = predicted_vals.mean()

return effect
return predicted_array
142 changes: 8 additions & 134 deletions tigramite/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ def get_general_fitted_model(self,
# Initialize the fit results
fit_results = {}
for y in self.Y:

# Construct array of shape (var, time) with first entry being
# a dummy, second is y followed by joint X and Z (ignore the notation in construct_array)
array, xyz = \
Expand All @@ -164,9 +165,6 @@ def get_general_fitted_model(self,
if self.data_transform is not None:
array = self.data_transform.fit_transform(X=array.T).T

# Cache array for use in prediction
self.observation_array = array
self.xyz = xyz
# Fit the model
# Copy and fit the model
a_model = deepcopy(self.model)
Expand All @@ -179,8 +177,11 @@ def get_general_fitted_model(self,
target_array = array[np.where(xyz==1)[0][0], :]

a_model.fit(X=predictor_array, y=target_array)

# Cache the results
fit_results[y] = {}
fit_results[y]['observation_array'] = array
fit_results[y]['xyz'] = xyz
fit_results[y]['model'] = a_model
# Cache the data transform
fit_results[y]['data_transform'] = deepcopy(self.data_transform)
Expand Down Expand Up @@ -254,13 +255,13 @@ def get_general_prediction(self,
conditions_data = a_transform.transform(X=conditions_data)

# Extract observational Z from stored array
z_indices = list(np.where(self.xyz==1)[0][1:])
z_array = self.observation_array[z_indices, :].T
z_indices = list(np.where(self.fit_results[y]['xyz']==1)[0][1:])
z_array = self.fit_results[y]['observation_array'][z_indices, :].T
Tobs = len(z_array)

if self.conditions is not None and conditions_data is not None:
s_indices = list(np.where(self.xyz==2)[0])
s_array = self.observation_array[s_indices, :].T
s_indices = list(np.where(self.fit_results[y]['xyz']==2)[0])
s_array = self.fit_results[y]['observation_array'][s_indices, :].T

# Now iterate through interventions (and potentially S)
for index, dox_vals in enumerate(intervention_data):
Expand Down Expand Up @@ -293,133 +294,6 @@ def get_general_prediction(self,
return predicted_array


# def get_general_prediction(self,
# intervention_data=None,
# conditions_data=None,
# pred_params=None,
# ):
# r"""Predict effect of intervention with fitted model.

# Uses the model.predict() function of the sklearn model.

# Parameters
# ----------
# intervention_data : data object, optional
# Tigramite dataframe object with optional new mask. Only the
# values for X will be extracted.
# conditions_data : data object, optional
# Tigramite dataframe object with optional new mask. Only the
# values for conditions will be extracted.
# pred_params : dict, optional
# Optional parameters passed on to sklearn prediction function.

# Returns
# -------
# Results from prediction.
# """

# pred_dict = {}
# for y in self.Y:
# # Print message
# if self.verbosity > 1:
# print("\n## Predicting target %s" % str(y))
# if pred_params is not None:
# for key in list(pred_params):
# print("%s = %s" % (key, pred_params[key]))
# # Default value for pred_params
# if pred_params is None:
# pred_params = {}
# # Check this is a valid target
# if y not in self.fit_results:
# raise ValueError("y = %s not yet fitted" % str(y))
# # Construct the array form of the data
# # Check if we've passed a new dataframe object
# observation_array, xyz = \
# self.dataframe.construct_array(X=self.X, Y=[y] + self.Z, Z=self.conditions,
# tau_max=self.tau_max,
# # mask=self.test_mask,
# mask_type=self.mask_type,
# cut_off=self.cut_off,
# verbosity=self.verbosity)

# intervention_array = np.copy(observation_array)
# if intervention_data is not None:
# tmp_array, _ = intervention_data.construct_array(X=self.X, Y=[y] + self.Z,
# Z=self.conditions,
# tau_max=self.tau_max,
# mask_type=self.mask_type,
# cut_off=self.cut_off,
# verbosity=self.verbosity)

# # Only replace X-variables in intervention_array (necessary if lags of
# # X are in Z...)
# for index in np.where(xyz==0)[0]:
# intervention_array[index] = tmp_array[index]

# if self.conditions is not None and conditions_data is not None:
# tmp_array, _ = conditions_data.construct_array(X=self.X, Y=[y] + self.Z,
# Z=self.conditions,
# tau_max=self.tau_max,
# mask_type=self.mask_type,
# cut_off=self.cut_off,
# verbosity=self.verbosity)

# # Only replace condition-variables in intervention_array
# # (necessary if lags of X are in Z...)
# for index in np.where(xyz==2)[0]:
# intervention_array[index] = tmp_array[index]

# # Transform the data if needed
# a_transform = self.fit_results[y]['data_transform']
# if a_transform is not None:
# intervention_array = a_transform.transform(X=intervention_array.T).T
# # Cache the test array
# self.intervention_array = intervention_array
# # Run the predictor, for Y only the Z-part is used, the first index is y
# predictor_indices = list(np.where(xyz==0)[0]) \
# + list(np.where(xyz==1)[0][1:]) \
# + list(np.where(xyz==2)[0])
# predictor_array = intervention_array[predictor_indices, :].T

# pred_dict[y] = self.fit_results[y]['model'].predict(
# X=predictor_array, **pred_params)

# # print(pred_dict[y])
# if self.conditions is not None and conditions_data is not None:

# a_conditional_model = deepcopy(self.conditional_model)

# # Fit Y|do(X) on S
# conditions_array = observation_array[list(np.where(xyz==2)[0])]
# target_array = pred_dict[y] # array[np.where(xyz==1)[0][0], :]

# if a_transform is not None:
# conditions_array = a_transform.transform(X=conditions_array.T).T
# target_array = a_transform.transform(X=target_array.T).T

# a_conditional_model.fit(X=conditions_array.T, y=target_array)
# self.fit_results[y]['conditional_model'] = a_conditional_model

# # Now predict conditional causal effect for new conditions
# tmp_array, _ = conditions_data.construct_array(X=self.X, Y=[y] + self.Z,
# Z=self.conditions,
# tau_max=self.tau_max,
# mask_type=self.mask_type,
# cut_off=self.cut_off,
# verbosity=self.verbosity)

# # Construct conditions array
# new_conditions_array = tmp_array[list(np.where(xyz==2)[0])]

# if a_transform is not None:
# new_conditions_array = a_transform.transform(X=new_conditions_array.T).T

# pred_dict[y] = a_conditional_model.predict(
# X=new_conditions_array.T, **pred_params)

# return pred_dict


def get_fit(self, all_parents,
selected_variables=None,
tau_max=None,
Expand Down
Loading

0 comments on commit b833174

Please sign in to comment.