more-effort.py

import numpy as np
from scipy import linalg
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from typing import Callable, Tuple, List
from dataclasses import dataclass

@dataclass
class FDCoefficients:
    """Dataclass to store finite difference coefficients and related information."""
    coeffs: np.ndarray
    order: int
    accuracy: int
    stencil: np.ndarray

class FiniteDifference1D:
    def __init__(self, order: int, accuracy: int):
        """
        Initialize the FiniteDifference1D class.
        
        Args:
            order (int): The order of the derivative to approximate.
            accuracy (int): The order of accuracy for the approximation.
        """
        self.order = order
        self.accuracy = accuracy
        self.fd_coeffs = self._compute_fd_coefficients()

    def _compute_fd_coefficients(self) -> FDCoefficients:
        """
        Compute finite difference coefficients for given order and accuracy.
        
        Returns:
            FDCoefficients: Dataclass containing coefficients and related information.
        """
        p = self.order
        n = self.accuracy + p - 1
        stencil = np.arange(-(n//2), n//2 + 1)
        
        A = np.zeros((n+1, n+1))
        b = np.zeros(n+1)
        
        for i in range(n+1):
            A[i, :] = stencil**i / np.math.factorial(i)
        
        b[p] = 1
        
        coeffs = linalg.solve(A, b)
        return FDCoefficients(coeffs, self.order, self.accuracy, stencil)

    def apply(self, f: np.ndarray, dx: float) -> np.ndarray:
        """
        Apply the finite difference scheme to a given function.
        
        Args:
            f (np.ndarray): The function values at grid points.
            dx (float): The grid spacing.
        
        Returns:
            np.ndarray: The approximated derivative values.
        """
        n = len(f)
        m = len(self.fd_coeffs.stencil)
        result = np.zeros(n)
        
        for i in range(n):
            for j, coeff in enumerate(self.fd_coeffs.coeffs):
                idx = i + self.fd_coeffs.stencil[j]
                if 0 <= idx < n:
                    result[i] += coeff * f[idx]
        
        return result / (dx ** self.order)

class DataFitter:
    def __init__(self, x: np.ndarray, y: np.ndarray, fd_scheme: FiniteDifference1D):
        """
        Initialize the DataFitter class.
        
        Args:
            x (np.ndarray): The x-coordinates of the data points.
            y (np.ndarray): The y-coordinates of the data points.
            fd_scheme (FiniteDifference1D): The finite difference scheme to use.
        """
        self.x = x
        self.y = y
        self.fd_scheme = fd_scheme
        self.fit_func = None

    def fit(self, num_points: int = 1000) -> Tuple[np.ndarray, np.ndarray]:
        """
        Fit the data using the finite difference scheme.
        
        Args:
            num_points (int): The number of points to use in the fitted function.
        
        Returns:
            Tuple[np.ndarray, np.ndarray]: The x and y values of the fitted function.
        """
        x_fine = np.linspace(self.x.min(), self.x.max(), num_points)
        dx = x_fine[1] - x_fine[0]
        
        # Interpolate the data to the fine grid
        interp = interp1d(self.x, self.y, kind='cubic')
        y_fine = interp(x_fine)
        
        # Apply the finite difference scheme
        derivative = self.fd_scheme.apply(y_fine, dx)
        
        # Integrate the derivative to get the fitted function
        y_fit = np.cumsum(derivative) * dx + y_fine[0]
        
        self.fit_func = interp1d(x_fine, y_fit, kind='cubic')
        return x_fine, y_fit

    def evaluate(self, x: np.ndarray) -> np.ndarray:
        """
        Evaluate the fitted function at given x values.
        
        Args:
            x (np.ndarray): The x values to evaluate the function at.
        
        Returns:
            np.ndarray: The y values of the fitted function at the given x values.
        """
        if self.fit_func is None:
            raise ValueError("Fit the data first using the 'fit' method.")
        return self.fit_func(x)

def visualize_results(x: np.ndarray, y: np.ndarray, x_fit: np.ndarray, y_fit: np.ndarray, title: str):
    """
    Visualize the original data and the fitted function.
    
    Args:
        x (np.ndarray): The x-coordinates of the original data.
        y (np.ndarray): The y-coordinates of the original data.
        x_fit (np.ndarray): The x-coordinates of the fitted function.
        y_fit (np.ndarray): The y-coordinates of the fitted function.
        title (str): The title of the plot.
    """
    plt.figure(figsize=(12, 6))
    plt.scatter(x, y, label='Original Data', color='red', alpha=0.7)
    plt.plot(x_fit, y_fit, label='Fitted Function', color='blue')
    plt.title(title)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

def calculate_error_metrics(y_true: np.ndarray, y_pred: np.ndarray) -> dict:
    """
    Calculate error metrics between true and predicted values.
    
    Args:
        y_true (np.ndarray): The true y values.
        y_pred (np.ndarray): The predicted y values.
    
    Returns:
        dict: A dictionary containing various error metrics.
    """
    mse = np.mean((y_true - y_pred)**2)
    rmse = np.sqrt(mse)
    mae = np.mean(np.abs(y_true - y_pred))
    mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
    
    return {
        'MSE': mse,
        'RMSE': rmse,
        'MAE': mae,
        'MAPE': mape
    }

def main():
    # Generate sample data
    x = np.linspace(0, 10, 50)
    y = np.sin(x) + 0.1 * np.random.randn(len(x))

    # Create finite difference scheme
    fd_scheme = FiniteDifference1D(order=1, accuracy=4)

    # Fit the data
    fitter = DataFitter(x, y, fd_scheme)
    x_fit, y_fit = fitter.fit(num_points=1000)

    # Visualize results
    visualize_results(x, y, x_fit, y_fit, "Finite Difference Data Fitting")

    # Calculate error metrics
    y_pred = fitter.evaluate(x)
    error_metrics = calculate_error_metrics(y, y_pred)

    print("Error Metrics:")
    for metric, value in error_metrics.items():
        print(f"{metric}: {value:.6f}")

if __name__ == "__main__":
    main()