# pade.py
# Andrew Davison, ad@coe.psu.ac.th, June 2026
"""
Obtain a Pade approximation for a function f(x), 
using its Maclaurin series coefficients.

A list of functions is presented, and the chosen one is 
given a Pade approximation based on the specified (n,m) 
order for that function.

The approx is evaluated against a series of values, and 
the results are compared to those for the actual function. 
Also, a graph plots the approximation against the 
original function.

The Pade polynomial's algebraic form is printed to stdout, 
and also printed using Latex on the plot.
"""

import math
from fractions import Fraction
import faUtils

import matplotlib.pyplot as plt
import sympy as sp


def errReport(fName, func, p, q, testVals, x0=0):
  print(f"\n{fName} Comparisons:")
  print("      x        Actual         Pade"
        "         Abs Err       Rel Err" )
  print( "-"*70)

  tot = 0.0
  validCount = 0
  for xVal in testVals:
    try:
      actual = func(xVal)
      approx = float( faUtils.evalPade(xVal, p, q, x0))
    except (ValueError, ZeroDivisionError, TypeError):
      print(f"{xVal:8.3f}   (domain error, skipped)")
      continue

    absErr = abs(approx - actual)
    relErr = abs((approx - actual) / actual) \
                           if actual != 0 else 0.0
    tot += relErr
    validCount += 1
    print(f"{xVal:8.3f}   {actual:12.8f}   {approx:12.8f}   "
            f"{absErr:12.8f}  {relErr:12.8f}")

  if validCount == 0:
    print("\nNo valid test points.")
    return None
  avg = tot / validCount
  print(f"\nAverage Relative Error = {avg:.8e}")
  return avg



def plotPade(fName, func, latexFnm,  p, q, xMin, xMax, nPts=400):
  """
  Plot func against its Padé approximation over [xMin, xMax].
  """
  step = (xMax-xMin)/(nPts-1)
  xs = [xMin + i*step for i in range(nPts)]

  ysFn = []
  ysPade  = []
  xsValid = []
  for x in xs:
    try:
      ye = func(x)
      yp = float( faUtils.evalPade(x, p, q))
    except (ValueError, ZeroDivisionError, TypeError):
      continue
    xsValid.append(x)
    ysFn.append(ye)
    ysPade.append(yp)
 
  if not xsValid:
    print("plotPade: no valid points to plot.")
    return

  plt.figure(figsize=(8, 5))
  plt.plot(xsValid, ysFn, label='Exact', lw=2, color='steelblue')
  plt.plot(xsValid, ysPade, lw=2,
           label=rf"Pade = $" + latexFnm + "$",
           color='tomato', linestyle='--')
  plt.xlabel('$x$')
  plt.ylabel(fName)
  plt.title(f'Padé Approx  (n={len(p)-1}, m={len(q)-1})'
            + (f' for {fName}'))

  leg = plt.legend()
  legendLabels = leg.get_texts()
  legendLabels[1].set_fontsize(14)  # to make the latex bigger

  plt.grid(True, alpha=0.3)
  plt.tight_layout()
  plt.show()


# ---------------------------------
x = sp.symbols('x')

fNames = ["sin(x)", "exp(x)", "exp(-x)",
          "1-0.5*ln(1+x^2)", "tan(x)" ]

spFns = [sp.sin, sp.exp, lambda x: sp.exp(-x), 
         lambda x: 1-0.5*sp.ln(1+x*x), sp.tan ]

# (n, m)
orders = [ (3,4), (2,1), (3,2), (8,8), (4,4)]

fns = [math.sin, math.exp, lambda x: math.exp(-x), 
       lambda x: 1-0.5*math.log(1+x*x), math.tan ]

testVals = [-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0]


print("Functions:", fNames)
i = int(input("Enter index of function name: "))
n, m = orders[i]
print(f"\nApproximating {fNames[i]} for order {orders[i]}")

p, q, latexFnm = faUtils.padeApprox(spFns[i], n, m)

errReport(fNames[i], fns[i], p, q, testVals)

plotPade(fNames[i], fns[i], latexFnm, p, q, 
            xMin=min(testVals), xMax=max(testVals))