|
|
@@ -0,0 +1,60 @@
|
|
|
+#!/usr/bin/env python3
|
|
|
+# -*- coding: UTF-8 -*-
|
|
|
+from scipy.optimize import differential_evolution
|
|
|
+import numpy as np
|
|
|
+voxel_num = 5
|
|
|
+phase_range = np.pi
|
|
|
+phase_init = 0.0
|
|
|
+
|
|
|
+total_periods = 50
|
|
|
+rf_samples_per_period = 10
|
|
|
+
|
|
|
+voxel_amplitudes = np.random.rand(voxel_num)
|
|
|
+test_amplitudes = np.zeros(voxel_num)
|
|
|
+test_amplitudes = voxel_amplitudes
|
|
|
+voxel_phases = np.linspace(0,phase_range, voxel_num)
|
|
|
+if len(voxel_amplitudes) != len(voxel_phases):
|
|
|
+ print("ERROR! Size of amplitude and phase arrays do not match!")
|
|
|
+ raise
|
|
|
+
|
|
|
+omega = 2*np.pi
|
|
|
+period = 2*np.pi/omega
|
|
|
+time_steps = np.linspace(0, period*total_periods, rf_samples_per_period*total_periods)
|
|
|
+
|
|
|
+def rf_signal(time):
|
|
|
+ signal = 0.0
|
|
|
+ for i in range(voxel_num):
|
|
|
+ signal += voxel_amplitudes[i]*np.sin(
|
|
|
+ omega*time + voxel_phases[i] + phase_init
|
|
|
+ )
|
|
|
+ return signal
|
|
|
+
|
|
|
+def assumed_signal(time, amplitudes):
|
|
|
+ signal = 0.0
|
|
|
+ for i in range(voxel_num):
|
|
|
+ signal += amplitudes[i]*np.sin(
|
|
|
+ omega*time + voxel_phases[i] + phase_init
|
|
|
+ )
|
|
|
+ return signal
|
|
|
+
|
|
|
+def fitness(amplitudes):
|
|
|
+ diff = rf_signal(time_steps) - assumed_signal(time_steps, amplitudes)
|
|
|
+ return np.sqrt(np.mean(np.square(diff)))
|
|
|
+
|
|
|
+#print(voxel_phases)
|
|
|
+print (voxel_amplitudes)
|
|
|
+# import matplotlib.pyplot as plt
|
|
|
+# plt.plot(time_steps, rf_signal(time_steps))
|
|
|
+# plt.show()
|
|
|
+#print(fitness(test_amplitudes))
|
|
|
+
|
|
|
+bounds = []
|
|
|
+minmax = (0,1)
|
|
|
+for i in range(voxel_num):
|
|
|
+ bounds.append(minmax)
|
|
|
+result = differential_evolution(fitness, bounds, polish=True)
|
|
|
+print(result.x, result.fun)
|
|
|
+print("Diff")
|
|
|
+print((voxel_amplitudes-result.x))
|
|
|
+print("percent")
|
|
|
+print((voxel_amplitudes-result.x)*100)
|
