Examples
========

This page contains detailed examples showing how to use pydelt for various tasks.

Example 1: Comparing Derivative Methods
---------------------------------------

Let's compare different derivative calculation methods on a known function:

.. code-block:: python

   import numpy as np
   import matplotlib.pyplot as plt
   from pydelt.derivatives import lla, fda, gold, glla
   
   # Generate test signal: f(x) = sin(2x) + 0.5*cos(5x)
   x = np.linspace(0, 2*np.pi, 200)
   signal = np.sin(2*x) + 0.5*np.cos(5*x)
   
   # Analytical derivative: f'(x) = 2*cos(2x) - 2.5*sin(5x)
   analytical = 2*np.cos(2*x) - 2.5*np.sin(5*x)
   
   # Compare different methods
   lla_result = lla(x.tolist(), signal.tolist(), window_size=5)
   fda_result = fda(signal, x)
   gold_result = gold(signal, x, embedding=3)
   glla_result = glla(signal, x, embedding=3)
   
   methods = {
       'LLA': lla_result[0],
       'FDA': fda_result['dsignal'][:, 1],
       'GOLD': gold_result['dsignal'][:, 0],
       'GLLA': glla_result['dsignal'][:, 0]
   }
   
   # Calculate errors for each method
   errors = {}
   
   for name, derivative in methods.items():
       try:
           errors[name] = np.mean(np.abs(derivative - analytical))
           print(f"{name}: Mean absolute error = {errors[name]:.6f}")
       except Exception as e:
           print(f"{name}: Failed with error {e}")
   
   # Plot results
   plt.figure(figsize=(12, 8))
   
   plt.subplot(2, 1, 1)
   plt.plot(x, signal, 'k-', label='Original signal', linewidth=2)
   plt.plot(x, analytical, 'r--', label='Analytical derivative', linewidth=2)
   plt.legend()
   plt.title('Original Signal and Analytical Derivative')
   
   plt.subplot(2, 1, 2)
   for name, derivative in methods.items():
       plt.plot(x, derivative, label=f'{name} (error: {errors[name]:.4f})')
   plt.plot(x, analytical, 'r--', label='Analytical', linewidth=2)
   plt.legend()
   plt.title('Derivative Comparison')
   plt.xlabel('x')
   
   plt.tight_layout()
   plt.show()

Example 2: Neural Network Derivatives
-------------------------------------

Using neural networks for derivative calculation:

.. code-block:: python

   import numpy as np
   from pydelt.autodiff import neural_network_derivative
   
   # Create training data
   x_train = np.linspace(0, 4*np.pi, 100)
   y_train = np.exp(-0.1*x_train) * np.sin(x_train)
   
   # Train neural network
   try:
       # Using PyTorch backend
       derivative_func = neural_network_derivative(
           x_train, y_train,
           framework='pytorch',
           hidden_layers=[128, 64, 32],
           epochs=1000,
           learning_rate=0.001,
           dropout=0.1
       )
       
       # Evaluate on test points
       x_test = np.linspace(0.5, 3.5*np.pi, 50)
       derivatives = derivative_func(x_test)
       
       # Analytical derivative for comparison
       analytical = -0.1*np.exp(-0.1*x_test)*np.sin(x_test) + np.exp(-0.1*x_test)*np.cos(x_test)
       
       error = np.mean(np.abs(derivatives.flatten() - analytical))
       print(f"Neural network derivative error: {error:.4f}")
       
       # Plot results
       import matplotlib.pyplot as plt
       plt.figure(figsize=(10, 6))
       plt.plot(x_test, analytical, 'r-', label='Analytical derivative', linewidth=2)
       plt.plot(x_test, derivatives.flatten(), 'b--', label='Neural network', linewidth=2)
       plt.legend()
       plt.title('Neural Network vs Analytical Derivative')
       plt.xlabel('x')
       plt.ylabel("f'(x)")
       plt.show()
       
   except ImportError:
       print("PyTorch not available - install with: pip install torch")

Example 3: Noisy Data Processing
--------------------------------

Handling noisy time series data:

