mix-finished

phase-encoding-method-of-moments.py

@@ -3,62 +3,60 @@
 #from scipy.special import gamma, binom
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.special import gamma
+
 from mpmath import mp, mpf
-mp.dps = 1000
+mp.dps = 2000
 
-voxel_num = 5
-phase_range = np.pi/2
-phase_init = np.pi/4 #(0.0)
-U_points = voxel_num * 1000
+voxel_num = 15
+phase_range = mp.pi/2
+phase_init = mp.pi/4 #mpf(0.0)
+U_points = voxel_num * 3000
 
-noise_ratio = (0.0*1e-38) #1e8
+noise_ratio = mpf(0.0*1e-38) #1e8
 
-total_periods = 10
-rf_samples_per_period = 20
+total_periods=10#DO NOT CHANGE to fit T2s_scale autoscale
+rf_samples_per_period = 40
 # max polynomial order equals  rf_samples_per_period * total_periods 
 
 # B0=1.5T freq=64Mhz, period = 15.6 ns
-period = (10) #ms
-omega = 2.0*np.pi/period
+period = mpf(10) #ms
+omega = 2.0*mp.pi/period
 #T2s_scale = 0.01 #ms # need to be 10ms
-T2s_scale = total_periods*period/15 #ms # need to be 10ms
-T2s_min = T2s_scale/10.0
+T2s_scale = total_periods*period/rf_samples_per_period/voxel_num*5 #ms # need to be 10ms
+T2s_min = T2s_scale/10.0 
 #print(period)
 #ms
-time_steps = np.array(np.linspace((0), (period*total_periods), rf_samples_per_period*total_periods))
-tmp = [np.random.rand() for n in range(voxel_num)]
+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 = [np.random.rand() for n in range(voxel_num)]
+tmp = [mp.rand() for n in range(voxel_num)]
 voxel_T2s_decay = np.array(tmp)*(T2s_scale-2*T2s_min) + T2s_min
 print(voxel_T2s_decay)
 voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
 if voxel_num == 5:
