scipy

What this skill does

# SciPy - Scientific Computing

Advanced scientific computing library built on NumPy, providing algorithms for optimization, integration, interpolation, and more.

## When to Use

- Integrating functions (numerical integration, ODEs)
- Optimizing functions (minimization, root finding, curve fitting)
- Interpolating data (1D, 2D, splines)
- Advanced linear algebra (sparse matrices, decompositions)
- Signal processing (filtering, Fourier transforms, wavelets)
- Statistical analysis (distributions, hypothesis tests)
- Image processing (filters, morphology, measurements)
- Spatial algorithms (distance matrices, clustering, Voronoi)
- Special mathematical functions (Bessel, gamma, error functions)
- Solving differential equations (ODEs, PDEs)

## Reference Documentation

**Official docs**: https://docs.scipy.org/  
**Search patterns**: `scipy.integrate.quad`, `scipy.optimize.minimize`, `scipy.interpolate`, `scipy.stats`, `scipy.signal`

## Core Principles

### Use SciPy For

| Task | Module | Example |
|------|--------|---------|
| Integration | `integrate` | `quad(f, 0, 1)` |
| Optimization | `optimize` | `minimize(f, x0)` |
| Interpolation | `interpolate` | `interp1d(x, y)` |
| Linear algebra | `linalg` | `linalg.solve(A, b)` |
| Signal processing | `signal` | `signal.butter(4, 0.5)` |
| Statistics | `stats` | `stats.norm.pdf(x)` |
| ODEs | `integrate` | `solve_ivp(f, t_span, y0)` |
| FFT | `fft` | `fft.fft(signal)` |

### Do NOT Use For

- Basic array operations (use NumPy)
- Machine learning (use scikit-learn)
- Deep learning (use PyTorch, TensorFlow)
- Symbolic mathematics (use SymPy)
- Data manipulation (use pandas)

## Quick Reference

### Installation

```bash
# pip
pip install scipy

# conda
conda install scipy

# With NumPy
pip install numpy scipy
```

### Standard Imports

```python
import numpy as np
from scipy import integrate, optimize, interpolate
from scipy import linalg, signal, stats
from scipy.integrate import odeint, solve_ivp
from scipy.optimize import minimize, root
from scipy.interpolate import interp1d, UnivariateSpline
```

### Basic Pattern - Integration

```python
from scipy import integrate
import numpy as np

# Define function
def f(x):
    return x**2

# Integrate from 0 to 1
result, error = integrate.quad(f, 0, 1)
print(f"Integral: {result:.6f} ± {error:.2e}")
```

### Basic Pattern - Optimization

```python
from scipy import optimize
import numpy as np

# Function to minimize
def f(x):
    return (x - 2)**2 + 1

# Minimize
result = optimize.minimize(f, x0=0)
print(f"Minimum at x = {result.x[0]:.6f}")
print(f"Minimum value = {result.fun:.6f}")
```

### Basic Pattern - Interpolation

```python
from scipy import interpolate
import numpy as np

# Data points
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

# Create interpolator
f = interpolate.interp1d(x, y, kind='cubic')

# Interpolate at new points
x_new = np.linspace(0, 4, 100)
y_new = f(x_new)
```

## Critical Rules

### ✅ DO

- **Check convergence** - Always verify optimization converged
- **Specify tolerances** - Set appropriate `rtol` and `atol`
- **Use appropriate methods** - Choose algorithm for problem type
- **Validate inputs** - Check array shapes and values
- **Handle edge cases** - Deal with singularities and discontinuities
- **Set integration limits carefully** - Watch for infinite limits
- **Use vectorization** - Functions should accept arrays
- **Check statistical assumptions** - Verify distribution assumptions
- **Specify degrees of freedom** - For interpolation and fitting
- **Use sparse matrices** - For large, sparse systems

### ❌ DON'T

- **Ignore convergence warnings** - They indicate problems
- **Use inappropriate tolerances** - Too loose or too tight
- **Apply wrong distribution** - Check data characteristics
- **Forget initial guesses** - Optimization needs good starting points
- **Integrate discontinuous functions** - Without special handling
- **Extrapolate beyond data** - Interpolation is not extrapolation
- **Mix incompatible units** - Keep consistent units
- **Ignore error estimates** - They provide confidence levels
- **Use wrong coordinate system** - Check Cartesian vs polar
- **Overfit with high-degree polynomials** - Causes oscillations

## Anti-Patterns (NEVER)

