OSCR

MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § MATERIALS AND METHODS › Magnetic resonance artificial intelligence velocimetry ↔ src/instant_aiv/models/NNpp.py, lines 201–233 · score 0.52 · activation function, weight normalization, connections, network, layers, models

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 824 lines · 31 KB · no license · 1 match

  1. # Libraries
  2. import numpy as np
  3. from jax import jit, vmap, grad
  4. import jax.numpy as jnp
  5. from jax.nn import sigmoid
  6. from typing import Tuple
  7. from typing import Tuple, List, Dict, Sequence
  8. import h5py
  9. from instant_aiv.models.metrics import *
  10. from instant_aiv.manage.dataloader import *
  11. import optax
  12. #Initialization
  13. from typing import Tuple
  14. def glorot_normal(in_dim: int, out_dim: int) -> jnp.ndarray:
  15. glorot_stddev = np.sqrt(2.0 / (in_dim + out_dim))
  16. return jnp.array(np.random.normal(loc=0.0, scale=glorot_stddev, size=(in_dim, out_dim)))
  17. def init_params(layers: List[int], initialization_type: str = 'xavier',Network_type: str='mlp',degree: int =5,Use_ResNet: bool =False) -> dict:
  18. def init_adaptive_params():
  19. F = 0.1 * jnp.ones(3 * len(layers) - 1)
  20. A = 0.1 * jnp.ones(3 * len(layers) - 1)
  21. return [{"a0": A[3*i], "a1": A[3*i + 1], "a2": A[3*i + 2],
  22. "f0": F[3*i], "f1": F[3*i + 1], "f2": F[3*i + 2]}
  23. for i in range(len(layers) - 1)]
  24. #Define Models:
  25. def init_layer_mlp(in_dim, out_dim):
  26. if initialization_type == 'xavier':
  27. W = glorot_normal(in_dim, out_dim)
  28. elif initialization_type == 'normal':
  29. W = jnp.array(np.random.normal(size=(in_dim, out_dim)))
  30. b = jnp.zeros(out_dim)
  31. g = jnp.ones(out_dim)
  32. return {"W": W, "b": b, "g": g}
  33. def init_layer_ResNet(in_dim, out_dim):
  34. if initialization_type == 'xavier':
  35. W1 = glorot_normal(in_dim, out_dim)
  36. W2 = glorot_normal(out_dim, out_dim)
  37. W3 = glorot_normal(out_dim, out_dim)
  38. elif initialization_type == 'normal':
  39. W1 = jnp.array(np.random.normal(size=(in_dim, out_dim)))
  40. W2 = jnp.array(np.random.normal(size=(out_dim, out_dim)))
  41. W3 = jnp.array(np.random.normal(size=(out_dim, out_dim)))
  42. b1 = jnp.zeros(out_dim)
  43. b2 = jnp.zeros(out_dim)
  44. b3 = jnp.zeros(out_dim)
  45. g1 = jnp.ones(out_dim)
  46. g2 = jnp.ones(out_dim)
  47. g3 = jnp.ones(out_dim)
  48. alpha=0.0
  49. return {"W": W1, "b": b1, "g": g1,
  50. "W2": W2, "b2": b2, "g2": g2,
  51. 'alpha':alpha}
  52. def init_layer_kan(in_dim, out_dim,degree=degree):
  53. std=1 / (in_dim * (degree + 1))
  54. W =jnp.array(np.random.normal(loc=0.0, scale=std, size=(in_dim, out_dim,degree+1)))
  55. b = jnp.zeros(out_dim)
  56. g = jnp.ones(out_dim)
  57. return {"W": W, "b": b, "g": g}
  58. def init_layer_kan_ResNet(in_dim, out_dim,degree=degree):
  59. std=1 / (in_dim * (degree + 1))
  60. W1 =jnp.array(np.random.normal(loc=0.0, scale=std, size=(in_dim, out_dim,degree+1)))
  61. W2 =jnp.array(np.random.normal(loc=0.0, scale=std, size=(in_dim, out_dim,degree+1)))
  62. b1 = jnp.zeros(out_dim)
  63. b2 = jnp.zeros(out_dim)
  64. g1 = jnp.ones(out_dim)
  65. g2 = jnp.ones(out_dim)
  66. alpha=0.0
  67. return {"W": W1, "b": b1, "g": g1,
  68. "W2": W2, "b2": b2, "g2": g2,
  69. 'alpha':alpha}
  70. #Select model
  71. if Network_type.lower()=='mlp':
  72. if Use_ResNet:
  73. init_layer_params=init_layer_ResNet
  74. else:
  75. init_layer_params=init_layer_mlp
  76. elif Network_type.lower()[:3]=='kan':
  77. if Use_ResNet:
  78. init_layer_params=init_layer_kan_ResNet
  79. else:
  80. init_layer_params=init_layer_kan
  81. else:
  82. print(f'Error: {Network_type.lower()} is not a valid option. The available options are:mlp and kan.')
  83. print(f'Initializing:{Network_type} parameters.')
  84. params = [init_layer_params(layers[i], layers[i + 1]) for i in range(len(layers) - 1)]
  85. U1, b1, g1 = glorot_normal(layers[0], layers[1]), jnp.zeros(layers[1]), jnp.ones(layers[1])
  86. U2, b2, g2 = glorot_normal(layers[0], layers[1]), jnp.zeros(layers[1]), jnp.ones(layers[1])
  87. mMLP_params = [{"U1": U1, "b1": b1, "g1": g1, "U2": U2, "b2": b2, "g2": g2}]
  88. return {
  89. 'params': params,
  90. 'AdaptiveAF': init_adaptive_params(),
  91. 'mMLP': mMLP_params
  92. }
  93. def init_params_res(layers: List[int], initialization_type: str = 'xavier',Network_type: str='mlp',degree: int =5,Use_ResNet: bool =False) -> dict:
  94. def init_adaptive_params():
  95. F = 0.1 * jnp.ones(3 * len(layers) - 1)
  96. A = 0.1 * jnp.ones(3 * len(layers) - 1)
  97. return [{"a0": A[3*i], "a1": A[3*i + 1], "a2": A[3*i + 2],
  98. "f0": F[3*i], "f1": F[3*i + 1], "f2": F[3*i + 2]}
  99. for i in range(len(layers) - 1)]
  100. def init_layer_params(in_dim, out_dim):
  101. if initialization_type == 'xavier':
  102. W1 = glorot_normal(in_dim, out_dim)
  103. W2 = glorot_normal(out_dim, out_dim)
  104. W3 = glorot_normal(out_dim, out_dim)
  105. elif initialization_type == 'normal':
  106. W1 = jnp.array(np.random.normal(size=(in_dim, out_dim)))
  107. W2 = jnp.array(np.random.normal(size=(out_dim, out_dim)))
  108. W3 = jnp.array(np.random.normal(size=(out_dim, out_dim)))
  109. b1 = jnp.zeros(out_dim)
  110. b2 = jnp.zeros(out_dim)
  111. b3 = jnp.zeros(out_dim)
  112. g1 = jnp.ones(out_dim)
  113. g2 = jnp.ones(out_dim)
  114. g3 = jnp.ones(out_dim)
  115. alpha=0.0
  116. return {"W": W1, "b": b1, "g": g1,
  117. "W2": W2, "b2": b2, "g2": g2,
  118. 'alpha':alpha}
  119. params = [init_layer_params(layers[i], layers[i + 1]) for i in range(len(layers) - 1)]
  120. U1, b1, g1 = glorot_normal(layers[0], layers[1]), jnp.zeros(layers[1]), jnp.ones(layers[1])
  121. U2, b2, g2 = glorot_normal(layers[0], layers[1]), jnp.zeros(layers[1]), jnp.ones(layers[1])
  122. mMLP_params = [{"U1": U1, "b1": b1, "g1": g1, "U2": U2, "b2": b2, "g2": g2}]
  123. return {
  124. 'params': params,
  125. 'AdaptiveAF': init_adaptive_params(),
  126. 'mMLP': mMLP_params
  127. }
  128. def init_params_dict(layer_dict, initialization,Use_ResNet=False,Network_type='mlp',degree=5):
  129. print(f'You selected: Network {Network_type} with degree(if KAN) {degree}, initialization {initialization},Use_ResNet {Use_ResNet}')
  130. if Network_type.lower()=='mlp':
  131. if Use_ResNet:
  132. init_function=init_params_res
  133. else:
  134. init_function=init_params
  135. elif Network_type[:3].lower()=='kan':
  136. init_function=init_params
  137. initialized_params = {}
  138. for key, layer_structure in layer_dict.items():
  139. # Initialize parameters for each key
  140. params = init_function(layer_structure,
  141. initialization_type=initialization.lower(),
  142. Network_type=Network_type,
  143. degree=degree,
  144. Use_ResNet=Use_ResNet)
  145. # Store in the dictionary
  146. initialized_params[key] = params
  147. return initialized_params
  148. def FCN(params, X_in, M1, M2, activation_fn, norm_fn):
  149. """
  150. Fully Connected Network (FCN) with a given normalization and activation function.
  151. Parameters:
  152. - params: Dictionary containing parameters of the neural network layers.
  153. - X_in: Input tensor to the network.
  154. - M1, M2: Parameters for normalization function.
  155. - activation_fn: Callable activation function.
  156. - norm_fn: Callable normalization function.
  157. Returns:
  158. - Output tensor from the network.
  159. """
  160. params_N = params["params"]
  161. inputs = norm_fn(X_in, M1, M2)
  162. for layer in params_N[:-1]:
  163. outputs = activation_fn(jnp.dot(inputs, layer["W"]) + layer["b"])
  164. inputs = outputs
  165. W = params_N[-1]["W"]
  166. b = params_N[-1]["b"]
  167. outputs = jnp.dot(inputs, W) + b
  168. return outputs
  169. def FCN_WN(params, X_in, M1, M2, activation, norm_fn):
  170. """
  171. Fully Connected Network (FCN) with Weight Normalization.
  172. Parameters:
  173. - params: Dictionary containing parameters of the neural network layers.
  174. - X_in: Input tensor to the network.
  175. - M1, M2: Parameters for normalization function.
  176. - activation: Callable activation function.
  177. - norm_fn: Callable normalization function.
  178. Returns:
  179. - Output tensor from the network.
  180. """
  181. # Normalize the input
  182. H = norm_fn(X_in, M1, M2)
  183. # Iterate through the layers
  184. for i, layer_params in enumerate(params["params"]):
  185. W, b, g = layer_params["W"], layer_params["b"], layer_params["g"]
  186. # Weight normalization
  187. V = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  188. # Linear transformation
  189. H = g * jnp.matmul(H, V) + b
  190. # Apply activation function for all layers except the last one
  191. if i != len(params["params"]) - 1:
  192. H = activation(H)
  193. return H
  194. def FCN_WN_MMLP(params, X_in, M1, M2, activation, norm_fn):
  195. """
  196. Fully Connected Network (FCN) with Weight Normalization and Modified MLP (MMLP) transformations.
  197. Parameters:
  198. - params: Dictionary containing parameters of the neural network layers and MMLP parameters.
  199. - X_in: Input tensor to the network.
  200. - M1, M2: Parameters for normalization function.
  201. - activation: Callable activation function.
  202. - norm_fn: Callable normalization function.
  203. Returns:
  204. - Output tensor from the network.
  205. """
  206. # Normalize the input
  207. H = norm_fn(X_in, M1, M2)
  208. # Unpack MMLP parameters and apply weight normalization
  209. mMLP_params = params["mMLP"][0]
  210. U1, U2, b1, b2, g1, g2 = (mMLP_params[key] for key in ["U1", "U2", "b1", "b2", "g1", "g2"])
  211. U1_norm, U2_norm = U1 / jnp.linalg.norm(U1, axis=0, keepdims=True), U2 / jnp.linalg.norm(U2, axis=0, keepdims=True)
  212. # Calculate U and V transformations
  213. U = activation(g1 * jnp.dot(H, U1_norm) + b1)
  214. V = activation(g2 * jnp.dot(H, U2_norm) + b2)
  215. # Iterate through the layers
  216. for idx, layer in enumerate(params["params"][:-1]):
  217. W, b, g = layer["W"], layer["b"], layer["g"]
  218. # Apply weight normalization
  219. W_norm = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  220. # Compute activations and apply MMLP combination step
  221. H = activation(g * jnp.dot(H, W_norm) + b)
  222. H = jnp.multiply(H, U) + jnp.multiply(1 - H, V)
  223. # Process the last layer
  224. W, b, g = params["params"][-1]["W"], params["params"][-1]["b"], params["params"][-1]["g"]
  225. W_norm = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  226. H = g * jnp.dot(H, W_norm) + b
  227. return H
  228. def FCN_MMLP(params, X_in, M1, M2, activation, norm_fn):
  229. """
  230. Fully Connected Network (FCN) with Modified MLP (MMLP) transformations.
  231. Parameters:
  232. - params: Dictionary containing parameters of the neural network layers and MMLP parameters.
  233. - X_in: Input tensor to the network.
  234. - M1, M2: Parameters for normalization function.
  235. - activation: Callable activation function.
  236. - norm_fn: Callable normalization function.
  237. Returns:
  238. - Output tensor from the network.
  239. """
  240. # Normalize the input
  241. inputs = norm_fn(X_in, M1, M2)
  242. # Unpack MMLP parameters
  243. mMLP_params = params["mMLP"][0]
  244. U1, U2, b1, b2 = mMLP_params["U1"], mMLP_params["U2"], mMLP_params["b1"], mMLP_params["b2"]
  245. # Calculate U and V transformations
  246. U = activation(jnp.dot(inputs, U1) + b1)
  247. V = activation(jnp.dot(inputs, U2) + b2)
  248. # Iterate through all layers except the last
  249. for layer in params["params"][:-1]:
  250. W, b = layer["W"], layer["b"]
  251. # Compute activations
  252. act_values = activation(jnp.dot(inputs, W) + b)
  253. # MMLP combination step
  254. inputs = jnp.multiply(act_values, U) + jnp.multiply(1 - act_values, V)
  255. # Compute output from the last layer
  256. W, b = params["params"][-1]["W"], params["params"][-1]["b"]
  257. outputs = jnp.dot(inputs, W) + b
  258. return outputs
  259. # Save Resuls
  260. def save_list(Loss,path,name='loss-'):
  261. filename=path+name+".npy"
  262. np.save(filename, np.array(Loss))
  263. def save_MLP_params(params: List[Dict[str, np.ndarray]], save_path,WN=False,Mod_MLP=False):
  264. with h5py.File(save_path, "w") as f:
  265. for layer_idx, layer_params in enumerate(params):
  266. layer_group = f.create_group(f"Layer_{layer_idx/2.0:.2f}")
  267. if WN:
  268. if Mod_MLP:
  269. W, b, g= layer_params.values()
  270. layer_group.create_dataset("W", shape=W.shape, dtype=np.float32, data=W)
  271. layer_group.create_dataset("b", shape=b.shape, dtype=np.float32, data=b)
  272. layer_group.create_dataset("g", shape=g.shape, dtype=np.float32, data=g)
  273. else:
  274. W, b, g= layer_params.values()
  275. layer_group.create_dataset("W", shape=W.shape, dtype=np.float32, data=W)
  276. layer_group.create_dataset("b", shape=b.shape, dtype=np.float32, data=b)
  277. layer_group.create_dataset("g", shape=g.shape, dtype=np.float32, data=g)
  278. else:
  279. W, b = layer_params.values()
  280. layer_group.create_dataset("W", shape=W.shape, dtype=np.float32, data=W)
  281. layer_group.create_dataset("b", shape=b.shape, dtype=np.float32, data=b)
  282. def read_params(filename,WN=False):
  283. data = h5py.File(filename, 'r')
  284. recover_params=[]
  285. for layer in data.keys() :
  286. if WN:
  287. stored={'W':data[layer]['W'][:],'b': data[layer]['b'][:],'g': data[layer]['g'][:]}
  288. else:
  289. stored={'W':data[layer]['W'][:],'b': data[layer]['b'][:]}
  290. recover_params.append(stored)
  291. return recover_params
  292. def select_model(WN=False, Mod_MLP=False,Use_ResNet=False,Adaptive=False,Light=False,Network_type='mlp',degree=5):
  293. if Network_type.lower()=='mlp':
  294. model_map = {
  295. (True, True,False, False, False): FCN_WN_MMLP,
  296. (True, False,False, False, False): FCN_WN,
  297. (False, True,False, False, False): FCN_MMLP,
  298. (False, False,False, False, False): FCN,
  299. (False, False,True, False, False): ResNet,
  300. (False, False,True, False, True): ResNet_light,
  301. (True, False,True, False, False): WN_ResNet,
  302. (True, False,True, False, True): WN_ResNet_light,
  303. (True, False,True, True, False): WN_ResNet_adaptive,
  304. }
  305. return model_map[(WN, Mod_MLP,Use_ResNet,Adaptive,Light)]
  306. elif Network_type.lower()=='kan':
  307. if degree==3:
  308. return KAN_Net3
  309. elif degree==5:
  310. if Use_ResNet:
  311. return KAN5_ResNet
  312. else:
  313. return KAN_Net5
  314. elif degree==7:
  315. return KAN_Net7
  316. elif degree==9:
  317. return KAN_Net9
  318. elif degree==11:
  319. return KAN_Net11
  320. else:
  321. return KAN_Net
  322. elif Network_type.lower()=='kan_theta':
  323. return KAN_Net_theta
  324. def initialize_optimizer(lr0, decay_rate, lrf, decay_step, T_e,optimizer_type='Adam',weight_decay=1e-5):
  325. print('Optimizer',optimizer_type.lower())
  326. if optimizer_type.lower()=='adam':
  327. if decay_rate == 0 or lrf == lr0:
  328. print('No decay')
  329. return optax.adam(lr0), decay_step
  330. else:
  331. if decay_step == 0:
  332. decay_step = T_e * np.log(decay_rate) / np.log(lrf / lr0)
  333. print(f'The decay step will be {decay_step}')
  334. return optax.adam(optax.exponential_decay(lr0, decay_step, decay_rate,)),decay_step
  335. elif optimizer_type.lower()=='adamw':
  336. print('Weight decay:',weight_decay)
  337. if decay_rate == 0 or lrf == lr0:
  338. print('No decay')
  339. return optax.adamw(learning_rate=lr0, weight_decay=weight_decay), decay_step
  340. else:
  341. if decay_step == 0:
  342. decay_step = T_e * np.log(decay_rate) / np.log(lrf / lr0)
  343. print(f'The decay step will be {decay_step}')
  344. # Use adamw with the specified learning rate schedule
  345. return optax.adamw(optax.exponential_decay(lr0, decay_step, decay_rate), weight_decay=weight_decay), decay_step
  346. elif optimizer_type.lower()=='lion':
  347. if decay_rate == 0 or lrf == lr0:
  348. weight_decay=weight_decay*3
  349. print('No decay')
  350. return optax.lion(learning_rate=lr0, weight_decay=weight_decay), decay_step
  351. else:
  352. if decay_step == 0:
  353. weight_decay=weight_decay*3
  354. decay_step = T_e * np.log(decay_rate) / np.log(lrf / lr0)
  355. print(f'The decay step will be {decay_step}')
  356. # Use adamw with the specified learning rate schedule
  357. return optax.lion(optax.exponential_decay(lr0, decay_step, decay_rate), weight_decay=weight_decay), decay_step
  358. def load_params_dict(result_path, dataset_name, layer_dict, initialization, type='Test',Use_ResNet=False):
  359. loaded_params = {}
  360. for key in layer_dict.keys():
  361. # Construct the file path
  362. file_path = f"{result_path}{dataset_name}-{type}_params_{key}.h5"
  363. # Read parameters from the file
  364. raw_params = read_all_params(file_path)[0]
  365. # Initialize a test parameter set for getting lengths
  366. test_params = init_params(layer_dict[key], initialization_type=initialization.lower())
  367. params_length = len(test_params['params'])
  368. params_length_AF = len(test_params['AdaptiveAF'])
  369. # Reconstruct the parameters
  370. reconstructed_params = reconstruct_params(raw_params, params_length, params_length_AF,Use_ResNet)
  371. # Store in the dictionary
  372. loaded_params[key] = reconstructed_params
  373. return loaded_params
  374. def save_params_dict(params,result_path,dataset_name,type='Test',Use_ResNet=False):
  375. # Assuming 'result_path' and 'dataset_name' are defined elsewhere in your code
  376. for key, params in params.items():
  377. # Construct the file name for each set of parameters
  378. output_path = f"{result_path}{dataset_name}-{type}_params_{key}.h5"
  379. # Extract arrays from params
  380. params_to_save = [extract_arrays_from_params(params,Use_ResNet)]
  381. # Save the parameters
  382. save_all_params(params_to_save, output_path)
  383. print(f'Params {key} have been saved!')
  384. def transfer_params(params_s,params_t,levels=[0,1,2]):
  385. for level in levels:
  386. params_t['params'][level]=params_s['params'][level]
  387. return params_t
  388. def L2_regularization(params,beta=0.001):
  389. Loss_L2=0
  390. for param in params:
  391. Loss_L2=jnp.mean(param['W']**2)+Loss_L2
  392. return beta*Loss_L2
  393. def L1_regularization(params,beta=0.001):
  394. Loss_L2=0
  395. for param in params:
  396. Loss_L2=jnp.abs(jnp.mean(param['W']))+Loss_L2
  397. return beta*Loss_L2
  398. def L2_regularization_ResNet(params,beta=0.001):
  399. Loss_L2=0
  400. for param in params:
  401. Loss_L2=jnp.sum(param['W']**2)+jnp.sum(param['W2']**2)+Loss_L2
  402. return beta*Loss_L2
  403. # RESNET ARCHITECTURES
  404. def ResNet(params, X_in, M1, M2, activation_fn, norm_fn):
  405. def linear_layer(H,layer_params,W_key="W",b_key="b"):
  406. W, b= layer_params[W_key], layer_params[b_key]
  407. H= jnp.dot(H,W) + b
  408. return H
  409. params_N = params["params"]
  410. inputs = norm_fn(X_in, M1, M2)
  411. layer_params=params_N[0]
  412. inputs=activation_fn(linear_layer(inputs,layer_params,W_key="W",b_key="b"))
  413. for ly in range(1,len(params_N[:-1])):
  414. layer_params=params_N[ly]
  415. g= activation_fn(linear_layer(inputs,layer_params,W_key="W",b_key="b"))
  416. h= linear_layer(g,layer_params,W_key="W2",b_key="b2")
  417. inputs= activation_fn(h+inputs)
  418. outputs = jnp.dot(inputs, params_N[-1]["W"]) + params_N[-1]["b"]
  419. return outputs
  420. def ResNet_light(params, X_in, M1, M2, activation_fn, norm_fn):
  421. def linear_layer(H,layer_params,W_key="W",b_key="b"):
  422. W, b= layer_params[W_key], layer_params[b_key]
  423. H= jnp.dot(H,W) + b
  424. return H
  425. params_N = params["params"]
  426. inputs = norm_fn(X_in, M1, M2)
  427. layer_params=params_N[0]
  428. inputs=activation_fn(linear_layer(inputs,layer_params,W_key="W",b_key="b"))
  429. inputs0=jnp.copy(inputs)
  430. for ly in range(1,len(params_N[:-2])):
  431. layer_params=params_N[ly]
  432. g= activation_fn(linear_layer(inputs,layer_params,W_key="W",b_key="b"))
  433. h= linear_layer(g,layer_params,W_key="W2",b_key="b2")
  434. inputs= activation_fn(h)
  435. layer_params=params_N[-2]
  436. inputs=activation_fn(linear_layer(inputs+inputs0,layer_params,W_key="W",b_key="b"))
  437. layer_params=params_N[-1]
  438. outputs=activation_fn(linear_layer(inputs,layer_params,W_key="W",b_key="b"))
  439. return outputs
  440. def WN_ResNet(params, X_in, M1, M2, activation, norm_fn):
  441. def linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"):
  442. W, b, g = layer_params[W_key], layer_params[b_key], layer_params[g_key]
  443. # Weight normalization
  444. V = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  445. # Linear transformation
  446. H = g * jnp.matmul(H, V) + b
  447. return H
  448. params_N = params["params"]
  449. # Normalize the input
  450. H = norm_fn(X_in, M1, M2)
  451. # Encoder Layer
  452. H = activation(linear_layer(H,params_N[0],W_key="W",b_key="b",g_key="g"))
  453. # Iterate through the layers
  454. for i in range(1,len(params_N[:-1])):
  455. layer_params=params_N[i]
  456. F=activation(linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"))
  457. G=linear_layer(F,layer_params,W_key="W2",b_key="b2",g_key="g2")
  458. H=activation(G+H)
  459. H=linear_layer(H,params_N[-1],W_key="W",b_key="b",g_key="g")
  460. return H
  461. def WN_ResNet_adaptive(params, X_in, M1, M2, activation, norm_fn):
  462. def linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"):
  463. W, b, g = layer_params[W_key], layer_params[b_key], layer_params[g_key]
  464. # Weight normalization
  465. V = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  466. # Linear transformation
  467. H = g * jnp.matmul(H, V) + b
  468. return H
  469. params_N = params["params"]
  470. # Normalize the input
  471. H = norm_fn(X_in, M1, M2)
  472. # Encoder Layer
  473. H = activation(linear_layer(H,params_N[0],W_key="W",b_key="b",g_key="g"))
  474. # Iterate through the layers
  475. for i in range(1,len(params_N[:-1])):
  476. layer_params=params_N[i]
  477. F=activation(linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"))
  478. G=linear_layer(F,layer_params,W_key="W2",b_key="b2",g_key="g2")
  479. H=activation(layer_params["alpha"]*G+(1-layer_params["alpha"])*H)
  480. H=linear_layer(H,params_N[-1],W_key="W",b_key="b",g_key="g")
  481. return H
  482. def WN_ResNet_light(params, X_in, M1, M2, activation, norm_fn):
  483. def linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"):
  484. W, b, g = layer_params[W_key], layer_params[b_key], layer_params[g_key]
  485. # Weight normalization
  486. V = W / jnp.linalg.norm(W, axis=0, keepdims=True)
  487. # Linear transformation
  488. H = g * jnp.matmul(H, V) + b
  489. return H
  490. params_N = params["params"]
  491. # Normalize the input
  492. H = norm_fn(X_in, M1, M2)
  493. # Encoder Layer
  494. H = activation(linear_layer(H,params_N[0],W_key="W",b_key="b",g_key="g"))
  495. H0=jnp.copy(H)
  496. # Iterate through the layers
  497. for i in range(1,len(params_N[:-2])):
  498. layer_params=params_N[i]
  499. F=activation(linear_layer(H,layer_params,W_key="W",b_key="b",g_key="g"))
  500. G=linear_layer(F,layer_params,W_key="W2",b_key="b2",g_key="g2")
  501. H=activation(G)
  502. H=activation(linear_layer(H+H0,params_N[-2],W_key="W",b_key="b",g_key="g"))
  503. H=linear_layer(H,params_N[-1],W_key="W",b_key="b",g_key="g")
  504. return H
  505. def KAN_Net3(params, X_in, M1, M2, activation, norm_fn):
  506. def Cheby_KAN_layer3(x,layer_params,expanded_arr=[]):
  507. # Read chebyshev coefficients:
  508. cheby_coeffs= layer_params["W"]
  509. inputdim=cheby_coeffs.shape[0]
  510. outdim=cheby_coeffs.shape[1]
  511. # Normalize
  512. x = activation(x)
  513. # Reshape
  514. x = x.reshape((-1, inputdim, 1))
  515. x=jnp.stack((T0(x),
  516. T1(x),
  517. T2(x),
  518. T3(x)),axis=2)# Discard dummy dimension
  519. # Compute the Chebyshev interpolation
  520. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  521. # Remove extra dimension
  522. x= jnp.reshape(x, (outdim,))
  523. return x
  524. #Define params
  525. params_N = params["params"]
  526. #Normalize inputs
  527. x = norm_fn(X_in, M1, M2)
  528. for ly in range(len(params_N)):
  529. layer_params=params_N[ly]
  530. x=Cheby_KAN_layer3(x,layer_params)
  531. return x
  532. def KAN_Net5(params, X_in, M1, M2, activation, norm_fn):
  533. def Cheby_KAN_layer5(x,layer_params,expanded_arr=[]):
  534. # Read chebyshev coefficients:
  535. cheby_coeffs= layer_params["W"]
  536. inputdim=cheby_coeffs.shape[0]
  537. outdim=cheby_coeffs.shape[1]
  538. # Normalize
  539. x = activation(x)
  540. # Reshape
  541. x = x.reshape((-1, inputdim, 1))
  542. x=jnp.stack((T0(x),
  543. T1(x),
  544. T2(x),
  545. T3(x),
  546. T4(x),
  547. T5(x)),axis=2)# Discard dummy dimension
  548. # Compute the Chebyshev interpolation
  549. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  550. # Remove extra dimension
  551. x= jnp.reshape(x, (outdim,))
  552. return x
  553. #Define params
  554. params_N = params["params"]
  555. #Normalize inputs
  556. x = norm_fn(X_in, M1, M2)
  557. for ly in range(len(params_N)):
  558. layer_params=params_N[ly]
  559. x=Cheby_KAN_layer5(x,layer_params)
  560. return x
  561. def KAN_Net7(params, X_in, M1, M2, activation, norm_fn):
  562. def Cheby_KAN_layer7(x,layer_params,expanded_arr=[]):
  563. # Read chebyshev coefficients:
  564. cheby_coeffs= layer_params["W"]
  565. inputdim=cheby_coeffs.shape[0]
  566. outdim=cheby_coeffs.shape[1]
  567. # Normalize
  568. x = activation(x)
  569. # Reshape
  570. x = x.reshape((-1, inputdim, 1))
  571. x=jnp.stack((T0(x),
  572. T1(x),
  573. T2(x),
  574. T3(x),
  575. T4(x),
  576. T5(x),
  577. T6(x),
  578. T7(x)),axis=2)# Discard dummy dimension
  579. # Compute the Chebyshev interpolation
  580. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  581. # Remove extra dimension
  582. x= jnp.reshape(x, (outdim,))
  583. return x
  584. #Define params
  585. params_N = params["params"]
  586. #Normalize inputs
  587. x = norm_fn(X_in, M1, M2)
  588. for ly in range(len(params_N)):
  589. layer_params=params_N[ly]
  590. x=Cheby_KAN_layer7(x,layer_params)
  591. return x
  592. def KAN_Net9(params, X_in, M1, M2, activation, norm_fn):
  593. def Cheby_KAN_layer9(x,layer_params,expanded_arr=[]):
  594. # Read chebyshev coefficients:
  595. cheby_coeffs= layer_params["W"]
  596. inputdim=cheby_coeffs.shape[0]
  597. outdim=cheby_coeffs.shape[1]
  598. # Normalize
  599. x = activation(x)
  600. # Reshape
  601. x = x.reshape((-1, inputdim, 1))
  602. x=jnp.stack((T0(x),
  603. T1(x),
  604. T2(x),
  605. T3(x),
  606. T4(x),
  607. T5(x),
  608. T6(x),
  609. T7(x),
  610. T8(x),
  611. T9(x)),axis=2)# Discard dummy dimension
  612. # Compute the Chebyshev interpolation
  613. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  614. # Remove extra dimension
  615. x= jnp.reshape(x, (outdim,))
  616. return x
  617. #Define params
  618. params_N = params["params"]
  619. #Normalize inputs
  620. x = norm_fn(X_in, M1, M2)
  621. for ly in range(len(params_N)):
  622. layer_params=params_N[ly]
  623. x=Cheby_KAN_layer9(x,layer_params)
  624. return x
  625. def KAN_Net11(params, X_in, M1, M2, activation, norm_fn):
  626. def Cheby_KAN_layer11(x,layer_params,expanded_arr=[]):
  627. # Read chebyshev coefficients:
  628. cheby_coeffs= layer_params["W"]
  629. inputdim=cheby_coeffs.shape[0]
  630. outdim=cheby_coeffs.shape[1]
  631. # Normalize
  632. x = activation(x)
  633. # Reshape
  634. x = x.reshape((-1, inputdim, 1))
  635. x=jnp.stack((T0(x),
  636. T1(x),
  637. T2(x),
  638. T3(x),
  639. T4(x),
  640. T5(x),
  641. T6(x),
  642. T7(x),
  643. T8(x),
  644. T9(x),
  645. T10(x),
  646. T11(x)),axis=2)# Discard dummy dimension
  647. # Compute the Chebyshev interpolation
  648. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  649. # Remove extra dimension
  650. x= jnp.reshape(x, (outdim,))
  651. return x
  652. #Define params
  653. params_N = params["params"]
  654. #Normalize inputs
  655. x = norm_fn(X_in, M1, M2)
  656. for ly in range(len(params_N)):
  657. layer_params=params_N[ly]
  658. x=Cheby_KAN_layer11(x,layer_params)
  659. return x
  660. # WN KAN
  661. def WN_KAN_Net5(params, X_in, M1, M2, activation, norm_fn):
  662. def Cheby_KAN_layer5(x,layer_params,expanded_arr=[]):
  663. # Read chebyshev coefficients and other params:
  664. cheby_coeffs= layer_params["W"]
  665. g= layer_params["g"]
  666. b= layer_params["b"]
  667. # Normalize
  668. cheby_coeffs= cheby_coeffs / jnp.linalg.norm(cheby_coeffs, axis=0, keepdims=True)
  669. inputdim=cheby_coeffs.shape[0]
  670. outdim=cheby_coeffs.shape[1]
  671. # Normalize
  672. x = activation(x)
  673. # Reshape
  674. x = x.reshape((-1, inputdim, 1))
  675. x=jnp.stack((T0(x),
  676. T1(x),
  677. T2(x),
  678. T3(x),
  679. T4(x),
  680. T5(x)),axis=2)# Discard dummy dimension
  681. # Compute the Chebyshev interpolation
  682. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  683. # Remove extra dimension
  684. x= jnp.reshape(x, (outdim,))
  685. return g * x+b
  686. #Define params
  687. params_N = params["params"]
  688. #Normalize inputs
  689. x = norm_fn(X_in, M1, M2)
  690. for ly in range(len(params_N)):
  691. layer_params=params_N[ly]
  692. x=Cheby_KAN_layer5(x,layer_params)
  693. return x
  694. def KAN5_ResNet(params, X_in, M1, M2, activation, norm_fn):
  695. def Cheby_KAN_layer5(x,layer_params,W_key):
  696. # Read chebyshev coefficients:
  697. cheby_coeffs= layer_params[W_key]
  698. inputdim=cheby_coeffs.shape[0]
  699. outdim=cheby_coeffs.shape[1]
  700. # Normalize
  701. x = activation(x)
  702. # Reshape
  703. x = x.reshape((-1, inputdim, 1))
  704. x=jnp.stack((T0(x),
  705. T1(x),
  706. T2(x),
  707. T3(x),
  708. T4(x),
  709. T5(x)),axis=2)# Discard dummy dimension
  710. # Compute the Chebyshev interpolation
  711. x = jnp.einsum("bid0,iod->bo", x, cheby_coeffs) # shape = (batch_size, output_dim)
  712. # Remove extra dimension
  713. x= jnp.reshape(x, (outdim,))
  714. return x
  715. params_N = params["params"]
  716. inputs = norm_fn(X_in, M1, M2)
  717. layer_params=params_N[0]
  718. inputs=Cheby_KAN_layer5(inputs,layer_params,W_key="W")
  719. for ly in range(1,len(params_N[:-1])):
  720. layer_params=params_N[ly]
  721. g= Cheby_KAN_layer5(inputs,layer_params,W_key="W")
  722. h= Cheby_KAN_layer5(g,layer_params,W_key="W2")
  723. inputs= h+inputs
  724. layer_params=params_N[-1]
  725. inputs=Cheby_KAN_layer5(inputs,layer_params,W_key="W")
  726. return inputs

NNpp.py at commit f242e08, no license · at the source

Overview

  1. Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
  2. Department of Mechanical Engineering, University of Rochester, Rochester, NY 14627, USA
  3. School of Engineering, Brown University, Providence, RI 02912, USA
  4. Center for Translational Neuromedicine, University of Copenhagen, Copenhagen 2200, Denmark
Institutions: Brown University (United States); University of Rochester (United States); University of Copenhagen (Denmark)
Journal: Science advances, volume 12, issue 22, article eaeb0404
Dates: received 30 July 2025; accepted 20 April 2026; published online 27 May 2026; in print May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.aeb0404 · PMID 42202031 · PMCID PMC13215179 · OpenAlex W7162502040
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), systems (subfield)
Methods: Connectivity, fMRI & imaging, Physiology & signal measures
MeSH: Brain*, Magnetic Resonance Imaging*, Rheology*, Animals, Humans, Hydrodynamics (* major topic)
Journal subjects: Physical and Materials Sciences, Computational Biology, Machine Learning
Topic: Cerebrospinal fluid and hydrocephalus (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: NIH (R01AT012312, U19NS128613); Army Research Office (W911NF1910280); Air Force Office of Scientific Research (FA9550-20-1-0358); ONR (N00014-22-1-2795)
Citations: cited by 2 papers (Europe PMC); 85 references in the paper

Abstract

The circulation of cerebrospinal and interstitial fluid plays a vital role in clearing metabolic waste from the brain, and its disruption has been linked to neurological disorders. However, directly measuring brain-wide fluid transport, especially in the deep brain, has remained elusive. Here, we introduce magnetic resonance artificial intelligence velocimetry (MR-AIV), a framework featuring a specialized physics-informed architecture and optimization method that reconstructs three-dimensional fluid velocity fields from dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI). MR-AIV unveils brain-wide velocity maps while providing estimates of tissue permeability and pressure fields, quantities inaccessible to other methods. Applied to the brain, MR-AIV reveals a functional landscape of interstitial and perivascular flow, quantitatively distinguishing slow diffusion-driven transport [∼0.1 micrometers per second (μm/s)] from rapid advective flow (∼3 μm/s). This approach enables new investigations into brain clearance mechanisms and fluid dynamics in health and disease, with broad potential applications to other porous medium systems, from geophysics to tissue mechanics.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

jdtoscano94/MR-AIV

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: f242e08e8c4f1fb82f395e363b258d4259600cba, 31 May 2026
Languages: Python (12), Jupyter (2)
Size: 60 files, 14 scripts
Software Heritage: not archived
Found in: “Data, code, and materials availability:”
Holds: README, environment (pyproject.toml, setup.cfg), 2 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (11 files), JAX (10 files), h5py (9 files), Matplotlib (9 files), SciPy (9 files), pandas (7 files), scikit-learn (7 files), OpenCV (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
15 files

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 14 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. The full dataset used in this study (including original DCE-MRI data, synthetic validation data, and preprocessed training data), along with the complete MR-AIV codebase, is publicly archived on Zenodo at https://zenodo.org/records/15345392. Additionally, to ease accessibility and provide detailed usage instructions, the codebase has also been made available via a public GitHub repository at https://github.com/jdtoscano94/MR-AIV. This study did not generate new materials.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 6 MeSH terms, 4 funders, 69 references.

Cite

This paper

Toscano, J. D., Guo, Y., Wang, Z., Vaezi, M., Mori, Y., Karniadakis, G. E., Boster, K. A. S., & Kelley, D. H. (2026). MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI. Science advances, 12(22), eaeb0404. https://doi.org/10.1126/sciadv.aeb0404

BibTeX

@article{toscano2026mr,
author = {Toscano, Juan Diego and Guo, Yisen and Wang, Zhibo and Vaezi, Mohammad and Mori, Yuki and Karniadakis, George Em and Boster, Kimberly A. S. and Kelley, Douglas H.},
title = {{MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI}},
journal = {Science advances},
year = {2026},
month = may,
volume = {12},
number = {22},
pages = {eaeb0404},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.aeb0404},
url = {https://doi.org/10.1126/sciadv.aeb0404},
pmid = {42202031},
pmcid = {PMC13215179}
}

RIS

TY - JOUR
AU - Toscano, Juan Diego
AU - Guo, Yisen
AU - Wang, Zhibo
AU - Vaezi, Mohammad
AU - Mori, Yuki
AU - Karniadakis, George Em
AU - Boster, Kimberly A. S.
AU - Kelley, Douglas H.
TI - MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/05/27
VL - 12
IS - 22
SP - eaeb0404
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.aeb0404
UR - https://doi.org/10.1126/sciadv.aeb0404
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.aeb0404",
"type": "article-journal",
"title": "MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI",
"container-title": "Science advances",
"author": [
{
"family": "Toscano",
"given": "Juan Diego"
},
{
"family": "Guo",
"given": "Yisen"
},
{
"family": "Wang",
"given": "Zhibo"
},
{
"family": "Vaezi",
"given": "Mohammad"
},
{
"family": "Mori",
"given": "Yuki"
},
{
"family": "Karniadakis",
"given": "George Em"
},
{
"family": "Boster",
"given": "Kimberly A. S."
},
{
"family": "Kelley",
"given": "Douglas H."
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "22",
"page": "eaeb0404",
"DOI": "10.1126/sciadv.aeb0404",
"PMID": "42202031",
"PMCID": "PMC13215179",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.aeb0404",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
27
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1186/s12987-026-00782-w
Quantifying cerebrospinal fluid flow in pial perivascular spaces of rats.
Journal: Fluids and barriers of the CNS
In common: 11 references, author Douglas H. Kelley
[2] doi:10.1073/pnas.2526239123 [code]
Quantitative assessment of flow between cerebrospinal and interstitial fluid compartments in humans.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: structural MRI / diffusion, 8 references
[3] doi:10.1038/s41467-026-76306-9 [code]
Non-invasive characterization of perivascular subarachnoid spaces.
Journal: Nature communications
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, structural MRI / diffusion, 5 references
[4] doi:10.1073/pnas.2517059123 [code]
Stretch and flow at the gliovascular interface: High-fidelity modeling of astrocyte endfeet.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: pandas, SciPy, Matplotlib, 1 other tool, 5 references
[5] doi:10.1002/hipo.70131 [code]
Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences.
Journal: Hippocampus
In common: JAX, OpenCV, h5py, 5 other tools, systems
[6] doi:10.1038/s41593-026-02262-8 [code]
Cheese3D enables sensitive detection and analysis of whole-face movement in mice.
Journal: Nature neuroscience
In common: JAX, OpenCV, h5py, 5 other tools, systems
[7] doi:10.1093/brain/awag080 [code]
Early glymphatic failure in AppNL-F knock-in mice is linked to parenchymal border macrophages loss.
Journal: Brain : a journal of neurology
In common: 6 references
[8] doi:10.1002/mrm.70514 [code]
Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone.
Journal: Magnetic resonance in medicine
In common: pandas, SciPy, Matplotlib, 1 other tool, structural MRI / diffusion, 4 references
[9] doi:10.1177/0271678x261455452 [code]
Meningeal CSF transport varies across parasagittal dura subregions with age in humans.
Journal: Journal of cerebral blood flow and metabolism : official journal of the International Society of Cerebral Blood Flow and Metabolism
In common: structural MRI / diffusion, 5 references
[10] doi:10.1002/alz.71745 [code]
Regional astrocyte dysregulation and altered glymphatic-related markers in Alzheimer's disease frontal cortex.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: OpenCV, pandas, SciPy, 2 other tools, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.