#!/usr/bin/env Rscript
# ==============================================================================
# GOMD Official R Package v2.0.0 — Full Live Function Verification Suite
# Executes all exported functions across client, oceanography, maritime
# decarbonization, non-parametric/causal econometrics, and spatial geometry.
# ==============================================================================

library(gomd)

cat("================================================================================\n")
cat("GOMD R PACKAGE v2.0.0 — FULL LIVE FUNCTION VERIFICATION SUITE\n")
cat("================================================================================\n\n")

# 1. Citations & Client
cat("[MODULE 1] Client, Provenance & Citations\n")
gomd_cite("apa")
prov <- gomd_verify_provenance('{"source":"copernicus_cmems","indicator":"marmara_sst","value":24.6}')
cat(sprintf("  [OK] gomd_verify_provenance() -> sha256=%s\n", substr(prov$sha256_digest, 1, 24)))

# 2. Physical & Biogeochemical Oceanography + Post-2020 Frontier Methods
cat("\n[MODULE 2] Physical & Biogeochemical Oceanography (9 Functions)\n")
eos <- gomd_seawater_density(salinity_psu = 18.2, temp_c = 15.4, pressure_dbar = 0)
cat(sprintf("  [OK] gomd_seawater_density()        -> rho=%.4f kg/m3\n", eos$rho_kg_m3))

teos <- gomd_teos10_stratification(salinity_psu = 38.5, temp_c = 14.6, pressure_dbar = 25.0,
                                   lat = 40.70, surface_salinity_psu = 21.4, surface_temp_c = 23.5)
cat(sprintf("  [OK] gomd_teos10_stratification()   -> SA=%.4f g/kg, N=%.2f cph\n",
            teos$absolute_salinity_sa_g_kg, teos$buoyancy_frequency_cph))

ox <- gomd_oxygen_solubility(temp_c = 14.5, salinity_psu = 38.5, observed_do_umol_kg = 65.0)
cat(sprintf("  [OK] gomd_oxygen_solubility()       -> DO_sat=%.3f mg/L (%s)\n", ox$do_sat_mg_l, ox$hypoxia_status))

carb <- gomd_carbonate_system(ph_total = 8.04, temp_c = 20.0, salinity_psu = 38.5, ta_umol_kg = 2610.0)
cat(sprintf("  [OK] gomd_carbonate_system()        -> pCO2=%.2f uatm, Omega_arag=%.3f\n",
            carb$pCO2_uatm, carb$omega_aragonite))

# Pull live time series from GOMD (https://gomd.pirireis.edu.tr) — zero synthetic vectors
sst_df <- gomd_get_series("copernicus_sst_med", entity = "TUR")
ph_df  <- gomd_get_series("socat_surface_ph_total_scale", entity = "TUR")
do_df  <- gomd_get_series("marmara_subhalocline_do_mg_l", entity = "OCN_MARMARA")
co2_df <- gomd_get_series("gml_mauna_loa_co2_ppm", entity = "WLD")
wav_df <- gomd_get_series("metocean_global_basin_wave_height_max_m", entity = "OCN_MED_AEGEAN")

sst_30d <- as.numeric(sst_df$value)
mhw <- gomd_detect_mhw(sst_30d, pct_threshold = 80, min_duration = 2)
cat(sprintf("  [OK] gomd_detect_mhw()              -> events=%d, threshold=%.2f C\n", mhw$event_count, mhw$threshold_sst_c))

n_common <- min(nrow(sst_df), nrow(ph_df), nrow(do_df))
comp <- gomd_compound_mhw_oax(
  as.numeric(sst_df$value[seq_len(n_common)]),
  as.numeric(ph_df$value[seq_len(n_common)]),
  as.numeric(do_df$value[seq_len(n_common)]),
  sst_pct = 75, ph_pct = 25, do_pct = 25
)
cat(sprintf("  [OK] gomd_compound_mhw_oax()        -> triple_steps=%d, mean_AGI=%.3f\n",
            comp$triple_compound_steps, comp$mean_aerobic_growth_index))

tsi <- gomd_trophic_index(chl_a_ug_l = 8.4, din_ug_l = 310, tp_ug_l = 46, do_saturation_pct = 62, secchi_m = 2.1)
cat(sprintf("  [OK] gomd_trophic_index()           -> TSI=%.2f, TRIX=%.3f\n", tsi$carlson_tsi_chla, tsi$vollenweider_trix))

muc <- gomd_mucilage_risk(sst_anomaly_c = 2.4, n_p_molar_ratio = 28.5, wind_speed_ms = 1.8, residence_time_days = 42, chl_a_ug_l = 9.2)
cat(sprintf("  [OK] gomd_mucilage_risk()           -> score=%.2f/100 (%s)\n", muc$mucilage_risk_score, muc$risk_level))

bcarb <- gomd_blue_carbon_valuation(seagrass_ha = 12500, saltmarsh_ha = 3200, mangrove_ha = 800)
cat(sprintf("  [OK] gomd_blue_carbon_valuation()   -> annual=%.0f tCO2e/yr, NPV=$%.2f\n",
            bcarb$annual_sequestration_tco2e_yr, bcarb$natural_capital_npv_scc_usd))

