-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtask_16_readiness_potential_lineplot.py
140 lines (125 loc) · 4.41 KB
/
task_16_readiness_potential_lineplot.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
"""Perform and save time frequency analysis of given files."""
from __future__ import annotations
from pathlib import Path
from typing import Literal
import pandas as pd
import pte_decode
from matplotlib import pyplot as plt
from numpy import True_
import motor_intention.plotting_settings
import motor_intention.project_constants as constants
PLOT_DIR = constants.PLOTS / "readiness_potential"
PLOT_DIR.mkdir(exist_ok=True, parents=True)
RP_DIR = constants.DERIVATIVES / "readiness_potential"
IN_PATHS = {
(ch, stim): RP_DIR / f"stim_{stim.lower()}" / ch / "readiness_potential.csv"
for ch in ("ecog", "dbs")
for stim in ("Off", "On")
}
def task_rp_lineplot(
in_paths: dict[
tuple[Literal["ecog", "dbs"], Literal["Off", "On"]],
Path,
] = IN_PATHS,
show_plots: bool = False,
) -> None:
"""Main function of this script."""
motor_intention.plotting_settings.activate()
motor_intention.plotting_settings.medoff_medon_stimon()
for channel_type in ("ecog", "dbs"):
rp_lineplot(channel_type, in_paths, show_plots)
def rp_lineplot(
channel_type: Literal["ecog", "dbs"],
in_paths: dict[
tuple[Literal["ecog", "dbs"], Literal["Off", "On"]],
Path,
],
show_plots: bool = False,
) -> None:
fig, axs = plt.subplots(3, 1, sharex=True, sharey=True, figsize=(2.3, 3.4))
i = 0
legend = True_
ch_str = "motorcortex" if channel_type == "ecog" else "stn"
outpath = PLOT_DIR / (f"rp_lineplot_{ch_str}.svg")
for (ch_type, stimulation), in_path in in_paths.items():
if ch_type != channel_type:
continue
PIPELINE = in_path.parents[0].name
IN_DIR = constants.DERIVATIVES / "readiness_potential" / PIPELINE / ch_type
X_LABEL = "Time [s]"
if ch_type == "ecog":
Y_LIMS = (-68, 10)
elif ch_type == "dbs":
Y_LIMS = (-26, 4)
else:
msg = f"Unknown ch_type: {ch_type}"
raise ValueError(msg)
Y_LABEL = "Voltage [µV]"
THRESHOLD = 0.0
CORRECTION_METHOD = "cluster_pvals"
ALPHA = 0.05
N_PERM = 10000
rp = pd.read_csv(
IN_DIR / "readiness_potential.csv",
dtype={
"Subject": str,
"Medication": str,
"Stimulation": str,
"Channels": str,
},
).replace({"MotorCortex": "Motor Cortex"})
if stimulation == "Off":
cond = "Medication"
else:
cond = "Stimulation"
rp = rp.query("Medication == 'OFF'")
rp = rp.set_index(["Subject", "Medication", "Stimulation", "Channels"])
times = rp.columns.to_numpy(dtype=float)
rps = {
"OFF": rp.query(f"{cond} == 'OFF'").to_numpy().T,
"ON": rp.query(f"{cond} == 'ON'").to_numpy().T,
}
conds = ("OFF", "ON") if cond == "Medication" else ("ON",)
for cond in conds:
color = plt.rcParams["axes.prop_cycle"].by_key()["color"][i]
x_label = X_LABEL if i == 2 else None
pte_decode.lineplot_single(
data=rps[cond],
times=times,
ax=axs[i],
# figsize=(6.4, 3.2),
outpath=None,
x_label=x_label,
y_label=Y_LABEL,
y_lims=None, # Y_LIMS,
title=None, # f"{COND_ABB} {cond}",
threshold=THRESHOLD,
color=color,
alpha=ALPHA,
n_perm=N_PERM,
correction_method=CORRECTION_METHOD,
two_tailed=False,
one_tailed_test="smaller",
add_vline=0.0,
print_n=True,
legend=legend,
show=False,
)
axs[i].set_xlim([-3, 2])
axs[i].set_xticks([-3, 0, 2])
axs[i].set_ylim([Y_LIMS[0], Y_LIMS[1]])
axs[i].set_yticks([Y_LIMS[0], 0, Y_LIMS[1]])
axs[i].spines["left"].set_position(("outward", 3))
axs[i].spines["bottom"].set_position(("outward", 3))
legend = False
i += 1
for ax in axs:
bottom_lim, top_lim = ax.get_ylim()
print(f"{bottom_lim = }, {top_lim = }")
motor_intention.plotting_settings.save_fig(fig, outpath)
if show_plots:
plt.show()
else:
plt.close()
if __name__ == "__main__":
task_rp_lineplot(show_plots=True)