From b88bd0dbd691a30b296995e34529d607de6540a6 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Wed, 12 Jun 2024 19:34:45 +0200 Subject: [PATCH] first drafts rj and tj --- python/rj_spectrum.py | 114 +++++++++++++++++++ python/rj_spectrum_chunk.py | 123 +++++++++++++++++++++ python/tetrahedral_spectrum.py | 145 ++++++++++++++++++++++++ python/tetrahedral_spectrum_chunk.py | 158 +++++++++++++++++++++++++++ 4 files changed, 540 insertions(+) create mode 100644 python/rj_spectrum.py create mode 100644 python/rj_spectrum_chunk.py create mode 100644 python/tetrahedral_spectrum.py create mode 100644 python/tetrahedral_spectrum_chunk.py diff --git a/python/rj_spectrum.py b/python/rj_spectrum.py new file mode 100644 index 0000000..68145bf --- /dev/null +++ b/python/rj_spectrum.py @@ -0,0 +1,114 @@ +from time import time + +import numpy as np +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +# spectral parameter +delta = 161e3 # in Hz +eta = 0 +lb = 2e3 # in Hz + +# correlation time +tau = [1e-7] # in s + +# acquisition parameter +acq_length = 4096 +dt = 1e-6 # in s +t_echo = [0, 5e-6, 10e-6, 20e-6, 50e-6, 100e-6, 200e-6] # all in s + +# derived parameter +t_acq = np.arange(acq_length) * dt +t_max = acq_length*dt + 2*max(t_echo) +dampening = np.exp(-lb * t_acq) + +# random number generator +seed = 1234 +rng = np.random.default_rng(seed) + +# number of random walkers +num_traj = 5000 + + +def omega_q(delta_: float, eta_: float, theta_: float, phi_: float) -> float: + cos_theta = np.cos(theta_) + sin_theta = np.sin(theta_) + return 2 * np.pi * delta_ * (3 * cos_theta * cos_theta - 1 + eta_ * sin_theta*sin_theta * np.cos(2*phi_)) + + +def new_orientation(delta_: float, eta_: float) -> float: + z_theta, z_phi = rng.random(2) + theta = np.arccos(1 - 2 * z_theta) + phi = 2 * np.pi * z_phi + + return omega_q(delta_, eta_, theta, phi) + + +for tau_i in tau: + print(f'\nStart for tau={tau_i}') + + timesignal = np.zeros((acq_length, len(t_echo))) + + start = time() + expected_jumps = round(t_max/tau_i) + if expected_jumps > 1e7: + print(f'Too many jumps to process, Skip {tau_i}s') + continue + + for i in range(num_traj): + t_passed = 0 + t = [0] + phase = [0] + accumulated_phase = 0 + + while t_passed < t_max: + # orientation until the next jump + current_omega = new_orientation(delta, eta) + + # time to next jump + t_wait = rng.exponential(tau_i) + + t_passed += t_wait + accumulated_phase += t_wait * current_omega + + t.append(t_passed) + phase.append(accumulated_phase) + + # convenient interpolation to get phase at arbitrary times + phase_interpol = interp1d(t, phase) + + for j, t_echo_j in enumerate(t_echo): + # effect of de-phasing and re-phasing + start_amp = -2 * phase_interpol(t_echo_j) + + # start of actual acquisition + timesignal[:, j] += np.cos(start_amp + phase_interpol(t_acq + 2*t_echo_j)) * dampening + + if (i+1) % 200 == 0: + elapsed = time()-start + print(f'Step {i+1} of {num_traj}', end=' - ') + total = num_traj * elapsed / (i+1) + print(f'elapsed: {elapsed:.2f}s - total: {total:.2f}s - remaining: {total-elapsed:.2f}s') + + timesignal /= num_traj + + # FT to spectrum + freq = np.fft.fftshift(np.fft.fftfreq(acq_length, dt)) + spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real + spec -= spec[0] + + t_echo_strings = list(map(str, t_echo)) + + # plot spectra + fig, ax = plt.subplots() + lines = ax.plot(freq, spec) + ax.set_title(f'RJ (tau = {tau_i}s)') + ax.legend(lines, t_echo_strings) + # plt.savefig(f'RJ_{tau_i}.png') + + # # save time signals and spectra + # np.savetxt(f'rj_spectrum_{tau_i}.dat', np.c_[freq, spec], header='f\t' + '\t'.join(t_echo_strings)) + # np.savetxt(f'rj_timesignal_{tau_i}.dat', np.c_[t_acq, timesignal], header='t\t' + '\t'.join(t_echo_strings)) + + plt.show() diff --git a/python/rj_spectrum_chunk.py b/python/rj_spectrum_chunk.py new file mode 100644 index 0000000..53077b7 --- /dev/null +++ b/python/rj_spectrum_chunk.py @@ -0,0 +1,123 @@ +from time import time + +import numpy as np +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +# spectral parameter +delta = 161e3 # in Hz +eta = 0 +lb = 5e3 # in Hz + +# correlation time +tau = [1e-5] # in s + +# acquisition parameter +acq_length = 4096 +dt = 1e-6 # in s +t_echo = [0, 5e-6, 10e-6, 20e-6, 50e-6, 100e-6, 200e-6] # all in s + +# derived parameter +t_acq = np.arange(acq_length) * dt +t_max = acq_length*dt + 2*max(t_echo) +dampening = np.exp(-lb * t_acq) + +# random number generator +seed = None +rng = np.random.default_rng(seed) + +# number of random walkers +num_traj = 50000 + + +def omega_q(delta_: float, eta_: float, theta_: float, phi_: float) -> np.ndarray: + cos_theta = np.cos(theta_) + sin_theta = np.sin(theta_) + return 2 * np.pi * delta_ * (3 * cos_theta * cos_theta - 1 + eta_ * sin_theta*sin_theta * np.cos(2*phi_)) + + +def new_orientation(delta_: float, eta_: float, size=1) -> np.ndarray: + z_theta, z_phi = rng.random((2, size)) + theta = np.arccos(1 - 2 * z_theta) + phi = 2 * np.pi * z_phi + + return omega_q(delta_, eta_, theta, phi) + + +def new_tau(size=1) -> np.ndarray: + return rng.exponential(tau_i, size=size) + + +for tau_i in tau: + print(f'\nStart for tau={tau_i}') + + timesignal = np.zeros((acq_length, len(t_echo))) + + start = time() + expected_jumps = round(t_max/tau_i) + if expected_jumps > 1e7: + print(f'Too many jumps to process, Skip {tau_i}s') + continue + + chunks = int(0.6 * t_max / tau_i) + 1 + print(f'Chunk size for trajectories: {chunks}') + + for i in range(num_traj): + t_passed = 0 + t = [0] + phase = [0] + accumulated_phase = 0 + + while t_passed < t_max: + # orientation until the next jump + current_omega = new_orientation(delta, eta, size=chunks) + + # time to next jump + t_wait = new_tau(size=chunks) + accumulated_phase = np.cumsum(t_wait*current_omega) + phase[-1] + + t_wait = np.cumsum(t_wait) + t_passed + t_passed = t_wait[-1] + t.extend(t_wait.tolist()) + + phase.extend(accumulated_phase.tolist()) + + # convenient interpolation to get phase at arbitrary times + phase_interpol = interp1d(t, phase) + + for j, t_echo_j in enumerate(t_echo): + # effect of de-phasing and re-phasing + start_amp = -2 * phase_interpol(t_echo_j) + + # start of actual acquisition + timesignal[:, j] += np.cos(start_amp + phase_interpol(t_acq + 2*t_echo_j)) * dampening + + if (i+1) % 200 == 0: + elapsed = time()-start + print(f'Step {i+1} of {num_traj}', end=' - ') + total = num_traj * elapsed / (i+1) + print(f'elapsed: {elapsed:.2f}s - total: {total:.2f}s - remaining: {total-elapsed:.2f}s') + + timesignal /= num_traj + + # FT to spectrum + freq = np.fft.fftshift(np.fft.fftfreq(acq_length, dt)) + spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real + spec -= spec[0] + # spec /= np.max(spec, axis=0) + + t_echo_strings = list(map(str, t_echo)) + + # plot spectra + fig, ax = plt.subplots() + lines = ax.plot(freq, spec) + ax.set_title(f'RJ (tau = {tau_i}s)') + ax.legend(lines, t_echo_strings) + # plt.savefig(f'RJ_{tau_i}.png') + + # # save time signals and spectra + # np.savetxt(f'rj_spectrum_{tau_i}_chunky.dat', np.c_[freq, spec], header='f\t' + '\t'.join(t_echo_strings)) + # np.savetxt(f'rj_timesignal_{tau_i}_chunky.dat', np.c_[t_acq, timesignal], header='t\t' + '\t'.join(t_echo_strings)) + + plt.show() diff --git a/python/tetrahedral_spectrum.py b/python/tetrahedral_spectrum.py new file mode 100644 index 0000000..dbe2599 --- /dev/null +++ b/python/tetrahedral_spectrum.py @@ -0,0 +1,145 @@ +from time import time + +import numpy as np +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +# spectral parameter +delta = 161e3 # in Hz +eta = 0 +lb = 10e3 # in Hz + +# correlation time +tau = [1e-5] # in s + +# acquisition parameter +acq_length = 4096 +dt = 1e-6 # in s +t_echo = [0, 5e-6, 10e-6, 20e-6, 50e-6, 100e-6, 200e-6] # all in s + +# derived parameter +t_acq = np.arange(acq_length) * dt +t_max = acq_length*dt + 2*max(t_echo) +dampening = np.exp(-lb * t_acq) + +# random number generator +seed = None +rng = np.random.default_rng(seed) + +# number of random walkers +num_traj = 10000 + + +def omega_q(delta_: float, eta_: float, theta_: float, phi_: float) -> np.ndarray: + cos_theta = np.cos(theta_) + sin_theta = np.sin(theta_) + return 2 * np.pi * delta_ * (3 * cos_theta * cos_theta - 1 + eta_ * sin_theta*sin_theta * np.cos(2*phi_)) + + +def rotate(x_in: float, y_in: float, z_in: float, a: float) -> tuple[float, float]: + # rotation by tetrahedral angle is a given, only second angle is free parameter + beta = 109.45 * np.pi / 180. + cos_beta = np.cos(beta) + sin_beta = np.sin(beta) + + cos_alpha = np.cos(a) + sin_alpha = np.sin(a) + + scale = np.sqrt(1 - z_in * z_in) + 1e-12 + + x = x_in * cos_beta + sin_beta / scale * (x_in * z_in * cos_alpha - y_in * sin_alpha) + y = y_in * cos_beta + sin_beta / scale * (y_in * z_in * cos_alpha + x_in * sin_alpha) + z = z_in * cos_beta - scale * sin_beta*cos_alpha + z = max(-1, min(1, z)) + + return np.arccos(z), np.arctan2(y, x) + + +for tau_i in tau: + print(f'\nStart for tau={tau_i}') + + timesignal = np.zeros((acq_length, len(t_echo))) + + start = time() + expected_jumps = round(t_max/tau_i) + + for i in range(num_traj): + # draw orientation + z_theta, z_phi, z_alpha = rng.random(3) + theta = np.arccos(1 - 2 * z_theta) + phi = 2 * np.pi * z_phi + alpha = 2 * np.pi * z_alpha + + # orientation in cartesian coordinates + x_start = np.sin(theta) * np.cos(phi) + y_start = np.sin(theta) * np.sin(phi) + z_start = np.cos(theta) + + # calculate orientation of tetrahedral edges and their frequencies + orientation = np.zeros(4) + orientation[0] = omega_q(delta, eta, theta, phi) + orientation[1] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha)) + orientation[2] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha + 2 * np.pi / 3)) + orientation[3] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha + 4 * np.pi / 3)) + + t_passed = 0 + t = [0] + phase = [0] + accumulated_phase = 0 + current_orientation = rng.choice([0, 1, 2, 3]) + + while t_passed < t_max: + # orientation until the next jump + # always jump to a different position i -> i + {1, 2, 3} + current_orientation += rng.choice([1, 2, 3]) + current_orientation %= 4 + + current_omega = orientation[current_orientation] + + # time to next jump + t_wait = rng.exponential(tau_i) + + t_passed += t_wait + accumulated_phase += t_wait * current_omega + + t.append(t_passed) + phase.append(accumulated_phase) + + # convenient interpolation to get phase at arbitrary times + phase_interpol = interp1d(t, phase) + + for j, t_echo_j in enumerate(t_echo): + # effect of de-phasing and re-phasing + start_amp = -2 * phase_interpol(t_echo_j) + + # start of actual acquisition + timesignal[:, j] += np.cos(start_amp + phase_interpol(t_acq + 2*t_echo_j)) * dampening + + if (i+1) % 200 == 0: + elapsed = time()-start + print(f'Step {i+1} of {num_traj}', end=' - ') + total = num_traj * elapsed / (i+1) + print(f'elapsed: {elapsed:.2f}s - total: {total:.2f}s - remaining: {total-elapsed:.2f}s') + + timesignal /= num_traj + + # FT to spectrum + freq = np.fft.fftshift(np.fft.fftfreq(acq_length, dt)) + spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real + spec -= spec[0] + + t_echo_strings = list(map(str, t_echo)) + + # plot spectra + fig, ax = plt.subplots() + lines = ax.plot(freq, spec) + ax.set_title(f'tau = {tau_i}s') + ax.legend(lines, t_echo_strings) + # plt.savefig(f'{tau_i}.png') + + # save time signals and spectra + # np.savetxt(f'spectrum_{tau_i}.dat', np.c_[freq, spec], header='f\t' + '\t'.join(t_echo_strings)) + # np.savetxt(f'timesignal_{tau_i}.dat', np.c_[t_acq, timesignal], header='t\t' + '\t'.join(t_echo_strings)) + + plt.show() diff --git a/python/tetrahedral_spectrum_chunk.py b/python/tetrahedral_spectrum_chunk.py new file mode 100644 index 0000000..9b9fafb --- /dev/null +++ b/python/tetrahedral_spectrum_chunk.py @@ -0,0 +1,158 @@ +from time import time +from numpy.typing import ArrayLike +import numpy as np +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +# spectral parameter +delta = 161e3 # in Hz +eta = 0 +lb = 5e3 # in Hz + +# correlation time +tau = [1e-6] # in s + +# acquisition parameter +acq_length = 4096 +dt = 1e-6 # in s +t_echo = [0, 5e-6, 10e-6, 20e-6, 50e-6, 100e-6, 200e-6] # all in s + +# derived parameter +t_acq = np.arange(acq_length) * dt +t_max = acq_length*dt + 2*max(t_echo) +dampening = np.exp(-lb * t_acq) + +# random number generator +seed = None +rng = np.random.default_rng(seed) + +# number of random walkers +num_traj = 50000 + + +def omega_q(delta_: float, eta_: float, theta_: ArrayLike, phi_: ArrayLike) -> np.ndarray: + cos_theta = np.cos(theta_) + sin_theta = np.sin(theta_) + return 2 * np.pi * delta_ * (3 * cos_theta * cos_theta - 1 + eta_ * sin_theta*sin_theta * np.cos(2*phi_)) + + +def rotate(x_in: float, y_in: float, z_in: float, a: float) -> tuple[float, float]: + # rotation by tetrahedral angle is a given, only second angle is free parameter + beta = 109.45 * np.pi / 180. + cos_beta = np.cos(beta) + sin_beta = np.sin(beta) + + cos_alpha = np.cos(a) + sin_alpha = np.sin(a) + + scale = np.sqrt(1 - z_in * z_in) + 1e-12 + + x = x_in * cos_beta + sin_beta / scale * (x_in * z_in * cos_alpha - y_in * sin_alpha) + y = y_in * cos_beta + sin_beta / scale * (y_in * z_in * cos_alpha + x_in * sin_alpha) + z = z_in * cos_beta - scale * sin_beta*cos_alpha + z = max(-1, min(1, z)) + + return np.arccos(z), np.arctan2(y, x) + + +def new_tau(size=1) -> np.ndarray: + return rng.exponential(tau_i, size=size) + + +for tau_i in tau: + print(f'\nStart for tau={tau_i}') + + timesignal = np.zeros((acq_length, len(t_echo))) + + start = time() + expected_jumps = round(t_max/tau_i) + if expected_jumps > 1e7: + print(f'Too many jumps to process, Skip {tau_i}s') + continue + + chunks = int(0.6 * t_max / tau_i) + 1 + print(f'Chunk size for trajectories: {chunks}') + + for i in range(num_traj): + # draw orientation + z_theta, z_phi, z_alpha = rng.random(3) + theta = np.arccos(1 - 2 * z_theta) + phi = 2 * np.pi * z_phi + alpha = 2 * np.pi * z_alpha + + # orientation in cartesian coordinates + x_start = np.sin(theta) * np.cos(phi) + y_start = np.sin(theta) * np.sin(phi) + z_start = np.cos(theta) + + # calculate orientation of tetrahedral edges and their frequencies + orientation = np.zeros(4) + orientation[0] = omega_q(delta, eta, theta, phi) + orientation[1] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha)) + orientation[2] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha + 2 * np.pi / 3)) + orientation[3] = omega_q(delta, eta, *rotate(x_start, y_start, z_start, alpha + 4 * np.pi / 3)) + + t_passed = 0 + t = [0] + phase = [0] + accumulated_phase = 0 + start_position = rng.choice([0, 1, 2, 3]) + + while t_passed < t_max: + # orientation until the next jump + jumps = rng.choice([1, 2, 3], size=chunks) + start_position + jumps = np.cumsum(jumps) + jumps %= 4 + + current_omega = orientation[jumps] + # current_omega = rng.choice(orientation, size=chunks) + + # time to next jump + t_wait = new_tau(size=chunks) + accumulated_phase = np.cumsum(t_wait*current_omega) + phase[-1] + + t_wait = np.cumsum(t_wait) + t_passed + t_passed = t_wait[-1] + t.extend(t_wait.tolist()) + + phase.extend(accumulated_phase.tolist()) + + # convenient interpolation to get phase at arbitrary times + phase_interpol = interp1d(t, phase) + + for j, t_echo_j in enumerate(t_echo): + # effect of de-phasing and re-phasing + start_amp = -2 * phase_interpol(t_echo_j) + + # start of actual acquisition + timesignal[:, j] += np.cos(start_amp + phase_interpol(t_acq + 2*t_echo_j)) * dampening + + if (i+1) % 200 == 0: + elapsed = time()-start + print(f'Step {i+1} of {num_traj}', end=' - ') + total = num_traj * elapsed / (i+1) + print(f'elapsed: {elapsed:.2f}s - total: {total:.2f}s - remaining: {total-elapsed:.2f}s') + + timesignal /= num_traj + + # FT to spectrum + freq = np.fft.fftshift(np.fft.fftfreq(acq_length, dt)) + spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real + spec -= spec[0] + # spec /= np.max(spec, axis=0) + + t_echo_strings = list(map(str, t_echo)) + + # plot spectra + fig, ax = plt.subplots() + lines = ax.plot(freq, spec) + ax.set_title(f'tau = {tau_i}s') + ax.legend(lines, t_echo_strings) + # plt.savefig(f'{tau_i}.png') + + # save time signals and spectra + # np.savetxt(f'spectrum_{tau_i}.dat', np.c_[freq, spec], header='f\t' + '\t'.join(t_echo_strings)) + # np.savetxt(f'timesignal_{tau_i}.dat', np.c_[t_acq, timesignal], header='t\t' + '\t'.join(t_echo_strings)) + + plt.show()