publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-07 23:14:03 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

code refactor

浏览代码
这个提交包含在：
majsylw
2021-08-16 10:37:43 +02:00
父节点 ba997a648f
当前提交 5ebbe7f3f0
共有 2 个文件被更改，包括 102 次插入 和 82 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

69
user_libs/DebyeFit/Debye_Fit.py
查看文件

@@ -25,7 +25,7 @@ import sys
import scipy.interpolate
import warnings
from optimization import *
from optimization import PSO_DLS, DA_DLS, DE_DLS
class Relaxation(object):
@@ -158,7 +158,7 @@ class Relaxation(object):
funckwargs={'rl': self.rl,
'im': self.im,
'freq': self.freq}
)
)
return tau, weights, ee, rl, im
def run(self):
@@ -186,7 +186,7 @@ class Relaxation(object):
print("\n#########",
"Try to automaticaly fit number of Debye poles, up to 20!",
"##########\n", sep="")
error = np.infty # artificial best error starting value
error = np.infty # artificial best error starting value
self.number_of_debye_poles = 1
iteration = 1
# stop increasing number of Debye poles if error is smaller then 5%
@@ -267,7 +267,7 @@ class Relaxation(object):
for given frequency points.
"""
plt.close("all")
fig = plt.figure(figsize=(16,8), tight_layout=True)
fig = plt.figure(figsize=(16, 8), tight_layout=True)
gs = gridspec.GridSpec(2, 1)
ax = fig.add_subplot(gs[0])
ax.grid(b=True, which="major", linewidth=0.2, linestyle="--")
@@ -366,7 +366,7 @@ class HavriliakNegami(Relaxation):
For beta = 1 and alpha !=0 & alpha !=1 Havriliak-Negami
transforms to Cole-Cole function.
:type beta: float
:param de: The difference of relative permittivity at infinite frequency
:param de: The difference of relative permittivity at infinite frequency
and the relative permittivity at zero frequency.
:type de: float
:param tau_0: Real positive float number, tau_0 is the relaxation time.
@@ -394,9 +394,9 @@ class HavriliakNegami(Relaxation):
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
self.set_freq(self.f_min, self.f_max, self.f_n)
self.e_inf, self.alpha, self.beta, self.de, self.tau_0 = e_inf, alpha, beta, de, tau_0
self.params = {'f_min':self.f_min, 'f_max':self.f_max,
'eps_inf':self.e_inf, 'Delta_eps':self.de, 'tau_0':self.tau_0,
'alpha':self.alpha, 'beta':self.beta}
self.params = {'f_min': self.f_min, 'f_max': self.f_max,
'eps_inf': self.e_inf, 'Delta_eps': self.de, 'tau_0': self.tau_0,
'alpha': self.alpha, 'beta': self.beta}
def check_inputs(self):
""" Check the validity of the Havriliak Negami model's inputs. """
@@ -422,6 +422,7 @@ class HavriliakNegami(Relaxation):
self.freq * self.tau_0)**self.alpha
)**self.beta
class Jonscher(Relaxation):
""" Approximate a given Jonsher function
Jonscher function = ε_∞ - ap * (-1j * 2πf / omegap)**n_p,
@@ -464,9 +465,9 @@ class Jonscher(Relaxation):
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
self.set_freq(self.f_min, self.f_max, self.f_n)
self.e_inf, self.a_p, self.omega_p, self.n_p = e_inf, a_p, omega_p, n_p
self.params = {'f_min':self.f_min, 'f_max':self.f_max,
'eps_inf':self.e_inf, 'n_p':self.n_p,
'omega_p':self.omega_p, 'a_p':self.a_p}
self.params = {'f_min': self.f_min, 'f_max': self.f_max,
'eps_inf': self.e_inf, 'n_p': self.n_p,
'omega_p': self.omega_p, 'a_p': self.a_p}
def check_inputs(self):
""" Check the validity of the inputs. """
@@ -484,7 +485,7 @@ class Jonscher(Relaxation):
def calculation(self):
"""Calculates the Q function for the given parameters"""
return self.e_inf + (self.a_p * (2 * np.pi *
return self.e_inf + (self.a_p * (2 * np.pi *
self.freq / self.omega_p)**(self.n_p-1)) * (
1 - 1j / np.tan(self.n_p * np.pi/2))
@@ -508,7 +509,7 @@ class Crim(Relaxation):
"""
def __init__(self, f_min, f_max, a, volumetric_fractions,
materials, sigma, mu, mu_sigma, material_name,
materials, sigma, mu, mu_sigma, material_name,
number_of_debye_poles=-1, f_n=50,
plot=False, save=False,
optimizer=PSO_DLS,
@@ -531,9 +532,9 @@ class Crim(Relaxation):
self.a = a
self.volumetric_fractions = np.array(volumetric_fractions)
self.materials = np.array(materials)
self.params = {'f_min':self.f_min, 'f_max':self.f_max,
'a':self.a, 'volumetric_fractions':self.volumetric_fractions,
'materials':self.materials}
self.params = {'f_min': self.f_min, 'f_max': self.f_max,
'a': self.a, 'volumetric_fractions': self.volumetric_fractions,
'materials': self.materials}
def check_inputs(self):
""" Check the validity of the inputs. """
@@ -600,7 +601,7 @@ class Rawdata(Relaxation):
def __init__(self, filename,
sigma, mu, mu_sigma,
material_name, number_of_debye_poles=-1,
f_n=50, delimiter =',',
f_n=50, delimiter=',',
plot=False, save=False,
optimizer=PSO_DLS,
optimizer_options={}):
@@ -613,7 +614,7 @@ class Rawdata(Relaxation):
optimizer_options=optimizer_options)
self.delimiter = delimiter
self.filename = Path(filename).absolute()
self.params = {'filename':self.filename}
self.params = {'filename': self.filename}
def check_inputs(self):
""" Check the validity of the inputs. """
@@ -645,7 +646,7 @@ class Rawdata(Relaxation):
if __name__ == "__main__":
### Kelley et al. parameters
# Kelley et al. parameters
setup = HavriliakNegami(f_min=1e7, f_max=1e11,
alpha=0.91, beta=0.45,
e_inf=2.7, de=8.6-2.7, tau_0=9.4e-10,
@@ -653,13 +654,13 @@ if __name__ == "__main__":
material_name="Kelley", f_n=100,
number_of_debye_poles=6,
plot=True, save=False,
optimizer_options={'swarmsize':30,
'maxiter':100,
'omega':0.5,
'phip':1.4,
'phig':1.4,
'minstep':1e-8,
'minfun':1e-8,
optimizer_options={'swarmsize': 30,
'maxiter': 100,
'omega': 0.5,
'phip': 1.4,
'phig': 1.4,
'minstep': 1e-8,
'minfun': 1e-8,
'seed': 111,
'pflag': True})
setup.run()
@@ -670,7 +671,7 @@ if __name__ == "__main__":
material_name="Kelley", f_n=100,
number_of_debye_poles=6,
plot=True, save=False,
optimizer=DA,
optimizer=DA_DLS,
optimizer_options={'seed': 111})
setup.run()
setup = HavriliakNegami(f_min=1e7, f_max=1e11,
@@ -680,24 +681,24 @@ if __name__ == "__main__":
material_name="Kelley", f_n=100,
number_of_debye_poles=6,
plot=True, save=False,
optimizer=DE,
optimizer=DE_DLS,
optimizer_options={'seed': 111})
setup.run()
### Testing setup
setup = Rawdata("examples/Test.txt", 0.1, 1, 0.1, "M1",
# Testing setup
setup = Rawdata("examples/Test.txt", 0.1, 1, 0.1, "M1",
number_of_debye_poles=3, plot=True,
optimizer_options={'seed': 111})
setup.run()
np.random.seed(111)
setup = HavriliakNegami(1e12, 1e-3, 0.5, 1, 10, 5,
1e-6, 0.1, 1, 0, "M2",
1e-6, 0.1, 1, 0, "M2",
number_of_debye_poles=6,
plot=True)
setup.run()
setup = Jonscher(1e6, 1e-5, 50, 1, 1e5, 0.7,
0.1, 1, 0.1, "M3",
number_of_debye_poles=4,
plot=True)
0.1, 1, 0.1, "M3",
number_of_debye_poles=4,
plot=True)
setup.run()
f = np.array([0.5, 0.5])
material1 = [3, 25, 1e6]

