d-params fit works

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/eval_spectra.py

@@ -1,33 +1,49 @@
+from ast import Constant
 from numba import njit, prange, complex128, float64
 import numpy as np
 
 
+# omega_init = True
+params_count = 4
+
+
 @njit(complex128(float64, float64[:]),
       fastmath=True, cache=True)
 def lorentzian(omega, xvec):
+    # global omega_init
+    pc = params_count
     res = 0
-    poles = len(xvec)//4
+    poles = len(xvec)//params_count
     factor = 1
     for i in range(poles):
-        gamma = xvec[i+0] / (factor**i)
-        omega_n = xvec[i+1] / (factor**i)
-        f = (xvec[i+2] + 1j*xvec[i+3]) / (factor**i)
+        gamma = xvec[pc*i+0] / (factor**i)
+        omega_n = xvec[pc*i+1] / (factor**i)
+        f = (xvec[pc*i+2] + 1j*xvec[pc*i+3]) / (factor**i)
+        # if (omega_init):
+        #     print('init', gamma, omega_n, f)
         if (np.abs(f) == 0):
             return res
         res = res + f / (omega * (omega + 1j*gamma) - omega_n**2)
+    # omega_init = False
     return res
 
 
-@njit(complex128[:](float64[:], float64[:]),
-      parallel=True, fastmath=True, cache=True)
+# @njit(complex128[:](float64[:], float64[:]),
+#       parallel=True, fastmath=True, cache=True)
 def spectra(omega, xvec):
-    # print('in eval', xvec)
-    poles = len(xvec)//4
+    # global omega_init
+    # omega_init = True
+    xvec_in = xvec
+    poles = len(xvec)//params_count
     # for i in range(poles):
-    #     xvec[i] = np.absolute(xvec[i])
-    #     xvec[i+1] = np.absolute(xvec[i+1])
+    #     for j in range(2):
+    #         if (xvec_in[i+j] < 0):
+    #             xvec_in[i+j] = 0
+    #     # if (xvec_in[i+2] > 0):
+    #         # xvec_in[i+2] = 0
 
     val = np.zeros(omega.size, dtype=np.cdouble)
-    for i in prange(omega.size):
-        val[i] = lorentzian(omega[i], xvec)
+    # print('in eval', xvec_in)
+    for i in range(omega.size):
+        val[i] = lorentzian(omega[i], xvec_in)
     return val

examples/NatureComm-2020-Plasmon–emitter_interactions_at_the_nanoscale/fit_d_param.py

@@ -1,78 +1,74 @@
 #!/usr/bin/env python3
 # -*- coding: UTF-8 -*-
-from functools import lru_cache
-from matplotlib import markers, pyplot as plt
-from scipy.optimize import curve_fit
+from matplotlib import pyplot as plt
 import numpy as np
 import cma
-# import pyfde
-from numba import njit, float64
-from eval_spectra import spectra
-
-from mealpy.physics_based.EO import AdaptiveEO
-
+from eval_spectra import spectra, params_count as pc
 
 from_disk = np.loadtxt('rs4-d_perp_interpolated.txt')
 step = 5
 omega_ratio = np.copy(from_disk[0, ::step])
 d_rs4 = from_disk[1, ::step] + 1j*from_disk[2, ::step]
 
-d_rms = d_rs4
-
 
 def rms(x0):
     d_fit = spectra(omega_ratio, x0)
     diff_re = np.real(d_rms - d_fit)
-    rms = np.sqrt(np.sum(np.abs(diff_re)**2))
     diff_im = np.imag(d_rms - d_fit)
-    rms += np.sqrt(np.sum(np.abs(diff_im)**2))
+    # rms = np.sqrt(np.sum(np.abs(diff_re)**2))
+    # rms += np.sqrt(np.sum(np.abs(diff_im)**2))
+    rms = (np.sum(np.abs(diff_re)**2))
+    rms += (np.sum(np.abs(diff_im)**2))
     return rms
 
 
