mpmath_special_functions_test_generator.py 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172
  1. #!/usr/bin/env python3
  2. import mpmath as mp
  3. import mpmath_riccati_bessel as mrb
  4. import mpmath_input_arguments as mia
  5. import os.path
  6. class TestData:
  7. def __init__(self, list_to_parse, filetype):
  8. self.filetype = filetype
  9. if self.filetype == 'c++':
  10. self.cpp_parse(list_to_parse)
  11. else:
  12. raise NotImplementedError("Only C++ files *.hpp parsing was implemented")
  13. def cpp_parse(self, list_to_parse):
  14. self.comment = list_to_parse[0]
  15. if self.comment[:2] != '//': raise ValueError('Not a comment')
  16. self.typeline = list_to_parse[1]
  17. if 'std::vector' not in self.typeline: raise ValueError('Unexpected C++ container')
  18. self.testname = list_to_parse[2]
  19. self.opening = list_to_parse[3]
  20. if self.opening != '= {': raise ValueError('For C++ we expect opeing with = {');
  21. self.ending = list_to_parse[-1]
  22. if self.ending != '};': raise ValueError('For C++ we expect closing };')
  23. self.evaluated_data = list_to_parse[4:-1]
  24. def get_string(self):
  25. out_sting = self.comment + '\n' + self.typeline + '\n' + self.testname + '\n' + self.opening + '\n'
  26. for result in self.evaluated_data:
  27. out_sting += result + '\n'
  28. out_sting += self.ending + '\n'
  29. return out_sting
  30. class UpdateSpecialFunctionsEvaluations:
  31. def __init__(self, filename='default_out.hpp', complex_arguments=[],
  32. output_dps=16, max_num_elements_of_nlist=51):
  33. self.evaluated_data = []
  34. self.test_setup = []
  35. self.filename = filename
  36. self.read_evaluated_data()
  37. self.complex_arguments = complex_arguments
  38. self.output_dps = output_dps
  39. self.max_num_elements_of_nlist = max_num_elements_of_nlist
  40. def read_evaluated_data(self):
  41. self.filetype = 'undefined'
  42. if self.filename.endswith('.hpp'):
  43. self.filetype = 'c++'
  44. if self.filename.endswith('.f90'):
  45. self.filetype = 'fortran'
  46. if not os.path.exists(self.filename):
  47. print("WARNING! Found no data file:", self.filename)
  48. return
  49. with open(self.filename, 'r') as in_file:
  50. content = in_file.readlines()
  51. content = [x.strip() for x in content]
  52. while '' in content:
  53. record_end_index = content.index('')
  54. new_record = content[:record_end_index]
  55. content = content[record_end_index + 1:]
  56. self.add_record(new_record)
  57. self.add_record(content)
  58. def add_record(self, new_record):
  59. if len(new_record) == 0: return
  60. if len(new_record) < 6: raise ValueError('Not enough lines in record:', new_record)
  61. self.evaluated_data.append(TestData(new_record, self.filetype))
  62. def get_file_content(self):
  63. self.evaluated_data.sort(key=lambda x: x.testname)#, reverse=True)
  64. out_string = ''
  65. for record in self.evaluated_data:
  66. out_string += record.get_string() + '\n'
  67. return out_string[:-1]
  68. def remove(self, testname):
  69. for i, result in enumerate(self.evaluated_data):
  70. if result.testname == testname:
  71. del self.evaluated_data[i]
  72. def get_n_list(self, z, max_number_of_elements=10):
  73. nmax = mrb.LeRu_cutoff(z)
  74. factor = nmax ** (1 / (max_number_of_elements - 2))
  75. n_list = [int(factor ** i) for i in range(max_number_of_elements - 1)]
  76. n_list.append(0)
  77. n_set = set(n_list)
  78. return sorted(n_set)
  79. def get_test_data_nlist(self, z_record, output_dps, n, func):
  80. x = str(z_record[0])
  81. mr = str(z_record[1][0])
  82. mi = str(z_record[1][1])
  83. z_str = ''
  84. try:
  85. z = mp.mpf(x) * mp.mpc(mr, mi)
  86. D1nz = func(n, z)
  87. z_str = ('{{' +
  88. mp.nstr(z.real, output_dps * 2) + ',' +
  89. mp.nstr(z.imag, output_dps * 2) + '},' +
  90. str(n) + ',{' +
  91. mp.nstr(D1nz.real, output_dps) + ',' +
  92. mp.nstr(D1nz.imag, output_dps) + '},' +
  93. mp.nstr(mp.fabs(D1nz.real * 10 ** -output_dps), 2) + ',' +
  94. mp.nstr(mp.fabs(D1nz.imag * 10 ** -output_dps), 2) +
  95. '},')
  96. except:
  97. pass
  98. return z_str
  99. def get_test_data(self, Du_test, output_dps, max_num_elements_of_n_list, func, funcname):
  100. output_list = ['// complex(z), n, complex(f(n,z)), abs_err_real, abs_err_imag',
  101. 'std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >',
  102. str(funcname)+'_test_' + str(output_dps) + 'digits','= {']
  103. for z_record in Du_test:
  104. x = str(z_record[0])
  105. mr = str(z_record[1][0])
  106. mi = str(z_record[1][1])
  107. mp.mp.dps = 20
  108. z = mp.mpf(x) * mp.mpc(mr, mi)
  109. n_list = self.get_n_list(z, max_num_elements_of_n_list)
  110. if z_record[4] == 'Yang': n_list = [0,1,30,50,60,70,75,80,85,90,99,116,130]
  111. print(z, n_list)
  112. failed_evaluations = 0
  113. for n in n_list:
  114. mp.mp.dps = 20
  115. old_z_string = self.get_test_data_nlist(z_record, output_dps, n, func)
  116. mp.mp.dps = 37
  117. new_z_string = self.get_test_data_nlist(z_record, output_dps, n, func)
  118. while old_z_string != new_z_string:
  119. new_dps = int(mp.mp.dps * 1.41)
  120. if new_dps > 300: break
  121. mp.mp.dps = new_dps
  122. print("New dps = ", mp.mp.dps, 'n =', n, ' (max ',n_list[-1],') for z =', z, ' ', end='')
  123. old_z_string = new_z_string
  124. new_z_string = self.get_test_data_nlist(z_record, output_dps, n, func)
  125. if new_z_string != '':
  126. output_list.append(new_z_string)
  127. else:
  128. failed_evaluations += 1
  129. # break
  130. print("\nFailed evaluations ", failed_evaluations, ' of ', len(n_list))
  131. output_list.append('};')
  132. return output_list
  133. def run_test(self, func, funcname):
  134. out_list_result = self.get_test_data(mia.complex_arguments, self.output_dps,
  135. self.max_num_elements_of_nlist,
  136. func, funcname)
  137. testname = str(funcname)+'_test_' + str(self.output_dps) + 'digits'
  138. self.remove(testname)
  139. self.add_record(out_list_result)
  140. def main():
  141. sf_evals = UpdateSpecialFunctionsEvaluations(filename='test_spec_functions_data.hpp',
  142. complex_arguments=mia.complex_arguments,
  143. output_dps=16, max_num_elements_of_nlist=51)
  144. # output_dps=5, max_num_elements_of_nlist=3)
  145. # sf_evals.run_test(mrb.D1, 'D1')
  146. sf_evals.run_test(mrb.D2, 'D2')
  147. # sf_evals.run_test(mrb.D3, 'D3')
  148. # sf_evals.run_test(mrb.psi_div_ksi, 'psi_div_ksi')
  149. # sf_evals.run_test(mrb.psi_div_xi, 'psi_div_xi')
  150. with open(sf_evals.filename, 'w') as out_file:
  151. out_file.write(sf_evals.get_file_content())
  152. main()