// Accept two trajectory slices, each a list of (latitude, longitude) tuples in decimal degrees.
fn dtw_distance(a: &[(f64, f64)], b: &[(f64, f64)]) -> f64 {
    let n = a.len();                                 // Number of position reports in trajectory a.
    let m = b.len();                                 // Number of position reports in trajectory b.
    let mut dtw = vec![vec![f64::INFINITY; m + 1]; n + 1]; // (n+1)×(m+1) matrix initialised to infinity; extra row/col for the base case.
    dtw[0][0] = 0.0;                                 // Distance between two empty prefixes is zero.

    for i in 1..=n {                                 // Row index: iterate over every prefix of trajectory a.
        for j in 1..=m {                             // Column index: iterate over every prefix of trajectory b.
            let cost = haversine(a[i - 1].0, a[i - 1].1, b[j - 1].0, b[j - 1].1); // Great-circle distance (km) between the i‑th AIS point and j‑th AIS point.
            let min_prev = dtw[i - 1][j - 1]         // Cost of aligning the two previous points (no stretch).
                .min(dtw[i - 1][j])                  // Cost of stretching trajectory a (insertion in b).
                .min(dtw[i][j - 1]);                 // Cost of stretching trajectory b (insertion in a).
            dtw[i][j] = cost + min_prev;             // Accumulated DTW distance for prefixes i and j.
        }
    }

    dtw[n][m]                                        // Bottom-right cell holds the full DTW distance between a and b.
}

// Haversine formula: angular distance on a sphere → great-circle distance in kilometres.
fn haversine(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
    let r = 6371.0;                                  // Earth mean radius in km (WGS‑84 approximation).
    let d_lat = (lat2 - lat1).to_radians();           // Difference in latitude, converted from degrees to radians.
    let d_lon = (lon2 - lon1).to_radians();           // Difference in longitude, converted from degrees to radians.
    let a = (d_lat / 2.0).sin().powi(2)               // Square of half‑the‑difference sine (central angle numerator, part 1).
        + lat1.to_radians().cos() * lat2.to_radians().cos() * (d_lon / 2.0).sin().powi(2); // Part 2: cos product × squared sine of half longitude difference.
    let c = 2.0 * a.sqrt().asin();                   // Central angle in radians via arcsin of sqrt(a).
    r * c                                             // Arc length on Earth's surface = radius × central angle.
}

use rayon::prelude::*;                               // Brings in the ParallelIterator trait that adds .par_iter() to iterables.

// Compute DTW distance for every unique unordered pair of vessel trajectories, returning (index_a, index_b, distance_km).
fn pairwise_dtw(trajectories: &[Vec<(f64, f64)>]) -> Vec<(usize, usize, f64)> {
    let n = trajectories.len();                      // Total number of vessel trajectories in the dataset.
    let mut pairs: Vec<(usize, usize)> = Vec::with_capacity(n * (n - 1) / 2); // Pre‑allocate for n‑choose‑2 unique pairs (e.g. ~50M for 10k trajectories).

    for i in 0..n {                                  // Outer loop: first trajectory index.
        for j in (i + 1)..n {                        // Inner loop: second trajectory index, always after i to avoid duplicates and self‑pairs.
            pairs.push((i, j));                      // Store the pair indices for later parallel processing.
        }
    }

    pairs
        .par_iter()                                  // Rayon parallel iterator: pairs are distributed across the thread pool via work stealing.
        .map(|&(i, j)| {                             // Each closure runs on one thread, receiving a single (i, j) pair.
            let dist = dtw_distance(&trajectories[i], &trajectories[j]); // Full sequential DTW on this pair — the DP itself stays single‑threaded.
            (i, j, dist)                             // Return the pair indices and their DTW distance for downstream analysis (e.g. clustering).
        })
        .collect()                                   // Gather all results back into a single Vec<(usize, usize, f64)>.
}

