import numpy as np from time import time import matplotlib.pyplot as plt from fractions import Fraction import sys sys.path.append('../.') import odeg as od # from matplotlib import rc plt.rcParams['figure.dpi'] = 300 plt.rcParams['savefig.dpi'] = 300 # rc('text', usetex=True) def app_rel_err(fill_ratio, **cfg): err_max = 0 for pmax in (2, 3, 4): cfg['p_max'] = pmax cfg['pos_charge'] = round(fill_ratio * 2 * (2*pmax+1), 5) app_f = od.freeEnergy(**cfg) cfg['p_max'] = 5 true_f = od.freeEnergy(**cfg) err_max = max(err_max, abs((app_f-true_f) / true_f)) return err_max def err_repr(val, rel_err): abs_err = abs(val * rel_err) # print("ae1:", abs_err) digit = int(np.floor(np.log10(abs_err))) abs_err = round(abs_err, -digit) # print("ae2:", abs_err) digit = int(np.floor(np.log10(abs_err))) # in case abs_err was rounded up from .xx98 to .x10 abs_err_leading_digit = round(abs_err * 10**(-digit)) val_rounded = round(val, -digit) # if digit < 0: # res = (f"%.{-digit}f" % val_rounded) + r" $\pm$ " + (f"%.{-digit}f" % abs_err) # else: # res = ("%.f" % val_rounded) + r" $\pm$ " + ("%f" % abs_err) if digit < 0: res = f"%.{-digit}f(%d)" % (val_rounded, abs_err_leading_digit) else: res = "%f(%d)" % (val_rounded, abs_err_leading_digit) # print("val, rel_err = ", val, rel_err) # print("output:", res) return "$" + res + "$" if __name__ == '__main__': r0_grid = np.array([0.5, 1.]) rs_grid = np.array([2., 4., 8.]) theta_grid = np.array([.125, .5, 1.]) cfg = { 'p_max': 5, 'pos_charge': 5, 'path': "../results", } fill_ratio = cfg['pos_charge'] / (2 * (2*cfg['p_max']+1)) free_energy_per_pc = np.zeros((r0_grid.size, rs_grid.size, theta_grid.size)) rel_error = np.zeros(free_energy_per_pc.shape) for r0_ind, rs_ind, theta_ind in np.ndindex(r0_grid.size, rs_grid.size, theta_grid.size): cfg['r0'] = r0_grid[r0_ind] cfg['rs'] = rs_grid[rs_ind] cfg['theta'] = theta_grid[theta_ind] free_energy_per_pc[r0_ind, rs_ind, theta_ind] = od.freeEnergy(**cfg) / cfg['pos_charge'] rel_error[r0_ind, rs_ind, theta_ind] = app_rel_err(fill_ratio, **cfg) print("\n\ncfg:", cfg) print("fe:", free_energy_per_pc[r0_ind, rs_ind, theta_ind]) print("err:", rel_error[r0_ind, rs_ind, theta_ind]) file = open("table_f.tex", "w") for theta_ind in range(theta_grid.size): line = r"$\Theta = " + str(Fraction(theta_grid[theta_ind])) + "$" for r0_ind in range(r0_grid.size): for rs_ind in range(rs_grid.size): line += " & " + err_repr(free_energy_per_pc[r0_ind, rs_ind, theta_ind], rel_error[r0_ind, rs_ind, theta_ind]) if theta_ind+1 < theta_grid.size: line += r"\\" + "\n" file.write(line) file.close()