.. code-block:: python

   import numpy as np
   from pydelt.derivatives import lla
   from pydelt.interpolation import lowess_interpolation, spline_interpolation
   
   # Create noisy data
   np.random.seed(42)
   t = np.linspace(0, 10, 100)
   clean_signal = np.sin(t) * np.exp(-0.1*t)
   noise = 0.2 * np.random.randn(len(t))
   noisy_signal = clean_signal + noise
   
   # Method 1: Direct derivative of noisy data
   direct_result = lla(t.tolist(), noisy_signal.tolist(), window_size=5)
   direct_derivative = direct_result[0]
   
   # Method 2: Smooth first, then differentiate
   smoother = lowess_interpolation(t, noisy_signal, frac=0.1)
   smoothed_signal = smoother(t)
   smooth_result = lla(t.tolist(), smoothed_signal.tolist(), window_size=5)
   smooth_derivative = smooth_result[0]
   
   # Method 3: Spline smoothing
   spline_smoother = spline_interpolation(t, noisy_signal, smoothing_factor=1.0)
   spline_signal = spline_smoother(t)
   spline_result = lla(t.tolist(), spline_signal.tolist(), window_size=5)
   spline_derivative = spline_result[0]
   
   # Analytical derivative for comparison
   analytical_derivative = np.cos(t)*np.exp(-0.1*t) - 0.1*np.sin(t)*np.exp(-0.1*t)
   
   # Calculate errors
   errors = {
       'Direct': np.mean(np.abs(direct_derivative - analytical_derivative)),
       'LOWESS': np.mean(np.abs(smooth_derivative - analytical_derivative)),
       'Spline': np.mean(np.abs(spline_derivative - analytical_derivative))
   }
   
   for method, error in errors.items():
       print(f"{method} derivative error: {error:.4f}")
   
   # Plotting
   import matplotlib.pyplot as plt
   fig, axes = plt.subplots(2, 2, figsize=(15, 10))
   
   # Original signals
   axes[0,0].plot(t, clean_signal, 'g-', label='Clean signal', linewidth=2)
   axes[0,0].plot(t, noisy_signal, 'k.', alpha=0.5, label='Noisy signal')
   axes[0,0].plot(t, smoothed_signal, 'b-', label='LOWESS smoothed')
   axes[0,0].plot(t, spline_signal, 'r-', label='Spline smoothed')
   axes[0,0].legend()
   axes[0,0].set_title('Signal Comparison')
   
   # Derivatives
   axes[0,1].plot(t, analytical_derivative, 'g-', label='Analytical', linewidth=3)
   axes[0,1].plot(t, direct_derivative, 'k-', alpha=0.7, label=f'Direct (err: {errors["Direct"]:.3f})')
   axes[0,1].plot(t, smooth_derivative, 'b-', label=f'LOWESS (err: {errors["LOWESS"]:.3f})')
   axes[0,1].plot(t, spline_derivative, 'r-', label=f'Spline (err: {errors["Spline"]:.3f})')
   axes[0,1].legend()
   axes[0,1].set_title('Derivative Comparison')
   
   # Error plots
   axes[1,0].plot(t, np.abs(direct_derivative - analytical_derivative), 'k-', label='Direct')
   axes[1,0].plot(t, np.abs(smooth_derivative - analytical_derivative), 'b-', label='LOWESS')
   axes[1,0].plot(t, np.abs(spline_derivative - analytical_derivative), 'r-', label='Spline')
   axes[1,0].set_yscale('log')
   axes[1,0].legend()
   axes[1,0].set_title('Absolute Errors (log scale)')
   axes[1,0].set_xlabel('Time')
   
   # Histogram of errors
   axes[1,1].hist(np.abs(direct_derivative - analytical_derivative), alpha=0.5, label='Direct', bins=20)
   axes[1,1].hist(np.abs(smooth_derivative - analytical_derivative), alpha=0.5, label='LOWESS', bins=20)
   axes[1,1].hist(np.abs(spline_derivative - analytical_derivative), alpha=0.5, label='Spline', bins=20)
   axes[1,1].legend()
   axes[1,1].set_title('Error Distribution')
   axes[1,1].set_xlabel('Absolute Error')
   
   plt.tight_layout()
   plt.show()

Example 4: Integration and Roundtrip Accuracy
---------------------------------------------

Testing integration accuracy by doing derivative → integral roundtrips:

