Coverage for /builds/tms-localization/papers/tmsloc_proto/scripts/05_merge_exp_data.py: 82%

1#!/usr/bin/env python 

2 

3""" 

4Merges experimental data and saves results in fn_exp specified in subject object file. 

5Also corrects coil <-> head distance in separate file '*_corrected.csv'. 

6 

7merge_exp_data.py -s <subjectobjfile> -e <exp_id> -m <mesh_idx> -o -d -r 

8 

9-s <str> ... filename (incl. path) of subject object file 

10-e <int> ... experiment ID 

11-m <int> ... mesh idx 

12-o ... coil position outlier correction (for conditions with multiple same stimulation locations) 

13-d ... coil <-> head distance correction 

14-r ... remove coil position with distance > +-5 mm from skin surface 

15-p ... plot MEPs 

16--start_mep ... Start of time frame after TMS pulse where p2p value is evaluated (in ms) 

17--end_mep ... End of time frame after TMS pulse where p2p value is evaluated (in ms) 

18""" 

19 

20import os 

21import h5py 

22import pynibs 

23import argparse 

24import numpy as np 

25import pandas as pd 

26 

27# Read input 

28parser = argparse.ArgumentParser(description='Run SimNIBS to determine electric fields\n') 

29parser.add_argument('-s', '--fn_subject', 

30 help='Filename (incl. path) of subject object file', required=True, type=str) 

31parser.add_argument('-e', '--exp_idx', 

32 help='Experiment ID', required=True, type=str) 

33parser.add_argument('-m', '--mesh_idx', 

34 help='Mesh idx', required=True) 

35parser.add_argument('-o', '--outlier_correction', 

36 help='Remove outliers of coil position.', required=False, 

37 action='store_true') 

38parser.add_argument('-d', '--distance_correction', 

39 help='Perform coil <-> head distance correction', required=False, 

40 action='store_true') 

41parser.add_argument('-r', '--remove_coil_skin_distance_outlier', 

42 help='Remove coil position with distance > -5/+2 mm from skin surface', required=False, 

43 action='store_true') 

44parser.add_argument('-p', '--plot', 

45 help='Plot MEPs', required=False, 

46 action='store_true') 

47parser.add_argument('--start_mep', 

48 help='Start of time frame after TMS pulse where p2p value is evaluated (in ms)', required=False, 

49 default=18) 

50parser.add_argument('--end_mep', 

51 help='End of time frame after TMS pulse where p2p value is evaluated (in ms)', required=False, 

52 default=35) 

53 

54args = parser.parse_args() 

55 

56fn_subject = os.path.abspath(args.fn_subject) 

57exp_id = args.exp_idx 

58mesh_idx = args.mesh_idx 

59start_mep = float(args.start_mep) 

60end_mep = float(args.end_mep) 

61 

62if args.outlier_correction: 

63 coil_outlier_corr = True 

64else: 

65 coil_outlier_corr = False 

66 

67if args.distance_correction: 

68 coil_distance_corr = True 

69else: 

70 coil_distance_corr = False 

71 

72if args.remove_coil_skin_distance_outlier: 

73 remove_coil_skin_distance_outlier = True 

74else: 

75 remove_coil_skin_distance_outlier = False 

76 

77if args.plot: 

78 plot = True 

79else: 

80 plot = False 

81 

82# load subject object 

83subject = pynibs.load_subject(fn_subject) 

84drop_mep_idx = None 

85mep_onsets = None 

86cfs_data_column = 0 

87channels = ["channel_0"] 

88 

89if 'mep_onsets' in subject.exp[exp_id].keys(): 

90 mep_onsets = subject.exp[exp_id]['mep_onsets'] 

91 

92if 'drop_mep_idx' in subject.exp[exp_id].keys(): 

93 drop_mep_idx = subject.exp[exp_id]['drop_mep_idx'] 

94 

95if "cfs_data_column" in subject.exp[exp_id].keys(): 

96 cfs_data_column = subject.exp[exp_id]["cfs_data_column"] 

97 channels = ["channel_" + str(i) for i in range(len(cfs_data_column))] 

98 

99if "channels" in subject.exp[exp_id].keys(): 

100 channels = subject.exp[exp_id]["channels"] 

101 cfs_data_column = range(len(channels)) 

102 

103if "fn_exp_hdf5" in subject.exp[exp_id].keys(): 

104 if os.path.exists(subject.exp[exp_id]["fn_exp_hdf5"][0]): 

105 os.remove(subject.exp[exp_id]["fn_exp_hdf5"][0]) 

106 

107if "fn_exp_csv" in subject.exp[exp_id].keys(): 

108 if os.path.exists(subject.exp[exp_id]["fn_exp_csv"][0]): 

109 os.remove(subject.exp[exp_id]["fn_exp_csv"][0]) 

110 

111# merged data will be saved in mesh specific subfolder 

112subject.exp[exp_id]["fn_exp_hdf5"] = [os.path.join(subject.subject_folder, "exp", exp_id, f"mesh_{mesh_idx}", "experiment.hdf5")] 

