Implement SIP

.gitignore

@@ -41,6 +41,7 @@ dist/
 build/
 pip
 tutorials/IC2S2-2025/tutorial-sample/*
+tutorials/SICSS_2026/data/
 
 # Development environment configurations (personal/local only)
 AGENTS.md

nomad/constants.py

@@ -22,6 +22,7 @@
     "location_id" : "location_id"}
 
 SEC_PER_UNIT = {'s': 1, 'min': 60, 'h': 3_600, 'd': 86_400, 'w': 604_800}
+EARTH_RADIUS_METERS = 6_371_000
 
 ALLOWED_BUILDINGS = {
     0: ['home'], 1: ['home'], 2: ['home'], 3: ['home'], 4: ['home'],

nomad/contact_estimation.py

@@ -6,6 +6,301 @@
 import nomad.io.base as loader
 import numpy as np
 import pandas as pd
+from nomad.stop_detection.utils import _haversine_distance
+from sklearn.neighbors import BallTree, KDTree
+from constants import _EARTH_RADIUS_M
+
+_TEMPORAL_BLOCK_SECONDS = 60 * 60
+
+def _temporal_blocks(start, end):
+    """Return stop-to-hour block memberships for non-empty intervals."""
+    valid = end > start
+    if not valid.any():
+        return pd.DataFrame(columns=["block", "stop"])
+
+    stop = np.flatnonzero(valid)
+    start_bin = np.floor_divide(start[valid], _TEMPORAL_BLOCK_SECONDS).astype(np.int64)
+    end_bin = np.floor_divide(end[valid] - 1, _TEMPORAL_BLOCK_SECONDS).astype(np.int64)
+    counts = end_bin - start_bin + 1
+    offsets = np.repeat(np.r_[0, counts.cumsum()[:-1]], counts)
+    block = np.repeat(start_bin, counts) + np.arange(counts.sum()) - offsets
+    return pd.DataFrame({"block": block, "stop": np.repeat(stop, counts)})
+
+
+def _radius_candidates(
+    stops,
+    traj_cols,
+    input_traj_cols,
+    kwargs,
+    distance_threshold,
+    start,
+    end,
+    users,
+):
+    """Return stop row pairs within distance_threshold and their distances."""
+    coord_key1, coord_key2, use_lon_lat = loader._fallback_spatial_cols(
+        stops.columns,
+        input_traj_cols,
+        kwargs,
+    )
+
+    if use_lon_lat:
+        # Haversine uses (lat, lon) in radians.
+        query_coords = np.radians(
+            stops[[traj_cols["latitude"], traj_cols["longitude"]]].to_numpy()
+        )
+        radius = distance_threshold / _EARTH_RADIUS_M
+        tree_class = BallTree
+        tree_kwargs = {"metric": "haversine"}
+    else:
+        query_coords = stops[[traj_cols[coord_key1], traj_cols[coord_key2]]].to_numpy()
+        radius = distance_threshold
+        tree_class = KDTree
+        tree_kwargs = {}
+
+    blocks = _temporal_blocks(start, end)
+    rows = []
+    cols = []
+    for _, block_stops in blocks.groupby("block", sort=False)["stop"]:
+        block_stops = block_stops.to_numpy()
+        if len(block_stops) < 2:
+            continue
+
+        tree = tree_class(query_coords[block_stops], **tree_kwargs)
+        indices = tree.query_radius(
+            query_coords[block_stops],
+            r=radius,
+            return_distance=False,
+            sort_results=False,
+        )
+        counts = np.array([len(idx) for idx in indices])
+        row = block_stops[np.repeat(np.arange(len(block_stops)), counts)]
+        col = block_stops[np.concatenate(indices).astype(int)]
+        keep = row < col
+        rows.append(row[keep])
+        cols.append(col[keep])
+
+    if not rows:
+        return np.empty(0, dtype=int), np.empty(0, dtype=int), np.empty(0, dtype=float)
+
+    pairs = pd.DataFrame(
+        {"stop_1": np.concatenate(rows), "stop_2": np.concatenate(cols)}
+    ).drop_duplicates()
+    stop_1 = pairs["stop_1"].to_numpy()
+    stop_2 = pairs["stop_2"].to_numpy()
+    keep = (
+        (users[stop_1] != users[stop_2])
+        & (start[stop_1] < end[stop_2])
+        & (start[stop_2] < end[stop_1])
+    )
+    stop_1, stop_2 = stop_1[keep], stop_2[keep]
+    if use_lon_lat:
+        distance = _haversine_distance(
+            query_coords[stop_1].T,
+            query_coords[stop_2].T,
+        )
+    else:
+        distance = np.linalg.norm(query_coords[stop_1] - query_coords[stop_2], axis=1)
+    keep = distance <= distance_threshold
+    return stop_1[keep], stop_2[keep], distance[keep]
+
+def estimate_contacts(stops, distance_threshold=None, traj_cols=None, **kwargs):
+    """
+    Estimate undirected co-location/contact events from a stop table.
+
+    With ``distance_threshold=None``, contacts are stops from different users at
+    the same ``location_id`` with strictly overlapping times. With a distance
+    threshold, contacts are stops from different users within that radius and
+    with strictly overlapping times. Latitude/longitude thresholds are meters;
+    projected x/y thresholds use the coordinate units.
+
+    Parameters
+    ----------
+    stops : pd.DataFrame
+        Stop table with user_id, a start time (datetime or timestamp), an end
+        time or duration, and either location_id (exact mode) or coordinates
+        (radius mode).
+    distance_threshold : float, optional
+        Contact radius. None selects exact-location mode. Meters for lat/lon,
+        coordinate units for projected x/y.
+    traj_cols : dict, optional
+        Mapping for user_id, time, duration, location_id, and coordinates.
+    **kwargs
+        Column-name overrides forwarded to the traj_cols resolver.
+
+    Returns
+    -------
+    pd.DataFrame
+        Columns = [user_id_1, user_id_2, contact_start, contact_end,
+        overlap_duration, location_id or distance]. overlap_duration is in
+        minutes. Contact times are returned as datetime for datetime input and
+        as Unix seconds for timestamp input; timezone metadata on datetime
+        input is not preserved (output is timezone-naive).
+    """
+    contact_col = "location_id" if distance_threshold is None else "distance"
+    output_cols = [
+        "user_id_1",
+        "user_id_2",
+        "contact_start",
+        "contact_end",
+        "overlap_duration",
+        contact_col,
+    ]
+    if len(stops) == 0:
+        return pd.DataFrame(columns=output_cols)
+
+    input_traj_cols = traj_cols
+    traj_cols = loader._parse_traj_cols(stops.columns, traj_cols, kwargs)
+
+    uid_col = traj_cols["user_id"]
+    if uid_col not in stops.columns:
+        raise ValueError(f"Missing required user_id column '{uid_col}'.")
+
+    t_key, use_datetime = loader._fallback_time_cols_dt(
+        stops.columns,
+        input_traj_cols,
+        kwargs,
+    )
+    start_col = traj_cols[t_key]
+    end_key = "end_datetime" if use_datetime else "end_timestamp"
+    end_col = traj_cols[end_key]
+    if end_col in stops.columns:
+        end_time = stops[end_col]
+    elif loader._has_duration_cols(stops.columns, traj_cols):
+        if use_datetime:
+            end_time = stops[start_col] + pd.to_timedelta(
+                stops[traj_cols["duration"]],
+                unit="m",
+            )
+        else:
+            end_time = stops[start_col] + stops[traj_cols["duration"]] * 60
+    else:
+        raise ValueError(
+            "Contact estimation requires an end time or duration column."
+        )
+
+    if use_datetime:
+        start = filters.to_timestamp(stops[start_col]).to_numpy()
+        end = filters.to_timestamp(end_time).to_numpy()
+    else:
+        start = stops[start_col].to_numpy()
+        end = end_time.to_numpy()
+
+    users = stops[uid_col].to_numpy()
+
+    if distance_threshold is None:
+        loc_col = traj_cols["location_id"]
+        if loc_col not in stops.columns:
+            raise ValueError(
+                "Exact-location contact estimation requires a location_id column."
+            )
+        if stops[loc_col].isna().any():
+            raise ValueError(
+                "Exact-location contact estimation requires non-missing location_id values."
+            )
+        location = stops[loc_col].to_numpy()
+        # Pair stops sharing a location, then keep distinct users with overlap.
+        candidates = pd.DataFrame({"stop": np.arange(len(stops)), "loc": location})
+        pairs = candidates.merge(candidates, on="loc", suffixes=("_1", "_2"))
+        stop_1 = pairs["stop_1"].to_numpy()
+        stop_2 = pairs["stop_2"].to_numpy()
+        keep = (
+            (stop_1 < stop_2)
+            & (users[stop_1] != users[stop_2])
+            & (start[stop_1] < end[stop_2])
+            & (start[stop_2] < end[stop_1])
+        )
+        stop_1, stop_2 = stop_1[keep], stop_2[keep]
+        contact_values = location[stop_1]
+    else:
+        stop_1, stop_2, distance = _radius_candidates(
+            stops,
+            traj_cols,
+            input_traj_cols,
+            kwargs,
+            distance_threshold,
+            start,
+            end,
+            users,
+        )
+        contact_values = distance
+
+    contact_start = np.maximum(start[stop_1], start[stop_2])
+    contact_end = np.minimum(end[stop_1], end[stop_2])
+    return pd.DataFrame(
+        {
+            "user_id_1": users[stop_1],
+            "user_id_2": users[stop_2],
+            "contact_start": (
+                pd.to_datetime(contact_start, unit="s") if use_datetime else contact_start
+            ),
+            "contact_end": (
+                pd.to_datetime(contact_end, unit="s") if use_datetime else contact_end
+            ),
+            "overlap_duration": ((contact_end - contact_start) // 60).astype(int),
+            contact_col: contact_values,
+        },
+        columns=output_cols,
+    )
+
+
+def compute_contact_weights(
+    contacts,
+    method="duration",
+    distance_threshold=None,
+    overlap_duration_col="overlap_duration",
+    distance_col="distance",
+):
+    """
+    Compute contact weights from a contact event table.
+
+    Supported methods are ``"duration"`` and ``"linear_distance"``. Linear
+    distance weighting uses ``overlap_duration * max(0, 1 - distance / threshold)``.
+
+    Parameters
+    ----------
+    contacts : pd.DataFrame
+        Contact event table from estimate_contacts.
+    method : str
+        Either 'duration' or 'linear_distance'.
+    distance_threshold : float, optional
+        Contact radius used by 'linear_distance'; required for that method.
+    overlap_duration_col : str
+        Column holding overlap duration in minutes.
+    distance_col : str
+        Column holding contact distance; required for 'linear_distance'.
+
+    Returns
+    -------
+    pandas.Series
+        Series named ``contact_weight`` and indexed like ``contacts``.
+    """
+    if overlap_duration_col not in contacts.columns:
+        raise ValueError(
+            f"Missing required overlap duration column '{overlap_duration_col}'."
+        )
+
+    method = method.lower()
+
+    if method == "duration":
+        return contacts[overlap_duration_col].rename("contact_weight")
+
+    if method != "linear_distance":
+        raise ValueError("method must be one of 'duration' or 'linear_distance'.")
+
+    if distance_col not in contacts.columns:
+        raise ValueError(
+            f"method='linear_distance' requires a '{distance_col}' column."
+        )
+    if distance_threshold is None:
+        raise ValueError(
+            "method='linear_distance' requires an explicit distance_threshold."
+        )
+
+    return (
+        contacts[overlap_duration_col]
+        * np.maximum(0, 1 - contacts[distance_col] / distance_threshold)
+    ).rename("contact_weight")
 
 
 def overlapping_visits(left, right, match_location=False, traj_cols=None, right_traj_cols=None, **kwargs):