clean up

phase-encoding-method-of-moments.py

@@ -1,154 +1,153 @@
 #!/usr/bin/env python3
 # -*- coding: UTF-8 -*-
-from scipy.optimize import differential_evolution
-from scipy.optimize import minimize
-from scipy.special import gamma, binom
+#from scipy.special import gamma, binom
 import numpy as np
+import matplotlib.pyplot as plt
+
+from mpmath import mp, mpf
+mp.dps = 1000
+
 voxel_num = 2
-phase_range = np.pi/2
-phase_init = 0.0
+phase_range = mp.pi/2
+phase_init = mp.pi/20 #mpf(0.0)
 U_points = voxel_num * 1000
 
-noise_ratio = 0.0 #1e8
+# noise_ratio = mpf(0.0) #1e8
 
-total_periods = 8
-rf_samples_per_period = 16
+total_periods = 10
+rf_samples_per_period = 5
+# max polynomial order equals  rf_samples_per_period * total_periods 
 
 # B0=1.5T freq=64Mhz, period = 15.6 ns
-period = 15.6/1000/1000 #ms
-omega = 2.0*np.pi/period
+period = mpf(15.6/1000/1000) #ms
+omega = 2.0*mp.pi/period
 #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
-time_steps = np.linspace(0, period*total_periods, rf_samples_per_period*total_periods)
-voxel_amplitudes = np.random.rand(voxel_num)
-voxel_T2s_decay = np.random.rand(voxel_num)*(T2s_scale-T2s_min) + T2s_min
+time_steps = np.array(mp.linspace(mpf(0), mpf(period*total_periods), rf_samples_per_period*total_periods))
+tmp = [mp.rand() for n in range(voxel_num)]
+voxel_amplitudes = np.array(tmp)
+tmp = [mp.rand() for n in range(voxel_num)]
+voxel_T2s_decay = np.array(tmp)*(T2s_scale-T2s_min) + T2s_min
 voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
 if voxel_num == 5:
-    voxel_all = np.array([0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353,0.662462,0.720518])
-    #voxel_all = [0.592606,0.135168,0.365712,0.667536,0.437378,0.918822,0.943879,0.590338,0.685997,0.658839]
+    voxel_all = np.array([mpf(0.822628),mpf(0.691376),mpf(0.282906),mpf(0.226013),mpf(0.90703),mpf(0.144985),mpf(0.328563),mpf(0.440353),mpf(0.662462),mpf(0.720518)])
+    #voxel_all = np.array([mpf(0.592606),mpf(0.135168),mpf(0.365712),mpf(0.667536),mpf(0.437378),mpf(0.918822),mpf(0.943879),mpf(0.590338),mpf(0.685997),mpf(0.658839)])
     voxel_amplitudes = voxel_all[:voxel_num]
     voxel_T2s_decay = voxel_all[voxel_num:]*T2s_scale
-#first estimate    0.551777 0.190833 0.271438 0.814036 0.347389 0.926153 0.908453 0.581414 0.666012 0.673226
 
