PROJECT 103 · PYTHON SOURCE

NYC mobility operations

Study-specific code. Shared modules, dependency versions, and reproduction instructions are included in all project files.

from pathlib import Path
import sys
sys.path.insert(0, str(Path(__file__).resolve().parents[2]))
"""A multi-million-row mobility warehouse and future-month trip-duration model."""
from capstones.common import *
from capstones.fetch import SOURCES
import duckdb
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error
from threadpoolctl import threadpool_limits

SLUG='103-nyc-mobility-operations'

def run():
    p=project_dir(SLUG)
    sources=[(download(u,n),u) for u,n in SOURCES[:4]]
    provenance(sources,p)
    con=duckdb.connect(str(p/'data/mobility.duckdb'))
    con.execute("SET threads=4; SET memory_limit='2GB'")
    paths=[f.as_posix() for f,_ in sources[:3]]
    con.read_parquet(paths,union_by_name=True,filename=True).create_view('source_trips',replace=True)
    con.execute('CREATE OR REPLACE TABLE zones AS SELECT * FROM read_csv_auto(?)',[sources[3][0].as_posix()])
    sql="""-- Grain: one reported trip. No passenger or driver identifiers are available.
CREATE OR REPLACE TABLE trip_ledger AS
WITH measured AS (
 SELECT tpep_pickup_datetime AS pickup, tpep_dropoff_datetime AS dropoff,
 PULocationID AS pickup_zone, DOLocationID AS dropoff_zone,
 trip_distance AS miles, total_amount AS passenger_charge, fare_amount AS fare,
 tip_amount AS recorded_tip, payment_type,
 date_diff('second',tpep_pickup_datetime,tpep_dropoff_datetime)/60.0 AS minutes,
 regexp_extract(filename,'2025-([0-9]{2})',1) AS source_month
 FROM source_trips
)
SELECT m.*, coalesce(z.Borough,'Unknown') AS borough,
 CASE WHEN pickup IS NULL OR strftime(pickup,'%m')<>source_month OR year(pickup)<>2025 THEN 'Outside source month'
 WHEN minutes IS NULL OR minutes<1 OR minutes>180 THEN 'Duration outside 1–180 min'
 WHEN miles IS NULL OR miles<0.1 OR miles>100 THEN 'Distance outside 0.1–100 miles'
 WHEN passenger_charge IS NULL OR passenger_charge<=0 OR passenger_charge>1000 OR fare<0 THEN 'Invalid charge'
 WHEN z.LocationID IS NULL OR pickup_zone>263 OR dropoff_zone NOT BETWEEN 1 AND 263 THEN 'Unknown zone'
 ELSE 'Eligible' END AS quality_class
FROM measured m LEFT JOIN zones z ON m.pickup_zone=z.LocationID;
CREATE OR REPLACE VIEW eligible_trips AS SELECT * FROM trip_ledger WHERE quality_class='Eligible';
CREATE OR REPLACE TABLE daily_borough AS
SELECT CAST(pickup AS DATE) AS date, borough, count(*) AS trips,
sum(passenger_charge) AS passenger_charges_usd, sum(fare) AS metered_fares_usd,
sum(minutes) AS occupied_minutes, sum(miles) AS miles,
sum(CASE WHEN payment_type=1 THEN recorded_tip ELSE 0 END) AS card_tips_usd,
count(*) FILTER (WHERE payment_type=1) AS card_trips
FROM eligible_trips GROUP BY 1,2;
CREATE OR REPLACE TABLE zone_month AS
SELECT strftime(pickup,'%Y-%m') AS month, pickup_zone, z.Zone AS zone, e.borough,
count(*) AS trips, avg(minutes) AS mean_minutes, quantile_cont(minutes,.9) AS p90_minutes,
sum(passenger_charge) AS passenger_charges_usd
FROM eligible_trips e LEFT JOIN zones z ON pickup_zone=z.LocationID GROUP BY 1,2,3,4;
"""
    (p/'sql/warehouse.sql').write_text(sql,encoding='utf-8')
    con.execute(sql)
    raw=con.execute('SELECT count(*) FROM source_trips').fetchone()[0]
    audit=con.execute('SELECT quality_class,count(*) AS rows FROM trip_ledger GROUP BY 1 ORDER BY 2 DESC').df()
    daily=con.execute('SELECT * FROM daily_borough ORDER BY date,borough').df()
    zones=con.execute('SELECT * FROM zone_month ORDER BY trips DESC').df()
    save_csv(audit,p/'outputs/data_quality.csv');save_csv(daily,p/'outputs/daily_borough.csv');save_csv(zones,p/'outputs/zone_month.csv')
    print(f'Taxi warehouse: {raw:,} reported rows; {int(daily.trips.sum()):,} eligible',flush=True)
    # Deterministic hash ordering within disjoint months bounds model training costs.
    frames=[]
    for month,size in [('04',120000),('05',60000),('06',60000)]:
        frame=con.execute(f"""SELECT * FROM (
          SELECT pickup, dropoff, pickup_zone, dropoff_zone, hour(pickup) AS hour,
          dayofweek(pickup) AS weekday, minutes, borough
          FROM eligible_trips WHERE source_month='{month}' AND strftime(dropoff,'%m')='{month}'
        ) ORDER BY hash(42, pickup, dropoff, pickup_zone, dropoff_zone, minutes) LIMIT {size}""").df()
        frames.append(frame)
    train,val,test=frames
    features=['pickup_zone','dropoff_zone','hour','weekday']
    baseline=float(train.minutes.median())
    lookup=train.groupby(['pickup_zone','dropoff_zone','hour']).minutes.median()
    def base(frame):
        keys=pd.MultiIndex.from_frame(frame[['pickup_zone','dropoff_zone','hour']])
        return lookup.reindex(keys).fillna(baseline).to_numpy()
    candidates=[];models=[]
    with threadpool_limits(limits=3):
        for leaves in [15,31]:
            model=HistGradientBoostingRegressor(loss='absolute_error',max_iter=110,max_leaf_nodes=leaves,learning_rate=.12,l2_regularization=2,early_stopping=False,random_state=42)
            model.fit(train[features],train.minutes)
            candidates.append({'model':f'HGB {leaves} leaves','validation_mae_minutes':mean_absolute_error(val.minutes,model.predict(val[features]))})
            models.append(model)
        chosen=int(np.argmin([row['validation_mae_minutes'] for row in candidates]))
        selected=models[chosen]
        test=test.copy();test['model']=selected.predict(test[features]);test['baseline']=base(test)
    # Fixed training fit; the validation month selects configuration only.
    test['day']=test.pickup.dt.strftime('%Y-%m-%d')
    test['absolute_error']=(test.model-test.minutes).abs()
    rows=[]
    for label,pred in [('Training median',np.repeat(baseline,len(test))),('Route-hour median',test.baseline),('Gradient boosting',test.model)]:
        rows.append({'model':label,'test_rows':len(test),'mae_minutes':mean_absolute_error(test.minutes,pred),'rmse_minutes':np.sqrt(mean_squared_error(test.minutes,pred))})
    results=pd.DataFrame(rows);save_csv(results,p/'outputs/results.csv')
    slices=test.groupby('borough').apply(lambda g: pd.Series({'test_trips':len(g),'model_mae_minutes':mean_absolute_error(g.minutes,g.model),'baseline_mae_minutes':mean_absolute_error(g.minutes,g.baseline)}),include_groups=False).reset_index()
    save_csv(slices,p/'outputs/borough_errors.csv')
    save_csv(test,p/'data/test_predictions.csv')
    interval=bootstrap_days(test,'baseline','model','minutes','day')
    improved=100*(1-results.iloc[2].mae_minutes/results.iloc[1].mae_minutes)
    total=int(daily.trips.sum())
    metrics={'source_rows':raw,'eligible_trips':total,'excluded_rows':raw-total,'date_start':'2025-04-01','date_end':'2025-06-30','model_mae_minutes':float(results.iloc[2].mae_minutes),'baseline_mae_minutes':float(results.iloc[1].mae_minutes),'improvement_pct':improved,'passenger_charges_usd':float(daily.passenger_charges_usd.sum()),'uncertainty':interval,'sampling':'120,000 April training trips, 60,000 May validation trips, 60,000 June test trips; deterministic hash ordering with salt 42; exact duplicate model rows are interchangeable'}
    js(p/'outputs/metrics.json',metrics)
    js(p/'outputs/evaluation.json',{'split':'April 2025 training, May validation, June untouched test; trips ending outside their source month excluded from modeling. Warehouse uses all eligible rows; models use deterministic within-month hash samples.','train_rows':len(train),'test_rows':len(test),'selection':f'Configuration chosen by May MAE: {candidates[chosen]["model"]}. No refitting or selection on June.','target':'Trip duration in minutes','features':features,'candidates':candidates,'uncertainty':interval,'limitations':'Zones are integer-encoded, so tree splits impose an ordering. Baseline route-hour groups unseen in April use the training median. No traffic, vehicle availability, or route geometry is observed.'})
    quality(p,{'source_count_reconciles':int(audit.rows.sum())==raw,'eligible_daily_reconciles':total==con.execute('SELECT count(*) FROM eligible_trips').fetchone()[0],'zone_month_reconciles':int(zones.trips.sum())==total,'dimension_unique':con.execute('SELECT count(*)=count(DISTINCT LocationID) FROM zones').fetchone()[0],'daily_grain_unique':not daily.duplicated(['date','borough']).any(),'train_before_validation':train.dropoff.max()<val.pickup.min(),'validation_before_test':val.dropoff.max()<test.pickup.min(),'no_post_trip_predictors':not set(features)&{'minutes','miles','fare','passenger_charge','recorded_tip','payment_type','dropoff'},'finite_predictions':bool(np.isfinite(test.model).all()),'expected_source_months':set(zones.month)=={'2025-04','2025-05','2025-06'}})
    fig,axes=plt.subplots(2,2,figsize=(12,7.8))
    trend=daily.groupby('date').trips.sum();axes[0,0].plot(trend.index,trend/1000,color=BLUE);axes[0,0].set(title='Recorded eligible trips by day',ylabel='Thousands of trips');axes[0,0].tick_params(axis='x',rotation=25)
    byborough=daily.groupby('borough').trips.sum().sort_values();axes[0,1].barh(byborough.index,byborough/1e6,color=TEAL);axes[0,1].set(title='Pickup borough',xlabel='Millions of eligible trips')
    axes[1,0].barh(results.model,results.mae_minutes,color=[ORANGE,TEAL,BLUE]);axes[1,0].set(title='Future-month duration benchmark',xlabel='June mean absolute error · minutes')
    x=np.arange(len(slices));axes[1,1].bar(x-.18,slices.baseline_mae_minutes,.36,label='Route-hour median',color=TEAL);axes[1,1].bar(x+.18,slices.model_mae_minutes,.36,label='Gradient boosting',color=BLUE);axes[1,1].set_xticks(x,slices.borough,rotation=25);axes[1,1].set(title='Error by pickup borough',ylabel='MAE · minutes');axes[1,1].legend(fontsize=8)
    fig.suptitle('NYC mobility operations | April–June 2025',fontsize=17,fontweight='bold');fig.tight_layout();fig.savefig(p/'outputs/chart.png',dpi=150);plt.close(fig)
    js(p/'outputs/dashboard.json',{'kind':'taxi','metrics':metrics,'daily':json.loads(daily.assign(date=daily.date.dt.strftime('%Y-%m-%d')).to_json(orient='records')),'zones':json.loads(zones.to_json(orient='records')),'scores':rows,'slices':json.loads(slices.to_json(orient='records'))})
    con.close()
    return metrics

if __name__=='__main__': print(json.dumps(run(),indent=2))