.. code-block:: python

   import numpy as np
   from pydelt.derivatives import lla
   from pydelt.integrals import integrate_derivative
   
   # Original function: f(x) = x^3 - 2x^2 + x + 5
   x = np.linspace(0, 5, 100)
   original_function = x**3 - 2*x**2 + x + 5
   
   # Analytical derivative: f'(x) = 3x^2 - 4x + 1
   analytical_derivative = 3*x**2 - 4*x + 1
   
   # Step 1: Calculate derivative numerically
   result = lla(x.tolist(), original_function.tolist(), window_size=5)
   numerical_derivative = result[0]
   
   # Step 2: Integrate the derivative back
   integrated_function, integration_error = integrate_derivative(
       x, numerical_derivative, 
       initial_value=original_function[0]
   )
   
   # Step 3: Calculate errors
   derivative_error = np.mean(np.abs(numerical_derivative - analytical_derivative))
   roundtrip_error = np.mean(np.abs(integrated_function - original_function))
   
   print(f"Derivative calculation error: {derivative_error:.6f}")
   print(f"Integration roundtrip error: {roundtrip_error:.6f}")
   print(f"Estimated integration error: {integration_error:.6f}")
   
   # Plotting
   import matplotlib.pyplot as plt
   fig, axes = plt.subplots(2, 2, figsize=(12, 10))
   
   # Original function
   axes[0,0].plot(x, original_function, 'b-', linewidth=2, label='Original f(x)')
   axes[0,0].plot(x, integrated_function, 'r--', linewidth=2, label='Integrated f(x)')
   axes[0,0].legend()
   axes[0,0].set_title('Function Comparison')
   axes[0,0].set_ylabel('f(x)')
   
   # Derivatives
   axes[0,1].plot(x, analytical_derivative, 'b-', linewidth=2, label='Analytical f\'(x)')
   axes[0,1].plot(x, numerical_derivative, 'r--', linewidth=2, label='Numerical f\'(x)')
   axes[0,1].legend()
   axes[0,1].set_title('Derivative Comparison')
   axes[0,1].set_ylabel('f\'(x)')
   
   # Function error
   axes[1,0].plot(x, np.abs(integrated_function - original_function), 'g-', linewidth=2)
   axes[1,0].set_title(f'Function Roundtrip Error (mean: {roundtrip_error:.4f})')
   axes[1,0].set_ylabel('|Error|')
   axes[1,0].set_xlabel('x')
   
   # Derivative error
   axes[1,1].plot(x, np.abs(numerical_derivative - analytical_derivative), 'orange', linewidth=2)
   axes[1,1].set_title(f'Derivative Error (mean: {derivative_error:.4f})')
   axes[1,1].set_ylabel('|Error|')
   axes[1,1].set_xlabel('x')
   
   plt.tight_layout()
   plt.show()

Example 5: Multivariate Time Series
-----------------------------------

Working with multi-dimensional data:

.. code-block:: python

   import numpy as np
   from pydelt.derivatives import lla
   
   # Create 3D trajectory data (e.g., particle motion)
   t = np.linspace(0, 4*np.pi, 200)
   
   # Parametric equations for a 3D spiral
   x = np.cos(t) * np.exp(-0.1*t)
   y = np.sin(t) * np.exp(-0.1*t)  
   z = 0.1 * t
   
   # Combine into multivariate signal
   trajectory = np.column_stack([x, y, z])
   
   # Calculate velocity (derivative of position) using multivariate support
   velocity_result = lla(t.tolist(), trajectory.tolist(), window_size=5)
   velocity = velocity_result[0]  # Shape: (N, 3) for 3D trajectory
   
   # Calculate speed (magnitude of velocity)
   speed = np.sqrt(np.sum(velocity**2, axis=1))
   
   # Analytical derivatives for comparison
   dx_dt = -np.sin(t)*np.exp(-0.1*t) - 0.1*np.cos(t)*np.exp(-0.1*t)
   dy_dt = np.cos(t)*np.exp(-0.1*t) - 0.1*np.sin(t)*np.exp(-0.1*t)
   dz_dt = 0.1 * np.ones_like(t)
   
   analytical_velocity = np.column_stack([dx_dt, dy_dt, dz_dt])
   analytical_speed = np.sqrt(np.sum(analytical_velocity**2, axis=1))
   
   # Calculate errors
   velocity_error = np.mean(np.abs(velocity - analytical_velocity))
   speed_error = np.mean(np.abs(speed - analytical_speed))
   
   print(f"Velocity calculation error: {velocity_error:.6f}")
   print(f"Speed calculation error: {speed_error:.6f}")
   
   # 3D plotting
   import matplotlib.pyplot as plt
   from mpl_toolkits.mplot3d import Axes3D
   
   fig = plt.figure(figsize=(15, 5))
   
   # 3D trajectory
   ax1 = fig.add_subplot(131, projection='3d')
   ax1.plot(x, y, z, 'b-', linewidth=2, label='Trajectory')
   ax1.set_xlabel('X')
   ax1.set_ylabel('Y')
   ax1.set_zlabel('Z')
   ax1.set_title('3D Trajectory')
   
   # Velocity components
   ax2 = fig.add_subplot(132)
   ax2.plot(t, velocity[:, 0], 'r-', label='vx (numerical)')
   ax2.plot(t, velocity[:, 1], 'g-', label='vy (numerical)')
   ax2.plot(t, velocity[:, 2], 'b-', label='vz (numerical)')
   ax2.plot(t, dx_dt, 'r--', alpha=0.7, label='vx (analytical)')
   ax2.plot(t, dy_dt, 'g--', alpha=0.7, label='vy (analytical)')
   ax2.plot(t, dz_dt, 'b--', alpha=0.7, label='vz (analytical)')
   ax2.legend()
   ax2.set_title('Velocity Components')
   ax2.set_xlabel('Time')
   ax2.set_ylabel('Velocity')
   
   # Speed comparison
   ax3 = fig.add_subplot(133)
   ax3.plot(t, speed, 'k-', linewidth=2, label='Numerical speed')
   ax3.plot(t, analytical_speed, 'r--', linewidth=2, label='Analytical speed')
   ax3.legend()
   ax3.set_title(f'Speed Comparison (error: {speed_error:.4f})')
   ax3.set_xlabel('Time')
   ax3.set_ylabel('Speed')
   
   plt.tight_layout()
   plt.show()

