"""
There are three primary ways to write and execute Python code, each with its
own strengths.

1. The Python Console
-----------------------------------------------------------
This is the interactive prompt you get by typing `python` in your terminal.
- WHEN TO USE: For quick, simple, one-off tasks. Think of it as a powerful
  calculator. 
- LIMITATIONS: It has no memory between sessions and is not suitable for
  writing programs of more than a few lines.

2. Jupyter Notebooks (.ipynb files)
-----------------------------------
Jupyter provides an interactive, cell-based environment in your web browser.
- WHEN TO USE: For exploration, data analysis, and (some) scientific computing.
  The ability to mix code, formatted text (Markdown), mathematical equations,
  and plots in one document makes it good for developing ideas, visualising
  data, and creating reports. Running code in independent cells can be useful
  for experimentation.
- LIMITATIONS: The out-of-order execution of cells can make notebooks
  difficult to debug and reproduce. They are not ideal for building complex,
  reusable applications.

3. Python Scripts (.py files)
-----------------------------
These are plain text files containing Python code, intended to be executed
from top to bottom.
- WHEN TO USE: This is the standard for building applications, tools, and
  reusable libraries. Scripts encourage organised, modular code. You write
  functions and classes in modules (`.py` files) that can be imported into
  other scripts. This is the most robust and shareable way to write code.
- LIMITATIONS: Less interactive than a notebook for initial exploration. You
  typically write the code, save it, and then run the entire script.
"""

# This script imports and uses the functions from myfunc.py

# The 'from <filename> import <function_name>' syntax is key.
# Python looks for a file named myfunc.py in the same directory.
# Importing runs the entire file, so any top-level code (like print statements)
# will execute immediately.
from myfuncs import myfunc, gaussian, polynomial
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize, integrate
from scipy.optimize import curve_fit

# Now we can use myfunc as if it were defined in this file.
print("Calculating myfunc for a few values:")

for i in range(5):
    x_value = -2 + i  # Generate some values from 0 to 1
    result = myfunc(x_value)
    print(f"  myfunc({x_value:.2f}) = {result:.4f}")


print("\n==========================")
print("SciPy Examples with myfunc")
print("==========================")


# 1) Find x such that myfunc(x)=0
print("\n1) Root finding with optimize.root_scalar")
root_res = optimize.root_scalar(myfunc, bracket=(-2.0, 1.0), method="brentq")
if root_res.converged:
    x_root = root_res.root
    print(f"   Root in (0, 1): x ≈ {x_root:.8f}, myfunc(x) ≈ {myfunc(x_root):.3e}")
else:
    print("   Root finding did not converge.")


# 2) Minimise myfunc on a bounded interval
print("\n2) Bounded minimisation with optimize.minimize_scalar")
min_res = optimize.minimize_scalar(myfunc, bounds=(0.0, 2.0), method="bounded")
if min_res.success:
    x_min = min_res.x
    print(f"   Minimum on [0,2]: x* ≈ {x_min:.8f}, myfunc(x*) ≈ {myfunc(x_min):.6f}")
else:
    print("   Minimisation did not succeed.")


# 3) Multivariate minimisation with optimize.minimize
def multi_myfunc(x):
    return myfunc(x[0] * x[1])

print("\n3) Multivariate minimisation with optimize.minimize")
res = optimize.minimize(multi_myfunc, x0=[0.5, -0.5], bounds=[(-2, 2), (-2, 2)])
print(f"   x* = {res.x},  myfunc(x*) = {res.fun:.6f}")


# 4) ∫ myfunc(x) dx from 0 to 1
print("\n4) Integration with integrate.quad")
I, I_err = integrate.quad(myfunc, 0.0, 1.0)
print(f"   ∫_0^1 myfunc(x) dx ≈ {I:.8f}  (estimated abs error ≈ {I_err:.2e})")

#    Another integration example using our polynomial helper
print("   Integrate a polynomial: p(x) = 1 + 3x + 5x^2  over [0, 2]")
poly = lambda x: polynomial(x, a0=1, a1=3, a2=5)  # 1 + 3x + 5x^2
Ip, Ip_err = integrate.quad(poly, 0.0, 2.0)
print(f"   ∫_0^2 p(x) dx ≈ {Ip:.8f}  (estimated abs error ≈ {Ip_err:.2e})")


# 5) Fit the gaussian mean and sigma from noisy data. Our gaussian in myfuncs
#    returns exp(-((x-mean)^2)/(2 sigma^2)).
print("\n5) Curve fitting with optimize.curve_fit (fit mean and sigma of gaussian)")

# Make a NumPy-friendly wrapper (curve_fit expects vectorized numpy functions)
def gaussian_np(x, mean, sigma):
    return np.exp(-((x - mean) ** 2) / (2.0 * sigma ** 2))

# Generate synthetic noisy data
rng = np.random.default_rng(0)
x_data = np.linspace(-3, 3, 80)
true_mean, true_sigma = 0.5, 0.8
y_true = gaussian_np(x_data, true_mean, true_sigma)
noise = 0.05 * rng.standard_normal(size=x_data.shape)
y_obs = y_true + noise

# Initial guesses (mean0, sigma0)
p0 = (0.0, 1.0)
popt, pcov = curve_fit(gaussian_np, x_data, y_obs, p0=p0, bounds=([-np.inf, 1e-3],[np.inf, np.inf]))
est_mean, est_sigma = popt
perr = np.sqrt(np.diag(pcov))

print(f"   True params: mean={true_mean:.3f}, sigma={true_sigma:.3f}")
print(f"   Fit  params: mean={est_mean:.3f} ± {perr[0]:.3f}, sigma={est_sigma:.3f} ± {perr[1]:.3f}")

plt.figure()
plt.scatter(x_data, y_obs, s=20, label="Observed data", color="tab:blue")
plt.plot(x_data, y_true, "k--", label="True curve")
plt.plot(x_data, gaussian_np(x_data, *popt), "r", label="Fitted curve")
plt.title("Gaussian Fit with curve_fit")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.show()

print("\nAll SciPy demos completed.")