113 

114if not os.path.exists(os.path.split(subject.exp[exp_id]["fn_exp_hdf5"][0])[0]): 

115 os.makedirs(os.path.split(subject.exp[exp_id]["fn_exp_hdf5"][0])[0]) 

116 

117# merge experimental data to experiment.hdf5 

118verbose = True 

119if subject.exp[exp_id]['nnav_system'].lower() == "localite": 

120 

121 # this is the csv version 

122 # pynibs.merge_exp_data_localite(subject=subject, exp_id=exp_id, mesh_idx=mesh_idx, 

123 # coil_outlier_corr=coil_outlier_corr, 

124 # verbose=verbose, version='csv') 

125 

126 # this is the hdf5 version 

127 pynibs.expio.localite.merge_exp_data_localite(subject=subject, 

128 exp_idx=exp_id, 

129 mesh_idx=mesh_idx, 

130 coil_outlier_corr_cond=coil_outlier_corr, 

131 remove_coil_skin_distance_outlier=remove_coil_skin_distance_outlier, 

132 coil_distance_corr=coil_distance_corr, 

133 verbose=verbose, 

134 drop_mep_idx=drop_mep_idx, 

135 mep_onsets=mep_onsets, 

136 channels=channels, 

137 cfs_data_column=cfs_data_column, 

138 plot=plot, 

139 start_mep=start_mep, 

140 end_mep=end_mep) 

141 

142elif subject.exp[exp_id]['nnav_system'].lower() == "visor": 

143 pynibs.expio.visor.merge_exp_data_visor(subject=subject, 

144 exp_id=exp_id, 

145 mesh_idx=mesh_idx, 

146 verbose=verbose) 

147 

148 

149elif subject.exp[exp_id]['nnav_system'].lower() == "brainsight": 

150 pynibs.expio.brainsight.merge_exp_data_brainsight(subject=subject, 

151 exp_idx=exp_id, 

152 mesh_idx=mesh_idx, 

153 coil_outlier_corr_cond=coil_outlier_corr, 

154 remove_coil_skin_distance_outlier=remove_coil_skin_distance_outlier, 

155 coil_distance_corr=coil_distance_corr, 

156 verbose=verbose, 

157 plot=plot, 

158 start_mep=start_mep, 

159 end_mep=end_mep) 

160else: 

161 raise NotImplementedError(f"Neuronavigation system {subject.exp[exp_id]['nnav_system']} not supported!" 

162 f"Use either 'Localite' or 'Visor'...") 

163 

164# Plot coil positions with MEP data 

165####################################################################################################################### 

166if 'plot_coil_ps_fn' in subject.exp[exp_id].keys(): 

167 fn_coil_pos = subject.exp[exp_id]['plot_coil_ps_fn'] 

168else: 

169 fn_coil_pos = "plot_coil_pos.hdf5" 

170fn_coil_pos = os.path.join(os.path.split(subject.exp[exp_id]["fn_exp_hdf5"][0])[0], fn_coil_pos) 

171df_exp = pd.read_hdf(subject.exp[exp_id]["fn_exp_hdf5"][0], "phys_data/postproc/EMG") 

172 

173mep = dict() 

174for channel in subject.exp[exp_id]["channels"]: 

175 mep[channel] = np.array(df_exp[f"p2p_{channel}"]) 

176 

177# plot coil positions 

178pynibs.create_stimsite_from_exp_hdf5(fn_exp=subject.exp[exp_id]["fn_exp_hdf5"][0], 

179 fn_hdf=fn_coil_pos, 

180 datanames=subject.exp[exp_id]["channels"], 

181 data=np.vstack([mep[channel] for channel in mep.keys()]).T, 

182 overwrite=True) 

183 

184# Create matsimnibs.hdf5 

185####################################################################################################################### 

186if 'fn_matsimnibs' in subject.exp[exp_id].keys(): 

187 fn_matsimnibs = subject.exp[exp_id]['fn_matsimnibs'] 

188else: 

189 fn_matsimnibs = "matsimnibs.hdf5" 

190fn_matsimnibs = os.path.join(os.path.split(subject.exp[exp_id]["fn_exp_hdf5"][0])[0], fn_matsimnibs) 

191fn_exp_hdf5 = subject.exp[exp_id]["fn_exp_hdf5"][0] 

192 

193print(f"Reading: {fn_exp_hdf5}") 

194df_stim_data = pd.read_hdf(fn_exp_hdf5, "stim_data") 

195 

196matsimnibs = np.zeros((4, 4, df_stim_data.shape[0])) 

197 

198for i in range(df_stim_data.shape[0]): 

199 matsimnibs[:, :, i] = df_stim_data["coil_mean"].iloc[i] 

200 

201print(f"Writing: {fn_matsimnibs}") 

202with h5py.File(fn_matsimnibs, "w") as f: 

203 f.create_dataset("matsimnibs", data=matsimnibs)