Example 6: Understanding Numerical Limitations
----------------------------------------------

Demonstrating critical point smoothing in multivariate derivatives:

.. code-block:: python

   import numpy as np
   from pydelt.multivariate import MultivariateDerivatives
   from pydelt.interpolation import SplineInterpolator
   
   # Function with known critical points: f(x,y) = (x-y)²
   x = np.linspace(-3, 3, 25)
   y = np.linspace(-3, 3, 25)
   X, Y = np.meshgrid(x, y)
   Z = (X - Y)**2  # Zero gradient along x=y line
   
   # Fit multivariate derivatives
   input_data = np.column_stack([X.flatten(), Y.flatten()])
   output_data = Z.flatten()
   
   mv = MultivariateDerivatives(SplineInterpolator, smoothing=0.1)
   mv.fit(input_data, output_data)
   
   # Test at critical points where gradient should be zero
   critical_points = np.array([[-3, -3], [3, 3]])  # Should have zero gradient
   
   gradient_func = mv.gradient()
   grad_critical = gradient_func(critical_points)
   
   print("Critical Points (should be ~[0, 0]):")
   for i, point in enumerate(critical_points):
       print(f"  Point {point}: Numerical gradient {grad_critical[i]:.3f}")
   
   print("\nKey Insight:")
   print("  Numerical methods smooth out zero gradients at critical points.")
   print("  This is a fundamental limitation of interpolation-based derivatives.")

**Key Takeaways:**

* Numerical methods cannot perfectly capture sharp mathematical features
* Critical points often show non-zero numerical gradients due to smoothing
* Always validate numerical results against analytical solutions when possible
* Consider neural network methods for exact derivatives at critical points

Example 7: Tensor Calculus
--------------------------

Directional derivative, divergence, curl, and strain/stress

.. code-block:: python

   import numpy as np
   from pydelt.multivariate import MultivariateDerivatives
   from pydelt.tensor_derivatives import TensorDerivatives
   from pydelt.interpolation import SplineInterpolator

   # 2D vector field: rotation F(x, y) = [-y, x]
   xs = ys = np.linspace(-2, 2, 25)
   X, Y = np.meshgrid(xs, ys)
   U = -Y
   V = X
   inputs = np.column_stack([X.flatten(), Y.flatten()])
   outputs = np.column_stack([U.flatten(), V.flatten()])

   mv = MultivariateDerivatives(SplineInterpolator, smoothing=0.0)
   mv.fit(inputs, outputs)

   q = np.array([[1.0, 0.5]])
   e = np.array([0.0, 1.0])
   dF_e = mv.directional_derivative(e)(q)  # shape (1, 2)
   divF = mv.divergence()(q)               # ≈ 0
   curlF = mv.curl()(q)                    # ≈ 2
   print("Directional derivative along e:", dF_e)
   print("div(F):", float(divF), " curl(F):", float(curlF))

   # 2D displacement field for strain/stress tensors
   U = 0.05 * X**2 - 0.02 * Y
   V = 0.10 * X * Y
   displacement = np.column_stack([U.flatten(), V.flatten()])

   td = TensorDerivatives(SplineInterpolator, smoothing=0.0)
   td.fit(inputs, displacement)

   strain = td.strain_tensor()(inputs)  # (N, 2, 2)
   stress = td.stress_tensor(lambda_param=1.0, mu_param=0.5)(inputs)  # (N, 2, 2)

   # Correct component extraction using three indices
   exx = strain[:, 0, 0].reshape(X.shape)
   eyy = strain[:, 1, 1].reshape(X.shape)
   exy = strain[:, 0, 1].reshape(X.shape)
   sxy = stress[:, 0, 1].reshape(X.shape)
   print("mean |ε_xy|:", float(np.mean(np.abs(exy))))
   print("mean |σ_xy|:", float(np.mean(np.abs(sxy))))