115
user_libs/DebyeFit/optimization.py
查看文件

@@ -54,8 +54,10 @@ class Optimizer(object):
Returns:
tau (ndarray): The the best relaxation times.
weights (ndarray): Resulting optimised weights for the given relaxation times.
ee (float): Average error between the actual and the approximated real part.
weights (ndarray): Resulting optimised weights for the given
relaxation times.
ee (float): Average error between the actual and the approximated
real part.
rl (ndarray): Real parts of chosen relaxation function
for given frequency points.
im (ndarray): Imaginary parts of chosen relaxation function
@@ -67,7 +69,8 @@ class Optimizer(object):
# find the weights using a calc_weights method
if self.calc_weights is None:
raise NotImplementedError()
_, _, weights, ee, rl_exp, im_exp = self.calc_weights(tau, **funckwargs)
_, _, weights, ee, rl_exp, im_exp = \
self.calc_weights(tau, **funckwargs)
return tau, weights, ee, rl_exp, im_exp
def calc_relaxation_times(self):
@@ -84,7 +87,8 @@ class Optimizer(object):
the actual and the approximated electric permittivity.
Args:
x (ndarray): The logarithm with base 10 of relaxation times of the Debyes poles.
x (ndarray): The logarithm with base 10 of relaxation times
of the Debyes poles.
rl (ndarray): Real parts of chosen relaxation function
for given frequency points.
im (ndarray): Imaginary parts of chosen relaxation function
@@ -92,7 +96,8 @@ class Optimizer(object):
freq (ndarray): The frequencies vector for defined grid.
Returns:
cost (float): Sum of mean absolute errors for real and imaginary parts.
cost (float): Sum of mean absolute errors for real and
imaginary parts.
"""
cost_i, cost_r, _, _, _, _ = DLS(x, rl, im, freq)
return cost_i + cost_r
@@ -162,11 +167,13 @@ class PSO_DLS(Optimizer):
"""
np.random.seed(self.seed)
# check input parameters
assert len(lb) == len(ub), 'Lower- and upper-bounds must be the same length'
assert len(lb) == len(ub), \
'Lower- and upper-bounds must be the same length'
assert hasattr(func, '__call__'), 'Invalid function handle'
lb = np.array(lb)
ub = np.array(ub)
assert np.all(ub > lb), 'All upper-bound values must be greater than lower-bound values'
assert np.all(ub > lb), \
'All upper-bound values must be greater than lower-bound values'
vhigh = np.abs(ub - lb)
vlow = -vhigh
@@ -321,26 +328,28 @@ class DA_DLS(Optimizer):
fun (float): The objective value at the best solution.
"""
np.random.seed(self.seed)
result = scipy.optimize.dual_annealing(func,
bounds=list(zip(lb, ub)),
args=funckwargs.values(),
maxiter=self.maxiter,
local_search_options=self.local_search_options,
initial_temp=self.initial_temp,
restart_temp_ratio=self.restart_temp_ratio,
visit=self.visit,
accept=self.accept,
maxfun=self.maxfun,
no_local_search=self.no_local_search,
callback=self.callback,
x0=self.x0)
result = scipy.optimize.dual_annealing(
func,
bounds=list(zip(lb, ub)),
args=funckwargs.values(),
maxiter=self.maxiter,
local_search_options=self.local_search_options,
initial_temp=self.initial_temp,
restart_temp_ratio=self.restart_temp_ratio,
visit=self.visit,
accept=self.accept,
maxfun=self.maxfun,
no_local_search=self.no_local_search,
callback=self.callback,
x0=self.x0)
print(result.message)
return result.x, result.fun
class DE_DLS(Optimizer):
"""
Create Differential Evolution-Damped Least Squares object with predefined parameters.
Create Differential Evolution-Damped Least Squares object
with predefined parameters.
The current class is a modified edition of the scipy.optimize
package which can be found at:
https://docs.scipy.org/doc/scipy/reference/generated/
@@ -349,7 +358,8 @@ class DE_DLS(Optimizer):
def __init__(self, maxiter=1000,
strategy='best1bin', popsize=15, tol=0.01, mutation=(0.5, 1),
recombination=0.7, callback=None, disp=False, polish=True,
init='latinhypercube', atol=0, updating='immediate', workers=1,
init='latinhypercube', atol=0,
updating='immediate', workers=1,
constraints=(), seed=None):
super(DE_DLS, self).__init__(maxiter, seed)
self.strategy = strategy
@@ -360,7 +370,7 @@ class DE_DLS(Optimizer):
self.callback = callback
self.disp = disp
self.polish = polish
self.init= init
self.init = init
self.atol = atol
self.updating = updating
self.workers = workers
@@ -388,22 +398,23 @@ class DE_DLS(Optimizer):
fun (float): The objective value at the best solution.
"""
np.random.seed(self.seed)
result = scipy.optimize.differential_evolution(func,
bounds=list(zip(lb, ub)),
args=funckwargs.values(),
strategy=self.strategy,
popsize=self.popsize,
tol=self.tol,
mutation=self.mutation,
recombination=self.recombination,
callback=self.callback,
disp=self.disp,
polish=self.polish,
init=self.init,
atol=self.atol,
updating=self.updating,
workers=self.workers,
constraints=self.constraints)
result = scipy.optimize.differential_evolution(
func,
bounds=list(zip(lb, ub)),
args=funckwargs.values(),
strategy=self.strategy,
popsize=self.popsize,
tol=self.tol,
mutation=self.mutation,
recombination=self.recombination,
callback=self.callback,
disp=self.disp,
polish=self.polish,
init=self.init,
atol=self.atol,
updating=self.updating,
workers=self.workers,
constraints=self.constraints)
print(result.message)
return result.x, result.fun
@@ -415,8 +426,10 @@ def DLS(logt, rl, im, freq):
also known as the damped least-squares (DLS) method.
Args:
logt (ndarray): The best known position form optimization module (optimal design),
the logarithm with base 10 of relaxation times of the Debyes poles.
logt (ndarray): The best known position form optimization module
(optimal design),
the logarithm with base 10 of relaxation times
of the Debyes poles.
rl (ndarray): Real parts of chosen relaxation function
for given frequency points.
im (ndarray): Imaginary parts of chosen relaxation function
@@ -428,12 +441,16 @@ def DLS(logt, rl, im, freq):
the approximated imaginary part.
cost_r (float): Mean absolute error between the actual and
the approximated real part (plus average error).
x (ndarray): Resulting optimised weights for the given relaxation times.
ee (float): Average error between the actual and the approximated real part.
rp (ndarray): The real part of the permittivity for the optimised relaxation
times and weights for the frequnecies included in freq.
x (ndarray): Resulting optimised weights for the given
relaxation times.
ee (float): Average error between the actual and the approximated
real part.
rp (ndarray): The real part of the permittivity for the optimised
relaxation times and weights for the frequnecies included
in freq.
ip (ndarray): The imaginary part of the permittivity for the optimised
relaxation times and weights for the frequnecies included in freq.
relaxation times and weights for the frequnecies included
in freq.
"""
# The relaxation time of the Debyes are given at as logarithms
# logt=log10(t0) for efficiency during the optimisation
@@ -444,8 +461,10 @@ def DLS(logt, rl, im, freq):
freq, len(tt)).reshape((-1, len(tt))) * tt)
# Adding dumping (Levenberg–Marquardt algorithm)
# Solving the overdetermined system y=Ax
x = np.abs(np.linalg.lstsq(d.imag, im, rcond=None)[0]) # absolute damped least-squares solution
rp, ip = np.matmul(d.real, x[np.newaxis].T).T[0], np.matmul(d.imag, x[np.newaxis].T).T[0]
x = np.abs(np.linalg.lstsq(d.imag, im, rcond=None)[0])
# x - absolute damped least-squares solution
rp, ip = np.matmul(d.real, x[np.newaxis].T).T[0], np.matmul(
d.imag, x[np.newaxis].T).T[0]
cost_i = np.sum(np.abs(ip-im))/len(im)
ee = np.mean(rl - rp)
if ee < 1: