bgfs python

phase-encoding-simple.py

@@ -1,20 +1,22 @@
 #!/usr/bin/env python3
 # -*- coding: UTF-8 -*-
 from scipy.optimize import differential_evolution
+from scipy.optimize import minimize
 import numpy as np
-voxel_num = 3
+voxel_num = 5
 phase_range = np.pi/2
 phase_init = 0.0
 
-noise_ratio = 1e8
+noise_ratio = 0.0 #1e8
 
-total_periods = 1000
-rf_samples_per_period = 10
+total_periods = 50
+rf_samples_per_period = 1
 
 # B0=1.5T freq=64Mhz, period = 15.6 ns
 period = 15.6/1000/1000 #ms
 omega = 2.0*np.pi/period
-T2s_scale = 0.01 #ms # need to be 10ms
+#T2s_scale = 0.01 #ms # need to be 10ms
+T2s_scale = total_periods*period #ms # need to be 10ms
 T2s_min = T2s_scale/1000.0
 #print(period)
 #ms
@@ -22,9 +24,9 @@ time_steps = np.linspace(0, period*total_periods, rf_samples_per_period*total_pe
 
 voxel_amplitudes = np.random.rand(voxel_num)
 voxel_T2s_decay = np.random.rand(voxel_num)*(T2s_scale-T2s_min) + T2s_min
-print("amp/decay", voxel_amplitudes,voxel_T2s_decay)
 voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
-print("all   ", voxel_all)
+if voxel_num == 5:
+    voxel_all = [0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353,0.662462,0.720518]
 #voxel_amplitudes = [0.4, 0.8, 0]
 #voxel_amplitudes = [0.9, 0.092893218813452, 0.5]
 #voxel_amplitudes = [0.6, 0.517157287525381, 0.4]
@@ -39,41 +41,65 @@ if len(voxel_amplitudes) != len(voxel_phases):
 
 
 def gen_rf_signal(time):
-    signal = 0.0
+    mag_sin = 0.0
+    mag_cos = 0.0
     for i in range(voxel_num):
-        signal += voxel_amplitudes[i]*np.sin(
-            omega*time + voxel_phases[i] + phase_init
-            ) * (
+        amp = voxel_amplitudes[i] * (
                 np.exp(-time/voxel_T2s_decay[i])
                 ) + ( 0.0 
-                           + np.random.rand()/noise_ratio
+                           # + np.random.rand()*noise_ratio
                     )
-    return signal
+        mag_sin += amp * np.sin(
+            voxel_phases[i] + phase_init
+            )
+        mag_cos += amp * np.cos(
+            voxel_phases[i] + phase_init
+            )
+    return mag_sin, mag_cos
 
 def assumed_signal(time, values):
     amplitudes = values[:voxel_num]
     T2s_decay = values[voxel_num:]*T2s_scale
-    signal = 0.0
+    mag_sin = 0.0
+    mag_cos = 0.0
     for i in range(voxel_num):
-        signal += amplitudes[i]*np.sin(
-            omega*time + voxel_phases[i] + phase_init
-            ) * (
+        amp = amplitudes[i] * (
                 np.exp(-time/T2s_decay[i])
-                )
-    return signal
-
-rf_signal_measured = gen_rf_signal(time_steps)
-
-def fitness(amplitudes):
-    diff = rf_signal_measured - assumed_signal(time_steps, amplitudes)
-    return np.sqrt(np.mean(np.square(diff)))
-
+                ) 
+        mag_sin += amp * np.sin(
+            voxel_phases[i] + phase_init
+            )
+        mag_cos += amp * np.cos(
+            voxel_phases[i] + phase_init
+            )
+    return mag_sin, mag_cos
+
+
+mag_sin, mag_cos = gen_rf_signal(time_steps)
+
+def fitness(amplitudes, *args):
+    assumed_sin, assumed_cos = assumed_signal(time_steps, amplitudes)
+    diff = 0
+    #amp = np.sqrt(mag_sin**2 + mag_cos**2)
+    # diff += np.sqrt(np.mean(np.square((mag_sin - assumed_sin)/amp)))
+    # diff += np.sqrt(np.mean(np.square((mag_cos - assumed_cos)/amp)))
+    diff += np.sqrt(np.mean(np.square((mag_sin - assumed_sin))))
+    diff += np.sqrt(np.mean(np.square((mag_cos - assumed_cos))))
+    diff += np.sqrt(np.mean(np.square(mag_sin/mag_cos - assumed_sin/assumed_cos)))
+    return diff
+
+def hyper_fit(amplitudes, *args):
+    result = minimize(fitness, amplitudes, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': False})
+    print (result.fun)
+    return result.fun
 #print(voxel_phases)
 #print (voxel_amplitudes)
 
-import matplotlib.pyplot as plt
-plt.plot(time_steps, rf_signal_measured)
-plt.show()
+# import matplotlib.pyplot as plt
+# plt.plot(time_steps, mag_sin)
+# plt.plot(time_steps, mag_cos)
+# plt.show()
+
 # #print(fitness(test_amplitudes))
 
 bounds = []
@@ -83,11 +109,33 @@ for i in range(voxel_num):
     bounds.append(amplitude_minmax)
 for i in range(voxel_num):
     bounds.append(T2s_minmax)
-result = differential_evolution(fitness, bounds, polish=True
-                                    #, maxiter = voxel_num*2*500
-                                    )
+
+x0 = np.full(2*voxel_num,0.5)
+
+result = minimize(hyper_fit, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': True})
+# result = differential_evolution(hyper_fit, bounds, polish=True
+#                                     , maxiter = voxel_num*2*500
+#                                 , disp = True
+# )
+
+# bestresult = minimize(fitness, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': True})
+# for x in range(50):
+#     print("Try to solve: ", x)
+#     x0 = np.random.rand(2*voxel_num)
+#     result = minimize(fitness, x0, bounds = bounds, method='L-BFGS-B', options={'ftol': 1e-16,'gtol': 1e-16, 'disp': False})
+#     if result.fun < bestresult.fun:
+#         bestresult = result        
+#         print("new best: ", bestresult.fun)
+# result = bestresult
+
+#result = minimize(fitness, x0, bounds = bounds, method='COBYLA', tol=1e-16, options={ 'maxiter':20000,'disp': True})
+#result = minimize(fitness, result.x, bounds = bounds, method='L-BFGS-B', options={'gtol': 1e-16, 'disp': True})
+#result = minimize(hyper_fit, x0, bounds = bounds, method='Nelder-Mead', tol=1e-16,  options={'maxiter':100,'disp': True})
+
 #result.x[voxel_num:] = result.x[voxel_num:]/T2s_scale
-print("eval  ",result.x, result.fun)
+print("amp/decay", voxel_amplitudes,voxel_T2s_decay)
+print("all   ", voxel_all)
+print("eval  ",result.x, "\n=====> fun=",result.fun)