Analysis of acute stroke-like lesions in MELAS: Distribution, potential boundaries and spreading pattern.

1. import numpy as np
2. from numpy.linalg import eig
3. from numpy.linalg import norm
4. from scipy.optimize import minimize_scalar
5. import matplotlib.pyplot as plt
6. import seaborn as sns
9. def whiten(X, use_np=True):
10. """Utility function to whiten data without zero-mean centering. Whitening means removing correlation between
11. features and making the individual features have unit variance.
13. :param
14. X : np.array
15. Data matrix. This is assumed to be in the (uncommon) format n_features x n_samples
16. use_np : bool, optional (default: `True`)
17. Whether to use numpy to compute the covariance matrix
19. :returns
20. Z : np.array
21. Data matrix, whitened to remove covariance between the features
22. """
24. if use_np:
25. C_X = np.cov(X, rowvar=True)
26. else:
27. C_X = (X - X.mean(1)) @ (X - X.mean(1)).T
28. D, E = eig(C_X)
29. V = E @ np.diag(1 / np.sqrt(D)) @ E.T
30. Z = V @ X
32. return Z
35. def rotation(phi):
36. """Create 2D rotation matrix
38. :param
39. phi : float
40. The angle by which we want to rotate.
42. :returns
43. A : np.array
44. A 2D rotation matrix
45. """
47. return np.array([[np.cos(phi), np.sin(phi)],
48. [-np.sin(phi), np.cos(phi)]])
51. def loss(Y):
52. """Compute the loss for a given reconstruction Y.
53. This will simply be the sum of squared elements in the restriction
54. of Y to it's negative elements
56. :param
57. Y : np.array
58. Data matrix, reconstrcution of the sources
60. :returns
61. l : np.float
62. The loss
63. """
64. # restrict Y to it's negative elements
65. n_samples = Y.shape[1]
66. Y_neg = np.where(Y < 0, Y, 0)
68. return 1 / (2 * n_samples) * norm(Y_neg, ord='fro') ** 2
71. def obj_fun(phi, Z):
72. """Objective to be used for finding the optimum rotation angle
74. :param
75. phi : float
76. Rotation angle for which we wish to compute the los
77. Z : np.matrix
78. Whitened data matrix. Must have two rows, each corresponding to one feature.
80. :returns
81. l : float
82. loss corresponding to a rotation of Z in 2D around phi
83. """
85. # check input
86. if Z.shape[0] != 2:
87. raise ValueError('Z has more than two features.')
89. # rotate the data
90. W = rotation(phi)
91. Y = W @ Z
93. return loss(Y)
96. def givens(n, i, j, phi):
97. """Compute n-dimensional givens rotation
99. :param
100. n : int
101. Dimension of the rotation matrix to be computed
102. i, j : int
103. Dimensions i and j define the surface we wish to rotate in
104. phi : float
105. Rotation angle
107. :returns
108. R : np.array
109. Given's rotation
111. """
112. R = np.eye(n)
113. R[i, i], R[j, j] = np.cos(phi), np.cos(phi)
114. R[i, j], R[j, i] = np.sin(phi), -np.sin(phi)
116. return R
119. def torque(Y):
120. """Compute torque values of Y.
122. These correspond to the gradient if different directions, where
123. each direction is a possible rotation in a surface defines by two axis. The resulting matrix of
124. torque values will have zeroes on the diagonal and will by symmetric.
126. :param
127. Y : np.matrix
128. Reconstruction of the sourdes
130. :returns
131. t_max : float
132. Maximum torque value found
133. ixs : tuple
134. i-j coordinates corresponding to the max. This defines a hyperplane in n-dimensional space.
135. G : np.array
136. Matrix of torque values. Symmetric and zero on the diagonal. Only the upper half is computed.
137. """
139. # compute the rectified parts of Y
140. Y_pos = np.where(Y > 0, Y, 0)
141. Y_neg = np.where(Y < 0, Y, 0)
143. # compute torque values
144. n = Y.shape[0]
145. G = np.zeros((n, n))
146. for i in range(n):
147. for j in range(i + 1, n):
148. G[i, j] = np.dot(Y_pos[i, :], Y_neg[j, :]) - np.dot(Y_neg[i, :], Y_pos[j, :])
150. # find max and corresponding indices
151. t_max = np.amax(np.abs(G))
152. result = np.where(np.abs(G) == t_max)
153. ixs = [result[i][0] for i in range(len(result))]
155. return t_max, ixs, G
158. def run_nn_ica(X, t_tol=1e-1, t_neg=None, verbose=1, i_max=1e3, print_all=100,
159. whiten_mat=True, keep='last', return_all=False):
160. """Algorithm to run non-negative indipendent component analysis.
162. Given some data X, find matrices A and S such that
163. X = A S,
164. where S is non-negative and has indipendent rows. We pose no
165. constraint on the matrix A.
167. This algorithm is implemented as described in Plumbley, 2003, and
168. relies upon whitening and rotating the data. This is guaranteed to
169. converge only if the sources are 'well grounded', i.e. have probability
170. down to zero. Note that this is the implementation for a square mixing
171. matrix A.
173. Parameters
174. --------
175. X : np.array,
176. The data matrix of shape (n_features, n_samples)
177. t_tol : float, optional (default: `1e-1`)
178. Stopping tolerance. If the maximum torque falls below this
179. value, stop.
180. t_neg: float, optional (default: `None`)
181. Stopping number of negative elements. If #negative elements crosses
182. this threshold, stop.
183. verbose : int, optional (default: `1`)
184. How much output to give
185. i_max : int, optional (default: `1e3`)
186. Maximum number of iterations
187. print_all : int, optional (default: `100`)
188. Print every print_all iterations
189. whiten : bool, optional (default: `True`)
190. whether to whiten the input matrix
191. keep: Str, optional (default: `'last'`)
192. which reconstruction to keep, possible options are:
193. `'last'`, `'best_neg'`, `'best_tol'`
194. and correspond to last, smallest #negative elements and smallest tolerance
195. respectively
196. return_all: bool, optional (default: `False`)
197. whether to return all of the progress of Y, W
198. returns Y_best, W_best, Z, t_max_arr, ys, ws
200. Returns
201. --------
202. Y : np.array
203. The reconstructed sources, up to scaling and permutation
204. W : np.array
205. The final rotation matrix
206. Z : np.array
207. The whitened data
208. t_max_arr : np.array
209. The maximum torque values for each iteration
210. """
212. assert keep in ('last', 'best_neg', 'best_tol'), f'Unknown selection criterion `{keep}`.'
215. def set_best(Y, W, tol):
216. nonlocal Y_best, W_best, tol_best
218. if return_all:
219. ys.append(Y)
220. ws.append(W)
222. if keep == 'last':
223. is_better = True
224. if keep == 'best_neg':
225. is_better = Y_best is None or np.sum(Y_best < 0) > np.sum(Y < 0)
226. else:
227. is_better = tol_best > tol
229. if is_better:
230. Y_best, W_best = Y, W
231. tol_best = tol # need to keep memory of this
233. # lists that records the progress of Y and W
234. ys = []
235. ws = []
237. # best values and tolerance so far
238. Y_best, W_best = None, None
239. tol_best = np.inf
241. t_neg = -np.inf if t_neg is None else t_neg
243. # initialise
244. n = X.shape[0]
245. W = np.eye(n)
246. t_max_arr = []
248. # whiten the data
249. Z = whiten(X) if whiten_mat else X
250. Y = W @ Z
252. i = 0
253. while True:
254. # compute the max torque of Y and corresponding indices
255. t_max, ixs, _ = torque(Y)
256. t_max_arr.append(t_max)
258. set_best(Y, W, t_max)
260. if t_max < t_tol or np.sum(Y_best < 0) < t_neg: # converged
261. print('=' * 10)
262. print('i = {}, t_max = {:.2f}, ixs = {}, #negative = {}'.format(i, t_max, ixs, np.sum(Y_best < 0)))
263. print('Converged. Returning the reconstruction.')
264. if return_all:
265. return Y_best, W_best, Z, t_max_arr, ys, ws
267. return Y_best, W_best, Z, t_max_arr
269. if i > i_max and t_max > t_max_arr[-2]: # failed to converge
270. print('=' * 10)
271. print('i = {}, t_max = {:.2f}, ixs = {}, #negative = {}'.format(i, t_max, ixs, np.sum(Y_best < 0)))
272. print(f'Error: Failed to converge. Returning current matrices.')
273. if return_all:
274. return Y_best, W_best, Z, t_max_arr, ys, ws
276. return Y_best, W_best, Z, t_max_arr
278. # print some information
279. if (verbose > 0) and (i % print_all == 0):
280. print('i = {}, t_max = {:.2f}, ixs = {}, #negative = {}'.format(i, t_max, ixs, np.sum(Y < 0)))
282. # reduce to axis pair, find rotation angle and construct givens matrix
283. Y_red = Y[ixs, :]
284. opt_res = minimize_scalar(fun=obj_fun, bounds=(0, 2 * np.pi), method='bounded', args=Y_red)
285. R = givens(n, ixs[0], ixs[1], opt_res['x'])
287. # update the rotation matrix W and the reconstruction matrix Y
288. W = R @ W
289. Y = R @ Y
291. i += 1
293. def plot_ica_reconstruction(data, labels=None, *args, **kwargs):
294. """Plotting function for ICA reconstruction
296. Parameters
297. --------
298. data: List[np.ndarray]
299. Matrices to be plotted.
300. titles: Union[List, Str, NoneType], optional (default: `None`)
301. Corresponding titles, if present
302. *args, **kwargs:
303. additional arguments for `sns.distplot`
305. Returns
306. --------
307. None
308. """
310. plt.close('all')
312. n_sources, n_mat = data[0].shape[0], len(data)
313. for i, d in enumerate(data):
314. assert d.shape[0] == n_sources, f'Wrong number of features for `data[{i}].shape[0]` == `{d.shape[0]}` != `{n_sources}`.'
316. if labels is not None and len(labels) != n_mat:
317. print(f'Inconsistent number of titles given `{len(labels)}` != `{n_mat}`, showing no labels.')
318. labels = None
320. fig, axes = plt.subplots(nrows=n_sources, ncols=n_mat, figsize=(4 * n_mat, 4 * n_sources))
321. if not isinstance(axes, np.ndarray):
322. axes = np.array([[axes]])
323. elif axes.ndim != 2:
324. # treat only 1 source specifically
325. axes = np.expand_dims(axes, axis=(n_sources) != 1)
327. for i, row in enumerate(axes):
328. for j, (A, ax) in enumerate(zip(data, row)):
329. sns.distplot(A[i, :], *args, ax=ax, kde=False,
330. axlabel='{}_{}'.format(labels[j], i) if labels is not None else None,
331. **kwargs)
333. fig.suptitle("Non-negative ICA")
334. fig.tight_layout(rect=[0, 0, 1, 0.97])
335. fig.show()

main.py at commit 27c801b, under BSD-3-Clause · at the source