Python for Scientific Computing — NumPy, SciPy, and Matplotlib

import numpy as np

# Creating arrays
arr = np.array([1, 2, 3, 4, 5])
zeros = np.zeros((3, 4))
ones = np.ones((2, 3))
identity = np.eye(5)
random = np.random.randn(1000)  # Normal distribution
range_arr = np.arange(0, 10, 2)  # [0, 2, 4, 6, 8]
linspace = np.linspace(0, 1, 5)  # [0.0, 0.25, 0.5, 0.75, 1.0]

# Array attributes
print(arr.shape)      # (5,)
print(arr.ndim)       # 1
print(arr.dtype)      # int64
print(arr.size)       # 5
print(arr.nbytes)     # 40 (5 * 8 bytes)

arr = np.array([[1, 2, 3, 4],
                [5, 6, 7, 8],
                [9, 10, 11, 12]])

print(arr[0, 0])       # 1
print(arr[:, 1])       # [2, 6, 10] (column 1)
print(arr[1:, :2])     # [[5, 6], [9, 10]]
print(arr[arr > 5])    # [6, 7, 8, 9, 10, 11, 12] (boolean indexing)

# Fancy indexing
rows = np.array([0, 2])
cols = np.array([1, 3])
print(arr[rows, cols])  # [2, 12]

# Avoid Python loops — use vectorized operations
data = np.random.randn(1000000)

# SLOW: Python loop
total = 0
for x in data:
    total += x ** 2 + np.sin(x)

# FAST: vectorized
result = np.sum(data ** 2 + np.sin(data))

# Universal functions
np.sqrt(data)
np.exp(data)
np.log(np.abs(data) + 1e-10)
np.sin(data)
np.abs(data)

# Aggregations
print(data.mean())     # ~0 (normal distribution)
print(data.std())      # ~1
print(data.min())
print(data.max())
print(np.percentile(data, [25, 50, 75]))

# Broadcasting: operate on arrays of different shapes
matrix = np.random.randn(4, 3)
row_mean = matrix.mean(axis=1, keepdims=True)  # (4, 1)
centered = matrix - row_mean  # Broadcasting! (4,3) - (4,1) → (4,3)

# Add a constant to each row
matrix += np.array([1, 2, 3])

# Outer product
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
outer = x[:, np.newaxis] * y  # (3,1) * (3,) → (3,3)

from scipy import optimize, stats, signal, ndimage
import numpy as np

### Optimization

# Find minimum of a function
def f(x):
    return x ** 2 + 10 * np.sin(x)

result = optimize.minimize_scalar(f)
print(f"Minimum at x={result.x:.3f}, f(x)={result.fun:.3f}")

# Curve fitting
x_data = np.linspace(0, 10, 100)
y_data = 3 * np.sin(2 * x_data) + 0.5 * np.random.randn(100)

def model(x, a, b, c):
    return a * np.sin(b * x) + c

params, _ = optimize.curve_fit(model, x_data, y_data)
print(f"Fitted: a={params[0]:.3f}, b={params[1]:.3f}, c={params[2]:.3f}")

# Descriptive statistics
data = np.random.normal(170, 10, 1000)  # Heights (cm)
print(f"Mean: {np.mean(data):.1f}")
print(f"Std: {np.std(data):.1f}")
print(f"Skewness: {stats.skew(data):.2f}")
print(f"Kurtosis: {stats.kurtosis(data):.2f}")

# Hypothesis testing
sample1 = np.random.normal(100, 15, 50)
sample2 = np.random.normal(105, 15, 50)

t_stat, p_value = stats.ttest_ind(sample1, sample2)
print(f"t-test: t={t_stat:.3f}, p={p_value:.4f}")
if p_value < 0.05:
    print("Significant difference!")

# Probability distributions
normal_rv = stats.norm(loc=0, scale=1)
print(normal_rv.pdf(0))     # 0.399 (PDF at 0)
print(normal_rv.cdf(1.96))  # 0.975 (CDF at 1.96)
print(normal_rv.ppf(0.975)) # 1.96 (inverse CDF)

# Moving average filter
signal_data = np.sin(np.linspace(0, 10, 1000)) + np.random.randn(1000) * 0.1
smoothed = signal.savgol_filter(signal_data, window_length=51, polyorder=3)

# FFT
freqs = np.fft.fftfreq(len(signal_data))
fft = np.fft.fft(signal_data)
magnitude = np.abs(fft)[:len(fft)//2]

import matplotlib.pyplot as plt

### Basic Plot

x = np.linspace(0, 10, 100)
y1 = np.sin(x)
y2 = np.cos(x)

plt.figure(figsize=(10, 6))
plt.plot(x, y1, label="sin(x)", linewidth=2)
plt.plot(x, y2, label="cos(x)", linewidth=2, linestyle="--")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Sine and Cosine")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

data = np.random.randn(10000)
plt.figure(figsize=(10, 6))
plt.hist(data, bins=50, density=True, alpha=0.7, color="steelblue")
plt.title("Normal Distribution (n=10000)")
plt.xlabel("Value")
plt.ylabel("Density")
plt.show()

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

x = np.linspace(0, 10, 100)

axes[0, 0].plot(x, np.sin(x))
axes[0, 0].set_title("sin(x)")

axes[0, 1].plot(x, np.cos(x), color="red")
axes[0, 1].set_title("cos(x)")

axes[1, 0].plot(x, np.sin(x) * np.cos(x), color="green")
axes[1, 0].set_title("sin(x) * cos(x)")

axes[1, 1].scatter(x, np.sin(x) + np.random.randn(100) * 0.1, s=5)
axes[1, 1].set_title("Noisy sin(x)")

plt.tight_layout()
plt.show()

# Heatmap
matrix = np.random.randn(10, 10)
plt.imshow(matrix, cmap="RdBu", aspect="auto")
plt.colorbar(label="Value")
plt.title("Heatmap")
plt.show()

# 3D plot
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
X, Y = np.meshgrid(np.linspace(-5, 5, 50), np.linspace(-5, 5, 50))
Z = np.sin(np.sqrt(X**2 + Y**2))
ax.plot_surface(X, Y, Z, cmap="viridis")
plt.show()