```python
from scipy import integrate, optimize
import numpy as np

# ❌ BAD: Ignoring convergence
result = optimize.minimize(f, x0=0)
optimal_x = result.x  # Didn't check if converged!

# ✅ GOOD: Check convergence
result = optimize.minimize(f, x0=0)
if result.success:
    optimal_x = result.x
else:
    print(f"Optimization failed: {result.message}")

# ❌ BAD: Non-vectorized function for integration
def bad_func(x):
    if x < 0.5:
        return x
    else:
        return 1 - x

# ✅ GOOD: Vectorized function
def good_func(x):
    return np.where(x < 0.5, x, 1 - x)

# ❌ BAD: Poor initial guess
result = optimize.minimize(complex_func, x0=[1000, 1000])
# May converge to local minimum or fail!

# ✅ GOOD: Reasonable initial guess
x0 = np.array([0.0, 0.0])  # Near expected minimum
result = optimize.minimize(complex_func, x0=x0)

# ❌ BAD: Extrapolation with interpolation
f = interpolate.interp1d(x_data, y_data)
y_new = f(100)  # x_data max is 10, this will crash!

# ✅ GOOD: Check bounds or use extrapolation
f = interpolate.interp1d(x_data, y_data, fill_value='extrapolate')
y_new = f(100)  # Now works (but be cautious!)

# ❌ BAD: Wrong statistical test
# Using t-test for non-normal data
stats.ttest_ind(non_normal_data1, non_normal_data2)

# ✅ GOOD: Use appropriate test
# Use Mann-Whitney U for non-normal data
stats.mannwhitneyu(non_normal_data1, non_normal_data2)
```

## Integration (scipy.integrate)

### Numerical Integration (Quadrature)

```python
from scipy import integrate
import numpy as np

# Single integral
def f(x):
    return np.exp(-x**2)

result, error = integrate.quad(f, 0, np.inf)
print(f"∫exp(-x²)dx from 0 to ∞ = {result:.6f}")
print(f"Error estimate: {error:.2e}")

# Integral with parameters
def g(x, a, b):
    return a * x**2 + b

result, error = integrate.quad(g, 0, 1, args=(2, 3))
print(f"Result: {result:.6f}")

# Integral with singularity
def h(x):
    return 1 / np.sqrt(x)

# Specify singularity points
result, error = integrate.quad(h, 0, 1, points=[0])
```

### Double and Triple Integrals

```python
from scipy import integrate
import numpy as np

# Double integral: ∫∫ x*y dx dy over [0,1]×[0,2]
def f(y, x):  # Note: y first, x second
    return x * y

result, error = integrate.dblquad(f, 0, 1, 0, 2)
print(f"Double integral: {result:.6f}")

# Triple integral
def g(z, y, x):
    return x * y * z

result, error = integrate.tplquad(g, 0, 1, 0, 1, 0, 1)
print(f"Triple integral: {result:.6f}")

# Variable limits
def lower(x):
    return 0

def upper(x):
    return x

result, error = integrate.dblquad(f, 0, 1, lower, upper)
```

### Solving ODEs

```python
from scipy.integrate import odeint, solve_ivp
import numpy as np

# Solve dy/dt = -k*y (exponential decay)
def exponential_decay(y, t, k):
    return -k * y

# Initial condition and time points
y0 = 100
t = np.linspace(0, 10, 100)
k = 0.5

# Solve with odeint (older interface)
solution = odeint(exponential_decay, y0, t, args=(k,))

# Solve with solve_ivp (modern interface)
def decay_ivp(t, y, k):
    return -k * y

sol = solve_ivp(decay_ivp, [0, 10], [y0], args=(k,), t_eval=t)

print(f"Final value (odeint): {solution[-1, 0]:.6f}")
print(f"Final value (solve_ivp): {sol.y[0, -1]:.6f}")
```

### System of ODEs

```python
from scipy.integrate import solve_ivp
import numpy as np

# Lotka-Volterra equations (predator-prey)
def lotka_volterra(t, z, a, b, c, d):
    x, y = z
    dxdt = a*x - b*x*y
    dydt = -c*y + d*x*y
    return [dxdt, dydt]

# Parameters
a, b, c, d = 1.5, 1.0, 3.0, 1.0

# Initial conditions
z0 = [10, 5]  # [prey, predator]

# Time span
t_span = (0, 15)
t_eval = np.linspace(0, 15, 1000)

# Solve
sol = solve_ivp(lotka_volterra, t_span, z0, args=(a, b, c, d), 
                t_eval=t_eval, method='