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

code refactor

浏览代码
这个提交包含在：
majsylw
2021-07-06 12:01:21 +02:00
父节点 0e12a19f9c
当前提交 7a1d9289d3
共有 3 个文件被更改，包括 548 次插入 和 387 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

674
user_libs/DebyeFit/Debye_Fit.py 普通文件
查看文件

@@ -0,0 +1,674 @@
# Author: Iraklis Giannakis, Sylwia Majchrowska
# E-mail: i.giannakis@ed.ac.uk
#
# Copyright (c) 2021 gprMax
# All rights reserved.
#
# Redistribution and use in source and binary forms are permitted
# provided that the above copyright notice and this paragraph are
# duplicated in all such forms and that any documentation,
# advertising materials, and other materials related to such
# distribution and use acknowledge that the software was developed
# as part of gprMax. The name of gprMax may not be used to
# endorse or promote products derived from this software without
# specific prior written permission.
# THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
import numpy as np
import os
from matplotlib import pylab as plt
import sys
import scipy.interpolate
from optimization import *
class Relaxation(object):
def __init__(self, number_of_debye_poles,
sigma, mu, mu_sigma,
material_name, plot=True, save=True,
optimizer=Particle_swarm,
optimizer_options={}):
"""
Create Relaxation function object for complex material.
Args:
number_of_debye_poles (int): Number of Debye functions used to
approximate the given electric
permittivity.
sigma (float): Conductivity.
mu (float): Relative permabillity.
mu_sigma (float): Magnetic looses.
material_name (str): A string containing the given name of
the material (e.g. "Clay").
plot (bool): if True will plot the actual and the approximated
permittivity (it can be neglected).
The argument is optional and if neglected plot=False.
save (bool): if True will save approximated material parameters
The argument is optional and if neglected save=False.
optimizer (Optimizer class): chosen optimization method:
Particle Swarm, Genetic or Dual Annealing.
optimizer_options (dict): Additional keyword arguments passed to
optimizer class (Default: empty dict).
"""
self.number_of_debye_poles = number_of_debye_poles
self.sigma = sigma
self.mu = mu
self.mu_sigma = mu_sigma
self.material_name = material_name
self.plot = plot
self.save = save
self.optimizer = optimizer(**optimizer_options)
def run(self):
"""
Solve the problem described by the given relaxation function
(Havriliak-Negami function, Crim, Jonscher)
or data given from a text file.
"""
# Check the validity of the inputs
self.check_inputs()
# Print information about chosen approximation settings
self.print_info()
# Calculate both real and imaginary parts
# for the frequencies included in the vector freq
q = self.calculation()
# Set the real and the imaginary part of the relaxation function
self.rl, self.im = q.real, q.imag
# Calling the main optimisation module
# if one of the weights is negative increase the stabiliser
# and repeat the optimisation
xmp, mx, ee, rp, ip = self.optimize()
# Print the results in gprMax format style
properties = self.print_output(xmp, mx, ee)
if self.save:
self.save_result(properties)
# Plot the actual and the approximate dielectric properties
if self.plot:
self.plot_result(rp + ee, ip)
def set_freq(self, f_min, f_max, n=50):
"""
Interpolate frequency vector using n
equally logarithmicaly spaced frequencies.
Args:
f_min (float): First bound of the frequency range
used to approximate the given function (Hz).
f_max (float): Second bound of the frequency range
used to approximate the given function (Hz).
n (int): Number of frequency points in frequency grid.
Note:
f_min and f_max must satisfied f_min < f_max
"""
# diff_freq = np.log10(f_max) - np.log10(f_min)
self.freq = np.logspace(np.log10(f_min) + 0.00001,
np.log10(f_max) - 0.00001,
int(n))
# int(n * diff_freq))
def check_inputs(self):
"""
Check the validity of the inputs.
"""
try:
d = [float(i) for i in
[self.number_of_debye_poles,
self.sigma, self.mu, self.mu_sigma]]
except ValueError:
sys.exit("The inputs should be numeric.")
if not isinstance(self.number_of_debye_poles, int):
sys.exit("The number of Debye poles must be integer.")
if (np.array(d) < 0).sum() != 0:
sys.exit("The inputs should be positive.")
def calculation(self):
"""
Approximate the given relaxation function
(Havriliak-Negami function, Crim, Jonscher) or based on raw data.
"""
raise NotImplementedError()
def print_info(self):
"""
Print information about chosen approximation settings.
"""
raise NotImplementedError()
def optimize(self):
"""
Calling the main optimisation module.
Returns:
xpm (array): The best known position form optimization module
(optimal design).
mx (array): Resulting optimised weights for the given relaxation times
ee (float): Average error between the actual and the approximated real part
rp (matrix): The real part of the permittivity for the optimised relaxation
times and weights for the frequnecies included in freq
ip (matrix): The imaginary part of the permittivity for the optimised
relaxation times and weights for the frequnecies included in freq
"""
# Define the lower and upper boundaries of search
lb = np.full(self.number_of_debye_poles,
-np.log10(np.max(self.freq)) - 3)
ub = np.full(self.number_of_debye_poles,
-np.log10(np.min(self.freq)) + 3)
# Call particle swarm optimisation to minimize the cost function.
xmp, _ = self.optimizer.fit(func=cost_function,
lb=lb, ub=ub,
funckwargs={'rl_g': self.rl,
'im_g': self.im,
'freq_g': self.freq}
)
_, _, mx, ee, rp, ip = linear(self.rl, self.im, xmp, self.freq)
return xmp, mx, ee, rp, ip
def print_output(self, xmp, mx, ee):
"""
Print out the resulting Debye parameters in a gprMax format.
Args:
xpm (array): The best known position form optimization module
(optimal design).
mx (): Resulting optimised weights for the given relaxation times.
ee (): Average error between the actual and the approximated real part.
"""
print("Debye expansion parameters: ")
print(f" |{'e_inf':^14s}|{'De':^14s}|{'log(t0)':^25s}|")
print("_" * 65)
for i in range(0, len(xmp)):
print("Debye {0:}:|{1:^14.5f}|{2:^14.5f}|{3:^25.5f}|"
.format(i + 1, ee/len(xmp), mx[i],
xmp[i]))
print("_" * 65)
# Print the Debye expnasion in a gprMax format
material_prop = "#material: {} {} {} {} {}".format(ee, self.sigma,
self.mu,
self.mu_sigma,
self.material_name)
print(material_prop)
material_prop = [material_prop + '\n']
dispersion_prop = "#add_dispersion_debye: {} {} {}".format(len(xmp),
mx[0],
10**xmp[0])
for i in range(1, len(xmp)):
dispersion_prop += " {} {}".format(mx[i], 10**xmp[i])
dispersion_prop += " {}".format(self.material_name)
print(dispersion_prop)
material_prop.append(dispersion_prop + '\n')
return material_prop
def plot_result(self, rl_exp, im_exp):
"""
Plot the actual and the approximated electric permittivity
using a semilogarithm X axes.
Args:
rl_exp (array): Real parts of optimised Debye expansion
for given frequency points (plus average error).
im_exp (array): Imaginary parts of optimised Debye expansion
for given frequency points.
"""
plt.close("all")
plt.rcParams["axes.facecolor"] = "black"
plt.semilogx(self.freq * 1e-6, rl_exp, "b-", linewidth=2.0,
label="Debye Expansion: Real")
plt.semilogx(self.freq * 1e-6, -im_exp, "w-", linewidth=2.0,
label="Debye Expansion: Imaginary")
plt.semilogx(self.freq * 1e-6, self.rl, "ro",
linewidth=2.0, label="Chosen Function: Real")
plt.semilogx(self.freq * 1e-6, -self.im, "go", linewidth=2.0,
label="Chosen Function: Imaginary")
plt.rcParams["axes.facecolor"] = "white"
plt.grid(b=True, which="major", color="w", linewidth=0.2,
linestyle="--")
axes = plt.gca()
axes.set_xlim([np.min(self.freq * 1e-6), np.max(self.freq * 1e-6)])
axes.set_ylim([-1, np.max(np.concatenate([self.rl, -self.im])) + 1])
plt.legend()
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative permittivity")
plt.show()
@staticmethod
def save_result(output, fdir="../materials"):
"""
Save the resulting Debye parameters in a gprMax format
Args:
output (str): Material and resulting Debye parameters
in a gprMax format.
fdir (str): Path to saving directory.
"""
if fdir != "../materials" and os.path.isdir(fdir):
file_path = os.path.join(fdir, "my_materials.txt")
elif os.path.isdir("../materials"):
file_path = os.path.join("../materials",
"my_materials.txt")
elif os.path.isdir("materials"):
file_path = os.path.join("materials",
"my_materials.txt")
elif os.path.isdir("user_libs/materials"):
file_path = os.path.join("user_libs", "materials",
"my_materials.txt")
else:
sys.exit("Cannot save material properties "
f"in {os.path.join(fdir, 'my_materials.txt')}!")
fileH = open(file_path, "a")
fileH.write(f"## {output[0].split(' ')[-1]}")
fileH.writelines(output)
fileH.write("\n")
fileH.close()
print(f"Material properties save at: {file_path}")
class HavriliakNegami(Relaxation):
def __init__(self, number_of_debye_poles,
f_min, f_max, alfa, bita, einf, de, t0,
sigma, mu, mu_sigma,
material_name, plot=False, save=True,
optimizer=Particle_swarm,
optimizer_options={}):
"""
Approximate a given Havriliak-Negami function
Havriliak-Negami function = einf + de / (1 + (1j * 2 * pi * f *t0)**alfa )**bita,
where f is the frequency in Hz
Args:
number_of_debye_poles (int): Number of Debye functions used to
approximate the given electric
permittivity.
f_min (float): Define the first bound of the frequency range
used to approximate the given function (Hz).
f_max (float): Define the second bound of the frequency range
used to approximate the given function (Hz).
Note: f_min and f_max can be either f_min > f_max
or f_min < f_max but not f_min = f_max.
einf (float): The real relative permittivity at infinity frequency
alfa (float): Havriliak-Negami parameter. Real positive float number
which varies 0 < alfa < 1. For alfa = 1 and bita !=0 & bita !=1
Havriliak-Negami transforms to Cole-Davidson function.
bita (float): Havriliak-Negami parameter. Real positive float number
which varies 0 < bita < 1. For bita = 1 and alfa !=0 & alfa !=1
Havriliak-Negami transforms to Cole-Cole function.
de (float): Havriliak-Negami parameter. Real positive float number.
de is the relative permittivity at infinite frequency
minus the relative permittivity at zero frequency.
t0 (float): Havriliak_Negami parameter. Real positive float number.
t0 is the relaxation time.
sigma (float): Conductivity.
mu (float): Relative permabillity.
mu_sigma (float): Magnetic looses.
material_name (str): A string containing the given name of
the material (e.g. "Clay").
plot (bool): if True will plot the actual and the approximated
permittivity. The argument is optional and
if neglected plot=False.
save (bool): if True will save approximated material parameters
The argument is optional and if neglected save=False.
optimizer (Optimizer class): chosen optimization method:
Particle Swarm, Genetic or Dual Annealing.
(Default: Partocle_swarm)
optimizer_options (dict): Additional keyword arguments passed to
optimizer class (Default: empty dict).
"""
super(HavriliakNegami, self).__init__(number_of_debye_poles,
sigma, mu, mu_sigma,
material_name, plot, save,
optimizer, optimizer_options)
# Place the lower frequency bound at f_min and the upper frequency bound at f_max
if f_min > f_max:
self.f_min, self.f_max = f_max, f_min
else:
self.f_min, self.f_max = f_min, f_max
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
self.set_freq(self.f_min, self.f_max)
self.einf, self.alfa, self.bita, self.de, self.t0 = einf, alfa, bita, de, t0
def check_inputs(self):
"""
Check the validity of the inputs.
"""
super(HavriliakNegami, self).check_inputs()
try:
d = [float(i) for i in
[self.f_min, self.f_max, self.alfa,
self.bita, self.einf, self.de, self.t0]]
except ValueError:
sys.exit("The inputs should be numeric.")
if (np.array(d) < 0).sum() != 0:
sys.exit("The inputs should be positive.")
if self.alfa > 1:
sys.exit("Alfa value must range between 0-1 (0 <= Alfa <= 1)")
if self.bita > 1:
sys.exit("Beta value must range between 0-1 (0 <= Beta <= 1)")
if self.f_min == self.f_max:
sys.exit("Null frequency range")
def print_info(self):
"""Print information about chosen approximation settings."""
print(f"Approximating Havriliak-Negami function"
f" using {self.number_of_debye_poles} Debye poles")
print("Havriliak-Negami parameters : ")
print("De = {} \ne_inf = {} \nt0 = {} \nalfa = {} \nbita = {} "
.format(self.de, self.einf, self.t0, self.alfa, self.bita))
def calculation(self):
"""Calculates the Havriliak-Negami function for
the given parameters."""
return self.einf + self.de / (np.array(
1 + np.array(1j * 2 * np.pi *
self.freq * self.t0
) ** self.alfa)**self.bita)
class Jonscher(Relaxation):
def __init__(self, number_of_debye_poles,
f_min, f_max, einf, ap, omegap, n_p,
sigma, mu, mu_sigma,
material_name, plot=False, save=True,
optimizer=Particle_swarm,
optimizer_options={}):
"""
Approximate a given Johnsher function
Jonscher function = einf - ap * ( -1j * 2 * pi * f / omegap)**n_p,
where f is the frequency in Hz
Args:
number_of_debye_poles (int): Number of Debye functions used to
approximate the given electric
permittivity.
f_min (float): Define the first bound of the frequency range
used to approximate the given function (Hz).
f_max (float): Define the second bound of the frequency range
used to approximate the given function (Hz).
f_min and f_max can be either f_min > f_max
or f_min < f_max but not f_min = f_max.
einf (float): The real relative permittivity at infinity frequency
ap (float): Jonscher parameter. Real positive float number.
omegap (float): Jonscher parameter. Real positive float number.
n_p (float): Jonscher parameter.
Real positive float number which varies 0 < n_p < 1.
sigma (float): Conductivity.
mu (float): Relative permabillity.
mu_sigma (float): Magnetic looses.
material_name (str): A string containing the given name of
the material (e.g. "Clay").
plot (bool): if True will plot the actual and the approximated
permittivity. The argument is optional and
if neglected plot=False.
save (bool): if True will save approximated material parameters
The argument is optional and if neglected save=False.
optimizer (Optimizer class): chosen optimization method:
Particle Swarm, Genetic or Dual Annealing.
optimizer_options (dict): Additional keyword arguments passed to
optimizer class (Default: empty dict).
"""
super(Jonscher, self).__init__(number_of_debye_poles,
sigma, mu, mu_sigma,
material_name, plot, save,
optimizer, optimizer_options)
# Place the lower frequency bound at f_min and the upper frequency bound at f_max
if f_min > f_max:
self.f_min, self.f_max = f_max, f_min
else:
self.f_min, self.f_max = f_min, f_max
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
self.set_freq(self.f_min, self.f_max)
self.einf, self.ap, self.omegap, self.n_p = einf, ap, omegap, n_p
def check_inputs(self):
"""
Check the validity of the inputs.
"""
super(Jonscher, self).check_inputs()
try:
d = [float(i) for i in
[self.f_min, self.f_max, self.n_p,
self.einf, self.omegap, self.ap]]
except ValueError:
sys.exit("The inputs should be numeric.")
if (np.array(d) < 0).sum() != 0:
sys.exit("The inputs should be positive.")
if self.n_p > 1:
sys.exit("n_p value must range between 0-1 (0 <= n_p <= 1)")
if self.f_min == self.f_max:
sys.exit("Error: Null frequency range")
def print_info(self):
"""
Print information about chosen approximation settings
"""
print(f"Approximating Jonsher function"
f" using {self.number_of_debye_poles} Debye poles")
print("Jonhser function parameters : ")
print(f"omega_p = {self.omegap}\n"
f"e_inf = {self.einf}\n"
f"n_p = {self.n_p}\n"
f"A_p = {self.ap}")
def calculation(self):
"""Calculates the Q function for the given parameters"""
return self.einf + (self.ap * np.array(
2 * np.pi * self.freq / self.omegap
)**(self.n_p-1)) * (
1 - 1j / np.tan(self.n_p * np.pi/2))
class Crim(Relaxation):
def __init__(self, number_of_debye_poles,
f_min, f_max, a, f1, e1, sigma,
mu, mu_sigma, material_name, plot=False, save=True,
optimizer=Particle_swarm,
optimizer_options={}):
"""
Approximate a given CRIM function
CRIM = (sum([volumetric_fraction[i]*(material[i][0] + material[i][1] /
(1 + (1j * 2 * pi * f *material[i][2])))**m_param
for i in range(0,len(material))]))**1/m_param
Args:
number_of_debye_poles (int): Number of Debye functions used to
approximate the given electric
permittivity.
f_min (float): Define the first bound of the frequency range
used to approximate the given function (Hz).
f_max (float): Define the second bound of the frequency range
used to approximate the given function (Hz).
f_min and f_max can be either f_min > f_max
or f_min < f_max but not f_min = f_max.
a (float): shape factor
f1 (list): volumetric fraction
e1 (list): materials
sigma (float): Conductivity.
mu (float): Relative permabillity.
mu_sigma (float): Magnetic looses.
material_name (str): A string containing the given name of
the material (e.g. "Clay").
plot (bool): if True will plot the actual and the approximated
permittivity. The argument is optional and
if neglected plot=False.
save (bool): if True will save approximated material parameters
The argument is optional and if neglected save=False.
optimizer (Optimizer class): chosen optimization method:
Particle Swarm, Genetic or Dual Annealing.
optimizer_options (dict): Additional keyword arguments passed to
optimizer class (Default: empty dict).
"""
super(Crim, self).__init__(number_of_debye_poles,
sigma, mu, mu_sigma,
material_name, plot, save,
optimizer, optimizer_options)
# Place the lower frequency bound at f_min and the upper frequency bound at f_max
if f_min > f_max:
self.f_min, self.f_max = f_max, f_min
else:
self.f_min, self.f_max = f_min, f_max
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
self.set_freq(self.f_min, self.f_max)
self.a, self.f1, self.e1 = a, f1, e1
def check_inputs(self):
"""
Check the validity of the inputs.
"""
super(Crim, self).check_inputs()
try:
d = [float(i) for i in
[self.f_min, self.f_max, self.a]]
except ValueError:
sys.exit("The inputs should be numeric.")
if (np.array(d) < 0).sum() != 0:
sys.exit("The inputs should be positive.")
if len(self.f1) != len(self.e1):
sys.exit("Number of volumetric volumes does not match the dielectric properties")
# Check if the materials are at least two
if len(self.f1) < 2:
sys.exit("The materials should be at least 2")
# Check if the frequency range is null
if self.f_min == self.f_max:
sys.exit("Null frequency range")
# Check if the inputs are positive
f = [i for i in self.f1 if i < 0]
if len(f) != 0:
sys.exit("Error: The inputs should be positive")
for i in range(0, len(self.f1)):
f = [i for i in self.e1[i][:] if i < 0]
if len(f) != 0:
sys.exit("Error: The inputs should be positive")
# Check if the summation of the volumetric fractions equal to one
if np.sum(self.f1) != 1:
sys.exit("Error: The summation of volumetric volumes should be equal to 1")
def print_info(self):
"""
Print information about chosen approximation settings
"""
print(f"Approximating Complex Refractive Index Model (CRIM)"
f" using {self.number_of_debye_poles} Debye poles")
print("CRIM parameters : ")
for i in range(0, len(self.f1)):
print("Material {} :".format(i+1))
print("---------------------------------")
print(" Vol. fraction = {}".format(self.f1[i]))
print(" e_inf = {}".format(self.e1[i][0]))
print(" De = {}".format(self.e1[i][1]))
print(" log(t0) = {}".format(np.log10(self.e1[i][2])))
def calculation(self):
"""Calculates the Crim function for the given parameters"""
q = np.zeros(len(self.freq))
for i in range(len(self.f1)):
q = q + self.f1[i]*np.array(
[self.e1[i][0] + self.e1[i][1] /
(np.array(1 + np.array(1j * 2 * np.pi * f * self.e1[i][2])))
for f in self.freq])**self.a
return q**(1 / self.a)
class Rawdata(Relaxation):
def __init__(self, number_of_debye_poles,
filename,
sigma, mu, mu_sigma,
material_name, plot=False, save=True,
optimizer=Particle_swarm,
optimizer_options={}):
"""
Interpolate data given from a text file.
Args:
number_of_debye_poles (int): Number of Debye functions used to
approximate the given electric
permittivity.
filename (str): text file which contains three columns:
frequency (Hz),Real,Imaginary (separated by comma).
sigma (float): Conductivity.
mu (float): Relative permabillity.
mu_sigma (float): Magnetic looses.
material_name (str): A string containing the given name of
the material (e.g. "Clay").
plot (bool): if True will plot the actual and the approximated
permittivity (it can be neglected).
The argument is optional and if neglected plot=False.
save (bool): if True will save approximated material parameters
The argument is optional and if neglected save=False.
optimizer (Optimizer class): chosen optimization method:
Particle Swarm, Genetic or Dual Annealing.
optimizer_options (dict): Additional keyword arguments passed to
optimizer class (Default: empty dict).
"""
super(Rawdata, self).__init__(number_of_debye_poles,
sigma, mu, mu_sigma,
material_name, plot, save,
optimizer, optimizer_options)
self.filename = filename
def check_inputs(self):
"""
Check the validity of the inputs.
"""
super(Rawdata, self).check_inputs()
if not os.path.isfile(self.filename):
sys.exit("File doesnt exists!")
def print_info(self):
"""
Print information about chosen approximation settings
"""
print(f"Approximating the function given"
f" from file name {self.filename}"
f" using {self.number_of_debye_poles} Debye poles")
def calculation(self):
"""
Interpolate real and imaginary part from data.
Column framework of the input file three columns comma-separated
Frequency(Hz),Real,Imaginary
"""
# Read the file
with open(self.filename) as f:
try:
array = np.array(
[[float(x) for x in line.split(",")] for line in f]
)
except ValueError:
sys.exit("Error: The inputs should be numeric")
self.set_freq(min(array[:, 0]), max(array[:, 0]))
rl_interp = scipy.interpolate.interp1d(array[:, 0], array[:, 1])
im_interp = scipy.interpolate.interp1d(array[:, 0], array[:, 2])
return rl_interp(self.freq) - 1j * im_interp(self.freq)
if __name__ == "__main__":
setup = Rawdata(3, "Test.txt", 0.1, 1, 0.1, "M1", plot=True,
optimizer_options={'seed':111,
'pflag':True})
setup.run()
setup = HavriliakNegami(6, 1e12, 1e-3, 0.5, 1, 10, 5,
1e-6, 0.1, 1, 0, "M2", plot=True,
optimizer=Dual_annealing,
optimizer_options={'seed':111,
'maxiter':50})
setup.run()
setup = Jonscher(4, 1e6, 1e-5, 50, 1, 1e5, 0.7,
0.1, 1, 0.1, "M3", plot=True,
optimizer_options={'seed':111})
setup.run()
f = [0.5, 0.5]
material1 = [3, 25, 1e6]
material2 = [3, 0, 1e3]
materials = [material1, material2]
setup = Crim(2, 1*1e-1, 1e-9, 0.5, f, materials, 0.1,
1, 0, "M4", plot=True,
optimizer_options={'seed':111})
setup.run()

