optimize can find the result (sometimes)

phase-encoding-method-of-moments.py

@@ -3,19 +3,23 @@
 #from scipy.special import gamma, binom
 import numpy as np
 import matplotlib.pyplot as plt
+import scipy.sparse as sparse
+from scipy.sparse.linalg import eigsh
+from numpy.linalg import eig
 
 from mpmath import mp, mpf
 mp.dps = 2000
 
-voxel_num = 5
+voxel_num = 10
 phase_range = np.pi/2 *0.0
 phase_init = np.pi/2  #(0.0)
 U_points = voxel_num * 3000
 
-noise_ratio = (1e-60) #1e8
+noise_ratio = (1.2e-16) #(1.2e-16) #1e8
 
 total_periods=10#DO NOT CHANGE to fit T2s_scale autoscale
-rf_samples_per_period = 10
+rf_samples_per_period = 4
+
 # max polynomial order equals  rf_samples_per_period * total_periods 
 
 # B0=1.5T freq=64Mhz, period = 15.6 ns
@@ -25,7 +29,9 @@ omega = 2.0*np.pi/period
 T2s_scale = total_periods*period/rf_samples_per_period/voxel_num*8 #ms # need to be 10ms
 T2s_min = T2s_scale/20.0 
 #ms
-time_steps = np.array(mp.linspace((0), (period*total_periods), rf_samples_per_period*total_periods))
+total_steps = rf_samples_per_period*total_periods
+if total_steps % 2 != 0: total_steps -= 1
+time_steps = np.array(mp.linspace((0), (period*total_periods), total_steps))
 tmp = [np.random.rand()*0.0+(1.0) for n in range(voxel_num)]
 voxel_amplitudes = np.array(tmp)
 tmp = [np.random.rand() for n in range(voxel_num)]
@@ -73,7 +79,7 @@ def gen_rf_signal(time):
                 #print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
             mag_sin[t] += amp * np.sin((voxel_phases[i] + phase_init))
             mag_cos[t] += amp * np.cos((voxel_phases[i] + phase_init))
-        mag_sin[t] += np.random.rand()*noise_ratio
+        mag_sin[t] *= (1 + (np.random.rand()-0.5)*noise_ratio)
         #mag_cos[t] += np.random.rand()*noise_ratio
     return mag_sin, mag_cos
 
@@ -83,12 +89,6 @@ def factorial(n):
 def binom(n,k):
     return factorial(n)/(factorial(k)*factorial(n-k))
 
-def shiftedLegendre(n):
-    coeffs = []
-    for k in range(n+1):
-        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)
 
 def K ( i, j):
     polyL = L[i] #precomputed shiftedLegendre(i)
@@ -123,31 +123,154 @@ def GetLambda(mag_rf):
     # all_lambda[10] = mpf(1.0)
     return all_lambda
 
+def GetAbsLastLambda(mag_rf):
+    M_cutoff = len(mag_rf)
+    all_lambda = []
+    i = M_cutoff - 1
+    lambd = mpf(0.0)
+    for j in range(i+1):
+        lambd += K(i,j)*mpf(mag_rf[j])
+    #    print("K({:d},{:d}) =".format(i,j), float(K(i,j)))
+    all_lambda.append(float(lambd))
+        # tmp = [mpf(0.0) for n in range(M_cutoff)]
+    # all_lambda = np.array(tmp)
+    # all_lambda[10] = mpf(1.0)
+    return np.abs(all_lambda[0])
+
+
+def shiftedLegendre(n):
+    coeffs = []
+    for k in range(n+1):
+        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)
 
 def GetGnk(c,n,k):
-    Gnk = 0
+    Gnk = mpf(0.0)
     for m in range(k+1):
-        Gnk+=(-1)**m * binom(k,m) * c[n+m]
+        Gnk+=(-1)**m * binom(mpf(k),mpf(m)) * c[n+m]
     return Gnk
 
-    
+def GetGnkSum(c):
+    total = 0
+    N = len(c)
+    for n in range(2):
+        for k in range (N-n-2,N-n):
+            Gnk = float(GetGnk(c , n, k))
+            if Gnk <0:
+                total += Gnk
+    return total
+
 
-mag_sin, mag_cos = gen_rf_signal(time_steps)
 