// Memory‑reuse variant: the caller provides a pre‑allocated scratch buffer so we avoid allocating an L×L matrix per call.
fn dtw_distance_scratch(
    a: &[(f64, f64)],                                // First vessel trajectory: slice of (lat, lon) position reports.
    b: &[(f64, f64)],                                // Second vessel trajectory: slice of (lat, lon) position reports.
    scratch: &mut [f64],                             // Flat scratch buffer, sized (max_len+1)×(max_len+1) so it can hold the largest matrix needed.
) -> f64 {
    let n = a.len();                                 // Length of trajectory a.
    let m = b.len();                                 // Length of trajectory b.
    let stride = m + 1;                              // Number of columns per row in the flat representation (m+1 to include the zero‑index column).

    scratch.fill(f64::INFINITY);                     // Reset every cell to +infinity — the DTW recurrence only writes visited cells.
    scratch[0] = 0.0;                                // Base case: DTW distance of two empty prefixes is zero.

    for i in 1..=n {                                 // Row i of the virtual matrix (prefix of trajectory a).
        for j in 1..=m {                             // Column j (prefix of trajectory b).
            let cost = haversine(a[i - 1].0, a[i - 1].1, b[j - 1].0, b[j - 1].1); // Haversine distance between i‑th and j‑th AIS position reports.
            let prev = scratch[(i - 1) * stride + (j - 1)] // Cell (i-1, j-1): aligning point i-1 with point j-1 (no warping).
                .min(scratch[(i - 1) * stride + j])        // Cell (i-1, j):   skipping a point in trajectory b (warping a).
                .min(scratch[i * stride + (j - 1)]);       // Cell (i, j-1):   skipping a point in trajectory a (warping b).
            scratch[i * stride + j] = cost + prev;   // Write the new accumulated distance into the flat buffer at position (i, j).
        }
    }

    scratch[n * stride + m]                          // Read the bottom‑right cell, which holds the full DTW distance between a and b.
}

// Parallel pairwise DTW with per‑thread scratch‑buffer reuse — avoids millions of heap allocations.
fn pairwise_dtw_fast(trajectories: &[Vec<(f64, f64)>]) -> Vec<(usize, usize, f64)> {
    let n = trajectories.len();                      // Number of vessel trajectories in the corpus.
    let max_len = trajectories.iter().map(|t| t.len()).max().unwrap_or(0); // Longest trajectory length, used to size the scratch buffer.
    let scratch_size = (max_len + 1) * (max_len + 1); // Square matrix size for the worst‑case (max_len×max_len) DTW, plus the zero row/col.

    let mut pairs: Vec<(usize, usize)> = Vec::with_capacity(n * (n - 1) / 2); // Pre‑allocate for all n‑choose‑2 unique pairs.
    for i in 0..n {                                  // Outer index over vessel trajectories.
        for j in (i + 1)..n {                        // Inner index: only j > i to avoid computing both (i,j) and (j,i).
            pairs.push((i, j));                      // Collect the pair into the work list.
        }
    }

    pairs
        .par_iter()                                  // Rayon distributes the work list across all available CPU cores.
        .map_init(                                   // Like .map() but with an init closure that runs once per Rayon worker thread.
            || vec![0.0_f64; scratch_size],          // Each thread allocates exactly one scratch buffer sized for the longest trajectory.
            |scratch, &(i, j)| {                     // Per‑pair closure receives the thread‑local scratch buffer and one (i,j) pair.
                let dist = dtw_distance_scratch(&trajectories[i], &trajectories[j], scratch); // DTW using the reused buffer — no allocation inside this hot loop.
                (i, j, dist)                         // Return the trajectory indices and their DTW distance in kilometres.
            },
        )
        .collect()                                   // Gather every thread's results into a single output vector.
}

use serde::Deserialize;                              // Serde's auto‑derive for deserialising JSON from the AISHub API response.

#[derive(Deserialize, Debug)]                        // Auto‑implement Deserialize (from JSON) and Debug (for logging).
struct AisReport {
    mmsi: u32,                                       // Maritime Mobile Service Identity — a unique 9‑digit vessel identifier (like a ship's "license plate").
    latitude: f64,                                   // Vessel latitude in decimal degrees at report time (e.g. 33.89075).
    longitude: f64,                                  // Vessel longitude in decimal degrees at report time (e.g. -128.26693).
    time: String,                                    // ISO‑like timestamp string (e.g. "2026-07-17 12:09:25 GMT") — used for chronological ordering.
}