41
user_libs/DebyeFit/Test.txt 普通文件
查看文件

@@ -0,0 +1,41 @@
10000000,29.9213536481435,1.25169556541143
11220184.5430196,29.9010908161542,1.40299702432563
12589254.1179417,29.8756399586587,1.57217537410568
14125375.4462276,29.8436916400072,1.76117429474009
15848931.9246111,29.8036168365262,1.972079339829
17782794.1003892,29.7533952203391,2.20709802135445
19952623.1496888,29.6905309907513,2.46852367160776
22387211.3856834,29.611956285414,2.75867656083315
25118864.3150958,29.5139238169428,3.0798138913592
28183829.3126445,29.3918930811742,3.43399830398795
31622776.6016838,29.2404187150832,3.82291275991739
35481338.9233575,29.0530558254938,4.24760876499865
39810717.0553497,28.8223057349563,4.70817602888384
44668359.2150963,28.5396365567361,5.20332657814755
50118723.3627272,28.1956253280391,5.729897510637
56234132.5190349,27.7802793753847,6.28229678071347
63095734.4480193,27.283598801889,6.85194778826586
70794578.4384137,26.6964309897057,7.42683054876774
79432823.4724282,26.0116311645087,7.99126340085894
89125093.8133746,25.2254715223637,8.52610387978976
100000000,24.339136006498,9.00954486736777
112201845.430197,23.3600192468154,9.41861299029905
125892541.179417,22.302462791172,9.73132391277098
141253754.462276,21.1875720188749,9.92923323827099
158489319.246111,20.0419098344834,9.99991217790304
177827941.003892,18.8951464050085,9.93877751706085
199526231.496888,17.7770654023544,9.74979803763085
223872113.856834,16.714546101342,9.44488182455412
251188643.150958,15.7291491674702,9.04211441899954
281838293.126445,14.8357276987863,8.56331078490773
316227766.016838,14.0421664572876,8.03145189099272
354813389.233576,13.3500622522389,7.46848900052093
398107170.553497,12.7559948799797,6.89379357256362
446683592.150963,12.2530185094144,6.32331224790652
501187233.627271,11.8320839792422,5.76932819987279
562341325.190349,11.4832204965472,5.24065519659195
630957344.480194,11.1964125953287,4.74308431385235
707945784.384137,10.9621824341848,4.27993617324392
794328234.724282,10.7719266477577,3.85261757843109
891250938.133744,10.6180694756136,3.46112404221322
1000000000,10.4940904606372,3.10446192269295