# 3. Maritime Decarbonization & Green Finance
cat("\n[MODULE 3] Maritime Decarbonization & Green Finance (7 Functions)\n")
cii <- gomd_calc_cii(fuel_tonnes = 11800, dwt = 82000, distance_nm = 58400, ship_type = "bulk_carrier", year = 2026)
cat(sprintf("  [OK] gomd_calc_cii()                -> Rating=%s (AER=%.4f)\n", cii$cii_rating, cii$attained_aer_g_dwt_nm))

ets <- gomd_calc_eu_ets(co2_tonnes = 28500, ch4_tonnes = 12.0, n2o_tonnes = 1.4, voyage_scope = "intra_eu", year = 2026, eua_price_eur = 74.5)
cat(sprintf("  [OK] gomd_calc_eu_ets()             -> EUAs=%.1f, Cost=EUR %.2f\n", ets$surrender_euas_required, ets$estimated_ets_cost_eur))

fueleu <- gomd_calc_fueleu(list(vlsfo = 6200, bio_lng = 1400, e_methanol = 450), year = 2026, voyage_scope_share = 0.75)
cat(sprintf("  [OK] gomd_calc_fueleu()             -> WtW=%.4f gCO2e/MJ (%s)\n", fueleu$attained_wtw_gco2e_mj, fueleu$status))

strive <- gomd_imo2023_trajectory(baseline_2008_aer = 6.12, current_aer = 4.85, annual_emissions_tco2e = 29400)
cat(sprintf("  [OK] gomd_imo2023_trajectory()      -> Stranded=%d, Carbon_VaR=$%.2f\n",
            strive$stranded_asset_year_strive, strive$maritime_carbon_var_10yr_usd))

ais <- gomd_ais_propulsion_emissions(mcr_kw = 11800, design_speed_kn = 14.5, observed_speed_kn = 12.2, voyage_hours = 360)
cat(sprintf("  [OK] gomd_ais_propulsion_emissions()-> Fuel=%.2f t, CO2=%.2f t\n", ais$fuel_consumption_tonnes, ais$co2_emissions_tonnes))

pos <- gomd_poseidon_alignment(attained_aer = 4.52, required_trajectory_aer = 4.68, loan_outstanding_usd = 28500000)
cat(sprintf("  [OK] gomd_poseidon_alignment()      -> delta=%.2f%% (Aligned=%s)\n", pos$alignment_delta_pct, pos$poseidon_aligned))

psc <- gomd_psc_risk(ship_age_years = 17, flag_tier = "grey", ro_performance = "medium", deficiencies_36m = 9, detentions_36m = 1)
cat(sprintf("  [OK] gomd_psc_risk()                -> NIR_pts=%d, P(detention)=%.1f%%\n", psc$nir_risk_points, psc$estimated_detention_risk_pct))

# 4. Non-Parametric Climate & Causal Econometrics (Live GOMD Series)
cat("\n[MODULE 4] Non-Parametric Climate & Causal Econometrics (8 Functions)\n")
sst_yr <- as.numeric(sst_df$value)
sst_periods <- as.numeric(substr(as.character(sst_df$period), 1, 4))
mk <- gomd_mann_kendall(sst_yr, sst_periods)
cat(sprintf("  [OK] gomd_mann_kendall()            -> tau=%.4f, Sen_slope=%.4f\n", mk$kendall_tau, mk$theil_sen_slope))

pet <- gomd_pettitt_test(sst_yr, sst_periods)
cat(sprintf("  [OK] gomd_pettitt_test()            -> break=%s, delta=+%.4f C\n", pet$changepoint_label, pet$regime_shift_delta))

evt <- gomd_gumbel_return_levels(as.numeric(wav_df$value))
cat(sprintf("  [OK] gomd_gumbel_return_levels()    -> RL_100yr=%.4f\n", evt$return_levels$return_level_100yr))

n_eg <- min(length(sst_yr), nrow(co2_df))
eg <- gomd_cointegration_eg(as.numeric(tail(co2_df$value, n_eg)), tail(sst_yr, n_eg))
cat(sprintf("  [OK] gomd_cointegration_eg()        -> ADF_t=%.4f\n", eg$engle_granger_adf_t_stat))

ipt <- gomd_innovative_polygon_trend(sst_yr)
cat(sprintf("  [OK] gomd_innovative_polygon_trend()-> ITA_slope=%.4f, tail_amp=%.3fx\n", ipt$ita_slope_per_step, ipt$tail_amplification_ratio))

n_te <- min(length(sst_yr), nrow(ph_df))
te <- gomd_transfer_entropy(tail(sst_yr, n_te), as.numeric(tail(ph_df$value, n_te)))
cat(sprintf("  [OK] gomd_transfer_entropy()        -> TE(X->Y)=%.4f bits\n", te$transfer_entropy_x_to_y_bits))

# 5. Spatial & Chokepoints
cat("\n[MODULE 5] Spatial & Chokepoint Analytics (3 Functions)\n")
hav <- gomd_haversine(41.0082, 28.9784, 36.8969, 30.7133)
cat(sprintf("  [OK] gomd_haversine()               -> %.3f km / %.2f NM\n", hav$distance_km, hav$distance_nm))

bas <- gomd_classify_basin(40.72, 28.15)
cat(sprintf("  [OK] gomd_classify_basin()          -> %s (%s)\n", bas$basin_name, bas$basin_code))

chk <- gomd_chokepoint_proximity(41.02, 29.00)
cat(sprintf("  [OK] gomd_chokepoint_proximity()    -> Nearest=%s (%.3f NM)\n",
            chk$nearest_chokepoint$name, chk$nearest_chokepoint$distance_nm))