-poles = 1
-dim = poles*4
-x0 = np.random.random(dim)
+x0 = np.random.random(pc)
 d_rms = d_rs4
-# x, es = cma.fmin2(rms, x0, sigma0=0.2)
-x = np.array([0.13421489, 0.82250415, -0.50359304, -0.0591722])
-x1 = x
+x, es = cma.fmin2(rms, x0, sigma0=0.2)
+# x = np.array([0.1332754793043937, 0.8248539955310855, -
+#              0.5024906405674647, -0.08076797734842008])
+x1 = np.copy(x)
+
 
-x0 = np.random.random(dim)
+x0 = np.random.random(pc)
 d_rms = d_rs4 - spectra(omega_ratio, x1)
 x, es = cma.fmin2(rms, x0, sigma0=2)
-# x = [0.00051888486267353, 0.9991499897780783,
-#      0.06926055304806265, -0.030608812209114836] # fitness = 4.7
-x2 = x
+# x = np.array([0.0019369902204593222, 0.9978752165162739,
+#              0.05801769075873917, -0.032110084386726336])
+x2 = np.copy(x)
+
+
+x0 = np.hstack((x1, x2))
+d_rms = d_rs4
+x, es = cma.fmin2(rms, x0, sigma0=0.02)
+# x = np.array([0.11878109939932953, 0.8142072049467617, -0.43466301805510305, -0.012772472816983446,
+#               0.012573279034397847, 1.0010563127596699, 0.07665592968844606, -0.07679477702750433])
+x12 = x
+
+
+x0 = np.random.random(pc)
+d_rms = d_rs4 - spectra(omega_ratio, x12)
+x, es = cma.fmin2(rms, x0, sigma0=0.2)
+# x = np.array([0.1434266344187499, 0.5188802822956783, -
+#              0.00950433613285183, 0.013585684987833985])
+x3 = np.copy(x)
 
-d_fit = spectra(omega_ratio, x)
 
+x0 = np.hstack((x1, x2, x3))
+d_rms = d_rs4
+x, es = cma.fmin2(rms, x0, sigma0=0.02)
+# x = np.array([0.09914803, 0.81158511, -0.34941712, 0.01388308,
+#               0.01560184, 1.00551237, 0.11006553, -0.09818891,
+#               0.34432793, 0.64182428, -0.0803532, 0.04641341])
+x123 = x
+
+
+d_fit = spectra(omega_ratio, x)
 plt.figure('rs4')
-# plt.title('rms = '+str(rms(x)/x.size))
 plt.plot(omega_ratio, np.real(d_rms), label='re d')
 plt.plot(omega_ratio, np.imag(d_rms), label='im d')
-
-# plt.plot(omega_ratio, np.real(
-#     d_rs4 - spectra(omega_ratio, x1)), label='diff re d')
-# plt.plot(omega_ratio, np.imag(
-#     d_rs4 - spectra(omega_ratio, x1)), label='diff im d')
-
 plt.plot(omega_ratio, np.real(d_fit), label='re d fit', alpha=0.2, lw=3)
 plt.plot(omega_ratio, np.imag(d_fit), label='im d fit', alpha=0.2, lw=3)
-
-# plt.plot(omega_ratio, func(xdata, *popt), 'r-',
-
-#          label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
+plt.axhline(y=0.0, color='black', linestyle='-', lw=1)
 plt.legend()
 plt.show()
-
-
-# from_disk = np.loadtxt('silver-d_perp_interpolated.txt')
-# plt.figure('silver')
-# plt.plot(from_disk[0, :], from_disk[1, :], label='re d perp')
-# plt.plot(from_disk[0, :], from_disk[2, :], label='im d perp')
-# from_disk = np.loadtxt('silver-d_parl_interpolated.txt')
-# plt.plot(from_disk[0, :], from_disk[1, :], label='re d parl')
-# plt.plot(from_disk[0, :], from_disk[2, :], label='im d parl')
-
-# plt.legend()
-# plt.show()