385
user_libs/DebyeFit/optimization.py 普通文件
查看文件

@@ -0,0 +1,385 @@
# Author: Iraklis Giannakis, Sylwia Majchrowska
# E-mail: i.giannakis@ed.ac.uk
#
# Copyright (c) 2021 gprMax
# All rights reserved.
#
# Redistribution and use in source and binary forms are permitted
# provided that the above copyright notice and this paragraph are
# duplicated in all such forms and that any documentation,
# advertising materials, and other materials related to such
# distribution and use acknowledge that the software was developed
# as part of gprMax. The name of gprMax may not be used to
# endorse or promote products derived from this software without
# specific prior written permission.
# THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
import numpy as np
from matplotlib import pylab as plt
import scipy.optimize
from tqdm import tqdm
class Optimizer(object):
def __init__(self, maxiter=1000, seed=None):
"""
Create particle swarm optimisation object.
Args:
maxiter (int): The maximum number of iterations for the swarm
to search (Default: 1000).
seed (int): Seed for RandomState.
Must be convertible to 32 bit unsigned integers.
"""
self.maxiter = maxiter
self.seed = seed
def fit(self):
"""
Call the optimization function that tries to find an optimal set
of relaxation times that minimise the error
between the actual and the approximated electric permittivity.
"""
raise NotImplementedError()
@staticmethod
def plot(x, y):
"""
Dynamically plot the error as the optimisation takes place.
Args:
x (array): The number of current iterations.
y (array): The objective value at for all x points.
"""
plt.rcParams["axes.facecolor"] = "black"
plt.plot(x, y, "b-", linewidth=3.0)
plt.ylim(min(y) - 0.1 * min(y),
max(y) + 0.1 * max(y))
plt.xlim(min(x), max(x))
plt.grid(b=True, which="major", color="w",
linewidth=0.2, linestyle="--")
plt.suptitle("Debye fitting process")
plt.xlabel("Iteration")
plt.ylabel("Average Error")
plt.pause(0.0001)
class Particle_swarm(Optimizer):
def __init__(self, swarmsize=40, maxiter=50,
omega=0.9, phip=0.9, phig=0.9,
minstep=1e-8, minfun=1e-8,
pflag=False, seed=None):
"""
Create particle swarm optimisation object with predefined parameters.
Args:
swarmsize (int): The number of particles in the swarm (Default: 40).
maxiter (int): The maximum number of iterations for the swarm
to search (Default: 50).
omega (float): Particle velocity scaling factor (Default: 0.9).
phip (float): Scaling factor to search away from the particle's
best known position (Default: 0.9).
phig (float): Scaling factor to search away from the swarm's
best known position (Default: 0.9).
minstep (float): The minimum stepsize of swarm's best position
before the search terminates (Default: 1e-8).
minfun (float): The minimum change of swarm's best objective value
before the search terminates (Default: 1e-8)
pflag (bool): if True will plot the actual and the approximated
value during optimization process (Default: False).
"""
super(Particle_swarm, self).__init__(maxiter, seed)
self.swarmsize = swarmsize
self.omega = omega
self.phip = phip
self.phig = phig
self.minstep = minstep
self.minfun = minfun
self.pflag = pflag
def fit(self, func, lb, ub, funckwargs={}):
"""
A particle swarm optimisation that tries to find an optimal set
of relaxation times that minimise the error
between the actual and the approximated electric permittivity.
The current class is a modified edition of the pyswarm package
which can be found at https://pythonhosted.org/pyswarm/
Args:
func (function): The function to be minimized.
lb (array): The lower bounds of the design variable(s).
ub (array): The upper bounds of the design variable(s).
funckwargs (dict): Additional keyword arguments passed to
objective and constraint function
(Default: empty dict).
Returns:
g (array): The swarm's best known position (optimal design).
fg (float): The objective value at ``g``.
"""
np.random.seed(self.seed)
# check input parameters
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'
vhigh = np.abs(ub - lb)
vlow = -vhigh
# Initialize objective function
obj = lambda x: func(x=x, **funckwargs)
# Initialize the particle swarm
d = len(lb) # the number of dimensions each particle has
x = np.random.rand(self.swarmsize, d) # particle positions
v = np.zeros_like(x) # particle velocities
p = np.zeros_like(x) # best particle positions
fp = np.zeros(self.swarmsize) # best particle function values
g = [] # best swarm position
fg = np.inf # artificial best swarm position starting value
for i in range(self.swarmsize):
# Initialize the particle's position
x[i, :] = lb + x[i, :] * (ub - lb)
# Initialize the particle's best known position
p[i, :] = x[i, :]
# Calculate the objective's value at the current particle's
fp[i] = obj(p[i, :])
# At the start, there may not be any feasible starting point,
# so just
# give it a temporary "best" point since it's likely to change
if i == 0:
g = p[0, :].copy()
# If the current particle's position is better than the swarm's,
# update the best swarm position
if fp[i] < fg:
fg = fp[i]
g = p[i, :].copy()
# Initialize the particle's velocity
v[i, :] = vlow + np.random.rand(d) * (vhigh - vlow)
# Iterate until termination criterion met
for it in tqdm(range(self.maxiter), desc='Debye fitting'):
rp = np.random.uniform(size=(self.swarmsize, d))
rg = np.random.uniform(size=(self.swarmsize, d))
for i in range(self.swarmsize):
# Update the particle's velocity
v[i, :] = self.omega * v[i, :] + self.phip * rp[i, :] * \
(p[i, :] - x[i, :]) + \
self.phig * rg[i, :] * (g - x[i, :])
# Update the particle's position,
# correcting lower and upper bound
# violations, then update the objective function value
x[i, :] = x[i, :] + v[i, :]
mark1 = x[i, :] < lb
mark2 = x[i, :] > ub
x[i, mark1] = lb[mark1]
x[i, mark2] = ub[mark2]
fx = obj(x[i, :])
# Compare particle's best position
if fx < fp[i]:
p[i, :] = x[i, :].copy()
fp[i] = fx
# Compare swarm's best position to current
# particle's position
if fx < fg:
tmp = x[i, :].copy()
stepsize = np.sqrt(np.sum((g - tmp) ** 2))
if np.abs(fg - fx) <= self.minfun:
print(f'Stopping search: Swarm best objective '
f'change less than {self.minfun}')
return tmp, fx
elif stepsize <= self.minstep:
print(f'Stopping search: Swarm best position '
f'change less than {self.minstep}')
return tmp, fx
else:
g = tmp.copy()
fg = fx
# Dynamically plot the error as the optimisation takes place
if self.pflag:
if it == 0:
xpp = [it]
ypp = [fg]
else:
xpp.append(it)
ypp.append(fg)
Particle_swarm.plot(xpp, ypp)
return g, fg
class Dual_annealing(Optimizer):
def __init__(self, maxiter=100,
local_search_options={}, initial_temp=5230.0,
restart_temp_ratio=2e-05, visit=2.62, accept=- 5.0,
maxfun=1e7, no_local_search=False,
callback=None, x0=None, seed=None):
"""
Create dual annealing object with predefined parameters.
Args:
maxiter (int): The maximum number of iterations for the swarm
to search (Default: 100).
local_search_options (dict): Extra keyword arguments to be passed
to the local minimizer, reffer to
scipy.optimize.minimize() function
(Default: empty dict).
initial_temp (float): The initial temperature, use higher values to
facilitates a wider search of the energy
landscape, allowing dual_annealing to escape
local minima that it is trapped in.
Range is (0.01, 5.e4] (Default: 5230).
restart_temp_ratio (float): During the annealing process,
temperature is decreasing, when it
reaches initial_temp * restart_temp_ratio,
the reannealing process is triggered.
Range is (0, 1) (Default: 2e-5).
visit (float): Parameter for visiting distribution. The value range is (1, 3]
(Default: 2.62).
accept (float): Parameter for acceptance distribution. It is used to control
the probability of acceptance. The lower the acceptance parameter,
the smaller the probability of acceptance. The value range (-1e4, -5]
(Default: -5.0).
no_local_search (bool):
maxfun (int): Soft limit for the number of objective function calls.
(Default: 1e7).
callback (callable): A callback function with signature callback(x, f, context),
which will be called for all minima found.
x and f are the coordinates and function value of
the latest minimum found, and context has value in [0, 1, 2],
with the following meaning:
0: minimum detected in the annealing process.
1: detection occurred in the local search process.
2: detection done in the dual annealing process.
If the callback implementation returns True,
the algorithm will stop.
x0 (ndarray): Coordinates of a single N-D starting point, shape(n,).
(Default: None).
seed (None, int): Specify seed for repeatable minimizations.
The random numbers generated with this seed only
affect the visiting distribution function and
new coordinates generation (Default: None).
pflag (bool): if True will plot the actual and the approximated
value during optimization process (Default: False).
"""
super(Dual_annealing, self).__init__(maxiter, seed)
self.local_search_options = local_search_options
self.initial_temp = initial_temp
self.restart_temp_ratio = restart_temp_ratio
self.visit = visit
self.accept = accept
self.maxfun = maxfun
self.no_local_search = no_local_search
self.callback = callback
self.x0 = x0
#self.pflag = pflag
def fit(self, func, lb, ub, funckwargs={}):
"""
Find the global minimum of a function using Dual Annealing.
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/
scipy.optimize.dual_annealing.html#scipy.optimize.dual_annealing
Args:
func (function): The function to be minimized
lb (array): The lower bounds of the design variable(s)
ub (array): The upper bounds of the design variable(s)
funckwargs (dict): Additional keyword arguments passed to
objective and constraint function
(Default: empty dict)
Returns:
g (array): The solution array (optimal design).
fg (float): The objective value at the 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)
return result.x, result.fun
def cost_function(x, rl_g, im_g, freq_g):
"""
The cost function is the average error between
the actual and the approximated electric permittivity.
Returns:
cost: The final error
"""
cost1, cost2, _, _, _, _ = linear(rl_g, im_g, x, freq_g)
cost = cost1 + cost2
return cost
def linear(rl, im, logt, freq):
"""
Returns:
cost1: Error (?)
cost2: Error (?)
x: Resulting optimised weights for the given relaxation times
ee: Average error between the actual and the approximated real part
rp: The real part of the permittivity for the optimised relaxation
times and weights for the frequnecies included in freq
ip: The imaginary part of the permittivity for the optimised
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
# Here they are transformed back t0=10**logt
tt = [10**logt[i] for i in range(0, len(logt))]
# y = Ax , here the A matrix for the real and the imaginary part is builded
d_r = np.array(
[[calc([1, 1, 0, 1, tt[i]], [freq[j]])[0]
for i in range(0, len(tt))] for j in
range(0, len(freq))])
d = np.array(
[[calc([1, 1, 0, 1, tt[i]], [freq[j]])[1]
for i in range(0, len(tt))] for j in
range(0, len(freq))])
# Adding dumping (Marquart least squares)
# Solving the overdetermined system y=Ax
x = np.abs(np.linalg.lstsq(d, im)[0])
mx, my, my2 = np.matrix(x), np.matrix(d), np.matrix(d_r)
rp, ip = my2 * np.transpose(mx), my * np.transpose(mx)
cost1 = np.sum([np.abs(ip[i]-im[i]) for i in range(0, len(im))])/len(im)
ee = (np.mean(rl - rp))
if ee < 1:
ee = 1
cost2 = np.sum([np.abs(rp[i] - rl[i] + ee)
for i in range(0, len(im))])/len(im)
return cost1, cost2, x, ee, rp, ip
def calc(cal_inputs, freq):
# Calculates the Havriliak-Negami function for the given cal_inputs
q = [cal_inputs[2] + cal_inputs[3] / (np.array(1 + np.array(
1j * 2 * np.pi * f * cal_inputs[4]) ** cal_inputs[0]
) ** cal_inputs[1]) for f in freq]
# Return the real and the imaginary part of the relaxation function
if len(q) > 1:
rl = [q[i].real for i in range(0, len(q))]
im = [q[i].imag for i in range(0, len(q))]
else:
rl = q[0].real
im = q[0].imag
return rl, im