1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#***********************************************************************
# This file is part of OpenMolcas. *
# *
# OpenMolcas is free software; you can redistribute it and/or modify *
# it under the terms of the GNU Lesser General Public License, v. 2.1. *
# OpenMolcas is distributed in the hope that it will be useful, but it *
# is provided "as is" and without any express or implied warranties. *
# For more details see the full text of the license in the file *
# LICENSE or in <http://www.gnu.org/licenses/>. *
# *
# Copyright (C) 2017, Ignacio Fdez. Galván *
#***********************************************************************
# Test the accuracy of the fitting polynomials generated by the genfit.py script.
# The roots and weights are computed at the "data" points and in the midpoints.
from mpmath import mp
from rys_aux import *
import sys
# Set floating point precision
# Single precision (32 bits): 24-bit significand, 6 decimal places
# Double precision (64 bits): 53-bit significand, 15 decimal places
# Quadruple precision (128 bits): 113-bit significand, 33 decimal places
mp.dps = 33
# Read the data from stdin, expected to be a file formated
# as a Molcas "rysrw" file
# Maximum number of roots and degree of fitting polynomials
nmax, fit_degree = (int(i) for i in sys.stdin.readline().split())
# Requested accuracy (not needed)
accuracy = []
while (len(accuracy) < nmax):
accuracy.extend([mp.mpf(i) for i in sys.stdin.readline().split()])
# Value of x from which the asymptotic approximation is accurate enough (not needed)
xmax = []
while (len(xmax) < nmax):
xmax.extend([mp.mpf(i) for i in sys.stdin.readline().split()])
# Step for the x values
dx = []
while (len(dx) < nmax):
dx.extend([mp.mpf(i) for i in sys.stdin.readline().split()])
# Number of x values
nx = []
while (len(nx) < nmax):
nx.extend([int(i) for i in sys.stdin.readline().split()])
# Number of fitting intervals
npol = []
while (len(npol) < nmax):
npol.extend([int(i) for i in sys.stdin.readline().split()])
# Map for which fitting interval to use for each x value
# (values are 1-based, so subtract 1)
polmap = [[] for i in range(nmax)]
# And center of each interval (offset value for x)
center = [[] for i in range(nmax)]
for n in range(nmax):
while (len(polmap[n]) < nx[n]):
polmap[n].extend([int(i)-1 for i in sys.stdin.readline().split()])
while (len(center[n]) < npol[n]):
center[n].extend([mp.mpf(i) for i in sys.stdin.readline().split()])
# Coefficients for all roots and weights and intervals
coeff = [[] for i in range(nmax)]
for n in range(nmax):
# First read a single list
tmp = []
while (len(tmp) < npol[n]*(n+1)*(fit_degree+1)*2):
tmp.extend([mp.mpf(i) for i in sys.stdin.readline().split()])
# Then arrange it, the nesting order in the file is:
# roots, weights
# 0th...6th order coefficients
# root 1...n
# interval 1...npol
m = 0
for j in range(2):
coeff[n].append([])
for k in range(fit_degree+1):
coeff[n][j].append([])
for l in range(n+1):
coeff[n][j][k].append([])
coeff[n][j][k][l] = tmp[m:m+npol[n]]
m += npol[n]
# Now it's time to actually test the fittings for each order of Rys polynomials
for n in range(nmax):
max_error = [zero, zero]
# Test all the "data points" (x values) and mid-points too
for ix in range(2*nx[n]-1):
x = ix*dx[n]/2
# Compute roots and weights
roots, weights = rysroots(n+1, x)
# Select approximating polynomial
ipol = polmap[n][int(mp.floor(x/dx[n]))]
xoff = x-center[n][ipol]
# For each root and weight calculate the value from the fitting polynomial
roots_fit = []
weights_fit = []
for i in range(n+1):
pol = []
for j in range(fit_degree+1):
pol.insert(0, coeff[n][0][j][i][ipol])
roots_fit.append(mp.polyval(pol, xoff))
pol = []
for j in range(fit_degree+1):
pol.insert(0, coeff[n][1][j][i][ipol])
weights_fit.append(mp.polyval(pol, xoff))
# Obtain the maximum relative error
error_r = max([abs((a-b)/b) for a,b in zip(roots_fit, roots)])
error_w = max([abs((a-b)/b) for a,b in zip(weights_fit, weights)])
error = max(error_r, error_w)
# Keep separate track of maximum error at data points and mid-points
i = ix%2
max_error[i] = max(max_error[i], error)
print(n+1, figs(max_error[0], exp=True), figs(max_error[1], exp=True))
|
