# this file is modified from CP2K
import numpy as np
#from scipy.optimize import lsq_linear
[docs]
def get_l_sq_wghts_cos_tf_t_to_w(tau_tj, omega_tj, E_min, E_max, num_points_per_magnitude=200):
def _calc_max_error_fit_tau_grid_with_cosine(omega, tau_tj, tau_wj_work, x_values, y_values):
func_val = tau_wj_work[None,:] * np.cos(omega * tau_tj[None,:]) * np.exp(-x_values[:,None] * tau_tj[None,:])
max_error_tmp = np.max(np.abs(y_values - np.sum(func_val, axis=1)))
return max_error_tmp
num_integ_points = tau_tj.shape[0]
num_x_nodes = (int(np.log10(E_max / E_min)) + 1) * num_points_per_magnitude
num_x_nodes = max(num_x_nodes, num_integ_points)
# set the x-values logarithmically in the interval [Emin,Emax]
multiplicator = np.power(E_max / E_min, 1.0 / (np.double(num_x_nodes) - 1.0))
x_values = E_min * np.power(multiplicator, np.arange(num_x_nodes))
# loop over all omega frequency points
weights_cos_tf_t_to_w = np.zeros(shape=[num_integ_points, num_integ_points], dtype=np.double)
max_errors = np.zeros(shape=[num_integ_points], dtype=np.double)
for jquad, omega in enumerate(omega_tj):
# y=2x/(x^2+omega_k^2)
y_values = 2.0 * x_values / (np.square(x_values) + np.square(omega))
# calculate mat_A
mat_A = np.cos(omega * tau_tj[None,:]) * np.exp(-x_values[:, None] * tau_tj[None, :])
# least square of Ax=y by svd
tau_wj_work = _get_least_square_svd(mat_A, y_values)
#tau_wj_work = lsq_linear(mat_A, y_values, lsq_solver="exact", verbose=0).x
weights_cos_tf_t_to_w[jquad, :] = tau_wj_work
max_errors[jquad] = _calc_max_error_fit_tau_grid_with_cosine(omega, tau_tj, tau_wj_work, x_values, y_values)
max_error = np.max(max_errors)
#print("wghts_cos_tf_t_to_w max error=%-15.8e" % max_error)
return weights_cos_tf_t_to_w, max_error
[docs]
def get_l_sq_wghts_cos_tf_w_to_t(tau_tj, omega_tj, E_min, E_max, num_points_per_magnitude=200):
num_integ_points = tau_tj.shape[0]
num_x_nodes = (int(np.log10(E_max / E_min)) + 1) * num_points_per_magnitude
num_x_nodes = max(num_x_nodes, num_integ_points)
# set the x-values logarithmically in the interval [Emin,Emax]
multiplicator = np.power(E_max / E_min, (1.0 / (np.double(num_x_nodes) - 1.0)))
x_values = E_min * np.power(multiplicator, np.arange(num_x_nodes))
def _calc_max_error_fit_omega_grid_with_cosine(tau, omega_tj, omega_wj_work, x_values, y_values):
func_val = omega_wj_work[None,:] * np.cos(tau*omega_tj[None,:]) * 2.0 * x_values[:,None] \
/ (np.square(x_values[:,None]) + np.square(omega_tj[None,:]))
max_error_tmp = np.max(np.abs(y_values - np.sum(func_val, axis=1)))
return max_error_tmp
weights_cos_tf_w_to_t = np.zeros(shape=[num_integ_points, num_integ_points], dtype=np.double)
max_errors = np.zeros(shape=[num_integ_points], dtype=np.double)
for jquad, tau in enumerate(tau_tj):
# y=exp(-x*|tau_k|)
y_values = np.exp(-x_values * tau)
# calculate mat_A
mat_A = np.cos(tau * omega_tj[None,:]) * 2.0 * x_values[:,None] \
/ (np.square(x_values[:,None]) + np.square(omega_tj[None,:]))
# least square of Ax=y by svd
omega_wj_work = _get_least_square_svd(mat_A, y_values)
#omega_wj_work = lsq_linear(mat_A, y_values, lsq_solver="exact", verbose=0).x
weights_cos_tf_w_to_t[jquad, :] = omega_wj_work
max_errors[jquad] = _calc_max_error_fit_omega_grid_with_cosine(tau, omega_tj, omega_wj_work, x_values, y_values)
max_error = np.max(max_errors)
#print("wghts_cos_tf_w_to_t max error=%-15.8e" % max_error)
return weights_cos_tf_w_to_t, max_error
[docs]
def get_l_sq_wghts_sin_tf_t_to_w(tau_tj, omega_tj, E_min, E_max, num_points_per_magnitude=200):
num_integ_points = tau_tj.shape[0]
num_x_nodes = (int(np.log10(E_max / E_min)) + 1) * num_points_per_magnitude
num_x_nodes = max(num_x_nodes, num_integ_points)
# set the x-values logarithmically in the interval [Emin,Emax]
multiplicator = np.power(E_max / E_min, 1.0 / (np.double(num_x_nodes) - 1.0))
x_values = E_min * np.power(multiplicator, np.arange(num_x_nodes))
def _calc_max_error_fit_tau_grid_with_sine(omega, tau_tj, tau_wj_work, x_values, y_values):
func_val = tau_wj_work[None,:] * np.sin(omega * tau_tj[None,:]) * np.exp(-x_values[:,None] * tau_tj[None,:])
max_error_tmp = np.max(np.abs(y_values - np.sum(func_val, axis=1)))
return max_error_tmp
weights_sin_tf_t_to_w = np.zeros(shape=[num_integ_points, num_integ_points], dtype=np.double)
max_errors = np.zeros(shape=[num_integ_points], dtype=np.double)
for jquad, omega in enumerate(omega_tj):
# y=2*omega/(x^2+omega_k^2)
y_values = 2.0 * omega / (np.square(x_values) + np.square(omega))
# calculate mat_A
mat_A = np.sin(omega * tau_tj[None,:]) * np.exp(-x_values[:,None] * tau_tj[None,:])
# least square of Ax=y by svd
tau_wj_work = _get_least_square_svd(mat_A, y_values)
#tau_wj_work = lsq_linear(mat_A, y_values, lsq_solver="exact", verbose=0).x
weights_sin_tf_t_to_w[jquad, :] = tau_wj_work
max_errors[jquad] = _calc_max_error_fit_tau_grid_with_sine(omega, tau_tj, tau_wj_work, x_values, y_values)
max_error = np.max(max_errors)
#print("wghts_sin_tf_t_to_w max error=%-15.8e" % max_error)
return weights_sin_tf_t_to_w, max_error
[docs]
def get_l_sq_wghts_sin_tf_w_to_t(tau_tj, omega_tj, E_min, E_max, num_points_per_magnitude=200):
num_integ_points = tau_tj.shape[0]
num_x_nodes = (int(np.log10(E_max / E_min)) + 1) * num_points_per_magnitude
num_x_nodes = max(num_x_nodes, num_integ_points)
# set the x-values logarithmically in the interval [Emin,Emax]
multiplicator = np.power(E_max / E_min, 1.0 / (num_x_nodes - 1.0))
x_values = E_min * np.power(multiplicator, np.arange(num_x_nodes))
def _calc_max_error_fit_omega_grid_with_sine(tau, omega_tj, omega_wj_work, x_values, y_values):
func_val = omega_wj_work[None,:] * np.sin(tau*omega_tj[None,:]) * 2.0 * omega_tj[None,:] \
/ (np.square(x_values[:,None]) + np.square(omega_tj[None,:]))
max_error_tmp = np.max(np.abs(y_values - np.sum(func_val, axis=1)))
return max_error_tmp
weights_sin_tf_w_to_t = np.zeros(shape=[num_integ_points, num_integ_points], dtype=np.double)
max_errors = np.zeros(shape=[num_integ_points], dtype=np.double)
for jquad, tau in enumerate(tau_tj):
# y=exp(-x*|tau_k|)
y_values = np.exp(-x_values * tau)
# calculate mat_A
mat_A = np.sin(tau*omega_tj[None,:]) * 2.0 * omega_tj[None,:] / (np.square(x_values[:,None]) + np.square(omega_tj[None,:]))
# least square of Ax=y by svd
omega_wj_work = _get_least_square_svd(mat_A, y_values)
#omega_wj_work = lsq_linear(mat_A, y_values, lsq_solver="exact", verbose=0).x
weights_sin_tf_w_to_t[jquad, :] = omega_wj_work
max_errors[jquad] = _calc_max_error_fit_omega_grid_with_sine(tau, omega_tj, omega_wj_work, x_values, y_values)
max_error = np.max(max_errors)
#print("wghts_sin_tf_w_to_t max error=%-15.8e" % max_error)
return weights_sin_tf_w_to_t, max_error
def _get_least_square_svd(mat_A, y_values):
num_x_nodes, num_integ_points = mat_A.shape
# Singular value decomposition of mat_A
mat_U, sing_values, mat_SinvVSinvT = np.linalg.svd(mat_A)
# integration weights = V Sigma U^T y
# 1) V*Sigma
mat_SinvVSinvSigma = np.zeros(shape=[num_integ_points, num_x_nodes], dtype=np.double)
mat_SinvVSinvSigma[:num_integ_points, :num_integ_points] = mat_SinvVSinvT.T / sing_values
# 2) U^T y
vec_UTy = np.matmul(mat_U.T, y_values)
# 3) (V*Sigma) * (U^T y)
wj_work = np.matmul(mat_SinvVSinvSigma, vec_UTy)
return wj_work