-#    voxel_all = np.array([(0.2),(0.6),(0.02),(0.1)])
-    voxel_all = np.array([(0.822628),(0.691376),(0.282906),(0.226013),(0.90703),(0.144985),(0.228563),(0.340353),(0.462462),(0.720518)])
-    #voxel_all = np.array([(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.2),mpf(0.6),mpf(0.02),mpf(0.1)])
+    voxel_all = np.array([mpf(0.822628),mpf(0.691376),mpf(0.282906),mpf(0.226013),mpf(0.90703),mpf(0.144985),mpf(0.228563),mpf(0.340353),mpf(0.462462),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
 
-# a_i = np.array([(0.3),(0.1),(0.15),(0.1)])
-# d_i = np.array([(0.7),(0.9),(0.2),(0.67)])
+# a_i = np.array([mpf(0.3),mpf(0.1),mpf(0.15),mpf(0.1)])
+# d_i = np.array([mpf(0.7),mpf(0.9),mpf(0.2),mpf(0.67)])
 # voxel_num = len(a_i)
 
 
-voxel_phases = np.array(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
 
 ampl = []
-
-
 def gen_rf_signal(time):
     '''Generates demodulated signal at radio frequence using voxels
     amplitudes, T2s decays, and phases.
 
     '''
-    tmp = [(0.0) for n in range(len(time))]
+    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)):
@@ -71,17 +69,17 @@ def gen_rf_signal(time):
             #       "d_{:d} =".format(i),float(d_i[i]),"+", np.random.rand()*noise_ratio)
 
             amp = voxel_amplitudes[i] * (
-                np.exp(-time[t]/voxel_T2s_decay[i])
+                mp.exp(-time[t]/voxel_T2s_decay[i])
                 ) + ( 0.0 
                           # + np.random.rand()*noise_ratio
                         )
             if t == 0:
-                #print("a_{:d}".format(i),float(voxel_amplitudes[i]* np.sin(voxel_phases[i] + phase_init)))
-                print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
-            mag_sin[t] += amp * np.sin(
+                #print("a_{:d}".format(i),float(voxel_amplitudes[i]* mp.sin(voxel_phases[i] + phase_init)))
+                print("d_{:d}".format(i),float( mp.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
+            mag_sin[t] += amp * mp.sin(
                 voxel_phases[i] + phase_init
                 )
-            mag_cos[t] += amp * np.cos(
+            mag_cos[t] += amp * mp.cos(
                 voxel_phases[i] + phase_init
                 )
     return mag_sin, mag_cos
@@ -95,7 +93,7 @@ def binom(n,k):
 def shiftedLegendre(n):
     coeffs = []
     for k in range(n+1):
-        val = mpf(-1)**mpf(n) * binom(mpf(n),mpf(k)) * binom(mpf(n+k),mpf(k)) * mpf(-1)**mpf(k)
+        val = mpf(-1)**n * binom(mpf(n),mpf(k)) * binom(n+k,k) * (-1)**k
         coeffs.insert(0,val*mp.sqrt(2*n+1))
     return np.poly1d(coeffs)
 
@@ -110,12 +108,12 @@ def GetU (lambdas):
     P[1] = 2*x-1
     for n in range(1,n_max):
         P[n+1] = ((2*n+1)/(n+1)) * (2*x-1) * P[n]  - (n/(n+1))*P[n-1]
-
     U = np.zeros(U_points)
     for i in range (len(lambdas)):
         if i%10 == 0: print(i)
-        #polyL = L[i] #shiftedLegendre(i)        
         U += lambdas[i]*P[i]*np.sqrt(2*i+1)
+        #polyL = L[i] #shiftedLegendre(i)
+        #U += lambdas[i]*polyL(x)
     return U
 
 def GetLambda(mag_rf):
@@ -125,12 +123,12 @@ def GetLambda(mag_rf):
 #        print("M = ", i)
         lambd = (0)
         for j in range(i+1):
-            lambd += float(K(i,j))*mag_rf[j]
+            lambd += K(i,j)*mag_rf[j]
 #            print("K({:d},{:d}) =".format(i,j), float(K(i,j)))
-        all_lambda.append(lambd)
-    # tmp = [(0.0) for n in range(M_cutoff)]
+        all_lambda.append(float(lambd))
+    # tmp = [mpf(0.0) for n in range(M_cutoff)]
     # all_lambda = np.array(tmp)
-    # all_lambda[10] = (1.0)
+    # all_lambda[10] = mpf(1.0)
     return all_lambda
 
 
@@ -142,7 +140,7 @@ sign = ""
 #         sign+="\n"
 #     sign = sign + '{:3.2g}'.format(float(a_i[i]))+"/"+'{:3.2g}'.format(float(d_i[i]))+", "
 
-# #    print(np.exp(-1.0/voxel_T2s_decay[i]))
+# #    print(mp.exp(-1.0/voxel_T2s_decay[i]))
                
  
 plt.plot(mag_sin, ls=' ', marker='o')
@@ -159,14 +157,12 @@ for i in range(len(mag_sin)):
     L += [polyL]
 
 x = np.linspace(0,1, U_points)
-polyL_val = np.array([float(L[-1](x[n])) for n in range(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.")
+print("Before lambda.")
 
 lambdas = GetLambda(mag_sin)
 print(lambdas)
@@ -174,7 +170,7 @@ U = GetU(lambdas)
 
 x = np.linspace(0,1, U_points)
 mag_x = np.linspace(0,1, len(mag_sin))
-plt.plot(x,U)
+plt.plot(x,U, lw=0.2)
 plt.title(r"$\sum^M \lambda_i L_i(x)$", y=1.02)
 plt.savefig("plt.pdf")
 plt.clf()