-for i in range(3):
-    print("c",i," = ",float(mag_sin[i]))
+def FixGnk(mag_sin):
+    total = GetGnkSum(mag_sin)
+    print(total)
+
+    
+def mat1mp (mag_sin):
+    N = len(mag_sin)
+    NN = int((N-1)/2)
+    tmp = np.array([mpf(0.0) for n in range((NN+1)**2)])
+    tmp2 = np.reshape(tmp,( NN+1 , NN+1 ))
+    mat = mp.matrix(tmp2)
+    # mat = np.zeros(( NN+1 , NN+1 ))
+    for i in range(NN+1):
+        for j in range(NN+1):
+            mat[i,j] = mag_sin[i+j+1]
+    # print(mat)
+    return mat
+def mat1 (mag_sin):
+    N = len(mag_sin)
+    NN = int((N-1)/2)
+    mat = np.zeros(( NN+1 , NN+1 ))
+    for i in range(NN+1):
+        for j in range(NN+1):
+            mat[i,j] = mag_sin[i+j+1]
+    # print(mat)
+    return mat
+
+def mat2 (mag_sin):
+    N = len(mag_sin)
+    NN = int((N-1)/2)
+    mat = np.zeros(( NN+1 , NN+1 ))
+    for i in range(NN+1):
+        for j in range(NN+1):
+            mat[i,j] = mag_sin[i+j] - mag_sin[i+j+1]
+    return mat
+
+
+def fitness(corrector):
+    mag_sin = mag_sin_0*(1 + mpf(1.1)*(corrector-0.5)*noise_ratio)
+    res = GetAbsLastLambda(mag_sin)
+    global globi
+    globi = globi + 1    
+    if globi % 100 == 0: print("REsult", res)
+    return res
+
+
+mag_sin_0, mag_cos_0 = gen_rf_signal(time_steps)
+mag_sin = np.copy(mag_sin_0)
+mag_cos = np.copy(mag_cos_0)
+
+L = [] # Shifted Legendre polinomials
+for i in range(len(mag_sin)):
+    polyL = shiftedLegendre(i)
+    L += [polyL]
 
-for n in range(3):
-    for k in range (3):
-        print("G",n,k,float(GetGnk(mag_sin, n, k)))
+globi = 1
+
+import time
+time1 = time.time()
+corr = np.zeros(len(mag_sin))
+fitness(corr)
+from scipy.optimize import minimize
+# = minimize(fitness,corr,method='L-BFGS-B')
+v = minimize(fitness,corr,method='TNC',options={'disp': False, 'minfev': 0, 'scale': None, 'rescale': -1, 'offset': None, 'gtol': -1, 'eps': 1e-08, 'eta': -1, 'maxiter': None, 'maxCGit': -1, 'mesg_num': None, 'ftol': -1, 'xtol': -1, 'stepmx': 0, 'accuracy': 0})
+
+# from scipy.optimize import differential_evolution
+# bounds = []
+# for i in range(len(mag_sin)):
+#     bounds.append((0,1))
+# v = differential_evolution(fitness,bounds)
+
+fitness(v.x)
+print(v.fun, "at", v.x)
+time2 = time.time()
+print( 'it took %0.5f s' % ((time2-time1)))
+
+print("N = ",len(mag_sin))
+print("Noise = ", noise_ratio )
+# print("c_n =",["%0.8f"%(x) for x in mag_sin])
+# print("mat1 =",mat1(mag_sin))
+# print("mat2 =",mat2(mag_sin))
+
+#print( 'eigsh took %0.5f s' % ((time2-time1)))
+#print("eigsh1 = ",v1,"eigsh2 = ",v2)
+# v1f = eig(mat1(mag_sin))
+# v2f = eig(mat2(mag_sin))
+# print("eig1 = ",np.sort(v1f[0]))
+# print("eig2 = ",np.sort(v2f[0]))
+
+# v = eigsh(mat2(mag_sin), k=6, return_eigenvectors=False, which='LA')
+# print(v)
+
+N = len(mag_sin)
+print("N = ", N)
+for n in range(2):
+    for k in range (N-n-2,N-n):
+        Gnk = float(GetGnk(mag_sin, n, k))
+        if Gnk <0:
+            print("G",n,k," = ", Gnk)
         
+#FixGnk(mag_sin)
 
 sign = ""
 # for i in range(voxel_num):
 #     if i%5 == 0:
 #         sign+="\n"
 #     sign = sign + '{:3.2g}'.format(float(a_i[i]))+"/"+'{:3.2g}'.format(float(d_i[i]))+", "
-
+# 
 # #    print(mp.exp(-1.0/voxel_T2s_decay[i]))
                
  
@@ -157,13 +280,10 @@ plt.savefig("signal.pdf")
 plt.clf()
 
 
-L = [] # Shifted Legendre polinomials
-for i in range(len(mag_sin)):
-    polyL = shiftedLegendre(i)
-    # print("i=",i,"  L_i:")
-    # print(polyL)
-    L += [polyL]
 
+recL = []
+
+    
 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)
@@ -172,8 +292,12 @@ x = np.linspace(0,1, U_points)
 # plt.clf()
 #print("Before lambda.")
 
+import time
+time1 = time.time()
 lambdas = GetLambda(mag_sin)
-#print(lambdas)
+time2 = time.time()
+print( 'GetLambda took %0.5f s' % ((time2-time1)))
+#print((lambdas))
 U = GetU(lambdas)
 
 sign ="" 
@@ -192,4 +316,8 @@ plt.xlabel(sign, fontsize=28)
 plt.savefig("plt.pdf")
 plt.clf()
 
-
+# rms = np.sqrt(np.mean(U**2))/np.max(U)
+# #print(rms)
+# status = "OK"
+# if rms < 0.2: status = "BAD"
+# print("rms =", rms, status)