-#voxel_amplitudes = [0.4, 0.8, 0]
-#voxel_amplitudes = [0.9, 0.092893218813452, 0.5]
-#voxel_amplitudes = [0.6, 0.517157287525381, 0.4]
-
-test_amplitudes = np.zeros(voxel_num)
-test_amplitudes = voxel_amplitudes
-voxel_phases = np.linspace(0,phase_range, voxel_num)
+voxel_phases = np.array(mp.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
 
-
 def gen_rf_signal(time):
-    mag_sin = 0.0
-    mag_cos = 0.0
-    for i in range(voxel_num):
-        amp = voxel_amplitudes[i] * (
-                np.exp(-time/voxel_T2s_decay[i])
+    '''Generates demodulated signal at radio frequence using voxels
+    amplitudes, T2s decays, and phases.
+
+    '''
+    tmp = [mpf(0.0) for n in range(len(time))]
+    mag_sin = np.array(tmp)
+    mag_cos = np.array(tmp)
+    for t in range(len(time)):
+        for i in range(voxel_num):
+            amp = voxel_amplitudes[i] * (
+                mp.exp(-time[t]/voxel_T2s_decay[i])
                 ) + ( 0.0 
-                           # + np.random.rand()*noise_ratio
-                    )
-        mag_sin += amp * np.sin(
-            voxel_phases[i] + phase_init
-            )
-        mag_cos += amp * np.cos(
-            voxel_phases[i] + phase_init
-            )
+                          # + np.random.rand()*noise_ratio
+                        )
+            mag_sin[t] += amp * mp.sin(
+                voxel_phases[i] + phase_init
+                )
+            mag_cos[t] += amp * mp.cos(
+                voxel_phases[i] + phase_init
+                )
     return mag_sin, mag_cos
 
 def factorial(n):
-    return gamma(n+1)
+    return mp.gamma(n+1)
+
+def binom(n,k):
+    return factorial(n)/(factorial(k)*factorial(n-k))
 
 def shiftedLegendre(n):
     coeffs = []
     for k in range(n+1):
-        val = (-1)**n * binom(n,k) * binom(n+k,k) * (-1)**k
+        val = mpf(-1)**n * binom(mpf(n),mpf(k)) * binom(n+k,k) * (-1)**k
         coeffs.insert(0,val)
     return np.poly1d(coeffs)
 
 def K ( i, j):
-    polyL = L[i] #shiftedLegendre(i)
+    polyL = L[i] #precomputed shiftedLegendre(i)
     return polyL.coeffs[-j-1]
 
+def GetU (lambdas):
+    x = np.array(mp.linspace(0,1, U_points))
+    tmp = [mpf(0.0) for n in range(len(lambdas))]
+    mag_sin = np.array(tmp)
+    tmp = [mpf(0.0) for n in range(U_points)]
+    U = np.array(tmp)
+    for i in range (len(lambdas)):
+        polyL = L[i] #shiftedLegendre(i)        
+        U += lambdas[i]*polyL(x)
+    return U
+
 def GetLambda(mag_rf):
     M_cutoff = len(mag_rf)
     all_lambda = []
     for i in range(M_cutoff):
-        lambd = 0
+        lambd = mpf(0)
         for j in range(i+1):
             lambd += K(i,j)*mag_rf[i]
         all_lambda.append(lambd)
-    all_lambda = np.zeros(M_cutoff)
-    all_lambda[19] = 1
+    # tmp = [mpf(0.0) for n in range(M_cutoff)]
+    # all_lambda = np.array(tmp)
+    # all_lambda[29] = mpf(1.0)
     return all_lambda
 
-def GetU (lambdas):
-    x = np.linspace(0,1, U_points)
-    U = np.zeros(U_points)
-    for i in range (len(lambdas)):
-        polyL = L[i] #shiftedLegendre(i)        
-        U += lambdas[i]*polyL(x)
-    return U
 
+mag_sin, mag_cos = gen_rf_signal(time_steps)
+
+sign = ""
+for i in range(voxel_num):
+    if i%5 == 0:
+        sign+="\n"
+    sign = sign + '{:3.2g}'.format(float(voxel_amplitudes[i] * mp.sin(
+        voxel_phases[i] + phase_init
+        )))+"/"+'{:3.2g}'.format(float(voxel_T2s_decay[i]))+", "
+
+
+plt.plot(mag_sin, ls='-', marker='o')
+plt.title("Signal to restore amp/decay_T:"+sign)
+plt.savefig("signal.pdf")
+plt.clf()
 
 
-mag_sin, mag_cos = gen_rf_signal(time_steps)
 L = [] # Shifted Legendre polinomials
 for i in range(len(mag_sin)):
     polyL = shiftedLegendre(i)        
     L += [polyL]
-#print(len(L))
-print(L[20])
+
+x = np.linspace(0,1, U_points)
+polyL_val = np.array([float(L[-1](x[n])) for n in range(U_points)])
+
+
+plt.plot(x,polyL_val)
+plt.title("Legendre polynom of order "+str(len(L)))
+plt.savefig("polyL.pdf")
+plt.clf()
+print("Output of last poly done.")
+
 lambdas = GetLambda(mag_sin)
+print(len(lambdas))
 U = GetU(lambdas)
 
-import matplotlib.pyplot as plt
 x = np.linspace(0,1, U_points)
 mag_x = np.linspace(0,1, len(mag_sin))
-# crop = 1
-# plt.plot(x[:-crop],U[:-crop])
 plt.plot(x,U)
-# plt.plot(mag_x,mag_sin)
-#plt.xlim(0.2,0.8)
-#plt.ylim(0.0,2.8)
 plt.savefig("plt.pdf")
-#plt.show()
-
-#print(voxel_phases)
-#print (voxel_amplitudes)
-
-# import matplotlib.pyplot as plt
-# plt.plot(time_steps, mag_sin)
-# plt.plot(time_steps, mag_cos)
-# plt.show()
-
-# #print(fitness(test_amplitudes))
-
-amplitude_minmax = (0,1)
-T2s_minmax = (T2s_min/T2s_scale,1)
-
-x0 = np.full(2*voxel_num,0.5)
-x0 = [0.551777,0.190833,0.271438,0.814036,0.347389,0.926153,0.908453,0.581414,0.666012,0.673226]
-
-
-# #result.x[voxel_num:] = result.x[voxel_num:]/T2s_scale
-# print("amp/decay", voxel_amplitudes,voxel_T2s_decay)
-# print("all   ", voxel_all)
-# print("eval  ",result.x, "\n=====> fun=",result.fun)
-
-
+plt.clf()
 
-# # print("Diff")
-# # print((voxel_amplitudes-result.x))
-# # print("percent")
-# print("percent",np.abs(voxel_all-result.x)*100)
-# if np.max(np.abs(voxel_all[:voxel_num]-result.x[:voxel_num])*100)>0.5:
-#     print ("==============  !!!LARGE!!! ===============")
 
-# print("\n")