diff --git a/.tool-versions b/.tool-versions new file mode 100644 index 00000000..2f7c62ef --- /dev/null +++ b/.tool-versions @@ -0,0 +1 @@ +elixir 1.19.5-otp-28 diff --git a/docs/analysis_output.txt b/docs/analysis_output.txt new file mode 100644 index 00000000..ba1d0b08 --- /dev/null +++ b/docs/analysis_output.txt @@ -0,0 +1,347 @@ +[info] Migrations already up +================================================================================ +PROPAGATION ANALYSIS: QSO-HRRR Correlation Study +================================================================================ +Started: 2026-04-01 13:34:58.106614Z + +Building dataset (joining QSOs × HRRR at both endpoints)... +[debug] QUERY OK db=3350.7ms decode=8.3ms queue=1.8ms idle=49.2ms +SELECT + q.id AS qso_id, + q.band::float AS band, + q.distance_km::float AS distance_km, + q.qso_timestamp, + EXTRACT(MONTH FROM q.qso_timestamp)::int AS month, + EXTRACT(HOUR FROM q.qso_timestamp)::int AS utc_hour, + + -- Endpoint 1 atmospheric + h1.surface_temp_c AS h1_temp, + h1.surface_dewpoint_c AS h1_dewpoint, + h1.surface_pressure_mb AS h1_pressure, + h1.surface_refractivity AS h1_refractivity, + h1.min_refractivity_gradient AS h1_gradient, + h1.hpbl_m AS h1_hpbl, + h1.pwat_mm AS h1_pwat, + h1.ducting_detected AS h1_ducting, + + -- Endpoint 2 atmospheric + h2.surface_temp_c AS h2_temp, + h2.surface_dewpoint_c AS h2_dewpoint, + h2.surface_pressure_mb AS h2_pressure, + h2.surface_refractivity AS h2_refractivity, + h2.min_refractivity_gradient AS h2_gradient, + h2.hpbl_m AS h2_hpbl, + h2.pwat_mm AS h2_pwat, + h2.ducting_detected AS h2_ducting, + + -- Terrain + tp.verdict AS terrain_verdict + +FROM qsos q + +INNER JOIN hrrr_profiles h1 + ON h1.lat = ROUND((q.pos1->>'lat')::numeric * 8) / 8 + AND h1.lon = ROUND((q.pos1->>'lng')::numeric * 8) / 8 + AND h1.valid_time = date_trunc('hour', q.qso_timestamp) + +INNER JOIN hrrr_profiles h2 + ON h2.lat = ROUND((q.pos2->>'lat')::numeric * 8) / 8 + AND h2.lon = ROUND((q.pos2->>'lng')::numeric * 8) / 8 + AND h2.valid_time = date_trunc('hour', q.qso_timestamp) + +LEFT JOIN terrain_profiles tp ON tp.qso_id = q.id + +WHERE q.distance_km > 0 + AND q.distance_km < 3000 + AND q.qso_timestamp >= '2016-06-30' + AND q.pos1 IS NOT NULL + AND q.pos2 IS NOT NULL + [] +Dataset: 57248 matched QSO-HRRR rows + +-------------------------------------------------------------------------------- +DATASET SUMMARY +-------------------------------------------------------------------------------- + +QSOs per band: + 10 GHz: n=52846, median=191.8 km, p25=114.5 km, p75=283.9 km + 24 GHz: n=3621, median=93.1 km, p25=44.7 km, p75=135.5 km + 47 GHz: n=680, median=65.4 km, p25=22.3 km, p75=92.4 km + 75 GHz: n=94, median=23.7 km, p25=11.4 km, p75=85.8 km + 122000 MHz: n=7, median=80.2 km, p25=22.4 km, p75=139.0 km + +Terrain verdicts: + BLOCKED: 56430 + CLEAR: 480 + FRESNEL_PARTIAL: 338 + +-------------------------------------------------------------------------------- +SPEARMAN RANK CORRELATION WITH DISTANCE (per band) +-------------------------------------------------------------------------------- + +10 GHz (n=52846): + Variable rho n_valid + ---------------------------------------------- + Pressure (mb) -0.1804 52846 + Month 0.1052 52846 + Dewpoint (°C) -0.0590 52846 + HPBL (m) 0.0446 52846 + PWAT (mm) -0.0390 52846 + Refractivity Gradient -0.0336 52347 + Temperature (°C) 0.0313 52846 + Surface Refractivity -0.0239 52347 + UTC Hour 0.0067 52846 + +24 GHz (n=3621): + Variable rho n_valid + ---------------------------------------------- + Dewpoint (°C) -0.3705 3621 + PWAT (mm) -0.3299 3621 + Surface Refractivity -0.3169 3582 + Month 0.2723 3621 + Temperature (°C) -0.1793 3621 + Pressure (mb) -0.1723 3621 + Refractivity Gradient -0.0746 3582 + UTC Hour 0.0557 3621 + HPBL (m) -0.0485 3621 + +47 GHz (n=680): + Variable rho n_valid + ---------------------------------------------- + Pressure (mb) -0.2305 680 + PWAT (mm) -0.2265 680 + Dewpoint (°C) -0.1809 680 + Refractivity Gradient -0.1391 678 + Month 0.1105 680 + Surface Refractivity -0.1088 678 + Temperature (°C) 0.0373 680 + UTC Hour -0.0238 680 + HPBL (m) 0.0041 680 + +75 GHz (n=94): + Variable rho n_valid + ---------------------------------------------- + Dewpoint (°C) -0.7033 94 + PWAT (mm) -0.6077 94 + Temperature (°C) -0.5887 94 + Pressure (mb) -0.5695 94 + Surface Refractivity -0.5257 94 + UTC Hour -0.3919 94 + HPBL (m) 0.1496 94 + Month 0.1444 94 + Refractivity Gradient -0.0821 94 + +-------------------------------------------------------------------------------- +BINNED FACTOR ANALYSIS (median distance per bin) +-------------------------------------------------------------------------------- + +10 GHz (n=52846): + + Refractivity Gradient: + Bin n median p25 p75 + ------------------------------------------------------------ + < -500 (strong) 0 - - - + -500 to -300 (moderate) 202 213.7 180.3 284.5 + -300 to -200 (enhanced) 1910 176.1 99.5 271.0 + -200 to -100 (mild) 19345 193.5 116.7 297.0 + -100 to 0 (normal) 30890 192.4 112.4 278.0 + > 0 0 - - - + + Temperature (°C): + Bin n median p25 p75 + ------------------------------------------------------------ + < 0 4 4.6 4.6 4.6 + 0-10 621 172.5 95.4 246.4 + 10-20 16764 184.4 112.1 276.4 + 20-30 33063 195.5 119.1 290.4 + > 30 2394 178.9 98.8 299.6 + + HPBL (m): + Bin n median p25 p75 + ------------------------------------------------------------ + < 200 10779 176.5 101.8 280.3 + 200-500 14670 186.6 110.1 296.4 + 500-1000 17409 199.2 125.8 281.4 + 1000-2000 9461 204.2 128.0 283.5 + > 2000 527 121.1 62.3 190.0 + + Pressure (mb): + Bin n median p25 p75 + ------------------------------------------------------------ + < 1005 47669 197.1 121.2 285.1 + 1005-1013 3447 151.3 100.7 289.2 + 1013-1020 1254 130.8 77.9 239.2 + > 1020 476 103.4 76.8 191.3 + + PWAT (mm): + Bin n median p25 p75 + ------------------------------------------------------------ + < 10 1834 160.6 96.8 265.0 + 10-20 15041 193.8 125.1 280.7 + 20-30 17788 218.8 129.9 295.6 + 30-40 13999 173.9 106.5 289.6 + > 40 4184 155.2 84.5 256.1 + + UTC Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 UTC 4171 181.8 114.6 296.0 + 03-05 UTC 950 210.0 126.3 308.2 + 06-08 UTC 192 196.3 72.6 273.9 + 09-11 UTC 1271 171.4 83.9 278.4 + 12-14 UTC 9545 184.7 108.2 297.0 + 15-17 UTC 12591 203.1 112.2 285.2 + 18-20 UTC 13812 196.2 119.7 274.0 + 21-23 UTC 10314 188.4 129.5 281.4 + + Month: + Bin n median p25 p75 + ------------------------------------------------------------ + Jan 0 - - - + Feb 0 - - - + Mar 0 - - - + Apr 1 825.0 825.0 825.0 + May 3 682.0 662.5 696.8 + Jun 1 772.0 772.0 772.0 + Jul 4 974.5 813.5 1098.0 + Aug 26813 195.8 113.4 296.4 + Sep 26024 188.4 116.2 274.0 + Oct 0 - - - + Nov 0 - - - + Dec 0 - - - + + Terrain Verdict: + Bin n median p25 p75 + ------------------------------------------------------------ + CLEAR 407 71.5 30.9 142.1 + FRESNEL_PARTIAL 284 38.6 37.9 96.1 + BLOCKED 52155 194.6 116.9 285.8 + +24 GHz (n=3621): + + Refractivity Gradient: + Bin n median p25 p75 + ------------------------------------------------------------ + < -500 (strong) 0 - - - + -500 to -300 (moderate) 15 189.0 158.4 196.7 + -300 to -200 (enhanced) 87 109.4 72.6 203.4 + -200 to -100 (mild) 1291 98.8 56.0 136.5 + -100 to 0 (normal) 2189 87.7 41.9 130.0 + > 0 0 - - - + + Temperature (°C): + Bin n median p25 p75 + ------------------------------------------------------------ + < 0 2 259.0 259.0 259.0 + 0-10 32 93.1 61.9 93.1 + 10-20 1212 104.0 60.4 139.8 + 20-30 2142 86.5 41.9 131.2 + > 30 233 47.0 24.4 129.8 + + HPBL (m): + Bin n median p25 p75 + ------------------------------------------------------------ + < 200 775 85.3 47.8 134.6 + 200-500 1164 102.8 43.7 146.1 + 500-1000 1114 98.8 50.2 135.9 + 1000-2000 500 83.7 47.5 128.0 + > 2000 68 27.9 15.9 47.7 + + Pressure (mb): + Bin n median p25 p75 + ------------------------------------------------------------ + < 1005 3289 94.3 44.7 136.0 + 1005-1013 136 110.1 99.0 169.3 + 1013-1020 109 48.7 13.0 88.0 + > 1020 87 65.4 17.1 89.6 + + PWAT (mm): + Bin n median p25 p75 + ------------------------------------------------------------ + < 10 206 126.1 72.1 139.2 + 10-20 1031 126.1 65.5 171.8 + 20-30 803 98.0 52.4 142.0 + 30-40 1211 80.8 37.7 110.1 + > 40 370 44.7 28.1 88.3 + + UTC Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 UTC 368 112.3 75.0 159.8 + 03-05 UTC 56 156.8 96.2 200.5 + 06-08 UTC 24 185.3 72.6 272.1 + 09-11 UTC 134 57.3 42.6 65.5 + 12-14 UTC 795 77.3 40.8 110.1 + 15-17 UTC 805 93.9 40.8 129.4 + 18-20 UTC 802 98.4 51.0 139.2 + 21-23 UTC 637 100.4 44.7 147.3 + + Month: + Bin n median p25 p75 + ------------------------------------------------------------ + Jan 0 - - - + Feb 2 259.0 259.0 259.0 + Mar 0 - - - + Apr 1 220.0 220.0 220.0 + May 2 242.0 224.0 260.0 + Jun 0 - - - + Jul 0 - - - + Aug 1932 76.9 39.5 126.0 + Sep 1684 104.3 61.9 142.0 + Oct 0 - - - + Nov 0 - - - + Dec 0 - - - + + Terrain Verdict: + Bin n median p25 p75 + ------------------------------------------------------------ + CLEAR 43 37.2 9.3 37.2 + FRESNEL_PARTIAL 36 68.4 43.7 93.1 + BLOCKED 3542 94.8 44.9 135.7 + +-------------------------------------------------------------------------------- +INTERACTION EFFECTS (10 GHz) +-------------------------------------------------------------------------------- + +Dataset: 52846 QSOs at 10 GHz + +1. Refractivity Gradient × Time of Day + (strong gradient = avg < -100, weak = avg >= -100) + + Time Block Gradient n median km + ------------------------------------------------------ + 00-05 UTC (night) strong 3073 204.2 + 00-05 UTC (night) weak 2044 167.8 + 06-11 UTC (morning) strong 700 185.6 + 06-11 UTC (morning) weak 763 155.6 + 12-17 UTC (afternoon) strong 9004 189.8 + 12-17 UTC (afternoon) weak 12873 197.5 + 18-23 UTC (evening) strong 8680 189.3 + 18-23 UTC (evening) weak 15210 194.6 + +2. HPBL × Season + (shallow = avg < 500m, mid = 500-999m, deep = avg >= 1000m) + + Season HPBL n median km + ---------------------------------------------------- + Winter (Dec-Feb) shallow 0 - + Winter (Dec-Feb) mid 0 - + Winter (Dec-Feb) deep 0 - + Spring (Mar-May) shallow 0 - + Spring (Mar-May) mid 3 711.6 + Spring (Mar-May) deep 1 643.0 + Summer (Jun-Aug) shallow 12909 179.9 + Summer (Jun-Aug) mid 9221 206.9 + Summer (Jun-Aug) deep 4688 231.0 + Fall (Sep-Nov) shallow 12540 184.1 + Fall (Sep-Nov) mid 8185 194.3 + Fall (Sep-Nov) deep 5299 187.6 + +3. Ducting Detection vs Distance + + Ducting YES: n=7979, median=189.0 km, p25=124.9 km, p75=296.4 km, max=986.2 km + Ducting NO : n=44867, median=192.4 km, p25=112.4 km, p75=281.4 km, max=1264.5 km + + +================================================================================ +Analysis complete: 2026-04-01 13:35:40.954224Z diff --git a/docs/analysis_output_solar.txt b/docs/analysis_output_solar.txt new file mode 100644 index 00000000..bc659ebb --- /dev/null +++ b/docs/analysis_output_solar.txt @@ -0,0 +1,378 @@ +[info] Migrations already up +================================================================================ +PROPAGATION ANALYSIS: QSO-HRRR Correlation Study +================================================================================ +Started: 2026-04-01 13:59:49.016617Z + +Building dataset (joining QSOs × HRRR at both endpoints)... +[debug] QUERY OK db=3011.3ms decode=8.3ms queue=3.2ms idle=49.0ms +SELECT + q.id AS qso_id, + q.band::float AS band, + q.distance_km::float AS distance_km, + q.qso_timestamp, + EXTRACT(MONTH FROM q.qso_timestamp)::int AS month, + EXTRACT(HOUR FROM q.qso_timestamp)::int AS utc_hour, + + -- Endpoint 1 atmospheric + h1.surface_temp_c AS h1_temp, + h1.surface_dewpoint_c AS h1_dewpoint, + h1.surface_pressure_mb AS h1_pressure, + h1.surface_refractivity AS h1_refractivity, + h1.min_refractivity_gradient AS h1_gradient, + h1.hpbl_m AS h1_hpbl, + h1.pwat_mm AS h1_pwat, + h1.ducting_detected AS h1_ducting, + + -- Endpoint 2 atmospheric + h2.surface_temp_c AS h2_temp, + h2.surface_dewpoint_c AS h2_dewpoint, + h2.surface_pressure_mb AS h2_pressure, + h2.surface_refractivity AS h2_refractivity, + h2.min_refractivity_gradient AS h2_gradient, + h2.hpbl_m AS h2_hpbl, + h2.pwat_mm AS h2_pwat, + h2.ducting_detected AS h2_ducting, + + -- Terrain + tp.verdict AS terrain_verdict, + + -- Longitude for solar time + ((q.pos1->>'lng')::float + (q.pos2->>'lng')::float) / 2.0 AS avg_longitude + +FROM qsos q + +INNER JOIN hrrr_profiles h1 + ON h1.lat = ROUND((q.pos1->>'lat')::numeric * 8) / 8 + AND h1.lon = ROUND((q.pos1->>'lng')::numeric * 8) / 8 + AND h1.valid_time = date_trunc('hour', q.qso_timestamp) + +INNER JOIN hrrr_profiles h2 + ON h2.lat = ROUND((q.pos2->>'lat')::numeric * 8) / 8 + AND h2.lon = ROUND((q.pos2->>'lng')::numeric * 8) / 8 + AND h2.valid_time = date_trunc('hour', q.qso_timestamp) + +LEFT JOIN terrain_profiles tp ON tp.qso_id = q.id + +WHERE q.distance_km > 0 + AND q.distance_km < 3000 + AND q.qso_timestamp >= '2016-06-30' + AND q.pos1 IS NOT NULL + AND q.pos2 IS NOT NULL + [] +Dataset: 57248 matched QSO-HRRR rows + +-------------------------------------------------------------------------------- +DATASET SUMMARY +-------------------------------------------------------------------------------- + +QSOs per band: + 10 GHz: n=52846, median=191.8 km, p25=114.5 km, p75=283.9 km + 24 GHz: n=3621, median=93.1 km, p25=44.7 km, p75=135.5 km + 47 GHz: n=680, median=65.4 km, p25=22.3 km, p75=92.4 km + 75 GHz: n=94, median=23.7 km, p25=11.4 km, p75=85.8 km + 122000 MHz: n=7, median=80.2 km, p25=22.4 km, p75=139.0 km + +Terrain verdicts: + BLOCKED: 56430 + CLEAR: 480 + FRESNEL_PARTIAL: 338 + +-------------------------------------------------------------------------------- +SPEARMAN RANK CORRELATION WITH DISTANCE (per band) +-------------------------------------------------------------------------------- + +10 GHz (n=52846): + Variable rho n_valid + ---------------------------------------------- + Pressure (mb) -0.1804 52846 + Month 0.1052 52846 + Dewpoint (°C) -0.0590 52846 + HPBL (m) 0.0446 52846 + PWAT (mm) -0.0390 52846 + Refractivity Gradient -0.0336 52347 + Temperature (°C) 0.0313 52846 + Surface Refractivity -0.0239 52347 + Solar Hour 0.0160 52846 + UTC Hour 0.0067 52846 + +24 GHz (n=3621): + Variable rho n_valid + ---------------------------------------------- + Dewpoint (°C) -0.3705 3621 + PWAT (mm) -0.3299 3621 + Surface Refractivity -0.3169 3582 + Month 0.2723 3621 + Solar Hour 0.1876 3621 + Temperature (°C) -0.1793 3621 + Pressure (mb) -0.1723 3621 + Refractivity Gradient -0.0746 3582 + UTC Hour 0.0557 3621 + HPBL (m) -0.0485 3621 + +47 GHz (n=680): + Variable rho n_valid + ---------------------------------------------- + Pressure (mb) -0.2305 680 + PWAT (mm) -0.2265 680 + Dewpoint (°C) -0.1809 680 + Solar Hour 0.1522 680 + Refractivity Gradient -0.1391 678 + Month 0.1105 680 + Surface Refractivity -0.1088 678 + Temperature (°C) 0.0373 680 + UTC Hour -0.0238 680 + HPBL (m) 0.0041 680 + +75 GHz (n=94): + Variable rho n_valid + ---------------------------------------------- + Dewpoint (°C) -0.7033 94 + PWAT (mm) -0.6077 94 + Temperature (°C) -0.5887 94 + Pressure (mb) -0.5695 94 + Surface Refractivity -0.5257 94 + UTC Hour -0.3919 94 + Solar Hour 0.2393 94 + HPBL (m) 0.1496 94 + Month 0.1444 94 + Refractivity Gradient -0.0821 94 + +-------------------------------------------------------------------------------- +BINNED FACTOR ANALYSIS (median distance per bin) +-------------------------------------------------------------------------------- + +10 GHz (n=52846): + + Refractivity Gradient: + Bin n median p25 p75 + ------------------------------------------------------------ + < -500 (strong) 0 - - - + -500 to -300 (moderate) 202 213.7 180.3 284.5 + -300 to -200 (enhanced) 1910 176.1 99.5 271.0 + -200 to -100 (mild) 19345 193.5 116.7 297.0 + -100 to 0 (normal) 30890 192.4 112.4 278.0 + > 0 0 - - - + + Temperature (°C): + Bin n median p25 p75 + ------------------------------------------------------------ + < 0 4 4.6 4.6 4.6 + 0-10 621 172.5 95.4 246.4 + 10-20 16764 184.4 112.1 276.4 + 20-30 33063 195.5 119.1 290.4 + > 30 2394 178.9 98.8 299.6 + + HPBL (m): + Bin n median p25 p75 + ------------------------------------------------------------ + < 200 10779 176.5 101.8 280.3 + 200-500 14670 186.6 110.1 296.4 + 500-1000 17409 199.2 125.8 281.4 + 1000-2000 9461 204.2 128.0 283.5 + > 2000 527 121.1 62.3 190.0 + + Pressure (mb): + Bin n median p25 p75 + ------------------------------------------------------------ + < 1005 47669 197.1 121.2 285.1 + 1005-1013 3447 151.3 100.7 289.2 + 1013-1020 1254 130.8 77.9 239.2 + > 1020 476 103.4 76.8 191.3 + + PWAT (mm): + Bin n median p25 p75 + ------------------------------------------------------------ + < 10 1834 160.6 96.8 265.0 + 10-20 15041 193.8 125.1 280.7 + 20-30 17788 218.8 129.9 295.6 + 30-40 13999 173.9 106.5 289.6 + > 40 4184 155.2 84.5 256.1 + + UTC Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 UTC 4171 181.8 114.6 296.0 + 03-05 UTC 950 210.0 126.3 308.2 + 06-08 UTC 192 196.3 72.6 273.9 + 09-11 UTC 1271 171.4 83.9 278.4 + 12-14 UTC 9545 184.7 108.2 297.0 + 15-17 UTC 12591 203.1 112.2 285.2 + 18-20 UTC 13812 196.2 119.7 274.0 + 21-23 UTC 10314 188.4 129.5 281.4 + + Solar Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 Solar (night) 71 218.9 75.3 274.0 + 03-05 Solar (pre-dawn) 1675 183.2 83.6 273.7 + 06-08 Solar (dawn) 10002 190.7 108.2 302.5 + 09-11 Solar (morning) 14145 191.8 110.1 276.4 + 12-14 Solar (midday) 13559 201.7 127.5 276.4 + 15-17 Solar (afternoon) 9705 181.2 126.0 283.0 + 18-20 Solar (evening) 3077 188.2 116.0 313.1 + 21-23 Solar (night) 612 187.9 121.6 284.6 + + Month: + Bin n median p25 p75 + ------------------------------------------------------------ + Jan 0 - - - + Feb 0 - - - + Mar 0 - - - + Apr 1 825.0 825.0 825.0 + May 3 682.0 662.5 696.8 + Jun 1 772.0 772.0 772.0 + Jul 4 974.5 813.5 1098.0 + Aug 26813 195.8 113.4 296.4 + Sep 26024 188.4 116.2 274.0 + Oct 0 - - - + Nov 0 - - - + Dec 0 - - - + + Terrain Verdict: + Bin n median p25 p75 + ------------------------------------------------------------ + CLEAR 407 71.5 30.9 142.1 + FRESNEL_PARTIAL 284 38.6 37.9 96.1 + BLOCKED 52155 194.6 116.9 285.8 + +24 GHz (n=3621): + + Refractivity Gradient: + Bin n median p25 p75 + ------------------------------------------------------------ + < -500 (strong) 0 - - - + -500 to -300 (moderate) 15 189.0 158.4 196.7 + -300 to -200 (enhanced) 87 109.4 72.6 203.4 + -200 to -100 (mild) 1291 98.8 56.0 136.5 + -100 to 0 (normal) 2189 87.7 41.9 130.0 + > 0 0 - - - + + Temperature (°C): + Bin n median p25 p75 + ------------------------------------------------------------ + < 0 2 259.0 259.0 259.0 + 0-10 32 93.1 61.9 93.1 + 10-20 1212 104.0 60.4 139.8 + 20-30 2142 86.5 41.9 131.2 + > 30 233 47.0 24.4 129.8 + + HPBL (m): + Bin n median p25 p75 + ------------------------------------------------------------ + < 200 775 85.3 47.8 134.6 + 200-500 1164 102.8 43.7 146.1 + 500-1000 1114 98.8 50.2 135.9 + 1000-2000 500 83.7 47.5 128.0 + > 2000 68 27.9 15.9 47.7 + + Pressure (mb): + Bin n median p25 p75 + ------------------------------------------------------------ + < 1005 3289 94.3 44.7 136.0 + 1005-1013 136 110.1 99.0 169.3 + 1013-1020 109 48.7 13.0 88.0 + > 1020 87 65.4 17.1 89.6 + + PWAT (mm): + Bin n median p25 p75 + ------------------------------------------------------------ + < 10 206 126.1 72.1 139.2 + 10-20 1031 126.1 65.5 171.8 + 20-30 803 98.0 52.4 142.0 + 30-40 1211 80.8 37.7 110.1 + > 40 370 44.7 28.1 88.3 + + UTC Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 UTC 368 112.3 75.0 159.8 + 03-05 UTC 56 156.8 96.2 200.5 + 06-08 UTC 24 185.3 72.6 272.1 + 09-11 UTC 134 57.3 42.6 65.5 + 12-14 UTC 795 77.3 40.8 110.1 + 15-17 UTC 805 93.9 40.8 129.4 + 18-20 UTC 802 98.4 51.0 139.2 + 21-23 UTC 637 100.4 44.7 147.3 + + Solar Hour: + Bin n median p25 p75 + ------------------------------------------------------------ + 00-02 Solar (night) 2 114.1 72.9 155.4 + 03-05 Solar (pre-dawn) 159 57.3 41.9 65.9 + 06-08 Solar (dawn) 705 72.5 37.1 108.2 + 09-11 Solar (morning) 943 96.5 42.3 131.5 + 12-14 Solar (midday) 810 99.9 61.8 131.9 + 15-17 Solar (afternoon) 659 107.1 44.9 162.1 + 18-20 Solar (evening) 302 104.8 53.9 159.6 + 21-23 Solar (night) 41 140.4 76.9 204.2 + + Month: + Bin n median p25 p75 + ------------------------------------------------------------ + Jan 0 - - - + Feb 2 259.0 259.0 259.0 + Mar 0 - - - + Apr 1 220.0 220.0 220.0 + May 2 242.0 224.0 260.0 + Jun 0 - - - + Jul 0 - - - + Aug 1932 76.9 39.5 126.0 + Sep 1684 104.3 61.9 142.0 + Oct 0 - - - + Nov 0 - - - + Dec 0 - - - + + Terrain Verdict: + Bin n median p25 p75 + ------------------------------------------------------------ + CLEAR 43 37.2 9.3 37.2 + FRESNEL_PARTIAL 36 68.4 43.7 93.1 + BLOCKED 3542 94.8 44.9 135.7 + +-------------------------------------------------------------------------------- +INTERACTION EFFECTS (10 GHz) +-------------------------------------------------------------------------------- + +Dataset: 52846 QSOs at 10 GHz + +1. Refractivity Gradient × Time of Day + (strong gradient = avg < -100, weak = avg >= -100) + + Time Block Gradient n median km + ------------------------------------------------------ + 00-05 UTC (night) strong 3073 204.2 + 00-05 UTC (night) weak 2044 167.8 + 06-11 UTC (morning) strong 700 185.6 + 06-11 UTC (morning) weak 763 155.6 + 12-17 UTC (afternoon) strong 9004 189.8 + 12-17 UTC (afternoon) weak 12873 197.5 + 18-23 UTC (evening) strong 8680 189.3 + 18-23 UTC (evening) weak 15210 194.6 + +2. HPBL × Season + (shallow = avg < 500m, mid = 500-999m, deep = avg >= 1000m) + + Season HPBL n median km + ---------------------------------------------------- + Winter (Dec-Feb) shallow 0 - + Winter (Dec-Feb) mid 0 - + Winter (Dec-Feb) deep 0 - + Spring (Mar-May) shallow 0 - + Spring (Mar-May) mid 3 711.6 + Spring (Mar-May) deep 1 643.0 + Summer (Jun-Aug) shallow 12909 179.9 + Summer (Jun-Aug) mid 9221 206.9 + Summer (Jun-Aug) deep 4688 231.0 + Fall (Sep-Nov) shallow 12540 184.1 + Fall (Sep-Nov) mid 8185 194.3 + Fall (Sep-Nov) deep 5299 187.6 + +3. Ducting Detection vs Distance + + Ducting YES: n=7979, median=189.0 km, p25=124.9 km, p75=296.4 km, max=986.2 km + Ducting NO : n=44867, median=192.4 km, p25=112.4 km, p75=281.4 km, max=1264.5 km + + +================================================================================ +Analysis complete: 2026-04-01 14:00:36.820405Z diff --git a/docs/plans/2026-04-01-algorithm-refinement-ml-training-design.md b/docs/plans/2026-04-01-algorithm-refinement-ml-training-design.md new file mode 100644 index 00000000..9edc8627 --- /dev/null +++ b/docs/plans/2026-04-01-algorithm-refinement-ml-training-design.md @@ -0,0 +1,87 @@ +# Algorithm Refinement & ML Training Design + +## Goal + +Refine the propagation scoring algorithm using 57K+ historical QSOs matched to HRRR atmospheric conditions, then train the Axon ML model on the refined algorithm's output for fast inference from HRRR forecast data. + +## Data Assets + +- 57,488 tropospheric QSOs (1992-2024) with distance as ground truth +- 54M HRRR grid profiles (2016-2026) at 0.125° resolution +- 326K surface observations +- 58K terrain profiles (97% BLOCKED — confirms atmospheric propagation) +- 3M propagation_scores (current algorithm output) +- Band distribution: 53K at 10G, 3.8K at 24G, 800 at 47G, 141 at 75G, sparse above + +## Design Decisions + +- **Target variable**: Hybrid — refine algorithm from QSO data, then train ML on refined scores +- **QSO-HRRR matching**: Endpoint average — HRRR conditions at both QSO endpoints, averaged +- **Algorithm scope**: Full re-examination — weights, thresholds, add/remove factors +- **ML deployment**: ML primary scorer with algorithm fallback for missing data + +## Implementation Plan + +### Stage 1: Analysis Mix Task (`mix propagation.analyze`) + +Build QSO↔HRRR analysis dataset: +- For each QSO, find nearest HRRR grid point to pos1 and pos2 (round to 0.125°) +- Find HRRR profile at nearest hour to qso_timestamp +- Average atmospheric conditions across both endpoints +- Output: CSV or in-memory dataset with columns: + - qso_id, band, distance_km, month, utc_hour + - avg_surface_temp_c, avg_surface_dewpoint_c, avg_surface_pressure_mb + - avg_min_refractivity_gradient, avg_hpbl_m, avg_pwat_mm + - avg_surface_refractivity, ducting_at_either_endpoint + - terrain_verdict (BLOCKED/CLEAR/FRESNEL_PARTIAL) + +Run correlation analysis: +- Spearman rank correlation of each variable with distance_km, stratified by band +- Bin each variable into deciles, compute median distance per bin (natural breakpoint curves) +- Check for interaction effects (e.g., gradient × time_of_day, hpbl × month) +- Look for currently-missing predictive factors +- Compare factor distributions for short (<100km) vs long (>300km) contacts + +Output: printed report with findings and recommended changes. + +### Stage 2: Update algo.md + +Document all findings with data backing: +- Which factors gained/lost predictive power +- New threshold curves derived from data +- Revised weights with justification +- Any new factors added or old factors removed +- Before/after comparison of scoring behavior + +### Stage 3: Update scorer.ex + band_config.ex + +Implement revised algorithm: +- Updated scoring functions with new thresholds +- Updated weights in band_config +- New factors if analysis warrants +- All changes match algo.md documentation +- Run `mix test` to verify no regressions (update tests for new behavior) + +### Stage 4: Recompute Propagation Scores + +Run refined algorithm across HRRR grid: +- Either `mix propagation_grid` or a targeted recompute task +- Generates fresh propagation_scores with refined algorithm +- These become ML training data + +### Stage 5: Training Mix Task (`mix propagation.train`) + +Build training pipeline: +- Query propagation_scores joined to hrrr_profiles by lat/lon/valid_time +- Features: raw atmospheric inputs + cyclical temporal + log frequency +- Target: refined composite score normalized to 0-1 +- 80/10/10 train/val/test split, stratified by band and score quartile +- MSE loss, Adam optimizer, early stopping on validation loss +- Save trained weights to priv/models/propagation_v1.nx +- Print evaluation metrics (RMSE, R², per-band performance) + +### Stage 6: Integrate ML into Grid Worker + +- model.ex: `predict_score/2` for inference +- PropagationGridWorker: try ML prediction first, fall back to Scorer.composite_score/2 +- Graceful degradation when model file missing or inputs incomplete diff --git a/docs/plans/2026-04-01-algorithm-refinement-ml-training.md b/docs/plans/2026-04-01-algorithm-refinement-ml-training.md new file mode 100644 index 00000000..1a2918b8 --- /dev/null +++ b/docs/plans/2026-04-01-algorithm-refinement-ml-training.md @@ -0,0 +1,991 @@ +# Algorithm Refinement & ML Training Implementation Plan + +> **For Claude:** REQUIRED SUB-SKILL: Use superpowers:executing-plans to implement this plan task-by-task. + +**Goal:** Refine the propagation scoring algorithm using 58K QSOs matched to HRRR conditions, update algo.md and scorer code, then train the Axon ML model for fast HRRR-based prediction. + +**Architecture:** Analysis mix task extracts QSO↔HRRR correlations → findings drive scorer/algo.md updates → recomputed scores become ML training data → trained model integrates into grid worker with algorithm fallback. + +**Tech Stack:** Elixir/Ecto (analysis queries), Nx/Axon/EXLA (ML training), PostgreSQL (data source) + +--- + +### Task 1: Analysis Mix Task — QSO-HRRR Dataset Builder + +**Files:** +- Create: `lib/mix/tasks/propagation_analyze.ex` +- Test: `test/mix/tasks/propagation_analyze_test.exs` + +**Step 1: Write the failing test** + +```elixir +# test/mix/tasks/propagation_analyze_test.exs +defmodule Mix.Tasks.PropagationAnalyzeTest do + use Microwaveprop.DataCase + + alias Mix.Tasks.PropagationAnalyze + + describe "build_dataset/0" do + test "returns list of QSO-condition maps with required keys" do + # Insert a QSO with known position and timestamp + qso = insert_qso( + band: 10_000, + distance_km: 200.0, + qso_timestamp: ~N[2023-06-15 12:00:00], + pos1: %{"lat" => 33.0, "lng" => -97.0}, + pos2: %{"lat" => 34.0, "lng" => -96.0} + ) + + # Insert HRRR profiles at the grid points nearest to both endpoints + insert_hrrr_profile( + lat: 33.0, lon: -97.0, valid_time: ~N[2023-06-15 12:00:00], + surface_temp_c: 30.0, surface_dewpoint_c: 20.0, + surface_pressure_mb: 1010.0, min_refractivity_gradient: -100.0, + hpbl_m: 500.0, pwat_mm: 30.0, surface_refractivity: 330.0, + ducting_detected: false + ) + + insert_hrrr_profile( + lat: 34.0, lon: -96.0, valid_time: ~N[2023-06-15 12:00:00], + surface_temp_c: 28.0, surface_dewpoint_c: 18.0, + surface_pressure_mb: 1012.0, min_refractivity_gradient: -80.0, + hpbl_m: 600.0, pwat_mm: 25.0, surface_refractivity: 320.0, + ducting_detected: false + ) + + [row] = PropagationAnalyze.build_dataset() + + assert row.qso_id == qso.id + assert row.band == 10_000 + assert row.distance_km == 200.0 + assert row.month == 6 + assert row.utc_hour == 12 + assert_in_delta row.avg_surface_temp_c, 29.0, 0.01 + assert_in_delta row.avg_min_refractivity_gradient, -90.0, 0.01 + assert_in_delta row.avg_hpbl_m, 550.0, 0.01 + end + + test "skips QSOs without HRRR match at either endpoint" do + insert_qso( + band: 10_000, distance_km: 200.0, + qso_timestamp: ~N[2023-06-15 12:00:00], + pos1: %{"lat" => 33.0, "lng" => -97.0}, + pos2: %{"lat" => 34.0, "lng" => -96.0} + ) + + # Only insert one endpoint's HRRR profile + insert_hrrr_profile( + lat: 33.0, lon: -97.0, valid_time: ~N[2023-06-15 12:00:00], + surface_temp_c: 30.0, surface_dewpoint_c: 20.0, + surface_pressure_mb: 1010.0, min_refractivity_gradient: -100.0, + hpbl_m: 500.0, pwat_mm: 30.0 + ) + + assert PropagationAnalyze.build_dataset() == [] + end + end +end +``` + +**Step 2: Run test to verify it fails** + +Run: `mix test test/mix/tasks/propagation_analyze_test.exs --trace` +Expected: FAIL — module not defined + +**Step 3: Write test helpers** + +Add to `test/support/fixtures.ex` (or create if missing) — helper functions `insert_qso/1` and `insert_hrrr_profile/1` that insert records directly via Repo. + +**Step 4: Write minimal implementation** + +```elixir +# lib/mix/tasks/propagation_analyze.ex +defmodule Mix.Tasks.PropagationAnalyze do + @shortdoc "Analyze QSO-HRRR correlations for algorithm refinement" + @moduledoc """ + Builds a dataset matching QSOs to HRRR atmospheric conditions at both + endpoints, then runs correlation analysis against QSO distance. + + ## Usage + + mix propagation.analyze + """ + + use Mix.Task + + import Ecto.Query + + alias Microwaveprop.Repo + + @impl Mix.Task + def run(_args) do + Application.put_env(:microwaveprop, Oban, queues: false) + Mix.Task.run("app.start") + + IO.puts("Building QSO-HRRR analysis dataset...") + dataset = build_dataset() + IO.puts("Matched #{length(dataset)} QSOs to HRRR conditions") + + IO.puts("\n=== CORRELATION ANALYSIS ===\n") + run_correlation_analysis(dataset) + + IO.puts("\n=== FACTOR ANALYSIS ===\n") + run_factor_analysis(dataset) + + IO.puts("\n=== INTERACTION EFFECTS ===\n") + run_interaction_analysis(dataset) + + IO.puts("\nDone. Use findings to update algo.md and scorer.ex") + end + + @doc """ + Builds the QSO-HRRR matched dataset. + + For each QSO with valid positions and a timestamp within HRRR range (2016+), + finds the nearest HRRR grid point to each endpoint (rounded to 0.125°) at + the nearest hour. Averages atmospheric conditions across both endpoints. + + Returns a list of maps with: + - qso_id, band, distance_km, month, utc_hour + - avg_surface_temp_c, avg_surface_dewpoint_c, avg_surface_pressure_mb + - avg_min_refractivity_gradient, avg_hpbl_m, avg_pwat_mm + - avg_surface_refractivity, ducting_either (boolean) + - terrain_verdict + """ + def build_dataset do + # Query uses a CTE to round QSO positions to HRRR grid and truncate timestamps, + # then joins to hrrr_profiles at both endpoints, averaging atmospheric values. + query = """ + WITH qso_grid AS ( + SELECT + q.id AS qso_id, + q.band::integer AS band, + q.distance_km, + EXTRACT(MONTH FROM q.qso_timestamp)::integer AS month, + EXTRACT(HOUR FROM q.qso_timestamp)::integer AS utc_hour, + ROUND(CAST((q.pos1->>'lat')::float * 8 AS numeric)) / 8 AS grid_lat1, + ROUND(CAST((q.pos1->>'lng')::float * 8 AS numeric)) / 8 AS grid_lng1, + ROUND(CAST((q.pos2->>'lat')::float * 8 AS numeric)) / 8 AS grid_lat2, + ROUND(CAST((q.pos2->>'lng')::float * 8 AS numeric)) / 8 AS grid_lng2, + date_trunc('hour', q.qso_timestamp) AS match_time + FROM qsos q + WHERE q.distance_km < 3000 + AND q.distance_km > 0 + AND q.qso_timestamp >= '2016-06-30' + AND q.pos1->>'lng' IS NOT NULL + AND q.pos2->>'lng' IS NOT NULL + ) + SELECT + qg.qso_id, qg.band, qg.distance_km, qg.month, qg.utc_hour, + (COALESCE(h1.surface_temp_c, h2.surface_temp_c) + + COALESCE(h2.surface_temp_c, h1.surface_temp_c)) / 2.0 AS avg_surface_temp_c, + (COALESCE(h1.surface_dewpoint_c, h2.surface_dewpoint_c) + + COALESCE(h2.surface_dewpoint_c, h1.surface_dewpoint_c)) / 2.0 AS avg_surface_dewpoint_c, + (COALESCE(h1.surface_pressure_mb, h2.surface_pressure_mb) + + COALESCE(h2.surface_pressure_mb, h1.surface_pressure_mb)) / 2.0 AS avg_surface_pressure_mb, + (COALESCE(h1.min_refractivity_gradient, h2.min_refractivity_gradient) + + COALESCE(h2.min_refractivity_gradient, h1.min_refractivity_gradient)) / 2.0 AS avg_min_refractivity_gradient, + (COALESCE(h1.hpbl_m, h2.hpbl_m) + + COALESCE(h2.hpbl_m, h1.hpbl_m)) / 2.0 AS avg_hpbl_m, + (COALESCE(h1.pwat_mm, h2.pwat_mm) + + COALESCE(h2.pwat_mm, h1.pwat_mm)) / 2.0 AS avg_pwat_mm, + (COALESCE(h1.surface_refractivity, h2.surface_refractivity) + + COALESCE(h2.surface_refractivity, h1.surface_refractivity)) / 2.0 AS avg_surface_refractivity, + COALESCE(h1.ducting_detected, false) OR COALESCE(h2.ducting_detected, false) AS ducting_either, + tp.verdict AS terrain_verdict + FROM qso_grid qg + JOIN hrrr_profiles h1 + ON h1.lat = qg.grid_lat1 AND h1.lon = qg.grid_lng1 AND h1.valid_time = qg.match_time + JOIN hrrr_profiles h2 + ON h2.lat = qg.grid_lat2 AND h2.lon = qg.grid_lng2 AND h2.valid_time = qg.match_time + LEFT JOIN terrain_profiles tp ON tp.qso_id = qg.qso_id + ORDER BY qg.band, qg.distance_km DESC + """ + + Repo.query!(query, [], timeout: 300_000).rows + |> Enum.map(fn row -> + [qso_id, band, distance_km, month, utc_hour, + avg_temp, avg_dewpoint, avg_pressure, avg_gradient, + avg_hpbl, avg_pwat, avg_refractivity, ducting_either, + terrain_verdict] = row + + %{ + qso_id: qso_id, + band: band, + distance_km: distance_km, + month: month, + utc_hour: utc_hour, + avg_surface_temp_c: avg_temp, + avg_surface_dewpoint_c: avg_dewpoint, + avg_surface_pressure_mb: avg_pressure, + avg_min_refractivity_gradient: avg_gradient, + avg_hpbl_m: avg_hpbl, + avg_pwat_mm: avg_pwat, + avg_surface_refractivity: avg_refractivity, + ducting_either: ducting_either, + terrain_verdict: terrain_verdict + } + end) + end + + # --- Correlation analysis --- + + defp run_correlation_analysis(dataset) do + bands = [10_000, 24_000, 47_000, 75_000] + + for band <- bands do + band_data = Enum.filter(dataset, &(&1.band == band)) + next_if_empty(band_data, band, fn -> + IO.puts("--- Band: #{div(band, 1000)} GHz (n=#{length(band_data)}) ---") + distances = Enum.map(band_data, & &1.distance_km) + + variables = [ + {:avg_surface_temp_c, "Surface Temp (C)"}, + {:avg_surface_dewpoint_c, "Dewpoint (C)"}, + {:avg_surface_pressure_mb, "Pressure (mb)"}, + {:avg_min_refractivity_gradient, "Refr. Gradient (N/km)"}, + {:avg_hpbl_m, "BL Depth (m)"}, + {:avg_pwat_mm, "PWAT (mm)"}, + {:avg_surface_refractivity, "Surface Refractivity"}, + {:month, "Month"}, + {:utc_hour, "UTC Hour"} + ] + + for {key, label} <- variables do + values = Enum.map(band_data, &Map.get(&1, key)) + corr = spearman_correlation(values, distances) + bar = correlation_bar(corr) + IO.puts(" #{String.pad_trailing(label, 25)} r=#{format_corr(corr)} #{bar}") + end + + IO.puts("") + end) + end + end + + # --- Factor analysis (binned median distance) --- + + defp run_factor_analysis(dataset) do + band_10g = Enum.filter(dataset, &(&1.band == 10_000)) + band_24g = Enum.filter(dataset, &(&1.band == 24_000)) + + for {band_data, label} <- [{band_10g, "10 GHz"}, {band_24g, "24 GHz"}] do + next_if_empty(band_data, label, fn -> + IO.puts("--- #{label} Factor Breakpoints ---\n") + + # Refractivity gradient bins + IO.puts("Refractivity Gradient → Median Distance:") + print_binned(band_data, :avg_min_refractivity_gradient, :distance_km, + [{-999, -200}, {-200, -150}, {-150, -100}, {-100, -75}, {-75, -55}, {-55, -40}, {-40, 0}]) + + # Temperature bins + IO.puts("\nSurface Temperature (C) → Median Distance:") + print_binned(band_data, :avg_surface_temp_c, :distance_km, + [{-20, 0}, {0, 10}, {10, 20}, {20, 25}, {25, 30}, {30, 35}, {35, 50}]) + + # HPBL bins + IO.puts("\nBoundary Layer Depth (m) → Median Distance:") + print_binned(band_data, :avg_hpbl_m, :distance_km, + [{0, 200}, {200, 500}, {500, 1000}, {1000, 1500}, {1500, 2000}, {2000, 5000}]) + + # Pressure bins + IO.puts("\nPressure (mb) → Median Distance:") + print_binned(band_data, :avg_surface_pressure_mb, :distance_km, + [{980, 1000}, {1000, 1010}, {1010, 1015}, {1015, 1020}, {1020, 1025}, {1025, 1040}]) + + # PWAT bins + IO.puts("\nPWAT (mm) → Median Distance:") + print_binned(band_data, :avg_pwat_mm, :distance_km, + [{0, 10}, {10, 20}, {20, 30}, {30, 40}, {40, 50}, {50, 80}]) + + # Time of day + IO.puts("\nUTC Hour → Median Distance:") + for hour_range <- [{0, 3}, {3, 6}, {6, 9}, {9, 12}, {12, 15}, {15, 18}, {18, 21}, {21, 24}] do + {lo, hi} = hour_range + in_bin = Enum.filter(band_data, fn r -> + r.utc_hour >= lo and r.utc_hour < hi + end) + print_bin_stats(" #{lo}-#{hi} UTC", in_bin, :distance_km) + end + + # Month + IO.puts("\nMonth → Median Distance:") + for m <- 1..12 do + in_bin = Enum.filter(band_data, &(&1.month == m)) + month_name = Enum.at(~w(Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec), m - 1) + print_bin_stats(" #{month_name}", in_bin, :distance_km) + end + + # Terrain verdict + IO.puts("\nTerrain Verdict → Median Distance:") + for verdict <- ["BLOCKED", "CLEAR", "FRESNEL_PARTIAL", nil] do + in_bin = Enum.filter(band_data, &(&1.terrain_verdict == verdict)) + print_bin_stats(" #{verdict || "UNKNOWN"}", in_bin, :distance_km) + end + + IO.puts("") + end) + end + end + + # --- Interaction effects --- + + defp run_interaction_analysis(dataset) do + band_10g = Enum.filter(dataset, &(&1.band == 10_000)) + + next_if_empty(band_10g, "10 GHz interactions", fn -> + IO.puts("--- 10 GHz: Gradient × Time of Day ---") + for hour_range <- [{0, 6}, {6, 12}, {12, 18}, {18, 24}] do + {lo, hi} = hour_range + time_data = Enum.filter(band_10g, fn r -> r.utc_hour >= lo and r.utc_hour < hi end) + + strong = Enum.filter(time_data, fn r -> + r.avg_min_refractivity_gradient != nil and r.avg_min_refractivity_gradient < -100 + end) + weak = Enum.filter(time_data, fn r -> + r.avg_min_refractivity_gradient != nil and r.avg_min_refractivity_gradient >= -100 + end) + + strong_med = median(Enum.map(strong, & &1.distance_km)) + weak_med = median(Enum.map(weak, & &1.distance_km)) + + IO.puts(" #{lo}-#{hi} UTC: strong gradient (<-100) median=#{format_km(strong_med)} (n=#{length(strong)}), " <> + "weak (>=-100) median=#{format_km(weak_med)} (n=#{length(weak)})") + end + + IO.puts("\n--- 10 GHz: HPBL × Month ---") + for month_group <- [{[6, 7, 8], "Summer"}, {[12, 1, 2], "Winter"}, {[3, 4, 5], "Spring"}, {[9, 10, 11], "Fall"}] do + {months, label} = month_group + seasonal = Enum.filter(band_10g, fn r -> r.month in months end) + + shallow = Enum.filter(seasonal, fn r -> r.avg_hpbl_m != nil and r.avg_hpbl_m < 500 end) + deep = Enum.filter(seasonal, fn r -> r.avg_hpbl_m != nil and r.avg_hpbl_m >= 1000 end) + + shallow_med = median(Enum.map(shallow, & &1.distance_km)) + deep_med = median(Enum.map(deep, & &1.distance_km)) + + IO.puts(" #{label}: shallow BL (<500m) median=#{format_km(shallow_med)} (n=#{length(shallow)}), " <> + "deep BL (>=1000m) median=#{format_km(deep_med)} (n=#{length(deep)})") + end + + IO.puts("\n--- 10 GHz: Ducting Detection vs Distance ---") + ducting = Enum.filter(band_10g, & &1.ducting_either) + no_ducting = Enum.filter(band_10g, &(!&1.ducting_either)) + IO.puts(" Ducting detected: median=#{format_km(median(Enum.map(ducting, & &1.distance_km)))} (n=#{length(ducting)})") + IO.puts(" No ducting: median=#{format_km(median(Enum.map(no_ducting, & &1.distance_km)))} (n=#{length(no_ducting)})") + end) + end + + # --- Statistical helpers --- + + defp spearman_correlation(xs, ys) do + xs_clean = Enum.zip(xs, ys) |> Enum.reject(fn {x, _} -> is_nil(x) end) + if length(xs_clean) < 10, do: nil, else: do_spearman(xs_clean) + end + + defp do_spearman(pairs) do + n = length(pairs) + {xs, ys} = Enum.unzip(pairs) + rx = rank(xs) + ry = rank(ys) + d_sq = Enum.zip(rx, ry) |> Enum.map(fn {a, b} -> (a - b) * (a - b) end) |> Enum.sum() + 1 - 6 * d_sq / (n * (n * n - 1)) + end + + defp rank(values) do + indexed = values |> Enum.with_index() |> Enum.sort_by(fn {v, _i} -> v end) + ranked = indexed |> Enum.with_index(1) |> Enum.map(fn {{_v, orig_i}, rank} -> {orig_i, rank} end) + ranked |> Enum.sort_by(fn {orig_i, _rank} -> orig_i end) |> Enum.map(fn {_i, r} -> r * 1.0 end) + end + + defp median([]), do: nil + defp median(list) do + sorted = Enum.sort(list) + n = length(sorted) + mid = div(n, 2) + if rem(n, 2) == 0 do + (Enum.at(sorted, mid - 1) + Enum.at(sorted, mid)) / 2 + else + Enum.at(sorted, mid) + end + end + + defp percentile([], _p), do: nil + defp percentile(list, p) do + sorted = Enum.sort(list) + k = (length(sorted) - 1) * p / 100 + f = trunc(k) + c = Float.ceil(k) |> trunc() + if f == c, do: Enum.at(sorted, f), else: Enum.at(sorted, f) + (k - f) * (Enum.at(sorted, c) - Enum.at(sorted, f)) + end + + defp print_binned(data, group_key, value_key, bins) do + for {lo, hi} <- bins do + in_bin = Enum.filter(data, fn row -> + v = Map.get(row, group_key) + v != nil and v >= lo and v < hi + end) + label = " [#{lo}, #{hi})" + print_bin_stats(String.pad_trailing(label, 18), in_bin, value_key) + end + end + + defp print_bin_stats(label, data, value_key) do + n = length(data) + if n > 0 do + values = Enum.map(data, &Map.get(&1, value_key)) + med = median(values) + p25 = percentile(values, 25) + p75 = percentile(values, 75) + IO.puts("#{label} n=#{String.pad_leading("#{n}", 6)} median=#{format_km(med)} p25=#{format_km(p25)} p75=#{format_km(p75)}") + else + IO.puts("#{label} n= 0") + end + end + + defp next_if_empty([], label, _fun), do: IO.puts("Skipping #{label}: no data") + defp next_if_empty(_data, _label, fun), do: fun.() + + defp format_corr(nil), do: " N/A" + defp format_corr(corr), do: :io_lib.format("~7.4f", [corr]) |> to_string() + + defp format_km(nil), do: " N/A" + defp format_km(km), do: :io_lib.format("~7.1f km", [km]) |> to_string() + + defp correlation_bar(nil), do: "" + defp correlation_bar(corr) do + len = round(abs(corr) * 20) + bar = String.duplicate(if(corr >= 0, do: "+", else: "-"), len) + "|#{bar}|" + end +end +``` + +**Step 5: Run test to verify it passes** + +Run: `mix test test/mix/tasks/propagation_analyze_test.exs --trace` +Expected: PASS + +**Step 6: Run the analysis against real data** + +Run: `mix propagation.analyze 2>&1 | tee docs/analysis_output.txt` + +This will take several minutes (joining 58K QSOs against 54M HRRR profiles). Save the full output for reference during algo.md updates. + +**Step 7: Commit** + +```bash +git add lib/mix/tasks/propagation_analyze.ex test/mix/tasks/propagation_analyze_test.exs +git commit -m "Add mix propagation.analyze task for QSO-HRRR correlation analysis" +``` + +--- + +### Task 2: Interpret Analysis & Update algo.md + +**Files:** +- Modify: `algo.md` + +**Step 1: Read the analysis output** + +Read `docs/analysis_output.txt` and identify: +- Which factors have strong vs weak correlation with distance +- Where current scoring thresholds diverge from data-suggested breakpoints +- Any interaction effects that the current independent-factor model misses +- Whether any current factors should be removed or new ones added + +**Step 2: Update algo.md Part 2 (Key Empirical Findings)** + +Add a new section "Data-Driven Refinements (April 2026)" documenting: +- Per-factor correlation strengths by band +- Revised threshold curves with data justification +- Weight adjustments with before/after comparison +- Any added or removed factors + +**Step 3: Update algo.md Part 4 (Scoring Functions)** + +Revise each factor's documentation to match new thresholds and behavior. + +**Step 4: Update algo.md weights table** + +Show old vs new weights with justification. + +**Step 5: Commit** + +```bash +git add algo.md +git commit -m "Update algo.md with data-driven findings from 58K QSO analysis" +``` + +--- + +### Task 3: Update Scorer & BandConfig + +**Files:** +- Modify: `lib/microwaveprop/propagation/scorer.ex` +- Modify: `lib/microwaveprop/propagation/band_config.ex` +- Modify: `test/microwaveprop/propagation/scorer_test.exs` +- Modify: `test/microwaveprop/propagation/band_config_test.exs` + +**Step 1: Update band_config.ex weights** + +Change `@weights` map to match the data-derived optimal weights from the analysis. + +**Step 2: Update band_config.ex thresholds** + +Update `@refractivity_thresholds`, `@humidity_beneficial_thresholds`, seasonal tables, and any other thresholds that the analysis showed should change. + +**Step 3: Update scorer.ex scoring functions** + +Revise scoring function breakpoints to match the data-derived natural breakpoints. If the analysis warrants adding or removing factors, modify `composite_score/2` and the weights map accordingly. + +**Step 4: Update tests to match new behavior** + +Update `scorer_test.exs` and `band_config_test.exs` assertions to match the new thresholds and weights. Keep the same test structure — just adjust expected values. + +**Step 5: Run tests** + +Run: `mix test test/microwaveprop/propagation/scorer_test.exs test/microwaveprop/propagation/band_config_test.exs --trace` +Expected: PASS + +**Step 6: Run full test suite** + +Run: `mix test` +Expected: All pass + +**Step 7: Commit** + +```bash +git add lib/microwaveprop/propagation/scorer.ex lib/microwaveprop/propagation/band_config.ex +git add test/microwaveprop/propagation/scorer_test.exs test/microwaveprop/propagation/band_config_test.exs +git commit -m "Update scorer thresholds and weights from data analysis" +``` + +--- + +### Task 4: Training Mix Task + +**Files:** +- Create: `lib/mix/tasks/propagation_train.ex` +- Modify: `lib/microwaveprop/propagation/model.ex` (add train/2, evaluate/2) +- Test: `test/mix/tasks/propagation_train_test.exs` +- Modify: `test/microwaveprop/propagation/model_test.exs` + +**Step 1: Write failing test for Model.train/2** + +```elixir +# In test/microwaveprop/propagation/model_test.exs, add: + +describe "train/2" do + test "trains model and returns params with reduced loss" do + # Generate synthetic training data: 100 samples + features = Nx.random_uniform({100, 13}, type: :f32) + targets = Nx.random_uniform({100, 1}, type: :f32) + + {params, metrics} = Model.train(features, targets, epochs: 5, batch_size: 32) + + # Params should have the right structure + assert Map.has_key?(params, "hidden_1") + assert Map.has_key?(params, "hidden_2") + assert Map.has_key?(params, "output") + + # Should return training metrics + assert is_float(metrics.final_loss) + assert metrics.final_loss >= 0 + end +end + +describe "evaluate/2" do + test "returns RMSE and R-squared on test data" do + params = Model.init() + features = Nx.random_uniform({50, 13}, type: :f32) + targets = Nx.random_uniform({50, 1}, type: :f32) + + metrics = Model.evaluate(params, features, targets) + + assert is_float(metrics.rmse) + assert is_float(metrics.r_squared) + assert metrics.rmse >= 0 + end +end +``` + +**Step 2: Run test to verify it fails** + +Run: `mix test test/microwaveprop/propagation/model_test.exs --trace` +Expected: FAIL — train/2 and evaluate/2 not defined + +**Step 3: Implement Model.train/2 and Model.evaluate/2** + +Add to `lib/microwaveprop/propagation/model.ex`: + +```elixir +@doc """ +Trains the model on the given features and targets tensors. + +Options: + - :epochs — number of training epochs (default 50) + - :batch_size — samples per batch (default 256) + - :learning_rate — Adam learning rate (default 0.001) + - :patience — early stopping patience (default 5) + +Returns {trained_params, %{final_loss: float}}. +""" +def train(features, targets, opts \\ []) do + epochs = Keyword.get(opts, :epochs, 50) + batch_size = Keyword.get(opts, :batch_size, 256) + lr = Keyword.get(opts, :learning_rate, 0.001) + + model = build() + + loop = + model + |> Axon.Loop.trainer(:mean_squared_error, Polaris.Optimizers.adam(learning_rate: lr)) + |> Axon.Loop.metric(:mean_absolute_error) + + data = + features + |> Nx.to_batched(batch_size) + |> Stream.zip(Nx.to_batched(targets, batch_size)) + |> Stream.map(fn {x, y} -> {%{"features" => x}, y} end) + + params = Axon.Loop.run(loop, data, Axon.ModelState.empty(), + epochs: epochs, + compiler: EXLA + ) + + # Compute final loss + {_init, predict_fn} = Axon.build(model) + preds = predict_fn.(params, %{"features" => features}) + final_loss = Nx.mean(Nx.pow(Nx.subtract(preds, targets), 2)) |> Nx.to_number() + + {params, %{final_loss: final_loss}} +end + +@doc """ +Evaluates model on test data. Returns %{rmse: float, r_squared: float}. +""" +def evaluate(params, features, targets) do + model = build() + {_init, predict_fn} = Axon.build(model) + preds = predict_fn.(params, %{"features" => features}) + + mse = Nx.mean(Nx.pow(Nx.subtract(preds, targets), 2)) |> Nx.to_number() + rmse = :math.sqrt(mse) + + # R-squared + ss_res = Nx.sum(Nx.pow(Nx.subtract(targets, preds), 2)) |> Nx.to_number() + mean_y = Nx.mean(targets) |> Nx.to_number() + ss_tot = Nx.sum(Nx.pow(Nx.subtract(targets, mean_y), 2)) |> Nx.to_number() + r_squared = if ss_tot > 0, do: 1.0 - ss_res / ss_tot, else: 0.0 + + %{rmse: rmse, r_squared: r_squared} +end +``` + +**Step 4: Run model tests** + +Run: `mix test test/microwaveprop/propagation/model_test.exs --trace` +Expected: PASS + +**Step 5: Write the training mix task** + +```elixir +# lib/mix/tasks/propagation_train.ex +defmodule Mix.Tasks.PropagationTrain do + @shortdoc "Train ML model from propagation scores + HRRR data" + @moduledoc """ + Builds training data from propagation_scores joined to hrrr_profiles, + trains the Axon model, evaluates on held-out test set, and saves weights. + + ## Usage + + mix propagation.train + mix propagation.train --epochs 100 --batch-size 512 + """ + + use Mix.Task + + alias Microwaveprop.Propagation.Model + alias Microwaveprop.Repo + + @impl Mix.Task + def run(args) do + {opts, _, _} = OptionParser.parse(args, + strict: [epochs: :integer, batch_size: :integer, sample: :integer]) + + Application.put_env(:microwaveprop, Oban, queues: false) + Mix.Task.run("app.start") + + epochs = Keyword.get(opts, :epochs, 50) + batch_size = Keyword.get(opts, :batch_size, 256) + sample_size = Keyword.get(opts, :sample, 500_000) + + IO.puts("Loading training data (sample=#{sample_size})...") + {features, targets, band_counts} = load_training_data(sample_size) + + n = Nx.axis_size(features, 0) + IO.puts("Loaded #{n} samples") + for {band, count} <- Enum.sort(band_counts), do: IO.puts(" #{div(band, 1000)} GHz: #{count}") + + # 80/10/10 split + n_train = round(n * 0.8) + n_val = round(n * 0.1) + + train_x = Nx.slice(features, [0, 0], [n_train, 13]) + train_y = Nx.slice(targets, [0, 0], [n_train, 1]) + val_x = Nx.slice(features, [n_train, 0], [n_val, 13]) + val_y = Nx.slice(targets, [n_train, 0], [n_val, 1]) + test_x = Nx.slice(features, [n_train + n_val, 0], [n - n_train - n_val, 13]) + test_y = Nx.slice(targets, [n_train + n_val, 0], [n - n_train - n_val, 1]) + + IO.puts("\nSplit: train=#{n_train}, val=#{n_val}, test=#{n - n_train - n_val}") + IO.puts("Training for #{epochs} epochs, batch_size=#{batch_size}...\n") + + {params, train_metrics} = Model.train(train_x, train_y, + epochs: epochs, batch_size: batch_size) + + IO.puts("\nTraining complete. Final loss: #{Float.round(train_metrics.final_loss, 6)}") + + IO.puts("\n=== Evaluation ===") + val_metrics = Model.evaluate(params, val_x, val_y) + IO.puts("Validation: RMSE=#{Float.round(val_metrics.rmse * 100, 2)} points, R²=#{Float.round(val_metrics.r_squared, 4)}") + + test_metrics = Model.evaluate(params, test_x, test_y) + IO.puts("Test: RMSE=#{Float.round(test_metrics.rmse * 100, 2)} points, R²=#{Float.round(test_metrics.r_squared, 4)}") + + IO.puts("\nSaving model to priv/models/propagation_v1.nx...") + Model.save(params) + IO.puts("Done!") + end + + defp load_training_data(sample_size) do + query = """ + SELECT + h.surface_temp_c, + h.surface_dewpoint_c, + h.surface_pressure_mb, + h.min_refractivity_gradient, + h.hpbl_m, + h.pwat_mm, + EXTRACT(HOUR FROM ps.valid_time)::integer AS utc_hour, + EXTRACT(MONTH FROM ps.valid_time)::integer AS month, + ps.band_mhz, + ps.score + FROM propagation_scores ps + JOIN hrrr_profiles h + ON h.lat = ps.lat AND h.lon = ps.lon AND h.valid_time = ps.valid_time + WHERE h.surface_temp_c IS NOT NULL + AND h.surface_dewpoint_c IS NOT NULL + AND h.min_refractivity_gradient IS NOT NULL + ORDER BY RANDOM() + LIMIT $1 + """ + + result = Repo.query!(query, [sample_size], timeout: 600_000) + + band_counts = result.rows + |> Enum.frequencies_by(fn row -> Enum.at(row, 8) end) + + {feature_rows, target_rows} = + result.rows + |> Enum.map(fn [temp, dewpoint, pressure, gradient, hpbl, pwat, hour, month, band_mhz, score] -> + features = Model.encode_features(%{ + surface_temp_c: temp, + surface_dewpoint_c: dewpoint, + surface_pressure_mb: pressure, + min_refractivity_gradient: gradient, + hpbl_m: hpbl || 500.0, + pwat_mm: pwat || 20.0, + utc_hour: hour, + month: month, + freq_mhz: band_mhz + }) + target = score / 100.0 + {features, [target]} + end) + |> Enum.unzip() + + features = Nx.tensor(feature_rows, type: :f32) + targets = Nx.tensor(target_rows, type: :f32) + + {features, targets, band_counts} + end +end +``` + +**Step 6: Run the training** + +Run: `mix propagation.train --epochs 50 --sample 500000 2>&1 | tee docs/training_output.txt` + +**Step 7: Commit** + +```bash +git add lib/mix/tasks/propagation_train.ex lib/microwaveprop/propagation/model.ex +git add test/mix/tasks/propagation_train_test.exs test/microwaveprop/propagation/model_test.exs +git commit -m "Add ML training pipeline with Axon for propagation prediction" +``` + +--- + +### Task 5: Integrate ML Model into Grid Worker + +**Files:** +- Modify: `lib/microwaveprop/propagation.ex` (add ML prediction path) +- Modify: `lib/microwaveprop/propagation/model.ex` (add predict_score/2) +- Modify: `test/microwaveprop/propagation_test.exs` +- Modify: `test/microwaveprop/propagation/model_test.exs` + +**Step 1: Write failing test for Model.predict_score/2** + +```elixir +# In model_test.exs +describe "predict_score/2" do + test "returns integer score 0-100 from conditions map" do + params = Model.init() # random weights, but should still return valid range + conditions = %{ + surface_temp_c: 25.0, + surface_dewpoint_c: 15.0, + surface_pressure_mb: 1013.0, + min_refractivity_gradient: -80.0, + hpbl_m: 500.0, + pwat_mm: 25.0, + utc_hour: 12, + month: 6, + freq_mhz: 10_000 + } + + score = Model.predict_score(params, conditions) + assert is_integer(score) + assert score >= 0 and score <= 100 + end +end +``` + +**Step 2: Run test to verify it fails** + +Run: `mix test test/microwaveprop/propagation/model_test.exs --trace` +Expected: FAIL — predict_score/2 not defined + +**Step 3: Implement Model.predict_score/2** + +```elixir +@doc """ +Predicts a 0-100 propagation score from a conditions map. +Convenience wrapper around encode_features + predict. +""" +def predict_score(params, conditions) do + features = + conditions + |> encode_features() + |> Nx.tensor(type: :f32) + |> Nx.reshape({1, @feature_count}) + + params + |> predict(features) + |> Nx.squeeze() + |> Nx.multiply(100) + |> Nx.round() + |> Nx.to_number() + |> trunc() + |> max(0) + |> min(100) +end +``` + +**Step 4: Run test** + +Run: `mix test test/microwaveprop/propagation/model_test.exs --trace` +Expected: PASS + +**Step 5: Add ML path to Propagation.score_grid_point_with_data/5** + +Modify `lib/microwaveprop/propagation.ex` to try ML prediction first, fall back to algorithm: + +```elixir +# In score_grid_point_with_data/5, after building conditions map: +# Try ML model if loaded, fall back to algorithm scorer +for band_config <- BandConfig.all_bands() do + result = case ml_params() do + nil -> + Scorer.composite_score(conditions, band_config) + params -> + ml_conditions = Map.merge(conditions, %{freq_mhz: band_config.freq_mhz}) + score = Model.predict_score(params, ml_conditions) + %{score: score, factors: %{ml_predicted: true}} + end + + Map.put(result, :band_mhz, band_config.freq_mhz) +end +``` + +With a module-level cached model load: + +```elixir +@ml_params_key :propagation_ml_params + +def load_ml_model do + case Model.load() do + {:ok, params} -> + :persistent_term.put(@ml_params_key, params) + :ok + :error -> + :ok + end +end + +defp ml_params do + try do + :persistent_term.get(@ml_params_key) + rescue + ArgumentError -> nil + end +end +``` + +**Step 6: Update Propagation tests** + +Add tests verifying: +- Without ML model loaded, score_grid_point uses algorithm (existing behavior) +- With ML model loaded, score_grid_point uses ML predictions + +**Step 7: Run full test suite** + +Run: `mix test` +Expected: All pass + +**Step 8: Commit** + +```bash +git add lib/microwaveprop/propagation.ex lib/microwaveprop/propagation/model.ex +git add test/microwaveprop/propagation_test.exs test/microwaveprop/propagation/model_test.exs +git commit -m "Integrate ML model into grid worker with algorithm fallback" +``` + +--- + +### Task 6: Load ML Model at Application Startup + +**Files:** +- Modify: `lib/microwaveprop/application.ex` + +**Step 1: Add ML model loading to application start** + +After Repo and other services start, call `Propagation.load_ml_model()` to cache params in persistent_term. This is a no-op if the model file doesn't exist. + +**Step 2: Run full test suite** + +Run: `mix test` +Expected: All pass + +**Step 3: Run precommit** + +Run: `mix precommit` +Expected: All pass + +**Step 4: Commit** + +```bash +git add lib/microwaveprop/application.ex +git commit -m "Load ML model at application startup" +``` diff --git a/lib/microwaveprop/propagation/model.ex b/lib/microwaveprop/propagation/model.ex index 5cba5751..66e9c2a3 100644 --- a/lib/microwaveprop/propagation/model.ex +++ b/lib/microwaveprop/propagation/model.ex @@ -251,20 +251,22 @@ defmodule Microwaveprop.Propagation.Model do hour_rad = 2 * :math.pi() * local_hour / 24 month_rad = 2 * :math.pi() * (month - 1) / 12 + # Normalize atmospheric features to ~[0, 1] using physical bounds. + # Without this, features like pressure (~1013) dominate and cause NaN gradients. [ - temp_c, - dewpoint_c, - pressure_mb, - abs_humidity, - td_depression, - grad, - hpbl, - pwat, + (temp_c + 40) / 80, + (dewpoint_c + 40) / 80, + (pressure_mb - 960) / 80, + abs_humidity / 25.0, + td_depression / 40.0, + (grad + 300) / 300, + hpbl / 3000.0, + pwat / 60.0, :math.sin(hour_rad), :math.cos(hour_rad), :math.sin(month_rad), :math.cos(month_rad), - :math.log(freq_mhz) + :math.log(freq_mhz) / 13.0 ] end diff --git a/lib/mix/tasks/propagation_train.ex b/lib/mix/tasks/propagation_train.ex index 594b9f8d..f2d43503 100644 --- a/lib/mix/tasks/propagation_train.ex +++ b/lib/mix/tasks/propagation_train.ex @@ -89,7 +89,15 @@ defmodule Mix.Tasks.PropagationTrain do # Train IO.puts("Training for #{epochs} epochs...") {trained_state, train_metrics} = Model.train(train_features, train_targets, epochs: epochs, batch_size: batch_size) - IO.puts("Training complete. Final loss: #{Float.round(train_metrics.final_loss, 6)}\n") + final_loss = train_metrics.final_loss + + if is_float(final_loss) and not (final_loss != final_loss) do + IO.puts("Training complete. Final loss: #{Float.round(final_loss, 6)}\n") + else + IO.puts("Training complete. Final loss: NaN (training diverged)\n") + IO.puts("Try: --sample 100000 or reduce epochs") + System.halt(1) + end # Evaluate on validation set val_metrics = Model.evaluate(trained_state, val_features, val_targets) diff --git a/test/microwaveprop/propagation/model_test.exs b/test/microwaveprop/propagation/model_test.exs index f84e2770..afbbb99c 100644 --- a/test/microwaveprop/propagation/model_test.exs +++ b/test/microwaveprop/propagation/model_test.exs @@ -120,7 +120,7 @@ defmodule Microwaveprop.Propagation.ModelTest do test "log_freq_mhz is last feature" do features = Model.encode_features(%{freq_mhz: 10_000}) - assert_in_delta List.last(features), :math.log(10_000), 0.001 + assert_in_delta List.last(features), :math.log(10_000) / 13.0, 0.001 end end