// Group raw AISHub reports by MMSI, sort each vessel's reports by time, and extract (lat, lon) sequences.
fn build_trajectories(reports: Vec<AisReport>) -> Vec<Vec<(f64, f64)>> {
    use std::collections::HashMap;                   // Standard hash map for grouping reports by vessel MMSI.

    let mut by_mmsi: HashMap<u32, Vec<AisReport>> = HashMap::new(); // Map from MMSI to its list of unsorted position reports.
    for r in reports {                               // Iterate over every parsed AIS report from the API response.
        by_mmsi.entry(r.mmsi).or_default().push(r);  // Insert into the per‑MMSI bucket, creating an empty vec if this is the first report for that vessel.
    }

    by_mmsi
        .into_values()                               // Take ownership of each vessel's report list, discarding the MMSI keys.
        .map(|mut traj| {                            // For each vessel's collection of reports:
            traj.sort_by(|a, b| a.time.cmp(&b.time)); // Sort chronologically by the timestamp string so the trajectory follows the vessel's actual path.
            traj.into_iter().map(|r| (r.latitude, r.longitude)).collect() // Keep only (lat, lon) pairs — discard MMSI and time now that ordering is established.
        })
        .collect()                                   // Collect all vessel trajectories into a Vec<Vec<(f64, f64)>> ready for DTW.
}

// DTW constrained by a Sakoe‑Chiba band of width w — limits the warping path to stay near the diagonal, giving O(n×w) instead of O(n×m).
fn dtw_band(a: &[(f64, f64)], b: &[(f64, f64)], w: usize) -> f64 {
    let n = a.len();                                 // Number of position reports in vessel trajectory a.
    let m = b.len();                                 // Number of position reports in vessel trajectory b.
    let mut dtw = vec![vec![f64::INFINITY; m + 1]; n + 1]; // Full matrix initialised to infinity; unvisited band cells remain ∞ and are never picked.
    dtw[0][0] = 0.0;                                 // Base case: distance between two empty prefixes.

    for i in 1..=n {                                 // Row index (prefix of trajectory a).
        let j_start = 1.max(i.saturating_sub(w));    // Left column boundary: stay at least 1 and at most w columns before the diagonal.
        let j_end = m.min(i + w);                    // Right column boundary: stay at most w columns after the diagonal, but don't exceed m.
        for j in j_start..=j_end {                   // Only visit cells inside the diagonal band — cells outside stay INFINITY.
            let cost = haversine(a[i - 1].0, a[i - 1].1, b[j - 1].0, b[j - 1].1); // Haversine distance between the i‑th and j‑th AIS position.
            dtw[i][j] = cost + dtw[i - 1][j - 1]     // Accumulated cost: this cell's haversine distance plus the cheapest predecessor.
                .min(dtw[i - 1][j])                  // Predecessor from the same column (warping trajectory a).
                .min(dtw[i][j - 1]);                 // Predecessor from the same row (warping trajectory b).
        }
    }

    dtw[n][m]                                        // Final cell — if INFINITY the band was too narrow and the trajectories cannot be aligned within it.
}

use clap::Parser;                                    // Clap's derive macro for declarative command‑line argument parsing.

#[derive(Parser)]                                    // Auto‑generate a CLI parser from struct field annotations.
struct Args {
    #[arg(short, long)]                              // Short flag `-a`, long flag `--api-key` (clap auto‑kebabifies `api_key`).
    api_key: String,                                 // AISHub API key required to authenticate and fetch vessel position data.
    #[arg(short, long, default_value = "10")]        // `-m` / `--min-positions`, defaults to 10; vessels with fewer reports are filtered out.
    min_positions: usize,                            // Minimum number of position reports needed to keep a vessel trajectory.
}

fn main() -> Result<(), Box<dyn std::error::Error>> { // Return a generic error type so we can use `?` for fallible calls.
    let args = Args::parse();                        // Parse command‑line arguments (exits with help text on error or `--help`).
    let raw = fetch_ais_data(&args.api_key)?;        // Fetch and deserialise AIS reports from the AISHub API using the provided key.
    let trajectories = build_trajectories(raw)       // Group by MMSI, sort by timestamp, extract (lat, lon) sequences.
        .into_iter()
        .filter(|t| t.len() >= args.min_positions)   // Drop vessels that don't have enough position reports for meaningful DTW.
        .collect::<Vec<_>>();                         // Collect the filtered trajectories into a Vec for indexing and parallel iteration.

    println!("Comparing {} trajectories", trajectories.len()); // Log how many vessel trajectories will be compared.

    let pairs = pairwise_dtw_fast(&trajectories);    // Compute all pairwise DTW distances in parallel, reusing per‑thread scratch buffers.

    for (i, j, dist) in pairs.iter().take(20) {      // Print the first 20 pairwise distances as a quick sanity check.
        println!("Trajectory {} <-> {}: {:.2} km", i, j, dist); // Format distance to two decimal places — typical output: "Trajectory 3 <-> 17: 42.15 km".
    }

    Ok(())                                           // Success exit code.
}