1.3) Scientific Computing with numpy

import numpy as np

py_list = [1, 2, 3, 4]
arr = np.array([1, 2, 3, 4])

print(py_list * 2)     # [1, 2, 3, 4, 1, 2, 3, 4] — the list is repeated
print(arr * 2)         # [2 4 6 8] — every element is doubled

# a list can hold mixed types; an array cannot
print(type(1).__name__, type("two").__name__, type(3.0).__name__)   # int str float
print(np.array([1, 2, 3]).dtype)                                    # int64 — one type for all

import numpy as np

# a 2D field: 4 latitudes (rows) x 6 longitudes (cols), daily mean temp (°C)
temp_celsius = np.array([
    [ 5.2,  4.8,  6.1,  3.9,  2.7,  4.4],
    [ 1.3,  0.5, -0.8, -1.2,  0.9,  2.1],
    [-2.6, -3.1, -1.9,  0.2, -0.5,  1.1],
    [ 6.4,  7.0,  5.5,  8.1,  4.2,  5.8],
])

print(temp_celsius)
print("shape:", temp_celsius.shape, "| ndim:", temp_celsius.ndim)

print(np.zeros((2, 3)))          # all zeros, given shape
print(np.ones(4))                # all ones
print(np.full((2, 2), 7.0))      # filled with a constant
print(np.arange(0, 10, 2))       # evenly spaced by step: [0 2 4 6 8]
print(np.linspace(0.0, 1.0, 5))  # n evenly spaced points: [0. 0.25 0.5 0.75 1. ]
print(np.random.random(3))       # 3 random floats in [0, 1)

lon = np.linspace(6.0, 9.0, 6)     # 6 longitudes
lat = np.linspace(46.0, 47.5, 4)   # 4 latitudes

lon2d, lat2d = np.meshgrid(lon, lat)
print(lon2d.shape, lat2d.shape)    # (4, 6) (4, 6) — one value per grid cell
print(lon2d)

print(temp_celsius[0, 0])      # one element
print(temp_celsius[0, :])      # first row, all longitudes
print(temp_celsius[:, -1])     # last column, all latitudes
print(temp_celsius[1:3, 2:4])  # a 2x2 sub-block

# a boolean mask and what it selects
freezing = temp_celsius < 0.0
print("n freezing cells:", int(freezing.sum()))    # True counts as 1
print("freezing values:", temp_celsius[freezing])  # 1D of matches

# vectorised: no python loop
temp_kelvin = temp_celsius + 273.15
print(temp_kelvin)  # shape (4,6)

lat_gradient_celsius = np.array([0.0, -1.5, -3.0, -4.5])   # colder toward the north

print(lat_gradient_celsius.shape)                  # (4,)   — a flat vector
print(lat_gradient_celsius[:, None].shape)         # (4, 1) — now a column
print(lat_gradient_celsius[:, np.newaxis].shape)   # (4, 1) — the same thing

# without the extra axis, the shapes cannot be aligned:
# temp_celsius + lat_gradient_celsius
# ValueError: operands could not be broadcast together with shapes (4,6) (4,)

# broadcasting: a (4,1) column stretches across the 6 longitudes
lat_gradient_celsius = np.array([0.0, -1.5, -3.0, -4.5])   # colder toward the north
adjusted = temp_celsius + lat_gradient_celsius[:, None]    # (4,1) + (4,6) -> (4,6)
print(adjusted.shape)
print(adjusted)

print(temp_celsius.mean(), np.mean(temp_celsius))   # method and function agree
print(temp_celsius.sum(), np.sum(temp_celsius))

# a reduction collapses an axis: axis=0 over rows, axis=1 over columns
print("overall mean:", temp_celsius.mean())
print("mean per column (over rows):", temp_celsius.mean(axis=0))
print("mean per row (over columns):", temp_celsius.mean(axis=1))
print("min, max:", temp_celsius.min(), temp_celsius.max())
print("argmin, argmax (flat index):", temp_celsius.argmin(), temp_celsius.argmax())

# reshape: same 24 values, new shape; -1 infers the missing length
flat = temp_celsius.reshape(-1)
print("reshaped:", flat.shape)

# stacking: combine arrays along a new row axis
col_means = temp_celsius.mean(axis=0)
index_row = np.arange(6, dtype=float)
print("vstacked:", np.vstack([index_row, col_means]).shape)

daily_rain_mm = np.array([0.0, 2.5, 0.0, 7.0, 1.5])
print(daily_rain_mm.cumsum())                  # [ 0.   2.5  2.5  9.5 11. ] — rain to date
print(daily_rain_mm.sum())                     # 11.0 — only the final total

# on a 2D field, axis works as for any reduction
print(temp_celsius.cumsum(axis=1).shape)       # (4, 6): a running total along each row

# np.where(condition, a, b): element-wise choice
category = np.where(temp_celsius < 0.0, "freezing", "above")
print(category)

# missing data as NaN; nan-aware vs ordinary reduction
temp_with_gaps = temp_celsius.copy()
temp_with_gaps[0, 0] = np.nan
print("nanmean (skips gaps):", np.nanmean(temp_with_gaps))
print("plain mean is contaminated:", np.mean(temp_with_gaps))  # nan
print("n missing:", np.isnan(temp_with_gaps).sum())      # np.isnan flags them; True counts as 1

# np.interp: fill a 1D gap by linear interpolation against an index
profile = temp_celsius[:, 0].copy()      # the first column, 4 latitudes
x = np.arange(profile.size)
known = np.array([0, 1, 3])              # pretend index 2 is missing
filled = np.interp(x, known, profile[known])
print(profile)
print(filled)                            # index 2 interpolated from its neighbours

present = ~np.isnan(temp_with_gaps)            # True where a value exists
print("n present:", present.sum(), "of", temp_with_gaps.size)
print(np.allclose(temp_with_gaps[present].mean(), np.nanmean(temp_with_gaps)))   # True

from pathlib import Path
Path("_files").mkdir(exist_ok=True)      # _files/ is not committed, so create it first

np.save("_files/temp_field.npy", temp_celsius)          # writes temp_field.npy

reloaded = np.load("_files/temp_field.npy")
print(reloaded.shape, reloaded.dtype)
print(np.allclose(temp_celsius, reloaded))       # True — an exact round trip

def to_anomaly(field):
    # subtract the mean into a preallocated array (as an assistant returned it)
    result = np.zeros_like(field)        # inherits field's dtype!
    result[:] = field - field.mean()
    return result

precip_mm = np.array([[0, 2, 5], [1, 0, 8], [3, 4, 2]])   # integer mm
print("input dtype:", precip_mm.dtype)
print(to_anomaly(precip_mm))

def to_anomaly(field):
    # let numpy promote to float; no wrong-dtype preallocation
    return field - field.mean()

print(to_anomaly(precip_mm))
print("output dtype:", to_anomaly(precip_mm).dtype)